Changeset 1010


Ignore:
Timestamp:
Jul 24, 2013 9:04:01 PM (10 years ago)
Author:
toby
Message:

more CIF work

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIIO.py

    r1006 r1010  
    17251725
    17261726    def UpdateParameters(self,Type=None,Wave=None):
    1727         'Revise the '
     1727        'Revise the instrument parameters'
    17281728        if Type is not None:
    17291729            self.Parameters[0]['Type'] = [Type,Type]

  • trunk/GSASIImath.py

    r1003 r1010  
    12981298            valoff = 0
    12991299    if esd != 0:
     1300        if valoff+esdoff < 0:
     1301            valoff = esdoff = 0
    13001302        out = ("{:."+str(valoff+esdoff)+"f}").format(value/10**valoff) # format the value
    13011303    elif valoff != 0: # esd = 0; exponential notation ==> esdoff decimal places
    13021304        out = ("{:."+str(esdoff)+"f}").format(value/10**valoff) # format the value
    13031305    else: # esd = 0; non-exponential notation ==> esdoff+1 significant digits
    1304         extra = -math.log10(abs(value))
     1306        if abs(value) > 0:           
     1307            extra = -math.log10(abs(value))
     1308        else:
     1309            extra = 0
    13051310        if extra > 0: extra += 1
    13061311        out = ("{:."+str(max(0,esdoff+int(extra)))+"f}").format(value) # format the value

  • trunk/exports/G2cif.py

    r990 r1010  
    33'''
    44
     5# TODO: need a mechanism for editing of instrument names, bond pub flags, templates,...
     6
    57import datetime as dt
    68import os.path
     9import numpy as np
    710import GSASIIIO as G2IO
    811#reload(G2IO)
    … …  
    110113
    111114            # include an overall profile r-factor, if there is more than one powder histogram
    112             if self.npowder > 1:
     115            if len(self.powderDict) > 1:
    113116                WriteCIFitem('\n# OVERALL POWDER R-FACTOR')
    114117                try:
    … …  
    192195                    if s != "": s += '\n'
    193196                    s += 'March-Dollase correction'
    194                     if self.npowder > 1:
     197                    if len(self.powderDict) > 1:
    195198                        s += ', histogram '+str(Histogram['hId']+1)
    196199                    s += ' coef. = ' + G2mth.ValEsd(self.parmDict[aname],sig)
    … …  
    199202                    if s != "": s += '\n'
    200203                    s += 'Simple spherical harmonic correction'
    201                     if self.npowder > 1:
     204                    if len(self.powderDict) > 1:
    202205                        s += ', histogram '+str(Histogram['hId']+1)
    203206                    s += ' Order = '+str(hapData['Pref.Ori.'][4])+'\n'
    … …  
    602605            # see wrpowdhist.for & wrreflist.for
    603606           
    604             #refprx = '_refln.' # mm
    605             refprx = '_refln_' # normal
    606 
    607607            if 'Lam1' in inst:
    608608                ratio = self.parmDict.get('I(L2)/I(L1)',inst['I(L2)/I(L1)'][1])
    … …  
    627627                             PutInCol('2',5))               
    628628            else:
    629                 lam1 = self.parmDict.get('Lam',inst['Lam'])
     629                lam1 = self.parmDict.get('Lam',inst['Lam'][1])
    630630                slam1 = self.sigDict.get('Lam',-0.00009)
    631631                WriteCIFitem('_diffrn_radiation_wavelength',G2mth.ValEsd(lam1,slam1))
    … …  
    669669            # WriteCIFitem('_refine_ls_R_Fsqd_factor','?')
    670670
    671             print histblk['Instrument Parameters'][0]['Type']
    672            
    673671            if histblk['Instrument Parameters'][0]['Type'][1][1] == 'X':
    674672                WriteCIFitem('_diffrn_radiation_probe','x-ray')
    … …  
    709707            #CALL WRVAL(IUCIF,'_gsas_exptl_extinct_corr_T_max',TEXT(11:20))
    710708
    711             if not oneblock:
    712                 # instrumental profile terms go here
    713                 WriteCIFitem('_pd_proc_ls_profile_function','?')
    714 
    715             raise Exception, "testing"
    716             phasenam = self.Phases.keys()[0]
    717             for key in self.Phases[phasenam]['Histograms']:
    718                 print key
    719                 print '------------'
    720                 print self.Phases[phasenam]['Histograms'][key]
    721             print histblk.keys()
    722             for key in histblk:
    723                 print key,histblk[key]
    724             print inst
    725             #print self.parmDict.keys()
    726             #print self.sigDict.keys()
    727 
    728             #print 'Data'
    729             #for item in histblk['Data']:
    730             #    print item
    731                 #try:
    732                 #    print key,histblk[key].keys()
    733                 #except:
    734                 #    print key
    735                 #    print histblk[key]
    736             #print 'Background'
    737             print histblk['Reflection Lists']['Garnet'][1]
    738             for i in range(0,80):
    739                 for j in [0,1,2,13]:
    740                     print histblk['Reflection Lists']['Garnet'][i][j],
    741                 print
    742                 #print histblk['Reflection Lists']['Garnet'][i][12].shape
    743                 #print histblk['Reflection Lists']['Garnet'][i][14]
    744             #print histblk['Background'][0]
    745             #print histblk['Background'][1]
    746             import numpy as np
    747             refList = np.array([refl[:11] for refl in histblk['Reflection Lists']['Garnet']])
    748             #refList = histblk['Reflection Lists']['Garnet']
    749             Icorr = np.array([refl[13] for refl in histblk['Reflection Lists']['Garnet']])
    750             FO2 = np.array([refl[8] for refl in histblk['Reflection Lists']['Garnet']])
    751             print Icorr
    752             I100 = refList.T[8]*Icorr
    753             print I100
    754             print I100/max(I100)
    755             Icorr = np.array([refl[13] for refl in histblk['Reflection Lists']['Garnet']]) * np.array([refl[8] for refl in histblk['Reflection Lists']['Garnet']])
    756             print I100/max(I100)
    757 
     709            if not oneblock:                 # instrumental profile terms go here
     710                WriteCIFitem('_pd_proc_ls_profile_function',
     711                             FormatInstProfile(histblk["Instrument Parameters"]))
     712
     713            #refprx = '_refln.' # mm
     714            refprx = '_refln_' # normal
    758715            WriteCIFitem('\n# STRUCTURE FACTOR TABLE')           
    759             WriteCIFitem('loop_' +
    760                          '\n\t' + refprx + 'index_h' +
     716            # compute maximum intensity reflection
     717            Imax = 0
     718            for phasenam in histblk['Reflection Lists']:
     719                scale = self.Phases[phasenam]['Histograms'][histlbl]['Scale'][0]
     720                Icorr = np.array([refl[13] for refl in histblk['Reflection Lists'][phasenam]])
     721                FO2 = np.array([refl[8] for refl in histblk['Reflection Lists'][phasenam]])
     722                I100 = scale*FO2*Icorr
     723                Imax = max(Imax,max(I100))
     724
     725            WriteCIFitem('loop_')
     726            if len(histblk['Reflection Lists'].keys()) > 1:
     727                WriteCIFitem('\t_pd_refln_phase_id')
     728            WriteCIFitem('\t' + refprx + 'index_h' +
    761729                         '\n\t' + refprx + 'index_k' +
    762730                         '\n\t' + refprx + 'index_l' +
    763                          '\n\t_pd_refln_wavelength_id' +
    764                          '\n\t_pd_refln_phase_id' +
    765                          '\n\t' + refprx + 'status' +
     731#                         '\n\t_pd_refln_wavelength_id' +
     732#                         '\n\t' + refprx + 'status' +
    766733                         '\n\t' + refprx + 'F_squared_meas' +
    767                          '\n\t' + refprx + 'F_squared_sigma' +
     734#                         '\n\t' + refprx + 'F_squared_sigma' +
    768735                         '\n\t' + refprx + 'F_squared_calc' +
    769736                         '\n\t' + refprx + 'phase_calc' +
    770                          '\n\t_pd_refln_d_spacing' +
    771                          '\n\t_gsas_i100_meas')
    772 
    773             WriteCIFitem('_reflns_number_total', text)
    774             WriteCIFitem('_reflns_limit_h_min', text)
    775             WriteCIFitem('_reflns_limit_h_max', text)
    776             WriteCIFitem('_reflns_limit_k_min', text)
    777             WriteCIFitem('_reflns_limit_k_max', text)
    778             WriteCIFitem('_reflns_limit_l_min', text)
    779             WriteCIFitem('_reflns_limit_l_max', text)
    780             WriteCIFitem('_reflns_d_resolution_high', text)
    781             WriteCIFitem('_reflns_d_resolution_low', text)
    782            
     737                         '\n\t_pd_refln_d_spacing')
     738            if Imax > 0:
     739                WriteCIFitem('\t_gsas_i100_meas')
     740
     741            refcount = 0
     742            hklmin = None
     743            hklmax = None
     744            dmax = None
     745            dmin = None
     746            for phasenam in histblk['Reflection Lists']:
     747                scale = self.Phases[phasenam]['Histograms'][histlbl]['Scale'][0]
     748                phaseid = self.Phases[phasenam]['pId']
     749                refcount += len(histblk['Reflection Lists'][phasenam])
     750                for ref in histblk['Reflection Lists'][phasenam]:
     751                    if hklmin is None:
     752                        hklmin = ref[0:3]
     753                        hklmax = ref[0:3]
     754                        dmax = dmin = ref[4]
     755                    if len(histblk['Reflection Lists'].keys()) > 1:
     756                        s = PutInCol(phaseid,2)
     757                    else:
     758                        s = ""
     759                    for i,hkl in enumerate(ref[0:3]):
     760                        hklmax[i] = max(hkl,hklmax[i])
     761                        hklmin[i] = min(hkl,hklmin[i])
     762                        s += PutInCol(int(hkl),4)
     763                    for I in ref[8:10]:
     764                        s += PutInCol(G2mth.ValEsd(I,-0.0009),10)
     765                    s += PutInCol(G2mth.ValEsd(ref[10],-0.9),7)
     766                    dmax = max(dmax,ref[4])
     767                    dmin = min(dmin,ref[4])
     768                    s += PutInCol(G2mth.ValEsd(ref[4],-0.009),8)
     769                    if Imax > 0:
     770                        I100 = 100.*scale*ref[8]*ref[13]/Imax
     771                        s += PutInCol(G2mth.ValEsd(I100,-0.09),6)
     772                    WriteCIFitem("  "+s)
     773
     774            WriteCIFitem('_reflns_number_total', str(refcount))
     775            if hklmin is not None and len(histblk['Reflection Lists']) == 1: # hkl range has no meaning with multiple phases
     776                WriteCIFitem('_reflns_limit_h_min', str(int(hklmin[0])))
     777                WriteCIFitem('_reflns_limit_h_max', str(int(hklmax[0])))
     778                WriteCIFitem('_reflns_limit_k_min', str(int(hklmin[1])))
     779                WriteCIFitem('_reflns_limit_k_max', str(int(hklmax[1])))
     780                WriteCIFitem('_reflns_limit_l_min', str(int(hklmin[2])))
     781                WriteCIFitem('_reflns_limit_l_max', str(int(hklmax[2])))
     782            if hklmin is not None:
     783                WriteCIFitem('_reflns_d_resolution_high', G2mth.ValEsd(dmin,-0.009))
     784                WriteCIFitem('_reflns_d_resolution_low', G2mth.ValEsd(dmax,-0.0009))
     785
    783786            WriteCIFitem('\n# POWDER DATA TABLE')
    784             # is data fixed step?
    785             fixedstep = False
    786             # are ESDs sqrt(I)
    787             countsdata = False
    788             zero = 0.01
    789             if fixedstep:
    790                 WriteCIFitem('_pd_meas_2theta_range_min', text)
    791                 WriteCIFitem('_pd_meas_2theta_range_max', text)
    792                 WriteCIFitem('_pd_meas_2theta_range_inc', text)
    793                 # zero correct
    794                 if zero != 0.0:
    795                     WriteCIFitem('_pd_proc_2theta_range_min', text)
    796                     WriteCIFitem('_pd_proc_2theta_range_max', text)
    797                     WriteCIFitem('_pd_proc_2theta_range_inc', text)
    798             WriteCIFitem('loop_' +
    799                          '\n\t_pd_proc_d_spacing')
    800                          #'_pd_meas_time_of_flight'
     787            # is data fixed step? If the step varies by <0.01% treat as fixed step
     788            steps = histblk['Data'][0][1:] - histblk['Data'][0][:-1]
     789            if abs(max(steps)-min(steps)) > abs(max(steps))/10000.:
     790                fixedstep = False
     791            else:
     792                fixedstep = True
     793
     794            if fixedstep: # and not TOF
     795                WriteCIFitem('_pd_meas_2theta_range_min', G2mth.ValEsd(histblk['Data'][0][0],-0.00009))
     796                WriteCIFitem('_pd_meas_2theta_range_max', G2mth.ValEsd(histblk['Data'][0][-1],-0.00009))
     797                WriteCIFitem('_pd_meas_2theta_range_inc', G2mth.ValEsd(steps.sum()/len(steps),-0.00009))
     798                # zero correct, if defined
     799                zero = None
     800                zerolst = histblk['Instrument Parameters'][0].get('Zero')
     801                if zerolst: zero = zerolst[1]
     802                zero = self.parmDict.get('Zero',zero)
     803                # TODO: Bob is zero added or subtracted?
     804                if zero:
     805                    WriteCIFitem('_pd_proc_2theta_range_min', G2mth.ValEsd(histblk['Data'][0][0]-zero,-0.00009))
     806                    WriteCIFitem('_pd_proc_2theta_range_max', G2mth.ValEsd(histblk['Data'][0][-1]-zero,-0.00009))
     807                    WriteCIFitem('_pd_proc_2theta_range_inc', G2mth.ValEsd(steps.sum()/len(steps),-0.00009))
     808               
     809            if zero:
     810                WriteCIFitem('_pd_proc_number_of_points', str(len(histblk['Data'][0])))
     811            else:
     812                WriteCIFitem('_pd_meas_number_of_points', str(len(histblk['Data'][0])))
     813            WriteCIFitem('\nloop_')
     814            #            WriteCIFitem('\t_pd_proc_d_spacing') # need easy way to get this
    801815            if not fixedstep:
    802                 if zero != 0.0:
     816                if zero:
    803817                    WriteCIFitem('\t_pd_proc_2theta_corrected')
    804818                else:
    805819                    WriteCIFitem('\t_pd_meas_2theta_scan')
    806             if countsdata:
    807                 WriteCIFitem('\t_pd_meas_counts_total')
    808             else:
    809                 WriteCIFitem('\t_pd_meas_intensity_total')
     820            # at least for now, always report weights. TODO: Should they be multiplied by weights?
     821            #if countsdata:
     822            #    WriteCIFitem('\t_pd_meas_counts_total')
     823            #else:
     824            WriteCIFitem('\t_pd_meas_intensity_total')
     825            WriteCIFitem('\t_pd_calc_intensity_total')
     826            WriteCIFitem('\t_pd_proc_intensity_bkg_calc')
    810827            WriteCIFitem('\t_pd_proc_ls_weight')
    811             WriteCIFitem('\t_pd_proc_intensity_bkg_calc')
    812             WriteCIFitem('\t_pd_calc_intensity_total')
    813             if zero != 0.0:
    814                 WriteCIFitem('_pd_proc_number_of_points', text)
    815             else:
    816                 WriteCIFitem('_pd_meas_number_of_points', text)
     828            # TODO: are data excluded?
     829            for x,yobs,yw,ycalc,ybkg in zip(histblk['Data'][0],
     830                                            histblk['Data'][1],
     831                                            histblk['Data'][2],
     832                                            histblk['Data'][3],
     833                                            histblk['Data'][4]):
     834                if fixedstep:
     835                    s = ""
     836                else:
     837                    s = PutInCol(G2mth.ValEsd(x-zero,-0.00009),10)
     838                s += PutInCol(G2mth.ValEsd(yobs,-0.00009),14)
     839                s += PutInCol(G2mth.ValEsd(ycalc,-0.00009),14)
     840                s += PutInCol(G2mth.ValEsd(ybkg,-0.00009),14)
     841                s += PutInCol(G2mth.ValEsd(yw,-0.000009),14)
     842                WriteCIFitem("  "+s)
    817843                         
    818844        def WriteSingleXtalData(histlbl):
    819845            histblk = self.Histograms[histlbl]
    820846            print 'TODO: single xtal here data for',histblk["Instrument Parameters"][0]['InstrName']
    821             # see wrreflist.for
    822             refprx = '_refln.' # mm
     847
     848            #refprx = '_refln.' # mm
    823849            refprx = '_refln_' # normal
     850
     851            print histblk.keys()
    824852           
    825853            WriteCIFitem('\n# STRUCTURE FACTOR TABLE')           
    … …  
    827855                         '\n\t' + refprx + 'index_h' +
    828856                         '\n\t' + refprx + 'index_k' +
    829                          '\n\t' + refprx + 'index_l' +
    830                          '\n\t' + refprx + 'status' +
     857                         '\n\t' + refprx + 'index_l' +
    831858                         '\n\t' + refprx + 'F_squared_meas' +
    832859                         '\n\t' + refprx + 'F_squared_sigma' +
    833860                         '\n\t' + refprx + 'F_squared_calc' +
    834                          '\n\t' + refprx + 'phase_calc')
    835             WriteCIFitem('_reflns_number_total', text)
    836             WriteCIFitem('_reflns_number_observed', text)
    837             WriteCIFitem('_reflns_limit_h_min', text)
    838             WriteCIFitem('_reflns_limit_h_max', text)
    839             WriteCIFitem('_reflns_limit_k_min', text)
    840             WriteCIFitem('_reflns_limit_k_max', text)
    841             WriteCIFitem('_reflns_limit_l_min', text)
    842             WriteCIFitem('_reflns_limit_l_max', text)
    843             WriteCIFitem('_reflns_d_resolution_high', text)
    844             WriteCIFitem('_reflns_d_resolution_low', text)
     861                         '\n\t' + refprx + 'phase_calc' +
     862                         '\n\t_pd_refln_d_spacing'
     863                         )
     864
     865            hklmin = None
     866            hklmax = None
     867            dmax = None
     868            dmin = None
     869            refcount = len(histblk['Data'])
     870            for ref in histblk['Data']:
     871                s = "  "
     872                if hklmin is None:
     873                    hklmin = ref[0:3]
     874                    hklmax = ref[0:3]
     875                    dmax = dmin = ref[4]
     876                for i,hkl in enumerate(ref[0:3]):
     877                    hklmax[i] = max(hkl,hklmax[i])
     878                    hklmin[i] = min(hkl,hklmin[i])
     879                    s += PutInCol(int(hkl),4)
     880                sig = ref[6] * ref[8] / ref[5]
     881                s += PutInCol(G2mth.ValEsd(ref[8],-abs(sig/10)),12)
     882                s += PutInCol(G2mth.ValEsd(sig,-abs(sig)/10.),10)
     883                s += PutInCol(G2mth.ValEsd(ref[9],-abs(sig/10)),12)
     884                s += PutInCol(G2mth.ValEsd(ref[10],-0.9),7)
     885                dmax = max(dmax,ref[4])
     886                dmin = min(dmin,ref[4])
     887                s += PutInCol(G2mth.ValEsd(ref[4],-0.009),8)
     888                WriteCIFitem(s)
     889            WriteCIFitem('_reflns_number_total', str(refcount))
     890            if hklmin is not None:
     891                WriteCIFitem('_reflns_limit_h_min', str(int(hklmin[0])))
     892                WriteCIFitem('_reflns_limit_h_max', str(int(hklmax[0])))
     893                WriteCIFitem('_reflns_limit_k_min', str(int(hklmin[1])))
     894                WriteCIFitem('_reflns_limit_k_max', str(int(hklmax[1])))
     895                WriteCIFitem('_reflns_limit_l_min', str(int(hklmin[2])))
     896                WriteCIFitem('_reflns_limit_l_max', str(int(hklmax[2])))
     897                WriteCIFitem('_reflns_d_resolution_high', G2mth.ValEsd(dmin,-0.009))
     898                WriteCIFitem('_reflns_d_resolution_low', G2mth.ValEsd(dmax,-0.0009))
    845899
    846900        #============================================================
    … …  
    864918        # count phases, powder and single crystal histograms
    865919        self.nphase = len(self.Phases)
    866         self.npowder = 0
    867         self.nsingle = 0
     920        self.powderDict = {}
     921        self.xtalDict = {}
    868922        for hist in self.Histograms:
     923            i = self.Histograms[hist]['hId']
    869924            if hist.startswith("PWDR"):
    870                 self.npowder += 1
     925                self.powderDict[i] = hist
    871926            elif hist.startswith("HKLF"):
    872                 self.nsingle += 1
     927                self.xtalDict[i] = hist
    873928        # is there anything to export?
    874         if self.nphase + self.npowder + self.nsingle == 0:
     929        if self.nphase + len(self.powderDict) + len(self.xtalDict) == 0:
    875930            self.G2frame.ErrorDialog(
    876931                'Empty project',
    … …  
    890945        self.quickmode = False
    891946        phasenam = phasenum = None # include all phases
    892         if mode != "full" or self.npowder + self.nsingle == 0:
     947        if mode != "full" or len(self.powderDict) + len(self.xtalDict) == 0:
    893948            self.quickmode = True
    894949            oneblock = True
    … …  
    906961        elif self.nphase > 1:
    907962            oneblock = False
    908         elif self.npowder + self.nsingle > 1:
     963        elif len(self.powderDict) + len(self.xtalDict) > 1:
    909964            oneblock = False
    910965        else: # one phase, one dataset, Full CIF
    … …  
    10481103                WriteCIFitem(loopprefix,datablockidDict[phasenam])
    10491104            # loop over data blocks
    1050             if self.npowder + self.nsingle > 1:
     1105            if len(self.powderDict) + len(self.xtalDict) > 1:
    10511106                loopprefix = ''
    10521107                WriteCIFitem('loop_   _pd_block_diffractogram_id')
    10531108            else:
    10541109                loopprefix = '_pd_block_diffractogram_id'
    1055             for hist in self.Histograms:
     1110            for i in sorted(self.powderDict.keys()):
     1111                hist = self.powderDict[i]
    10561112                histblk = self.Histograms[hist]
    1057                 if hist.startswith("PWDR"):
    1058                     instnam = histblk["Sample Parameters"]['InstrName']
    1059                 elif hist.startswith("HKLF"):
    1060                     instnam = histblk["Instrument Parameters"][0]['InstrName']
     1113                instnam = histblk["Sample Parameters"]['InstrName']
     1114                instnam = instnam.replace(' ','')
     1115                i = histblk['hId']
     1116                datablockidDict[hist] = (str(self.CIFdate) + "|" + str(self.CIFname) + "|" +
     1117                                         str(self.shortauthorname) + "|" +
     1118                                         instnam + "_hist_"+str(i))
     1119                WriteCIFitem(loopprefix,datablockidDict[hist])
     1120            for i in sorted(self.xtalDict.keys()):
     1121                hist = self.xtalDict[i]
     1122                histblk = self.Histograms[hist]
     1123                instnam = histblk["Instrument Parameters"][0]['InstrName']
    10611124                instnam = instnam.replace(' ','')
    10621125                i = histblk['hId']
    … …  
    10931156            #============================================================
    10941157            # loop over histograms, exporting them
    1095             for hist in self.Histograms:
     1158            for i in sorted(self.powderDict.keys()):
     1159                hist = self.powderDict[i]
    10961160                histblk = self.Histograms[hist]
    1097                 i = histblk['hId']
    10981161                if hist.startswith("PWDR"):
    10991162                    WriteCIFitem('\ndata_'+self.CIFname+"_pwd_"+str(i))
    … …  
    11071170                    WritePowderTemplate()
    11081171                    WritePowderData(hist)
    1109                 elif hist.startswith("HKLF"):
     1172            for i in sorted(self.xtalDict.keys()):
     1173                hist = self.xtalDict[i]
     1174                histblk = self.Histograms[hist]
     1175                if hist.startswith("HKLF"):
    11101176                    WriteCIFitem('\ndata_'+self.CIFname+"_sx_"+str(i))
    11111177                    #instnam = histblk["Instrument Parameters"][0]['InstrName']
Note: See TracChangeset for help on using the changeset viewer.