source: trunk/GSASIIIO.py @ 397

Last change on this file since 397 was 397, checked in by vondreele, 10 years ago

collect default settings for Sample in one routine
add recalibrate routine for images
azimuths from image integration are now the center angle of each azimuth bin
put in 1/2 pixel offset in image calibration/integration calcs

  • Property svn:keywords set to Date Author Revision URL Id
File size: 55.5 KB
Line 
1"""GSASIIIO: functions for IO of data
2   Copyright: 2008, Robert B. Von Dreele (Argonne National Laboratory)
3"""
4########### SVN repository information ###################
5# $Date: 2011-10-21 17:15:50 +0000 (Fri, 21 Oct 2011) $
6# $Author: vondreele $
7# $Revision: 397 $
8# $URL: trunk/GSASIIIO.py $
9# $Id: GSASIIIO.py 397 2011-10-21 17:15:50Z vondreele $
10########### SVN repository information ###################
11import wx
12import math
13import numpy as np
14import cPickle
15import sys
16import random as ran
17import GSASIIpath
18import GSASIIgrid as G2gd
19import GSASIIspc as G2spc
20import GSASIIlattice as G2lat
21import GSASIIpwdGUI as G2pdG
22import GSASIIElem as G2el
23import os.path as ospath
24
25def sfloat(S):
26    if S.strip():
27        return float(S)
28    else:
29        return 0.0
30
31def sint(S):
32    if S.strip():
33        return int(S)
34    else:
35        return 0
36
37def SelectPowderData(self, filename):
38    """Selects banks of data from a filename of any GSAS powder data format
39    Input - filename: any GSAS powder data formatted file (currently STD, FXYE, FXY & ESD)
40    Returns - a list of banks to be read; each entry in list is a tuple containing:
41    filename: same as input filename
42    Pos: position for start of data; record just after BANK record
43    Bank: the BANK record
44    """
45    File = open(filename,'Ur')
46    Title = '''
47First line of this file:
48'''+File.readline()
49    dlg = wx.MessageDialog(self, Title, 'Is this the file you want?', 
50        wx.YES_NO | wx.ICON_QUESTION)
51    try:
52        result = dlg.ShowModal()
53    finally:
54        dlg.Destroy()
55    if result == wx.ID_NO: return (0,0)
56    Temperature = 300
57   
58    if '.xye' in filename:      #Topas style xye file (e.g. 2-th, I, sig) - no iparm file/no BANK record
59        dlg = wx.MessageDialog(self,'''Is this laboratory Cu Ka1/Ka2 data?
60(No = 0.6A wavelength synchrotron data)
61Change wavelength in Instrument Parameters if needed''','Data type?',
62            wx.YES_NO | wx.ICON_QUESTION)
63        try:
64            result = dlg.ShowModal()
65        finally:
66            dlg.Destroy()
67        print result
68        if result == wx.ID_YES:
69            Iparm = {}                                               #Assume CuKa lab data
70            Iparm['INS   HTYPE '] = 'PXC '
71            Iparm['INS  1 ICONS'] = '  1.540500  1.544300       0.0         0       0.7    0       0.5   '
72            Iparm['INS  1PRCF1 '] = '    3    8      0.01                                                '
73            Iparm['INS  1PRCF11'] = '   2.000000E+00  -2.000000E+00   5.000000E+00   0.000000E+00        '
74            Iparm['INS  1PRCF12'] = '   0.000000E+00   0.000000E+00   0.150000E-01   0.150000E-01        '
75        else:
76            Iparm = {}                                               #Assume 0.6A synchrotron data
77            Iparm['INS   HTYPE '] = 'PXC '
78            Iparm['INS  1 ICONS'] = '  0.600000  0.000000       0.0         0      0.99    0       0.5   '
79            Iparm['INS  1PRCF1 '] = '    3    8      0.01                                                '
80            Iparm['INS  1PRCF11'] = '   1.000000E+00  -1.000000E+00   0.300000E+00   0.000000E+00        '
81            Iparm['INS  1PRCF12'] = '   0.000000E+00   0.000000E+00   0.100000E-01   0.100000E-01        '
82                       
83       
84    else:                       #GSAS style fxye or fxy file (e.g. 100*2-th, I, sig)
85        self.IparmName = GetInstrumentFile(self,filename)
86        if self.IparmName:
87            Iparm = GetInstrumentData(self.IparmName)
88        else:
89            Iparm = {}                                               #Assume CuKa lab data if no iparm file
90            Iparm['INS   HTYPE '] = 'PXC '
91            Iparm['INS  1 ICONS'] = '  1.540500  1.544300       0.0         0       0.7    0       0.5   '
92            Iparm['INS  1PRCF1 '] = '    3    8      0.01                                                '
93            Iparm['INS  1PRCF11'] = '   2.000000E+00  -2.000000E+00   5.000000E+00   0.000000E+00        '
94            Iparm['INS  1PRCF12'] = '   0.000000E+00   0.000000E+00   0.150000E-01   0.150000E-01        '
95    S = 1
96    Banks = []
97    Pos = []
98    FoundData = []
99    Comments = []
100    wx.BeginBusyCursor()
101    try:
102        while S:
103            S = File.readline()
104            if S[:1] != '#':
105                if S[:4] == 'BANK':
106                    Banks.append(S)
107                    Pos.append(File.tell())
108                elif '.xye' in filename:    #No BANK in a xye file
109                    Banks.append('BANK 1 XYE')
110                    Pos.append(File.tell())
111                    break
112            else:
113                Comments.append(S[:-1])
114                if 'Temp' in S.split('=')[0]:
115                    Temperature = float(S.split('=')[1])
116        File.close()
117    finally:
118        wx.EndBusyCursor()
119    if Comments:
120       print 'Comments on file:'
121       for Comment in Comments: print Comment
122    if Banks:
123        result = [0]
124        if len(Banks) >= 2:
125            dlg = wx.MultiChoiceDialog(self, 'Which scans do you want?', 'Select scans', Banks, wx.CHOICEDLG_STYLE)
126            try:
127                if dlg.ShowModal() == wx.ID_OK:
128                    result = dlg.GetSelections()
129                else:
130                    result = []
131            finally:
132                dlg.Destroy()
133        for i in result:
134            FoundData.append((filename,Pos[i],Banks[i]))
135    else:
136        dlg = wx.MessageDialog(self, 'ERROR - this is not a GSAS powder data file', 'No BANK records', wx.OK | wx.ICON_ERROR)
137        try:
138            result = dlg.ShowModal()
139        finally:
140            dlg.Destroy()
141    return FoundData,Iparm,Comments,Temperature
142
143def GetInstrumentFile(self,filename):
144    import os.path as op
145    dlg = wx.FileDialog(self,'Choose an instrument file','.', '', 'GSAS iparm file (*.prm)|*.prm|All files(*.*)|*.*', wx.OPEN)
146    if self.dirname: 
147        dlg.SetDirectory(self.dirname)
148        Tname = filename[:filename.index('.')]+'.prm'
149        if op.exists(Tname):
150            self.IparmName = Tname       
151    if self.IparmName: dlg.SetFilename(self.IparmName)
152    filename = ''
153    try:
154        if dlg.ShowModal() == wx.ID_OK:
155            filename = dlg.GetPath()
156    finally:
157        dlg.Destroy()
158    return filename
159
160def GetInstrumentData(IparmName):
161    file = open(IparmName, 'Ur')
162    S = 1
163    Iparm = {}
164    while S:
165        S = file.readline()
166        Iparm[S[:12]] = S[12:-1]
167    return Iparm
168   
169def GetPowderPeaks(fileName):
170    sind = lambda x: math.sin(x*math.pi/180.)
171    asind = lambda x: 180.*math.asin(x)/math.pi
172    Cuka = 1.54052
173    File = open(fileName,'Ur')
174    Comments = []
175    peaks = []
176    S = File.readline()
177    while S:
178        if S[:1] == '#':
179            Comments.append(S[:-1])
180        else:
181            item = S.split()
182            if len(item) == 1:
183                peaks.append([float(item[0]),1.0])
184            elif len(item) > 1:
185                peaks.append([float(item[0]),float(item[0])])
186        S = File.readline()
187    File.close()
188    if Comments:
189       print 'Comments on file:'
190       for Comment in Comments: print Comment
191    Peaks = []
192    if peaks[0][0] > peaks[-1][0]:          # d-spacings - assume CuKa
193        for peak in peaks:
194            dsp = peak[0]
195            sth = Cuka/(2.0*dsp)
196            if sth < 1.0:
197                tth = 2.0*asind(sth)
198            else:
199                break
200            Peaks.append([tth,peak[1],True,False,0,0,0,dsp,0.0])
201    else:                                   #2-thetas - assume Cuka (for now)
202        for peak in peaks:
203            tth = peak[0]
204            dsp = Cuka/(2.0*sind(tth/2.0))
205            Peaks.append([tth,peak[1],True,False,0,0,0,dsp,0.0])
206    return Comments,Peaks
207
208def GetPawleyPeaks(filename):
209    rt2ln2x2 = 2.35482
210    File = open(filename,'Ur')
211    PawleyPeaks = []
212    S = File.readline()         #skip header
213    S = File.readline()
214    item = S.split()
215    while S:
216        h,k,l = int(item[0]),int(item[1]),int(item[2])
217        mult = int(item[3])
218        tth = float(item[5])
219        sig = float(item[6])/rt2ln2x2
220        Iobs = float(item[7])*mult
221        PawleyPeaks.append([h,k,l,mult,tth,False,Iobs,0.0])
222        S = File.readline()
223        item = S.split()
224        if item[3] == '-100.0000':       #find trailer
225            break
226    File.close()
227    return PawleyPeaks
228   
229def GetHKLData(filename):
230    print 'Reading: '+filename
231    File = open(filename,'Ur')
232    HKLref = []
233    HKLmin = [1000,1000,1000]
234    HKLmax = [0,0,0]
235    FoMax = 0
236    ifFc = False
237    S = File.readline()
238    while '#' in S[0]:        #get past comments if any
239        S = File.readline()       
240    if '_' in S:         #cif style .hkl file
241        while 'loop_' not in S:         #skip preliminaries if any - can't have 'loop_' in them!
242            S = File.readline()       
243        S = File.readline()             #get past 'loop_' line
244        pos = 0
245        hpos = kpos = lpos = Fosqpos = Fcsqpos = sigpos = -1
246        while S:
247            if '_' in S:
248                if 'index_h' in S:
249                    hpos = pos
250                elif 'index_k' in S:
251                    kpos = pos
252                elif 'index_l' in S:
253                    lpos = pos
254                elif 'F_squared_meas' in S:
255                    Fosqpos = pos
256                elif 'F_squared_calc' in S:
257                    Fcsqpos = pos
258                elif 'F_squared_sigma' in S:
259                    sigpos = pos
260                pos += 1
261            else:
262                data = S.split()
263                if data:                    #avoid blank lines
264                    HKL = np.array([int(data[hpos]),int(data[kpos]),int(data[lpos])])
265                    h,k,l = HKL
266                    Fosq = float(data[Fosqpos])
267                    if sigpos != -1:
268                        sigFosq = float(data[sigpos])
269                    else:
270                        sigFosq = 1.
271                    if Fcsqpos != -1:
272                        Fcsq = float(data[Fcsqpos])
273                        if Fcsq:
274                            ifFc = True
275                    else:
276                        Fcsq = 0.
277                       
278                    HKLmin = [min(h,HKLmin[0]),min(k,HKLmin[1]),min(l,HKLmin[2])]
279                    HKLmax = [max(h,HKLmax[0]),max(k,HKLmax[1]),max(l,HKLmax[2])]
280                    FoMax = max(FoMax,Fosq)
281                    HKLref.append([HKL,Fosq,sigFosq,Fcsq,0,0,0])                 #room for Fcp, Fcpp & phase
282            S = File.readline()
283    else:                   #dumb h,k,l,Fo,sigFo .hkl file
284        while S:
285            h,k,l,Fo,sigFo = S.split()
286            HKL = np.array([int(h),int(k),int(l)])
287            h,k,l = HKL
288            Fo = float(Fo)
289            sigFo = float(sigFo)
290            HKLmin = [min(h,HKLmin[0]),min(k,HKLmin[1]),min(l,HKLmin[2])]
291            HKLmax = [max(h,HKLmax[0]),max(k,HKLmax[1]),max(l,HKLmax[2])]
292            FoMax = max(FoMax,Fo)
293            HKLref.append([HKL,Fo**2,2.*Fo*sigFo,0,0,0,0])                 #room for Fc, Fcp, Fcpp & phase
294            S = File.readline()
295    File.close()
296    return HKLref,HKLmin,HKLmax,FoMax,ifFc
297
298def GetPowderData(filename,Pos,Bank,DataType):
299    '''Reads one BANK of data from GSAS raw powder data file
300    input:
301    filename: GSAS raw powder file dataname
302    Pos: start of data in file just after BANK record
303    Bank: the BANK record
304    DataType: powder data type, e.g. "PXC" for Powder X-ray CW data
305    returns: list [x,y,e,yc,yb]
306    x: np.array of x-axis values
307    y: np.array of powder pattern intensities
308    w: np.array of w=sig(intensity)^2 values
309    yc: np.array of calc. intensities (zero)
310    yb: np.array of calc. background (zero)
311    yd: np.array of obs-calc profiles
312    '''
313    print 'Reading: '+filename
314    print 'Bank:    '+Bank[:-1]
315    if 'FXYE' in Bank:
316        return GetFXYEdata(filename,Pos,Bank,DataType)
317    elif ' XYE' in Bank:
318        return GetXYEdata(filename,Pos,Bank,DataType)
319    elif 'FXY' in Bank:
320        return GetFXYdata(filename,Pos,Bank,DataType)
321    elif 'ESD' in Bank:
322        return GetESDdata(filename,Pos,Bank,DataType)
323    elif 'STD' in Bank:
324        return GetSTDdata(filename,Pos,Bank,DataType)
325    else:
326        return GetSTDdata(filename,Pos,Bank,DataType)
327    return []
328
329def GetFXYEdata(filename,Pos,Bank,DataType):
330    File = open(filename,'Ur')
331    File.seek(Pos)
332    x = []
333    y = []
334    w = []
335    S = File.readline()
336    while S and S[:4] != 'BANK':
337        vals = S.split()
338        if DataType[2] == 'C':
339            x.append(float(vals[0])/100.)               #CW: from centidegrees to degrees
340        elif DataType[2] == 'T':
341            x.append(float(vals[0])/1000.0)             #TOF: from musec to millisec
342        f = float(vals[1])
343        if f <= 0.0:
344            y.append(0.0)
345            w.append(1.0)
346        else:
347            y.append(float(vals[1]))
348            w.append(1.0/float(vals[2])**2)
349        S = File.readline()
350    File.close()
351    N = len(x)
352    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
353   
354def GetXYEdata(filename,Pos,Bank,DataType):
355    File = open(filename,'Ur')
356    File.seek(Pos)
357    x = []
358    y = []
359    w = []
360    S = File.readline()
361    while S:
362        vals = S.split()
363        try:
364            x.append(float(vals[0]))
365            f = float(vals[1])
366            if f <= 0.0:
367                y.append(0.0)
368                w.append(1.0)
369            else:
370                y.append(float(vals[1]))
371                w.append(1.0/float(vals[2])**2)
372            S = File.readline()
373        except ValueError:
374            break
375    File.close()
376    N = len(x)
377    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
378   
379   
380def GetFXYdata(filename,Pos,Bank,DataType):
381    File = open(filename,'Ur')
382    File.seek(Pos)
383    x = []
384    y = []
385    w = []
386    S = File.readline()
387    while S and S[:4] != 'BANK':
388        vals = S.split()
389        if DataType[2] == 'C':
390            x.append(float(vals[0])/100.)               #CW: from centidegrees to degrees
391        elif DataType[2] == 'T':
392            x.append(float(vals[0])/1000.0)             #TOF: from musec to millisec
393        f = float(vals[1])
394        if f > 0.0:
395            y.append(f)
396            w.append(1.0/f)
397        else:             
398            y.append(0.0)
399            w.append(1.0)
400        S = File.readline()
401    File.close()
402    N = len(x)
403    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
404   
405def GetESDdata(filename,Pos,Bank,DataType):
406    File = open(filename,'Ur')
407    cons = Bank.split()
408    if DataType[2] == 'C':
409        start = float(cons[5])/100.0               #CW: from centidegrees to degrees
410        step = float(cons[6])/100.0
411    elif DataType[2] == 'T':
412        start = float(cons[5])/1000.0              #TOF: from musec to millisec
413        step = float(cons[6])/1000.0
414    File.seek(Pos)
415    x = []
416    y = []
417    w = []
418    S = File.readline()
419    j = 0
420    while S and S[:4] != 'BANK':
421        for i in range(0,80,16):
422            xi = start+step*j
423            yi = sfloat(S[i:i+8])
424            ei = sfloat(S[i+8:i+16])
425            x.append(xi)
426            if yi > 0.0:
427                y.append(yi)
428                w.append(1.0/ei**2)
429            else:             
430                y.append(0.0)
431                w.append(1.0)
432            j += 1
433        S = File.readline()
434    File.close()
435    N = len(x)
436    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
437
438def GetSTDdata(filename,Pos,Bank,DataType):
439    File = open(filename,'Ur')
440    cons = Bank.split()
441    Nch = cons[2]
442    if DataType[2] == 'C':
443        start = float(cons[5])/100.0               #CW: from centidegrees to degrees
444        step = float(cons[6])/100.0
445    elif DataType[2] == 'T':
446        start = float(cons[5])/1000.0              #TOF: from musec to millisec - not likely!
447        step = float(cons[6])/1000.0
448    File.seek(Pos)
449    x = []
450    y = []
451    w = []
452    S = File.readline()
453    j = 0
454    while S and S[:4] != 'BANK':
455        for i in range(0,80,8):
456            xi = start+step*j
457            ni = max(sint(S[i:i+2]),1)
458            yi = max(sfloat(S[i+2:i+8]),0.0)
459            if yi:
460                ei = math.sqrt(yi*ni)
461            else:
462                yi = 0.0
463                ei = 1.0
464            j += 1
465            if j < Nch:
466                x.append(xi)
467                y.append(yi)
468                w.append(1.0/ei**2)
469        S = File.readline()
470    File.close()
471    N = len(x)
472    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
473   
474def CheckImageFile(self,imagefile):
475    if not ospath.exists(imagefile):
476        dlg = wx.FileDialog(self, 'Bad image file name; choose name', '.', '',\
477        'Any image file (*.tif;*.tiff;*.mar*;*.avg;*.sum;*.img)\
478        |*.tif;*.tiff;*.mar*;*.avg;*.sum;*.img|\
479        Any detector tif (*.tif;*.tiff)|*.tif;*.tiff|\
480        MAR file (*.mar*)|*.mar*|\
481        GE Image (*.avg;*.sum)|*.avg;*.sum|\
482        ADSC Image (*.img)|*.img|\
483        All files (*.*)|*.*',wx.OPEN)
484        if self.dirname:
485            dlg.SetDirectory(self.dirname)
486        try:
487            dlg.SetFilename(ospath.split(imagefile)[1])
488            if dlg.ShowModal() == wx.ID_OK:
489                self.dirname = dlg.GetDirectory()
490                imagefile = dlg.GetPath()
491            else:
492                imagefile = False
493        finally:
494            dlg.Destroy()
495    return imagefile
496       
497def GetImageData(self,imagefile,imageOnly=False):       
498    ext = ospath.splitext(imagefile)[1]
499    Comments = []
500    if ext == '.tif' or ext == '.tiff':
501        Comments,Data,Npix,Image = GetTifData(imagefile)
502    elif ext == '.img':
503        Comments,Data,Npix,Image = GetImgData(imagefile)
504        Image[0][0] = 0
505    elif ext == '.mar3450' or ext == '.mar2300':
506        Comments,Data,Npix,Image = GetMAR345Data(imagefile)
507    elif ext in ['.sum','.avg','']:
508        Comments,Data,Npix,Image = GetGEsumData(imagefile)
509    elif ext == '.G2img':
510        Comments,Data,Npix,Image = GetG2Image(imagefile)
511    if imageOnly:
512        return Image
513    else:
514        return Comments,Data,Npix,Image
515       
516def PutG2Image(filename,Comments,Data,Npix,image):
517    File = open(filename,'wb')
518    cPickle.dump([Comments,Data,Npix,image],File,1)
519    File.close()
520    return
521   
522def GetG2Image(filename):
523    File = open(filename,'rb')
524    Comments,Data,Npix,image = cPickle.load(File)
525    File.close()
526    return Comments,Data,Npix,image
527   
528def GetGEsumData(filename,imageOnly=False):
529    import struct as st
530    import array as ar
531    if not imageOnly:
532        print 'Read GE sum file: ',filename   
533    File = open(filename,'rb')
534    if '.sum' in filename:
535        head = ['GE detector sum data from APS 1-ID',]
536        sizexy = [2048,2048]
537    elif '.avg' in filename:
538        head = ['GE detector avg data from APS 1-ID',]
539        sizexy = [2048,2048]
540    else:
541        head = ['GE detector raw data from APS 1-ID',]
542        File.seek(18)
543        size,nframes = st.unpack('<ih',File.read(6))
544        sizexy = [2048,2048]
545        pos = 8192
546        File.seek(pos)
547    Npix = sizexy[0]*sizexy[1]
548    if '.sum' in filename:
549        image = np.array(ar.array('f',File.read(4*Npix)),dtype=np.int32)
550    elif '.avg' in filename:
551        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
552    else:
553        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
554        while nframes > 1:
555            image += np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
556            nframes -= 1
557    image = np.reshape(image,(sizexy[1],sizexy[0]))
558    data = {'pixelSize':(200,200),'wavelength':0.15,'distance':250.0,'center':[204.8,204.8],'size':sizexy} 
559    File.close()   
560    if imageOnly:
561        return image
562    else:
563        return head,data,Npix,image
564       
565def GetImgData(filename,imageOnly=False):
566    import struct as st
567    import array as ar
568    if not imageOnly:
569        print 'Read ADSC img file: ',filename
570    File = open(filename,'rb')
571    head = File.read(511)
572    lines = head.split('\n')
573    head = []
574    center = [0,0]
575    for line in lines[1:-2]:
576        line = line.strip()[:-1]
577        if line:
578            if 'SIZE1' in line:
579                size = int(line.split('=')[1])
580                Npix = size*size
581            elif 'WAVELENGTH' in line:
582                wave = float(line.split('=')[1])
583            elif 'BIN' in line:
584                if line.split('=')[1] == '2x2':
585                    pixel=(102,102)
586                else:
587                    pixel = (51,51)
588            elif 'DISTANCE' in line:
589                distance = float(line.split('=')[1])
590            elif 'CENTER_X' in line:
591                center[0] = float(line.split('=')[1])
592            elif 'CENTER_Y' in line:
593                center[1] = float(line.split('=')[1])
594            head.append(line)
595    data = {'pixelSize':pixel,'wavelength':wave,'distance':distance,'center':center,'size':[size,size]}
596    image = []
597    row = 0
598    pos = 512
599    File.seek(pos)
600    image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
601    image = np.reshape(image,(sizexy[1],sizexy[0]))
602#    image = np.zeros(shape=(size,size),dtype=np.int32)   
603#    while row < size:
604#        File.seek(pos)
605#        line = ar.array('H',File.read(2*size))
606#        image[row] = np.asarray(line)
607#        row += 1
608#        pos += 2*size
609    File.close()
610    if imageOnly:
611        return image
612    else:
613        return lines[1:-2],data,Npix,image
614       
615def GetMAR345Data(filename,imageOnly=False):
616    import array as ar
617    import struct as st
618    try:
619        import pack_f as pf
620    except:
621        msg = wx.MessageDialog(None, message="Unable to load the GSAS MAR image decompression, pack_f",
622                               caption="Import Error",
623                               style=wx.ICON_ERROR | wx.OK | wx.STAY_ON_TOP)
624        msg.ShowModal()
625        return None,None,None,None
626
627    if not imageOnly:
628        print 'Read Mar345 file: ',filename
629    File = open(filename,'rb')
630    head = File.read(4095)
631    numbers = st.unpack('<iiiiiiiiii',head[:40])
632    lines = head[128:].split('\n')
633    head = []
634    for line in lines:
635        line = line.strip()
636        if 'PIXEL' in line:
637            values = line.split()
638            pixel = (int(values[2]),int(values[4]))     #in microns
639        elif 'WAVELENGTH' in line:
640            wave = float(line.split()[1])
641        elif 'DISTANCE' in line:
642            distance = float(line.split()[1])           #in mm
643        elif 'CENTER' in line:
644            values = line.split()
645            center = [float(values[2])/10.,float(values[4])/10.]    #make in mm from pixels
646        if line: 
647            head.append(line)
648    data = {'pixelSize':pixel,'wavelength':wave,'distance':distance,'center':center}
649    for line in head:
650        if 'FORMAT' in line[0:6]:
651            items = line.split()
652            size = int(items[1])
653            Npix = size*size
654    pos = 4096
655    data['size'] = [size,size]
656    File.seek(pos)
657    line = File.read(8)
658    while 'CCP4' not in line:       #get past overflow list for now
659        line = File.read(8)
660        pos += 8
661    pos += 37
662    File.seek(pos)
663    raw = File.read()
664    File.close()
665    image = np.zeros(shape=(size,size),dtype=np.int32)
666    image = pf.pack_f(len(raw),raw,size,image)
667    if imageOnly:
668        return image.T              #transpose to get it right way around
669    else:
670        return head,data,Npix,image.T
671       
672def GetTifData(filename,imageOnly=False):
673    import struct as st
674    import array as ar
675    File = open(filename,'rb')
676    dataType = 5
677    try:
678        Meta = open(filename+'.metadata','Ur')
679        head = Meta.readlines()
680        for line in head:
681            line = line.strip()
682            if 'dataType' in line:
683                dataType = int(line.split('=')[1])
684        Meta.close()
685    except IOError:
686        print 'no metadata file found - will try to read file anyway'
687        head = ['no metadata file found',]
688       
689    tag = File.read(2)
690    byteOrd = '<'
691    if tag == 'II' and int(st.unpack('<h',File.read(2))[0]) == 42:     #little endian
692        IFD = int(st.unpack(byteOrd+'i',File.read(4))[0])
693    elif tag == 'MM' and int(st.unpack('>h',File.read(2))[0]) == 42:   #big endian
694        byteOrd = '>'
695        IFD = int(st.unpack(byteOrd+'i',File.read(4))[0])       
696    else:
697        lines = ['not a detector tiff file',]
698        return lines,0,0,0
699    File.seek(IFD)                                                  #get number of directory entries
700    NED = int(st.unpack(byteOrd+'h',File.read(2))[0])
701    print IFD,NED
702    IFD = {}
703    for ied in range(NED):
704        Tag,Type = st.unpack(byteOrd+'Hh',File.read(4))
705        nVal = st.unpack(byteOrd+'i',File.read(4))[0]
706        if Type == 1:
707            Value = st.unpack(byteOrd+nVal*'b',File.read(nVal))
708        elif Type == 2:
709            Value = st.unpack(byteOrd+'i',File.read(4))
710        elif Type == 3:
711            Value = st.unpack(byteOrd+nVal*'h',File.read(nVal*2))
712            x = st.unpack(byteOrd+nVal*'h',File.read(nVal*2))
713        elif Type == 4:
714            Value = st.unpack(byteOrd+nVal*'i',File.read(nVal*4))
715        elif Type == 5:
716            Value = st.unpack(byteOrd+nVal*'i',File.read(nVal*4))
717        elif Type == 11:
718            Value = st.unpack(byteOrd+nVal*'f',File.read(nVal*4))
719        IFD[Tag] = [Type,nVal,Value]
720#    for key in IFD:
721#        print key,IFD[key]
722    sizexy = [IFD[256][2][0],IFD[257][2][0]]
723    [nx,ny] = sizexy
724    Npix = nx*ny
725    if 272 in IFD:
726        ifd = IFD[272]
727        File.seek(ifd[2][0])
728        S = File.read(ifd[1])
729        if 'PILATUS' in S:
730            tifType = 'Pilatus'
731            dataType = 0
732            pixy = (172,172)
733            File.seek(4096)
734            if not imageOnly:
735                print 'Read Pilatus tiff file: ',filename
736            image = ar.array('L',File.read(4*Npix))
737            image = np.array(np.asarray(image),dtype=np.int32)
738    elif 262 in IFD and IFD[262][2][0] > 4:
739        tifType = 'DND'
740        pixy = (158,158)
741        File.seek(512)
742        if not imageOnly:
743            print 'Read DND SAX/WAX-detector tiff file: ',filename
744        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
745    elif sizexy == [1536,1536]:
746        tifType = 'APS Gold'
747        pixy = (150,150)
748        File.seek(64)
749        if not imageOnly:
750            print 'Read Gold tiff file:',filename
751        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
752    elif sizexy == [2048,2048] or sizexy == [1024,1024]:
753        if IFD[273][2][0] == 8:
754            if IFD[258][2][0] == 32:
755                tifType = 'PE'
756                pixy = (200,200)
757                File.seek(8)
758                if not imageOnly:
759                    print 'Read APS PE-detector tiff file: ',filename
760                if dataType == 5:
761                    image = np.array(ar.array('f',File.read(4*Npix)),dtype=np.float32)
762                else:
763                    image = np.array(ar.array('I',File.read(4*Npix)),dtype=np.int32)
764        elif IFD[273][2][0] == 4096:
765            tifType = 'MAR'
766            pixy = (158,158)
767            File.seek(4096)
768            if not imageOnly:
769                print 'Read MAR CCD tiff file: ',filename
770            image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
771#    elif sizexy == [960,960]:
772#        tiftype = 'PE-BE'
773#        pixy = (200,200)
774#        File.seek(8)
775#        if not imageOnly:
776#            print 'Read Gold tiff file:',filename
777#        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
778           
779    else:
780        lines = ['not a known detector tiff file',]
781        return lines,0,0,0
782       
783    image = np.reshape(image,(sizexy[1],sizexy[0]))
784    center = [pixy[0]*sizexy[0]/2000,pixy[1]*sizexy[1]/2000]
785    data = {'pixelSize':pixy,'wavelength':0.10,'distance':100.0,'center':center,'size':sizexy}
786    File.close()   
787    if imageOnly:
788        return image
789    else:
790        return head,data,Npix,image
791   
792def ProjFileOpen(self):
793    file = open(self.GSASprojectfile,'rb')
794    print 'load from file: ',self.GSASprojectfile
795    wx.BeginBusyCursor()
796    try:
797        while True:
798            try:
799                data = cPickle.load(file)
800            except EOFError:
801                break
802            datum = data[0]
803            print 'load: ',datum[0]
804           
805            Id = self.PatternTree.AppendItem(parent=self.root,text=datum[0])
806            if 'PWDR' in datum[0]:               
807                self.PatternTree.SetItemPyData(Id,datum[1][:3])     #temp. trim off junk
808            else:
809                self.PatternTree.SetItemPyData(Id,datum[1])
810            for datus in data[1:]:
811                print '    load: ',datus[0]
812                sub = self.PatternTree.AppendItem(Id,datus[0])
813                self.PatternTree.SetItemPyData(sub,datus[1])
814            if 'IMG' in datum[0]:                   #retreive image default flag & data if set
815                Data = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Image Controls'))
816                if Data['setDefault']:
817                    self.imageDefault = Data               
818        file.close()
819       
820    finally:
821        wx.EndBusyCursor()
822    print 'project load successful'
823    self.NewPlot = True
824   
825def ProjFileSave(self):
826    if not self.PatternTree.IsEmpty():
827        file = open(self.GSASprojectfile,'wb')
828        print 'save to file: ',self.GSASprojectfile
829        wx.BeginBusyCursor()
830        try:
831            item, cookie = self.PatternTree.GetFirstChild(self.root)
832            while item:
833                data = []
834                name = self.PatternTree.GetItemText(item)
835                print 'save: ',name
836                data.append([name,self.PatternTree.GetItemPyData(item)])
837                item2, cookie2 = self.PatternTree.GetFirstChild(item)
838                while item2:
839                    name = self.PatternTree.GetItemText(item2)
840                    print '    save: ',name
841                    data.append([name,self.PatternTree.GetItemPyData(item2)])
842                    item2, cookie2 = self.PatternTree.GetNextChild(item, cookie2)                           
843                item, cookie = self.PatternTree.GetNextChild(self.root, cookie)                           
844                cPickle.dump(data,file,1)
845            file.close()
846        finally:
847            wx.EndBusyCursor()
848        print 'project save successful'
849       
850def SaveIntegration(self,PickId,data):
851    azms = self.Integrate[1]
852    X = self.Integrate[2][:-1]
853    Xminmax = [X[0],X[-1]]
854    N = len(X)
855    Id = self.PatternTree.GetItemParent(PickId)
856    name = self.PatternTree.GetItemText(Id)
857    name = name.replace('IMG ','PWDR ')
858    Comments = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Comments'))
859    names = ['Type','Lam','Zero','Polariz.','U','V','W','X','Y','SH/L','Azimuth'] 
860    codes = [0 for i in range(11)]
861    LRazm = data['LRazimuth']
862    Azms = []
863    if data['fullIntegrate'] and data['outAzimuths'] == 1:
864        Azms = [45.0,]                              #a poor man's average?
865    else:
866        for i,azm in enumerate(azms[:-1]):
867            Azms.append((azms[i+1]+azm)/2.)
868    for i,azm in enumerate(azms[:-1]):
869        item, cookie = self.PatternTree.GetFirstChild(self.root)
870        Id = 0
871        while item:
872            Name = self.PatternTree.GetItemText(item)
873            if name == Name:
874                Id = item
875            item, cookie = self.PatternTree.GetNextChild(self.root, cookie)
876        parms = ['PXC',data['wavelength'],0.0,0.99,1.0,-0.10,0.4,0.30,1.0,0.0001,Azms[i]]    #set polarization for synchrotron radiation!
877        Y = self.Integrate[0][i]
878        W = 1./Y                    #probably not true
879        Sample = G2pdG.SetDefaultSample()
880        if Id:
881            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Comments'),Comments)                   
882            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Limits'),[tuple(Xminmax),Xminmax])
883            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Background'),[['chebyschev',1,3,1.0,0.0,0.0]])
884            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Instrument Parameters'),[tuple(parms),parms[:],codes,names])
885            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Peak List'),[])
886            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Index Peak List'),[])
887            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Unit Cells List'),[])             
888            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Reflection Lists'),{})             
889        else:
890            Id = self.PatternTree.AppendItem(parent=self.root,text=name+" Azm= %.2f"%(Azms[i]))
891            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Comments'),Comments)                   
892            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Limits'),[tuple(Xminmax),Xminmax])
893            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Background'),[['chebyschev',1,3,1.0,0.0,0.0]])
894            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Instrument Parameters'),[tuple(parms),parms[:],codes,names])
895            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Sample Parameters'),Sample)
896            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Peak List'),[])
897            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Index Peak List'),[])
898            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Unit Cells List'),[])
899            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Reflection Lists'),{})             
900        self.PatternTree.SetItemPyData(Id,[[''],[np.array(X),np.array(Y),np.array(W),np.zeros(N),np.zeros(N),np.zeros(N)]])
901    self.PatternTree.SelectItem(Id)
902    self.PatternTree.Expand(Id)
903    self.PatternId = Id
904           
905def powderFxyeSave(self,exports,powderfile):
906    head,tail = ospath.split(powderfile)
907    name,ext = tail.split('.')
908    wx.BeginBusyCursor()
909    for i,export in enumerate(exports):
910        filename = ospath.join(head,name+'-%03d.'%(i)+ext)
911        prmname = filename.strip(ext)+'prm'
912        prm = open(prmname,'w')      #old style GSAS parm file
913        PickId = G2gd.GetPatternTreeItemId(self, self.root, export)
914        Values,Names = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self, \
915            PickId, 'Instrument Parameters'))[1::2]     #get values & names
916        Inst = dict(zip(Names,Values))
917        print Inst['Type']
918        prm.write( '            123456789012345678901234567890123456789012345678901234567890        '+'\n')
919        prm.write( 'INS   BANK      1                                                               '+'\n')
920        prm.write(('INS   HTYPE   %sR                                                              '+'\n')%(Inst['Type']))
921        if 'Lam1' in Inst:              #Ka1 & Ka2
922            prm.write(('INS  1 ICONS%10.7f%10.7f    0.0000               0.990    0     0.500   '+'\n')%(Inst['Lam1'],Inst['Lam2']))
923        elif 'Lam' in Inst:             #single wavelength
924            prm.write(('INS  1 ICONS%10.7f%10.7f    0.0000               0.990    0     0.500   '+'\n')%(Inst['Lam'],0.0))
925        prm.write( 'INS  1 IRAD     0                                                               '+'\n')
926        prm.write( 'INS  1I HEAD                                                                    '+'\n')
927        prm.write( 'INS  1I ITYP    0    0.0000  180.0000         1                                 '+'\n')
928        prm.write(('INS  1DETAZM%10.3f                                                          '+'\n')%(Inst['Azimuth']))
929        prm.write( 'INS  1PRCF1     3    8   0.00100                                                '+'\n')
930        prm.write(('INS  1PRCF11     %15.6g%15.6g%15.6g%15.6g   '+'\n')%(Inst['U'],Inst['V'],Inst['W'],0.0))
931        prm.write(('INS  1PRCF12     %15.6g%15.6g%15.6g%15.6g   '+'\n')%(Inst['X'],Inst['Y'],Inst['SH/L']/2.,Inst['SH/L']/2.))
932        prm.close()
933        file = open(filename,'w')
934        print 'save powder pattern to file: ',filename
935        try:
936            x,y,w,yc,yb,yd = self.PatternTree.GetItemPyData(PickId)[1]
937            file.write(powderfile+'\n')
938            file.write('BANK 1 %d %d CONS %.2f %.2f 0 0 FXYE\n'%(len(x),len(x),\
939                100.*x[0],100.*(x[1]-x[0])))
940            s = list(np.sqrt(1./np.array(w)))       
941            XYW = zip(x,y,s)
942            for X,Y,S in XYW:
943                file.write("%15.6g %15.6g %15.6g\n" % (100.*X,Y,max(S,1.0)))
944            file.close()
945        finally:
946            wx.EndBusyCursor()
947        print 'powder pattern file written'
948       
949def powderXyeSave(self,exports,powderfile):
950    head,tail = ospath.split(powderfile)
951    name,ext = tail.split('.')
952    for i,export in enumerate(exports):
953        filename = ospath.join(head,name+'-%03d.'%(i)+ext)
954        PickId = G2gd.GetPatternTreeItemId(self, self.root, export)
955        file = open(filename,'w')
956        file.write('#%s\n'%(export))
957        print 'save powder pattern to file: ',filename
958        wx.BeginBusyCursor()
959        try:
960            x,y,w,yc,yb,yd = self.PatternTree.GetItemPyData(PickId)[1]
961            s = list(np.sqrt(1./np.array(w)))       
962            XYW = zip(x,y,s)
963            for X,Y,W in XYW:
964                file.write("%15.6g %15.6g %15.6g\n" % (X,Y,W))
965            file.close()
966        finally:
967            wx.EndBusyCursor()
968        print 'powder pattern file written'
969       
970def PDFSave(self,exports):   
971    for export in exports:
972        PickId = G2gd.GetPatternTreeItemId(self, self.root, export)
973        SQname = 'S(Q)'+export[4:]
974        GRname = 'G(R)'+export[4:]
975        sqfilename = ospath.join(self.dirname,export.replace(' ','_')[5:]+'.sq')
976        grfilename = ospath.join(self.dirname,export.replace(' ','_')[5:]+'.gr')
977        sqId = G2gd.GetPatternTreeItemId(self, PickId, SQname)
978        grId = G2gd.GetPatternTreeItemId(self, PickId, GRname)
979        sqdata = np.array(self.PatternTree.GetItemPyData(sqId)[1][:2]).T
980        grdata = np.array(self.PatternTree.GetItemPyData(grId)[1][:2]).T
981        sqfile = open(sqfilename,'w')
982        grfile = open(grfilename,'w')
983        sqfile.write('#T S(Q) %s\n'%(export))
984        grfile.write('#T G(R) %s\n'%(export))
985        sqfile.write('#L Q     S(Q)\n')
986        grfile.write('#L R     G(R)\n')
987        for q,sq in sqdata:
988            sqfile.write("%15.6g %15.6g\n" % (q,sq))
989        sqfile.close()
990        for r,gr in grdata:
991            grfile.write("%15.6g %15.6g\n" % (r,gr))
992        grfile.close()
993   
994def PeakListSave(self,file,peaks):
995    print 'save peak list to file: ',self.peaklistfile
996    if not peaks:
997        dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
998        try:
999            result = dlg.ShowModal()
1000        finally:
1001            dlg.Destroy()
1002        return
1003    for peak in peaks:
1004        file.write("%10.4f %12.2f %10.3f %10.3f \n" % \
1005            (peak[0],peak[2],peak[4],peak[6]))
1006    print 'peak list saved'
1007             
1008def IndexPeakListSave(self,peaks):
1009    file = open(self.peaklistfile,'wa')
1010    print 'save index peak list to file: ',self.peaklistfile
1011    wx.BeginBusyCursor()
1012    try:
1013        if not peaks:
1014            dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
1015            try:
1016                result = dlg.ShowModal()
1017            finally:
1018                dlg.Destroy()
1019            return
1020        for peak in peaks:
1021            file.write("%12.6f\n" % (peak[7]))
1022        file.close()
1023    finally:
1024        wx.EndBusyCursor()
1025    print 'index peak list saved'
1026   
1027def ReadEXPPhase(self,filename):
1028    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1029    textureData = {'Order':0,'Model':'cylindrical','Sample omega':[False,0.0],
1030        'Sample chi':[False,0.0],'Sample phi':[False,0.0],'SH Coeff':[False,{}],
1031        'SHShow':False,'PFhkl':[0,0,1],'PFxyz':[0,0,1],'PlotType':'Pole figure'}
1032    shNcof = 0
1033    file = open(filename, 'Ur')
1034    S = 1
1035    Expr = [{},{},{},{},{},{},{},{},{}]
1036    while S:
1037        S = file.readline()
1038        if 'EXPR NPHAS' in S[:12]:
1039            Num = S[12:-1].count('0')
1040            NPhas = S[12:-1].split()
1041        if 'CRS' in S[:3]:
1042            N = int(S[3:4])-1
1043            Expr[N][S[:12]] = S[12:-1]
1044    file.close()
1045    PNames = []
1046    for n,N in enumerate(NPhas):
1047        if N != '0':
1048            result = n
1049            key = 'CRS'+str(n+1)+'    PNAM'
1050            PNames.append(Expr[n][key])
1051    if Num < 8:
1052        dlg = wx.SingleChoiceDialog(self, 'Which phase to read?', 'Read phase data', PNames, wx.CHOICEDLG_STYLE)
1053        try:
1054            if dlg.ShowModal() == wx.ID_OK:
1055                result = dlg.GetSelection()
1056        finally:
1057            dlg.Destroy()       
1058    EXPphase = Expr[result]
1059    keyList = EXPphase.keys()
1060    keyList.sort()
1061    SGData = {}
1062    if NPhas[result] == '1':
1063        Ptype = 'nuclear'
1064    elif NPhas[result] in ['2','3']:
1065        Ptype = 'magnetic'
1066    elif NPhas[result] == '4':
1067        Ptype = 'macromolecular'
1068    elif NPhas[result] == '10':
1069        Ptype = 'Pawley'
1070    for key in keyList:
1071        if 'PNAM' in key:
1072           PhaseName = EXPphase[key].strip()
1073        elif 'ABC   ' in key:
1074            abc = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                       
1075        elif 'ANGLES' in key:
1076            angles = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                                               
1077        elif 'SG SYM' in key:
1078            SpGrp = EXPphase[key][:15].strip()
1079            E,SGData = G2spc.SpcGroup(SpGrp)
1080        elif 'OD    ' in key:
1081            SHdata = EXPphase[key].split()
1082            textureData['Order'] = int(SHdata[0])
1083            textureData['Model'] = shModels[int(SHdata[2])]
1084            textureData['Sample omega'] = [False,float(SHdata[6])]
1085            textureData['Sample chi'] = [False,float(SHdata[7])]
1086            textureData['Sample phi'] = [False,float(SHdata[8])]
1087            shNcof = int(SHdata[1])
1088    Atoms = []
1089    if Ptype == 'nuclear':
1090        for key in keyList:
1091            if 'AT' in key:
1092                if key[11:] == 'A':
1093                    S = EXPphase[key]
1094                elif key[11:] == 'B':
1095                    S += EXPphase[key]
1096                    Atom = [S[50:58].strip(),S[:10].strip(),'',
1097                        float(S[10:20]),float(S[20:30]),float(S[30:40]),
1098                        float(S[40:50]),'',int(S[60:62]),S[130:131]]
1099                    if Atom[9] == 'I':
1100                        Atom += [float(S[68:78]),0.,0.,0.,0.,0.,0.]
1101                    elif Atom[9] == 'A':
1102                        Atom += [0.0,float(S[68:78]),float(S[78:88]),
1103                            float(S[88:98]),float(S[98:108]),
1104                            float(S[108:118]),float(S[118:128])]
1105                    XYZ = Atom[3:6]
1106                    Atom[7],Atom[8] = G2spc.SytSym(XYZ,SGData)
1107                    Atom.append(ran.randint(0,sys.maxint))
1108                    Atoms.append(Atom)
1109    elif Ptype == 'macromolecular':
1110        for key in keyList:
1111            if 'AT' in key[6:8]:
1112                S = EXPphase[key]
1113                Atom = [int(S[56:60]),S[50:54].strip().upper(),S[54:56],
1114                    S[46:51].strip(),S[:8].strip(),'',
1115                    float(S[16:24]),float(S[24:32]),float(S[32:40]),
1116                    float(S[8:16]),'1',1,'I',float(S[40:46]),0,0,0,0,0,0]
1117                XYZ = Atom[6:9]
1118                Atom[10],Atom[11] = G2spc.SytSym(XYZ,SGData)
1119                Atom.append(ran.randint(0,sys.maxint))
1120                Atoms.append(Atom)
1121    Volume = G2lat.calc_V(G2lat.cell2A(abc+angles))
1122    if shNcof:
1123        shCoef = {}
1124        nRec = [i+1 for i in range((shNcof-1)/6+1)]
1125        for irec in nRec:
1126            ODkey = keyList[0][:6]+'OD'+'%3dA'%(irec)
1127            indx = EXPphase[ODkey].split()
1128            ODkey = ODkey[:-1]+'B'
1129            vals = EXPphase[ODkey].split()
1130            for i,val in enumerate(vals):
1131                key = 'C(%s,%s,%s)'%(indx[3*i],indx[3*i+1],indx[3*i+2])
1132                shCoef[key] = float(val)
1133        textureData['SH Coeff'] = [False,shCoef]
1134       
1135    Phase = {
1136            'General':{
1137                'Name':PhaseName,
1138                'Type':Ptype,
1139                'SGData':SGData,
1140                'Cell':[False,]+abc+angles+[Volume,],
1141                'Pawley dmin':1.0},
1142            'Atoms':Atoms,
1143            'Drawing':{},
1144            'Histograms':{},
1145            'Pawley ref':[],
1146            'Models':{},
1147            'SH Texture':textureData
1148            }
1149           
1150    return Phase
1151       
1152def ReadPDBPhase(filename):
1153    EightPiSq = 8.*math.pi**2
1154    file = open(filename, 'Ur')
1155    Phase = {}
1156    Title = ''
1157    Compnd = ''
1158    Atoms = []
1159    A = np.zeros(shape=(3,3))
1160    S = file.readline()
1161    while S:
1162        Atom = []
1163        if 'TITLE' in S[:5]:
1164            Title = S[10:72].strip()
1165            S = file.readline()
1166        elif 'COMPND    ' in S[:10]:
1167            Compnd = S[10:72].strip()
1168            S = file.readline()
1169        elif 'CRYST' in S[:5]:
1170            abc = S[7:34].split()
1171            angles = S[34:55].split()
1172            cell=[float(abc[0]),float(abc[1]),float(abc[2]),
1173                float(angles[0]),float(angles[1]),float(angles[2])]
1174            Volume = G2lat.calc_V(G2lat.cell2A(cell))
1175            AA,AB = G2lat.cell2AB(cell)
1176            SpGrp = S[55:65]
1177            E,SGData = G2spc.SpcGroup(SpGrp)
1178            if E: 
1179                print ' ERROR in space group symbol ',SpGrp,' in file ',filename
1180                print ' N.B.: make sure spaces separate axial fields in symbol' 
1181                print G2spc.SGErrors(E)
1182                return None
1183            SGlines = G2spc.SGPrint(SGData)
1184            for line in SGlines: print line
1185            S = file.readline()
1186        elif 'SCALE' in S[:5]:
1187            V = (S[10:41].split())
1188            A[int(S[5])-1] = [float(V[0]),float(V[1]),float(V[2])]
1189            S = file.readline()
1190        elif 'ATOM' in S[:4] or 'HETATM' in S[:6]:
1191            XYZ = [float(S[31:39]),float(S[39:47]),float(S[47:55])]
1192            XYZ = np.inner(AB,XYZ)
1193            SytSym,Mult = G2spc.SytSym(XYZ,SGData)
1194            Uiso = float(S[61:67])/EightPiSq
1195            Type = S[12:14].upper()
1196            if Type[0] in '123456789':
1197                Type = Type[1:]
1198            Atom = [S[22:27].strip(),S[17:20].upper(),S[20:22],
1199                S[12:17].strip(),Type.strip(),'',XYZ[0],XYZ[1],XYZ[2],
1200                float(S[55:61]),SytSym,Mult,'I',Uiso,0,0,0,0,0,0]
1201            S = file.readline()
1202            if 'ANISOU' in S[:6]:
1203                Uij = S[30:72].split()
1204                Uij = [float(Uij[0])/10000.,float(Uij[1])/10000.,float(Uij[2])/10000.,
1205                    float(Uij[3])/10000.,float(Uij[4])/10000.,float(Uij[5])/10000.]
1206                Atom = Atom[:14]+Uij
1207                Atom[12] = 'A'
1208                S = file.readline()
1209            Atom.append(ran.randint(0,sys.maxint))
1210            Atoms.append(Atom)
1211        else:           
1212            S = file.readline()
1213    file.close()
1214    if Title:
1215        PhaseName = Title
1216    elif Compnd:
1217        PhaseName = Compnd
1218    else:
1219        PhaseName = 'None'
1220    Phase['General'] = {'Name':PhaseName,'Type':'macromolecular','SGData':SGData,
1221        'Cell':[False,]+cell+[Volume,]}
1222    Phase['Atoms'] = Atoms
1223    Phase['Drawing'] = {}
1224    Phase['Histograms'] = {}
1225   
1226    return Phase
1227   
1228def ReadCIFPhase(filename):
1229    anisoNames = ['aniso_u_11','aniso_u_22','aniso_u_33','aniso_u_12','aniso_u_13','aniso_u_23']
1230    file = open(filename, 'Ur')
1231    Phase = {}
1232    Title = ospath.split(filename)[-1]
1233    print '\n Reading cif file: ',Title
1234    Compnd = ''
1235    Atoms = []
1236    A = np.zeros(shape=(3,3))
1237    S = file.readline()
1238    while S:
1239        if '_symmetry_space_group_name_H-M' in S:
1240            SpGrp = S.split("_symmetry_space_group_name_H-M")[1].strip().strip('"').strip("'")
1241            E,SGData = G2spc.SpcGroup(SpGrp)
1242            if E:
1243                print ' ERROR in space group symbol ',SpGrp,' in file ',filename
1244                print ' N.B.: make sure spaces separate axial fields in symbol' 
1245                print G2spc.SGErrors(E)
1246                return None
1247            S = file.readline()
1248        elif '_cell' in S:
1249            if '_cell_length_a' in S:
1250                a = S.split('_cell_length_a')[1].strip().strip('"').strip("'").split('(')[0]
1251            elif '_cell_length_b' in S:
1252                b = S.split('_cell_length_b')[1].strip().strip('"').strip("'").split('(')[0]
1253            elif '_cell_length_c' in S:
1254                c = S.split('_cell_length_c')[1].strip().strip('"').strip("'").split('(')[0]
1255            elif '_cell_angle_alpha' in S:
1256                alp = S.split('_cell_angle_alpha')[1].strip().strip('"').strip("'").split('(')[0]
1257            elif '_cell_angle_beta' in S:
1258                bet = S.split('_cell_angle_beta')[1].strip().strip('"').strip("'").split('(')[0]
1259            elif '_cell_angle_gamma' in S:
1260                gam = S.split('_cell_angle_gamma')[1].strip().strip('"').strip("'").split('(')[0]
1261            S = file.readline()
1262        elif 'loop_' in S:
1263            labels = {}
1264            i = 0
1265            while S:
1266                S = file.readline()
1267                if '_atom_site' in S.strip()[:10]:
1268                    labels[S.strip().split('_atom_site_')[1].lower()] = i
1269                    i += 1
1270                else:
1271                    break
1272            if labels:
1273                if 'aniso_label' not in labels:
1274                    while S:
1275                        atom = ['','','',0,0,0,1.0,'','','I',0.01,0,0,0,0,0,0]
1276                        S.strip()
1277                        if len(S.split()) != len(labels):
1278                            if 'loop_' in S:
1279                                break
1280                            S += file.readline().strip()
1281                        data = S.split()
1282                        if len(data) != len(labels):
1283                            break
1284                        for key in labels:
1285                            if key == 'type_symbol':
1286                                atom[1] = data[labels[key]]
1287                            elif key == 'label':
1288                                atom[0] = data[labels[key]]
1289                            elif key == 'fract_x':
1290                                atom[3] = float(data[labels[key]].split('(')[0])
1291                            elif key == 'fract_y':
1292                                atom[4] = float(data[labels[key]].split('(')[0])
1293                            elif key == 'fract_z':
1294                                atom[5] = float(data[labels[key]].split('(')[0])
1295                            elif key == 'occupancy':
1296                                atom[6] = float(data[labels[key]].split('(')[0])
1297                            elif key == 'thermal_displace_type':
1298                                if data[labels[key]].lower() == 'uiso':
1299                                    atom[9] = 'I'
1300                                    atom[10] = float(data[labels['u_iso_or_equiv']].split('(')[0])
1301                                else:
1302                                    atom[9] = 'A'
1303                                    atom[10] = 0.0
1304                                   
1305                        atom[7],atom[8] = G2spc.SytSym(atom[3:6],SGData)
1306                        atom.append(ran.randint(0,sys.maxint))
1307                        Atoms.append(atom)
1308                        S = file.readline()
1309                else:
1310                    while S:
1311                        S.strip()
1312                        data = S.split()
1313                        if len(data) != len(labels):
1314                            break
1315                        name = data[labels['aniso_label']]
1316                        for atom in Atoms:
1317                            if name == atom[0]:
1318                                for i,uname in enumerate(anisoNames):
1319                                    atom[i+11] = float(data[labels[uname]].split('(')[0])
1320                        S = file.readline()
1321                                                                       
1322        else:           
1323            S = file.readline()
1324    file.close()
1325    if Title:
1326        PhaseName = Title
1327    else:
1328        PhaseName = 'None'
1329    SGlines = G2spc.SGPrint(SGData)
1330    for line in SGlines: print line
1331    cell = [float(a),float(b),float(c),float(alp),float(bet),float(gam)]
1332    Volume = G2lat.calc_V(G2lat.cell2A(cell))
1333    Phase['General'] = {'Name':PhaseName,'Type':'nuclear','SGData':SGData,
1334        'Cell':[False,]+cell+[Volume,]}
1335    Phase['Atoms'] = Atoms
1336    Phase['Drawing'] = {}
1337    Phase['Histograms'] = {}
1338   
1339    return Phase
1340
1341def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
1342    # returns value(esd) string; nTZ=True for no trailing zeros
1343    # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
1344    #get the 2 significant digits in the esd
1345    edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
1346    #get the number of digits to represent them
1347    epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
1348   
1349    mdec = lambda esd: -int(math.log10(abs(esd)))
1350    ndec = lambda esd: int(1.545-math.log10(abs(esd)))
1351    if esd > 0:
1352        fmt = '"%.'+str(ndec(esd))+'f(%d)"'
1353        print fmt,ndec(esd),esd*10**(mdec(esd)+1)
1354        return fmt%(value,int(esd*10**(mdec(esd)+1)))
1355    elif esd < 0:
1356         return str(round(value,mdec(esd)))
1357    else:
1358        text = "%F"%(value)
1359        if nTZ:
1360            return text.rstrip('0')
1361        else:
1362            return text
1363
1364def Fesd(value,esd=0,nTZ=False):
1365#pythonized version of fortran routine in GSAS cifsubs directory - doesn't work correctly
1366    nint = lambda x: int(round(x))
1367    iExp = 0
1368    if value == 0. and esd == 0.:
1369        iDec = 1
1370        iFld = 5
1371    elif value == 0.:
1372        iDec = max(0.,1.545-math.log10(abs(esd)))
1373        iFld = 4+iDec
1374    elif esd == 0.:
1375        iDec = 5
1376        iFld = max(1.,math.log10(abs(value)))+3+iDec
1377    else:
1378        iFld = math.log10(max(abs(esd),abs(value)))
1379        if iFld < -4:
1380            iExp = 1-iFld
1381            iFld -= iExp
1382        elif iFld > 8:
1383            iExp = -iFld
1384            iFld += iExp
1385        if iExp:
1386            value *= 10.**iExp
1387            esd *= 10.**iExp
1388        iDec = min(7,int(max(0.,1.545-math.log10(max(0.000001*abs(value),abs(esd))))))
1389        iFld = max(1,iFld)+3+iDec
1390    if esd <= 0.:
1391        iSigw = 0
1392    else:
1393        iSig = nint(esd*(10.**iDec))
1394        iSigw = 1
1395        if iSig > 0:
1396            iSigw = int(1.+math.log10(1.*iSig))
1397    if iSigw > 2:
1398        xmult = 10.**(iSigw-2)
1399        value = xmult*nint(value/xmult)
1400        iSig = xmult*nint(iSig/xmult)           
1401        iSigw = int(1.+math.log10(1.*iSig))
1402    if iSigw == 0:
1403        fmt = '%.'+str(iDec)+'f'
1404        string = fmt%(value) 
1405    elif iDec > 0:
1406        fmt = '%.'+str(iDec)+'f(%d)'
1407        string = fmt%(value,iSig)
1408    else:
1409        fmt = '%'+str(iFld)+'d(%d)'
1410        string = fmt%(nint(value),iSig)
1411    if iExp:
1412        iDig = 1+math.log10(abs(1.*iExp))
1413        if iExp > 0:
1414            iDig += 1
1415            string += str(-iExp)
1416    return string
1417   
1418
Note: See TracBrowser for help on using the repository browser.