remanence_anisotropy_magic.py

#!/usr/bin/env python

import matplotlib
#matplotlib.use('WXAgg')


import sys,pylab,scipy,os
from pylab import *

import pmag

def main():
    """
    NAME
        remanence_aniso_magic.py

    DESCRIPTION
        This program is similar to aarm_magic.py and atrm_magic.py with minor modifications.

        Converts magic measurement file with ATRM/AARM data to best-fit tensor (6 elements plus sigma)
        following Hext (1963), and calculates F-test statistics.
        
        Comments:
        - infield steps are marked with method codes LT-T-I:LP-AN-TRM; LT-AF-I:LP-AN-ARM
        - zerofield steps are marked with method codes LT-T-Z:LP-AN-TRM; LT-AF-Z:LP-AN-ARM
        - alteration check is marked with method codes LT-PTRM-I:LP-AN-TRM
        please notice;
        - ATRM: The program uses treatment_dc_field_phi/treatment_dc_field_theta columns to infer the direction of the applied field
                 (this is a change from atrm_magic.py)
        - ATRM: zerofield (baseline) magnetization is subtructed from all infield measurements
        - AARM: The program uses measurement number (running number) to to infer the direction of the applied field
                assuming the SIO protocol for 6,9,15 measurements scheme.
                See cookbook for diagram and details.
        - AARM: zerofield (baseline) are assumed to be before any infield, and the baseline is subtructed from the 
                subsequent infield magnetization.
      
    SYNTAX 
        remanence_aniso_magic.py [-h] [command line options]

    OPTIONS
        -h prints help message and quits
        -f FILE: specify input file, default is magic_measurements.txt

    INPUT  
        magic measurement file with ATRM and/or AARM data.
        if both types of measurements exist then the program calculates both.
        
    OUTPUT  
        rmag_anisotropy.log
                -I- information
                -W- Warning
                -E- Error
        rmag_anistropy.txt:
            This file contains in addition to some some magic information the following:
                - anistropy tensor s1 to s6 normalized by the trace:
                    |Mx|   |s1 s4 s6|   |Bx|
                    |My| = |s4 s2 s5| . |By|
                    |Mz|   |s6 s5 s3|   |Bz|
                - anisotropy_sigma (Hext, 1963)
                - anisotropy_alt (altertion check for ATRM in units of %):
                    100* [abs(M_first-Mlast)/max(M_first,M_last)]
                -                     
        rmag_results.txt:
            This file contains in addition to some  magic information the follow(ing:
                - anisotropy_t1,anisotropy_t2,anisotropy_t3 : eigenvalues
                - anisotropy_v*_dec,anisotropy_v*_inc: declination/inclination of the eigenvectors
                - anisotropy_ftest,anisotropy_ftest12,anisotropy_ftest13 
                - (the crtical F for 95% confidence level of anistropy is given in result_description column).
                
                
    """


    #==================================================================================
    
    meas_file="magic_measurements.txt"
    args=sys.argv
    dir_path='.'
    #
    # get name of file from command line
    #
    if '-WD' in args:
        ind=args.index('-WD')
        dir_path=args[ind+1]
    if "-h" in args:
        print main.__doc__
        sys.exit()
    if "-f" in args:
        ind=args.index("-f")
        meas_file=sys.argv[ind+1]
    else:
        meas_file=dir_path+'/'+meas_file
    WD=dir_path

    
    #======================================
    # functions
    #======================================

    def get_Data(magic_file):
    
        #------------------------------------------------
        # Read magic measurement file and sort to blocks
        #------------------------------------------------
        Data={}
        try:
            meas_data,file_type=pmag.magic_read(magic_file)
        except:
            print "-E- ERROR: Cant read magic_measurement.txt file. File is corrupted."
            return Data
            
        # get list of unique specimen names
        
        #sids=pmag.get_specs(meas_data) # samples ID's
        
        for rec in meas_data:
            s=rec["er_specimen_name"]
            method_codes= rec["magic_method_codes"].strip('\n')
            method_codes.replace(" ","")
            methods=method_codes.split(":")
            if "LP-AN-TRM" in methods:
                if s not in Data.keys():
                    Data[s]={}
                if 'atrmblock' not in Data[s].keys():
                    Data[s]['atrmblock']=[]
                Data[s]['atrmblock'].append(rec)
            
            
            if "LP-AN-ARM" in methods:
                if s not in Data.keys():
                    Data[s]={}
                if 'aarmblock' not in Data[s].keys():
                    Data[s]['aarmblock']=[]
                Data[s]['aarmblock'].append(rec)
        return (Data)        

    #======================================
    # better to put this one in pmagpy
    #======================================
        
    def calculate_aniso_parameters(B,K):
    
        aniso_parameters={}
        S_bs=dot(B,K)
        
        # normalize by trace
        trace=S_bs[0]+S_bs[1]+S_bs[2]
        S_bs=S_bs/trace
        s1,s2,s3,s4,s5,s6=S_bs[0],S_bs[1],S_bs[2],S_bs[3],S_bs[4],S_bs[5]
        s_matrix=[[s1,s4,s6],[s4,s2,s5],[s6,s5,s3]]
        
        # calculate eigen vector,
        t,evectors=eig(s_matrix)
        # sort vectors
        t=list(t)
        t1=max(t)
        ix_1=t.index(t1)
        t3=min(t)
        ix_3=t.index(t3)
        for tt in range(3):
            if t[tt]!=t1 and t[tt]!=t3:
                t2=t[tt]
                ix_2=t.index(t2)
                
        v1=[evectors[0][ix_1],evectors[1][ix_1],evectors[2][ix_1]]
        v2=[evectors[0][ix_2],evectors[1][ix_2],evectors[2][ix_2]]
        v3=[evectors[0][ix_3],evectors[1][ix_3],evectors[2][ix_3]]
    
    
        DIR_v1=pmag.cart2dir(v1)
        DIR_v2=pmag.cart2dir(v2)
        DIR_v3=pmag.cart2dir(v3)
    
                            
        aniso_parameters['anisotropy_s1']="%f"%s1
        aniso_parameters['anisotropy_s2']="%f"%s2
        aniso_parameters['anisotropy_s3']="%f"%s3
        aniso_parameters['anisotropy_s4']="%f"%s4
        aniso_parameters['anisotropy_s5']="%f"%s5
        aniso_parameters['anisotropy_s6']="%f"%s6
        aniso_parameters['anisotropy_degree']="%f"%(t1/t3)
        aniso_parameters['anisotropy_t1']="%f"%t1
        aniso_parameters['anisotropy_t2']="%f"%t2
        aniso_parameters['anisotropy_t3']="%f"%t3
        aniso_parameters['anisotropy_v1_dec']="%.1f"%DIR_v1[0]
        aniso_parameters['anisotropy_v1_inc']="%.1f"%DIR_v1[1]
        aniso_parameters['anisotropy_v2_dec']="%.1f"%DIR_v2[0]
        aniso_parameters['anisotropy_v2_inc']="%.1f"%DIR_v2[1]
        aniso_parameters['anisotropy_v3_dec']="%.1f"%DIR_v3[0]
        aniso_parameters['anisotropy_v3_inc']="%.1f"%DIR_v3[1]
    
        # modified from pmagpy:
        if len(K)/3==9 or len(K)/3==6 or len(K)/3==15:
            n_pos=len(K)/3
            tmpH = Matrices[n_pos]['tmpH']
            a=s_matrix
            S=0.
            comp=zeros((n_pos*3),'f')
            for i in range(n_pos):
                for j in range(3):
                    index=i*3+j
                    compare=a[j][0]*tmpH[i][0]+a[j][1]*tmpH[i][1]+a[j][2]*tmpH[i][2]
                    comp[index]=compare
            for i in range(n_pos*3):
                d=K[i]/trace - comp[i] # del values
                S+=d*d
            nf=float(n_pos*3-6) # number of degrees of freedom
            if S >0: 
                sigma=math.sqrt(S/nf)
            hpars=pmag.dohext(nf,sigma,[s1,s2,s3,s4,s5,s6])
            
            aniso_parameters['anisotropy_sigma']="%f"%sigma
            aniso_parameters['anisotropy_ftest']="%f"%hpars["F"]
            aniso_parameters['anisotropy_ftest12']="%f"%hpars["F12"]
            aniso_parameters['anisotropy_ftest23']="%f"%hpars["F23"]
            aniso_parameters['result_description']="Critical F: %s"%(hpars['F_crit'])
            aniso_parameters['anisotropy_F_crit']="%f"%float(hpars['F_crit'])
            aniso_parameters['anisotropy_n']=n_pos
            
        return(aniso_parameters)
    
    
    #======================================
    # Main
    #======================================
          
    
        
    aniso_logfile=open(WD+"/rmag_anisotropy.log",'w')
    aniso_logfile.write("------------------------\n")
    aniso_logfile.write( "-I- Start rmag anisrotropy script\n")
    aniso_logfile.write( "------------------------\n")

    Data=get_Data(meas_file)
    #try:    
    #    Data=get_Data(meas_file)
    #except:
    #    aniso_logfile.write( "-E- Cant open measurement file %s\n" %meas_file)
    #    print "-E- Cant open measurement file %s\n exiting" %meas_file
    #    exit()
        
    aniso_logfile.write( "-I-  Open measurement file %s\n" %meas_file)
    
        
    Data_anisotropy={}                
    specimens=Data.keys()
    specimens.sort()
    
    
    
    #-----------------------------------
    # Prepare rmag_anisotropy.txt file for writing
    #-----------------------------------
    
    rmag_anisotropy_file =open(WD+"/rmag_anisotropy.txt",'w')
    rmag_anisotropy_file.write("tab\trmag_anisotropy\n")
    
    rmag_results_file =open(WD+"/rmag_results.txt",'w')
    rmag_results_file.write("tab\trmag_results\n")
    
    rmag_anistropy_header=['er_specimen_name','er_sample_name','er_site_name','anisotropy_type','anisotropy_n','anisotropy_description','anisotropy_s1','anisotropy_s2','anisotropy_s3','anisotropy_s4','anisotropy_s5','anisotropy_s6','anisotropy_sigma','anisotropy_alt','magic_experiment_names','magic_method_codes','rmag_anisotropy_name']
    
    String=""
    for i in range (len(rmag_anistropy_header)):
        String=String+rmag_anistropy_header[i]+'\t'
    rmag_anisotropy_file.write(String[:-1]+"\n")
    
    
    
    rmag_results_header=['er_specimen_names','er_sample_names','er_site_names','anisotropy_type','magic_method_codes','magic_experiment_names','result_description','anisotropy_t1','anisotropy_t2','anisotropy_t3','anisotropy_ftest','anisotropy_ftest12','anisotropy_ftest23',\
                            'anisotropy_v1_dec','anisotropy_v1_inc','anisotropy_v2_dec','anisotropy_v2_inc','anisotropy_v3_dec','anisotropy_v3_inc']
    
    
    String=""
    for i in range (len(rmag_results_header)):
        String=String+rmag_results_header[i]+'\t'
    rmag_results_file.write(String[:-1]+"\n")
    
    #-----------------------------------
    # Matrices definitions:
    # A design matrix
    # B dot(inv(dot(A.transpose(),A)),A.transpose())
    # tmpH is used for sigma calculation (9,15 measurements only)
    # 
    #  Anisotropy tensor:
    #
    # |Mx|   |s1 s4 s6|   |Bx|
    # |My| = |s4 s2 s5| . |By|
    # |Mz|   |s6 s5 s3|   |Bz|
    #
    # A matrix (measurement matrix):
    # Each mesurement yields three lines in "A" matrix
    #
    # |Mi  |   |Bx 0  0   By  0   Bz|   |s1|
    # |Mi+1| = |0  By 0   Bx  Bz  0 | . |s2|
    # |Mi+2|   |0  0  Bz  0   By  Bx|   |s3|
    #                                   |s4|
    #                                   |s5|
    #
    #-----------------------------------
    
    Matrices={}
    
    for n_pos in [6,9,15]:
    
        Matrices[n_pos]={}
        
        A=zeros((n_pos*3,6),'f')
    
        if n_pos==6:
            positions=[[0.,0.,1.],[90.,0.,1.],[0.,90.,1.],\
                    [180.,0.,1.],[270.,0.,1.],[0.,-90.,1.]]
    
        if n_pos==15:
            positions=[[315.,0.,1.],[225.,0.,1.],[180.,0.,1.],[135.,0.,1.],[45.,0.,1.],\
                    [90.,-45.,1.],[270.,-45.,1.],[270.,0.,1.],[270.,45.,1.],[90.,45.,1.],\
                    [180.,45.,1.],[180.,-45.,1.],[0.,-90.,1.],[0,-45.,1.],[0,45.,1.]]
        if n_pos==9:
            positions=[[315.,0.,1.],[225.,0.,1.],[180.,0.,1.],\
                    [90.,-45.,1.],[270.,-45.,1.],[270.,0.,1.],\
                    [180.,45.,1.],[180.,-45.,1.],[0.,-90.,1.]]
    
        
        tmpH=zeros((n_pos,3),'f') # define tmpH
        for i in range(len(positions)):
            CART=pmag.dir2cart(positions[i])
            a=CART[0];b=CART[1];c=CART[2]
            A[3*i][0]=a
            A[3*i][3]=b
            A[3*i][5]=c
    
            A[3*i+1][1]=b
            A[3*i+1][3]=a
            A[3*i+1][4]=c
    
            A[3*i+2][2]=c
            A[3*i+2][4]=b
            A[3*i+2][5]=a
            
            tmpH[i][0]=CART[0]
            tmpH[i][1]=CART[1]
            tmpH[i][2]=CART[2]
    
        B=dot(inv(dot(A.transpose(),A)),A.transpose())
    
        Matrices[n_pos]['A']=A
        Matrices[n_pos]['B']=B
        Matrices[n_pos]['tmpH']=tmpH
    
    
    
    
    
    for specimen in specimens:
    
        if 'atrmblock' in Data[specimen].keys():
            
            #-----------------------------------
            # aTRM 6 positions
            #-----------------------------------
                
            aniso_logfile.write("-I- Start calculating ATRM tensor for specimen %s\n "%specimen)
            atrmblock=Data[specimen]['atrmblock']
            if len(atrmblock)<6:
                aniso_logfile.write("-W- specimen %s has not enough measurementf for ATRM calculation\n"%specimen)
                continue
            
            B=Matrices[6]['B']
                                
            Reject_specimen = False
    
            # The zero field step is a "baseline"
    
            # Search the baseline in the ATRM measurement
            
            baseline=""
            Alteration_check=""
            Alteration_check_index=""
            baselines=[]
    
            # search for baseline in atrm blocks
            for rec in atrmblock:
                dec=float(rec['measurement_dec'])
                inc=float(rec['measurement_inc'])
                moment=float(rec['measurement_magn_moment'])
                # find the temperature of the atrm
                if float(rec['treatment_dc_field'])!=0 and float(rec['treatment_temp'])!=273:
                    atrm_temperature=float(rec['treatment_temp'])
                # find baseline
                if float(rec['treatment_dc_field'])==0 and float(rec['treatment_temp'])!=273:
                    baselines.append(array(pmag.dir2cart([dec,inc,moment])))
                # Find alteration check
                #print rec['measurement_number']
            
            if len(baselines)!=0:
                aniso_logfile.write( "-I- found ATRM baseline for specimen %s\n"%specimen)
                baselines=array(baselines)
                baseline=array([mean(baselines[:,0]),mean(baselines[:,1]),mean(baselines[:,2])])                                 
                
            else:
                baseline=zeros(3,'f')
                aniso_logfile.write( "-I- No aTRM baseline for specimen %s\n"%specimen)
                        
            # sort measurements
            
            M=zeros([6,3],'f')
            
            for rec in atrmblock:
    
                dec=float(rec['measurement_dec'])
                inc=float(rec['measurement_inc'])
                moment=float(rec['measurement_magn_moment'])
                CART=array(pmag.dir2cart([dec,inc,moment]))-baseline
                
                if float(rec['treatment_dc_field'])==0: # Ignore zero field steps
                    continue
                if  "LT-PTRM-I" in rec['magic_method_codes'].split(":"): #  alteration check
                    Alteration_check=CART
                    Alteration_check_dc_field_phi=float(rec['treatment_dc_field_phi'])
                    Alteration_check_dc_field_theta=float(rec['treatment_dc_field_theta'])
                    if Alteration_check_dc_field_phi==0 and Alteration_check_dc_field_theta==0 :
                        Alteration_check_index=0
                    if Alteration_check_dc_field_phi==90 and Alteration_check_dc_field_theta==0 :
                        Alteration_check_index=1
                    if Alteration_check_dc_field_phi==0 and Alteration_check_dc_field_theta==90 :
                        Alteration_check_index=2
                    if Alteration_check_dc_field_phi==180 and Alteration_check_dc_field_theta==0 :
                        Alteration_check_index=3
                    if Alteration_check_dc_field_phi==270 and Alteration_check_dc_field_theta==0 :
                        Alteration_check_index=4
                    if Alteration_check_dc_field_phi==0 and Alteration_check_dc_field_theta==-90 :
                        Alteration_check_index=5
                    aniso_logfile.write(  "-I- found alteration check  for specimen %s\n"%specimen)
                    continue
                
                treatment_dc_field_phi=float(rec['treatment_dc_field_phi'])
                treatment_dc_field_theta=float(rec['treatment_dc_field_theta'])
                treatment_dc_field=float(rec['treatment_dc_field'])
                
                #+x, M[0]
                if treatment_dc_field_phi==0 and treatment_dc_field_theta==0 :
                    M[0]=CART
                #+Y , M[1]
                if treatment_dc_field_phi==90 and treatment_dc_field_theta==0 :
                    M[1]=CART
                #+Z , M[2]
                if treatment_dc_field_phi==0 and treatment_dc_field_theta==90 :
                    M[2]=CART
                #-x, M[3]
                if treatment_dc_field_phi==180 and treatment_dc_field_theta==0 :
                    M[3]=CART
                #-Y , M[4]
                if treatment_dc_field_phi==270 and treatment_dc_field_theta==0 :
                    M[4]=CART
                #-Z , M[5]
                if treatment_dc_field_phi==0 and treatment_dc_field_theta==-90 :
                    M[5]=CART
        
            # check if at least one measurement in missing
            for i in range(len(M)):
                if M[i][0]==0 and M[i][1]==0 and M[i][2]==0: 
                    aniso_logfile.write( "-E- ERROR: missing atrm data for specimen %s\n"%(specimen))
                    Reject_specimen=True
    
            # alteration check        
    
            anisotropy_alt=0
            if Alteration_check!="":
                for i in range(len(M)):
                    if Alteration_check_index==i:
                        M_1=sqrt(sum((array(M[i])**2)))
                        M_2=sqrt(sum(Alteration_check**2))
                        diff=abs(M_1-M_2)
                        diff_ratio=diff/max(M_1,M_2)
                        diff_ratio_perc=100*diff_ratio
                        if diff_ratio_perc > anisotropy_alt:
                            anisotropy_alt=diff_ratio_perc
            else:
                aniso_logfile.write( "-W- Warning: no alteration check for specimen %s \n "%specimen )
    
            # Check for maximum difference in anti parallel directions.
            # if the difference between the two measurements is more than maximum_diff
            # The specimen is rejected
            
            # i.e. +x versus -x, +y versus -y, etc.s
    
            for i in range(3):
                M_1=sqrt(sum(array(M[i])**2))
                M_2=sqrt(sum(array(M[i+3])**2))
                
                diff=abs(M_1-M_2)
                diff_ratio=diff/max(M_1,M_2)
                diff_ratio_perc=100*diff_ratio
                
                if diff_ratio_perc>anisotropy_alt:
                    anisotropy_alt=diff_ratio_perc
                    
            if not Reject_specimen:
            
                # K vector (18 elements, M1[x], M1[y], M1[z], ... etc.) 
                K=zeros(18,'f')
                K[0],K[1],K[2]=M[0][0],M[0][1],M[0][2]
                K[3],K[4],K[5]=M[1][0],M[1][1],M[1][2]
                K[6],K[7],K[8]=M[2][0],M[2][1],M[2][2]
                K[9],K[10],K[11]=M[3][0],M[3][1],M[3][2]
                K[12],K[13],K[14]=M[4][0],M[4][1],M[4][2]
                K[15],K[16],K[17]=M[5][0],M[5][1],M[5][2]
    
                if specimen not in Data_anisotropy.keys():
                    Data_anisotropy[specimen]={}
                aniso_parameters=calculate_aniso_parameters(B,K)
                Data_anisotropy[specimen]['ATRM']=aniso_parameters
                Data_anisotropy[specimen]['ATRM']['anisotropy_alt']="%.2f"%anisotropy_alt               
                Data_anisotropy[specimen]['ATRM']['anisotropy_type']="ATRM"
                Data_anisotropy[specimen]['ATRM']['er_sample_name']=atrmblock[0]['er_sample_name']
                Data_anisotropy[specimen]['ATRM']['er_specimen_name']=specimen
                Data_anisotropy[specimen]['ATRM']['er_site_name']=atrmblock[0]['er_site_name']
                Data_anisotropy[specimen]['ATRM']['anisotropy_description']='Hext statistics adapted to ATRM'
                Data_anisotropy[specimen]['ATRM']['magic_experiment_names']=specimen+";ATRM"
                Data_anisotropy[specimen]['ATRM']['magic_method_codes']="LP-AN-TRM:AE-H"
                Data_anisotropy[specimen]['ATRM']['rmag_anisotropy_name']=specimen
    
    
        if 'aarmblock' in Data[specimen].keys():    
    
            #-----------------------------------
            # AARM - 6, 9 or 15 positions
            #-----------------------------------
                
            aniso_logfile.write( "-I- Start calculating AARM tensors specimen %s\n"%specimen)
    
            aarmblock=Data[specimen]['aarmblock']
            if len(aarmblock)<12:
                aniso_logfile.write( "-W- WARNING: not enough aarm measurement for specimen %s\n"%specimen)
                continue
            elif len(aarmblock)==12:
                n_pos=6
                B=Matrices[6]['B']
                M=zeros([6,3],'f')
            elif len(aarmblock)==18:
                n_pos=9
                B=Matrices[9]['B']
                M=zeros([9,3],'f')
            # 15 positions
            elif len(aarmblock)==30:
                n_pos=15
                B=Matrices[15]['B']
                M=zeros([15,3],'f')
            else:
                aniso_logfile.write( "-E- ERROR: number of measurements in aarm block is incorrect sample %s\n"%specimen)
                continue
                
            Reject_specimen = False
    
            for i in range(n_pos):
                for rec in aarmblock:
                    if float(rec['measurement_number'])==i*2+1:
                        dec=float(rec['measurement_dec'])
                        inc=float(rec['measurement_inc'])
                        moment=float(rec['measurement_magn_moment'])                    
                        M_baseline=array(pmag.dir2cart([dec,inc,moment]))
                        
                    if float(rec['measurement_number'])==i*2+2:
                        dec=float(rec['measurement_dec'])
                        inc=float(rec['measurement_inc'])
                        moment=float(rec['measurement_magn_moment'])                    
                        M_arm=array(pmag.dir2cart([dec,inc,moment]))
                M[i]=M_arm-M_baseline
    
                
            K=zeros(3*n_pos,'f')
            for i in range(n_pos):
                K[i*3]=M[i][0]
                K[i*3+1]=M[i][1]
                K[i*3+2]=M[i][2]            
    
            if specimen not in Data_anisotropy.keys():
                Data_anisotropy[specimen]={}
            aniso_parameters=calculate_aniso_parameters(B,K)
            Data_anisotropy[specimen]['AARM']=aniso_parameters
            Data_anisotropy[specimen]['AARM']['anisotropy_alt']=""               
            Data_anisotropy[specimen]['AARM']['anisotropy_type']="AARM"
            Data_anisotropy[specimen]['AARM']['er_sample_name']=aarmblock[0]['er_sample_name']
            Data_anisotropy[specimen]['AARM']['er_site_name']=aarmblock[0]['er_site_name']
            Data_anisotropy[specimen]['AARM']['er_specimen_name']=specimen
            Data_anisotropy[specimen]['AARM']['anisotropy_description']='Hext statistics adapted to AARM'
            Data_anisotropy[specimen]['AARM']['magic_experiment_names']=specimen+";AARM"
            Data_anisotropy[specimen]['AARM']['magic_method_codes']="LP-AN-ARM:AE-H"
            Data_anisotropy[specimen]['AARM']['rmag_anisotropy_name']=specimen
            
    
    #-----------------------------------   
    
    specimens=Data_anisotropy.keys()
    specimens.sort
    
    # remove previous anistropy data, and replace with the new one:
    s_list=Data.keys()
    for sp in s_list:
        if 'AniSpec' in Data[sp].keys():
            del  Data[sp]['AniSpec']
    for specimen in specimens:
        # if both AARM and ATRM axist prefer the AARM !!
        if 'AARM' in Data_anisotropy[specimen].keys():
            TYPES=['AARM']
        if 'ATRM' in Data_anisotropy[specimen].keys():
            TYPES=['ATRM']
        if  'AARM' in Data_anisotropy[specimen].keys() and 'ATRM' in Data_anisotropy[specimen].keys():
            TYPES=['ATRM','AARM']
            aniso_logfile.write( "-W- WARNING: both aarm and atrm data exist for specimen %s. using AARM by default. If you prefer using one of them, delete the other!\n"%specimen)
        for TYPE in TYPES:
            String=""
            for i in range (len(rmag_anistropy_header)):
                try:
                    String=String+Data_anisotropy[specimen][TYPE][rmag_anistropy_header[i]]+'\t'
                except:
                    String=String+"%f"%(Data_anisotropy[specimen][TYPE][rmag_anistropy_header[i]])+'\t'
            rmag_anisotropy_file.write(String[:-1]+"\n")
    
            String=""
            Data_anisotropy[specimen][TYPE]['er_specimen_names']=Data_anisotropy[specimen][TYPE]['er_specimen_name']
            Data_anisotropy[specimen][TYPE]['er_sample_names']=Data_anisotropy[specimen][TYPE]['er_sample_name']
            Data_anisotropy[specimen][TYPE]['er_site_names']=Data_anisotropy[specimen][TYPE]['er_site_name']
            for i in range (len(rmag_results_header)):
                try:
                    String=String+Data_anisotropy[specimen][TYPE][rmag_results_header[i]]+'\t'
                except:
                    String=String+"%f"%(Data_anisotropy[specimen][TYPE][rmag_results_header[i]])+'\t'
            rmag_results_file.write(String[:-1]+"\n")
    
            if 'AniSpec' not in Data[specimen]:
                Data[specimen]['AniSpec']={}
            Data[specimen]['AniSpec'][TYPE]=Data_anisotropy[specimen][TYPE]
    
    aniso_logfile.write("------------------------\n")
    aniso_logfile.write("-I-  remanence_aniso_magic script finished sucsessfuly\n")
    aniso_logfile.write( "------------------------\n")
    
    rmag_anisotropy_file.close()
    print "Anisotropy tensors elements are saved in rmag_anistropy.txt"
    print "Other anisotropy statistics are saved in rmag_results.txt"
    print "log file is  in rmag_anisotropy.log"


    
main()