source: trunk/GSASIIIO.py @ 399

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

Implement sequential refinement
remove print "load" & "save" for each item in Tree
revise application of azimuth offset - azimuths are now all "true" with correction

  • Property svn:keywords set to Date Author Revision URL Id
File size: 55.4 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-25 21:53:33 +0000 (Tue, 25 Oct 2011) $
6# $Author: vondreele $
7# $Revision: 399 $
8# $URL: trunk/GSASIIIO.py $
9# $Id: GSASIIIO.py 399 2011-10-25 21:53:33Z 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           
804            Id = self.PatternTree.AppendItem(parent=self.root,text=datum[0])
805            if 'PWDR' in datum[0]:               
806                self.PatternTree.SetItemPyData(Id,datum[1][:3])     #temp. trim off junk
807            else:
808                self.PatternTree.SetItemPyData(Id,datum[1])
809            for datus in data[1:]:
810                sub = self.PatternTree.AppendItem(Id,datus[0])
811                self.PatternTree.SetItemPyData(sub,datus[1])
812            if 'IMG' in datum[0]:                   #retreive image default flag & data if set
813                Data = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Image Controls'))
814                if Data['setDefault']:
815                    self.imageDefault = Data               
816        file.close()
817       
818    finally:
819        wx.EndBusyCursor()
820    print 'project load successful'
821    self.NewPlot = True
822   
823def ProjFileSave(self):
824    if not self.PatternTree.IsEmpty():
825        file = open(self.GSASprojectfile,'wb')
826        print 'save to file: ',self.GSASprojectfile
827        wx.BeginBusyCursor()
828        try:
829            item, cookie = self.PatternTree.GetFirstChild(self.root)
830            while item:
831                data = []
832                name = self.PatternTree.GetItemText(item)
833                data.append([name,self.PatternTree.GetItemPyData(item)])
834                item2, cookie2 = self.PatternTree.GetFirstChild(item)
835                while item2:
836                    name = self.PatternTree.GetItemText(item2)
837                    data.append([name,self.PatternTree.GetItemPyData(item2)])
838                    item2, cookie2 = self.PatternTree.GetNextChild(item, cookie2)                           
839                item, cookie = self.PatternTree.GetNextChild(self.root, cookie)                           
840                cPickle.dump(data,file,1)
841            file.close()
842        finally:
843            wx.EndBusyCursor()
844        print 'project save successful'
845       
846def SaveIntegration(self,PickId,data):
847    azms = self.Integrate[1]
848    X = self.Integrate[2][:-1]
849    Xminmax = [X[0],X[-1]]
850    N = len(X)
851    Id = self.PatternTree.GetItemParent(PickId)
852    name = self.PatternTree.GetItemText(Id)
853    name = name.replace('IMG ','PWDR ')
854    Comments = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Comments'))
855    names = ['Type','Lam','Zero','Polariz.','U','V','W','X','Y','SH/L','Azimuth'] 
856    codes = [0 for i in range(11)]
857    LRazm = data['LRazimuth']
858    Azms = []
859    if data['fullIntegrate'] and data['outAzimuths'] == 1:
860        Azms = [45.0,]                              #a poor man's average?
861    else:
862        for i,azm in enumerate(azms[:-1]):
863            Azms.append((azms[i+1]+azm)/2.)
864    for i,azm in enumerate(azms[:-1]):
865        item, cookie = self.PatternTree.GetFirstChild(self.root)
866        Id = 0
867        while item:
868            Name = self.PatternTree.GetItemText(item)
869            if name == Name:
870                Id = item
871            item, cookie = self.PatternTree.GetNextChild(self.root, cookie)
872        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!
873        Y = self.Integrate[0][i]
874        W = 1./Y                    #probably not true
875        Sample = G2pdG.SetDefaultSample()
876        if Id:
877            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Comments'),Comments)                   
878            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Limits'),[tuple(Xminmax),Xminmax])
879            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Background'),[['chebyschev',1,3,1.0,0.0,0.0]])
880            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Instrument Parameters'),[tuple(parms),parms[:],codes,names])
881            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Peak List'),[])
882            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Index Peak List'),[])
883            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Unit Cells List'),[])             
884            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Reflection Lists'),{})             
885        else:
886            Id = self.PatternTree.AppendItem(parent=self.root,text=name+" Azm= %.2f"%(Azms[i]))
887            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Comments'),Comments)                   
888            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Limits'),[tuple(Xminmax),Xminmax])
889            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Background'),[['chebyschev',1,3,1.0,0.0,0.0]])
890            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Instrument Parameters'),[tuple(parms),parms[:],codes,names])
891            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Sample Parameters'),Sample)
892            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Peak List'),[])
893            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Index Peak List'),[])
894            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Unit Cells List'),[])
895            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Reflection Lists'),{})             
896        self.PatternTree.SetItemPyData(Id,[[''],[np.array(X),np.array(Y),np.array(W),np.zeros(N),np.zeros(N),np.zeros(N)]])
897    self.PatternTree.SelectItem(Id)
898    self.PatternTree.Expand(Id)
899    self.PatternId = Id
900           
901def powderFxyeSave(self,exports,powderfile):
902    head,tail = ospath.split(powderfile)
903    name,ext = tail.split('.')
904    wx.BeginBusyCursor()
905    for i,export in enumerate(exports):
906        filename = ospath.join(head,name+'-%03d.'%(i)+ext)
907        prmname = filename.strip(ext)+'prm'
908        prm = open(prmname,'w')      #old style GSAS parm file
909        PickId = G2gd.GetPatternTreeItemId(self, self.root, export)
910        Values,Names = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self, \
911            PickId, 'Instrument Parameters'))[1::2]     #get values & names
912        Inst = dict(zip(Names,Values))
913        print Inst['Type']
914        prm.write( '            123456789012345678901234567890123456789012345678901234567890        '+'\n')
915        prm.write( 'INS   BANK      1                                                               '+'\n')
916        prm.write(('INS   HTYPE   %sR                                                              '+'\n')%(Inst['Type']))
917        if 'Lam1' in Inst:              #Ka1 & Ka2
918            prm.write(('INS  1 ICONS%10.7f%10.7f    0.0000               0.990    0     0.500   '+'\n')%(Inst['Lam1'],Inst['Lam2']))
919        elif 'Lam' in Inst:             #single wavelength
920            prm.write(('INS  1 ICONS%10.7f%10.7f    0.0000               0.990    0     0.500   '+'\n')%(Inst['Lam'],0.0))
921        prm.write( 'INS  1 IRAD     0                                                               '+'\n')
922        prm.write( 'INS  1I HEAD                                                                    '+'\n')
923        prm.write( 'INS  1I ITYP    0    0.0000  180.0000         1                                 '+'\n')
924        prm.write(('INS  1DETAZM%10.3f                                                          '+'\n')%(Inst['Azimuth']))
925        prm.write( 'INS  1PRCF1     3    8   0.00100                                                '+'\n')
926        prm.write(('INS  1PRCF11     %15.6g%15.6g%15.6g%15.6g   '+'\n')%(Inst['U'],Inst['V'],Inst['W'],0.0))
927        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.))
928        prm.close()
929        file = open(filename,'w')
930        print 'save powder pattern to file: ',filename
931        try:
932            x,y,w,yc,yb,yd = self.PatternTree.GetItemPyData(PickId)[1]
933            file.write(powderfile+'\n')
934            file.write('BANK 1 %d %d CONS %.2f %.2f 0 0 FXYE\n'%(len(x),len(x),\
935                100.*x[0],100.*(x[1]-x[0])))
936            s = list(np.sqrt(1./np.array(w)))       
937            XYW = zip(x,y,s)
938            for X,Y,S in XYW:
939                file.write("%15.6g %15.6g %15.6g\n" % (100.*X,Y,max(S,1.0)))
940            file.close()
941        finally:
942            wx.EndBusyCursor()
943        print 'powder pattern file written'
944       
945def powderXyeSave(self,exports,powderfile):
946    head,tail = ospath.split(powderfile)
947    name,ext = tail.split('.')
948    for i,export in enumerate(exports):
949        filename = ospath.join(head,name+'-%03d.'%(i)+ext)
950        PickId = G2gd.GetPatternTreeItemId(self, self.root, export)
951        file = open(filename,'w')
952        file.write('#%s\n'%(export))
953        print 'save powder pattern to file: ',filename
954        wx.BeginBusyCursor()
955        try:
956            x,y,w,yc,yb,yd = self.PatternTree.GetItemPyData(PickId)[1]
957            s = list(np.sqrt(1./np.array(w)))       
958            XYW = zip(x,y,s)
959            for X,Y,W in XYW:
960                file.write("%15.6g %15.6g %15.6g\n" % (X,Y,W))
961            file.close()
962        finally:
963            wx.EndBusyCursor()
964        print 'powder pattern file written'
965       
966def PDFSave(self,exports):   
967    for export in exports:
968        PickId = G2gd.GetPatternTreeItemId(self, self.root, export)
969        SQname = 'S(Q)'+export[4:]
970        GRname = 'G(R)'+export[4:]
971        sqfilename = ospath.join(self.dirname,export.replace(' ','_')[5:]+'.sq')
972        grfilename = ospath.join(self.dirname,export.replace(' ','_')[5:]+'.gr')
973        sqId = G2gd.GetPatternTreeItemId(self, PickId, SQname)
974        grId = G2gd.GetPatternTreeItemId(self, PickId, GRname)
975        sqdata = np.array(self.PatternTree.GetItemPyData(sqId)[1][:2]).T
976        grdata = np.array(self.PatternTree.GetItemPyData(grId)[1][:2]).T
977        sqfile = open(sqfilename,'w')
978        grfile = open(grfilename,'w')
979        sqfile.write('#T S(Q) %s\n'%(export))
980        grfile.write('#T G(R) %s\n'%(export))
981        sqfile.write('#L Q     S(Q)\n')
982        grfile.write('#L R     G(R)\n')
983        for q,sq in sqdata:
984            sqfile.write("%15.6g %15.6g\n" % (q,sq))
985        sqfile.close()
986        for r,gr in grdata:
987            grfile.write("%15.6g %15.6g\n" % (r,gr))
988        grfile.close()
989   
990def PeakListSave(self,file,peaks):
991    print 'save peak list to file: ',self.peaklistfile
992    if not peaks:
993        dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
994        try:
995            result = dlg.ShowModal()
996        finally:
997            dlg.Destroy()
998        return
999    for peak in peaks:
1000        file.write("%10.4f %12.2f %10.3f %10.3f \n" % \
1001            (peak[0],peak[2],peak[4],peak[6]))
1002    print 'peak list saved'
1003             
1004def IndexPeakListSave(self,peaks):
1005    file = open(self.peaklistfile,'wa')
1006    print 'save index peak list to file: ',self.peaklistfile
1007    wx.BeginBusyCursor()
1008    try:
1009        if not peaks:
1010            dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
1011            try:
1012                result = dlg.ShowModal()
1013            finally:
1014                dlg.Destroy()
1015            return
1016        for peak in peaks:
1017            file.write("%12.6f\n" % (peak[7]))
1018        file.close()
1019    finally:
1020        wx.EndBusyCursor()
1021    print 'index peak list saved'
1022   
1023def ReadEXPPhase(self,filename):
1024    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1025    textureData = {'Order':0,'Model':'cylindrical','Sample omega':[False,0.0],
1026        'Sample chi':[False,0.0],'Sample phi':[False,0.0],'SH Coeff':[False,{}],
1027        'SHShow':False,'PFhkl':[0,0,1],'PFxyz':[0,0,1],'PlotType':'Pole figure'}
1028    shNcof = 0
1029    file = open(filename, 'Ur')
1030    S = 1
1031    Expr = [{},{},{},{},{},{},{},{},{}]
1032    while S:
1033        S = file.readline()
1034        if 'EXPR NPHAS' in S[:12]:
1035            Num = S[12:-1].count('0')
1036            NPhas = S[12:-1].split()
1037        if 'CRS' in S[:3]:
1038            N = int(S[3:4])-1
1039            Expr[N][S[:12]] = S[12:-1]
1040    file.close()
1041    PNames = []
1042    for n,N in enumerate(NPhas):
1043        if N != '0':
1044            result = n
1045            key = 'CRS'+str(n+1)+'    PNAM'
1046            PNames.append(Expr[n][key])
1047    if Num < 8:
1048        dlg = wx.SingleChoiceDialog(self, 'Which phase to read?', 'Read phase data', PNames, wx.CHOICEDLG_STYLE)
1049        try:
1050            if dlg.ShowModal() == wx.ID_OK:
1051                result = dlg.GetSelection()
1052        finally:
1053            dlg.Destroy()       
1054    EXPphase = Expr[result]
1055    keyList = EXPphase.keys()
1056    keyList.sort()
1057    SGData = {}
1058    if NPhas[result] == '1':
1059        Ptype = 'nuclear'
1060    elif NPhas[result] in ['2','3']:
1061        Ptype = 'magnetic'
1062    elif NPhas[result] == '4':
1063        Ptype = 'macromolecular'
1064    elif NPhas[result] == '10':
1065        Ptype = 'Pawley'
1066    for key in keyList:
1067        if 'PNAM' in key:
1068           PhaseName = EXPphase[key].strip()
1069        elif 'ABC   ' in key:
1070            abc = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                       
1071        elif 'ANGLES' in key:
1072            angles = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                                               
1073        elif 'SG SYM' in key:
1074            SpGrp = EXPphase[key][:15].strip()
1075            E,SGData = G2spc.SpcGroup(SpGrp)
1076        elif 'OD    ' in key:
1077            SHdata = EXPphase[key].split()
1078            textureData['Order'] = int(SHdata[0])
1079            textureData['Model'] = shModels[int(SHdata[2])]
1080            textureData['Sample omega'] = [False,float(SHdata[6])]
1081            textureData['Sample chi'] = [False,float(SHdata[7])]
1082            textureData['Sample phi'] = [False,float(SHdata[8])]
1083            shNcof = int(SHdata[1])
1084    Atoms = []
1085    if Ptype == 'nuclear':
1086        for key in keyList:
1087            if 'AT' in key:
1088                if key[11:] == 'A':
1089                    S = EXPphase[key]
1090                elif key[11:] == 'B':
1091                    S += EXPphase[key]
1092                    Atom = [S[50:58].strip(),S[:10].strip(),'',
1093                        float(S[10:20]),float(S[20:30]),float(S[30:40]),
1094                        float(S[40:50]),'',int(S[60:62]),S[130:131]]
1095                    if Atom[9] == 'I':
1096                        Atom += [float(S[68:78]),0.,0.,0.,0.,0.,0.]
1097                    elif Atom[9] == 'A':
1098                        Atom += [0.0,float(S[68:78]),float(S[78:88]),
1099                            float(S[88:98]),float(S[98:108]),
1100                            float(S[108:118]),float(S[118:128])]
1101                    XYZ = Atom[3:6]
1102                    Atom[7],Atom[8] = G2spc.SytSym(XYZ,SGData)
1103                    Atom.append(ran.randint(0,sys.maxint))
1104                    Atoms.append(Atom)
1105    elif Ptype == 'macromolecular':
1106        for key in keyList:
1107            if 'AT' in key[6:8]:
1108                S = EXPphase[key]
1109                Atom = [int(S[56:60]),S[50:54].strip().upper(),S[54:56],
1110                    S[46:51].strip(),S[:8].strip(),'',
1111                    float(S[16:24]),float(S[24:32]),float(S[32:40]),
1112                    float(S[8:16]),'1',1,'I',float(S[40:46]),0,0,0,0,0,0]
1113                XYZ = Atom[6:9]
1114                Atom[10],Atom[11] = G2spc.SytSym(XYZ,SGData)
1115                Atom.append(ran.randint(0,sys.maxint))
1116                Atoms.append(Atom)
1117    Volume = G2lat.calc_V(G2lat.cell2A(abc+angles))
1118    if shNcof:
1119        shCoef = {}
1120        nRec = [i+1 for i in range((shNcof-1)/6+1)]
1121        for irec in nRec:
1122            ODkey = keyList[0][:6]+'OD'+'%3dA'%(irec)
1123            indx = EXPphase[ODkey].split()
1124            ODkey = ODkey[:-1]+'B'
1125            vals = EXPphase[ODkey].split()
1126            for i,val in enumerate(vals):
1127                key = 'C(%s,%s,%s)'%(indx[3*i],indx[3*i+1],indx[3*i+2])
1128                shCoef[key] = float(val)
1129        textureData['SH Coeff'] = [False,shCoef]
1130       
1131    Phase = {
1132            'General':{
1133                'Name':PhaseName,
1134                'Type':Ptype,
1135                'SGData':SGData,
1136                'Cell':[False,]+abc+angles+[Volume,],
1137                'Pawley dmin':1.0},
1138            'Atoms':Atoms,
1139            'Drawing':{},
1140            'Histograms':{},
1141            'Pawley ref':[],
1142            'Models':{},
1143            'SH Texture':textureData
1144            }
1145           
1146    return Phase
1147       
1148def ReadPDBPhase(filename):
1149    EightPiSq = 8.*math.pi**2
1150    file = open(filename, 'Ur')
1151    Phase = {}
1152    Title = ''
1153    Compnd = ''
1154    Atoms = []
1155    A = np.zeros(shape=(3,3))
1156    S = file.readline()
1157    while S:
1158        Atom = []
1159        if 'TITLE' in S[:5]:
1160            Title = S[10:72].strip()
1161            S = file.readline()
1162        elif 'COMPND    ' in S[:10]:
1163            Compnd = S[10:72].strip()
1164            S = file.readline()
1165        elif 'CRYST' in S[:5]:
1166            abc = S[7:34].split()
1167            angles = S[34:55].split()
1168            cell=[float(abc[0]),float(abc[1]),float(abc[2]),
1169                float(angles[0]),float(angles[1]),float(angles[2])]
1170            Volume = G2lat.calc_V(G2lat.cell2A(cell))
1171            AA,AB = G2lat.cell2AB(cell)
1172            SpGrp = S[55:65]
1173            E,SGData = G2spc.SpcGroup(SpGrp)
1174            if E: 
1175                print ' ERROR in space group symbol ',SpGrp,' in file ',filename
1176                print ' N.B.: make sure spaces separate axial fields in symbol' 
1177                print G2spc.SGErrors(E)
1178                return None
1179            SGlines = G2spc.SGPrint(SGData)
1180            for line in SGlines: print line
1181            S = file.readline()
1182        elif 'SCALE' in S[:5]:
1183            V = (S[10:41].split())
1184            A[int(S[5])-1] = [float(V[0]),float(V[1]),float(V[2])]
1185            S = file.readline()
1186        elif 'ATOM' in S[:4] or 'HETATM' in S[:6]:
1187            XYZ = [float(S[31:39]),float(S[39:47]),float(S[47:55])]
1188            XYZ = np.inner(AB,XYZ)
1189            SytSym,Mult = G2spc.SytSym(XYZ,SGData)
1190            Uiso = float(S[61:67])/EightPiSq
1191            Type = S[12:14].upper()
1192            if Type[0] in '123456789':
1193                Type = Type[1:]
1194            Atom = [S[22:27].strip(),S[17:20].upper(),S[20:22],
1195                S[12:17].strip(),Type.strip(),'',XYZ[0],XYZ[1],XYZ[2],
1196                float(S[55:61]),SytSym,Mult,'I',Uiso,0,0,0,0,0,0]
1197            S = file.readline()
1198            if 'ANISOU' in S[:6]:
1199                Uij = S[30:72].split()
1200                Uij = [float(Uij[0])/10000.,float(Uij[1])/10000.,float(Uij[2])/10000.,
1201                    float(Uij[3])/10000.,float(Uij[4])/10000.,float(Uij[5])/10000.]
1202                Atom = Atom[:14]+Uij
1203                Atom[12] = 'A'
1204                S = file.readline()
1205            Atom.append(ran.randint(0,sys.maxint))
1206            Atoms.append(Atom)
1207        else:           
1208            S = file.readline()
1209    file.close()
1210    if Title:
1211        PhaseName = Title
1212    elif Compnd:
1213        PhaseName = Compnd
1214    else:
1215        PhaseName = 'None'
1216    Phase['General'] = {'Name':PhaseName,'Type':'macromolecular','SGData':SGData,
1217        'Cell':[False,]+cell+[Volume,]}
1218    Phase['Atoms'] = Atoms
1219    Phase['Drawing'] = {}
1220    Phase['Histograms'] = {}
1221   
1222    return Phase
1223   
1224def ReadCIFPhase(filename):
1225    anisoNames = ['aniso_u_11','aniso_u_22','aniso_u_33','aniso_u_12','aniso_u_13','aniso_u_23']
1226    file = open(filename, 'Ur')
1227    Phase = {}
1228    Title = ospath.split(filename)[-1]
1229    print '\n Reading cif file: ',Title
1230    Compnd = ''
1231    Atoms = []
1232    A = np.zeros(shape=(3,3))
1233    S = file.readline()
1234    while S:
1235        if '_symmetry_space_group_name_H-M' in S:
1236            SpGrp = S.split("_symmetry_space_group_name_H-M")[1].strip().strip('"').strip("'")
1237            E,SGData = G2spc.SpcGroup(SpGrp)
1238            if E:
1239                print ' ERROR in space group symbol ',SpGrp,' in file ',filename
1240                print ' N.B.: make sure spaces separate axial fields in symbol' 
1241                print G2spc.SGErrors(E)
1242                return None
1243            S = file.readline()
1244        elif '_cell' in S:
1245            if '_cell_length_a' in S:
1246                a = S.split('_cell_length_a')[1].strip().strip('"').strip("'").split('(')[0]
1247            elif '_cell_length_b' in S:
1248                b = S.split('_cell_length_b')[1].strip().strip('"').strip("'").split('(')[0]
1249            elif '_cell_length_c' in S:
1250                c = S.split('_cell_length_c')[1].strip().strip('"').strip("'").split('(')[0]
1251            elif '_cell_angle_alpha' in S:
1252                alp = S.split('_cell_angle_alpha')[1].strip().strip('"').strip("'").split('(')[0]
1253            elif '_cell_angle_beta' in S:
1254                bet = S.split('_cell_angle_beta')[1].strip().strip('"').strip("'").split('(')[0]
1255            elif '_cell_angle_gamma' in S:
1256                gam = S.split('_cell_angle_gamma')[1].strip().strip('"').strip("'").split('(')[0]
1257            S = file.readline()
1258        elif 'loop_' in S:
1259            labels = {}
1260            i = 0
1261            while S:
1262                S = file.readline()
1263                if '_atom_site' in S.strip()[:10]:
1264                    labels[S.strip().split('_atom_site_')[1].lower()] = i
1265                    i += 1
1266                else:
1267                    break
1268            if labels:
1269                if 'aniso_label' not in labels:
1270                    while S:
1271                        atom = ['','','',0,0,0,1.0,'','','I',0.01,0,0,0,0,0,0]
1272                        S.strip()
1273                        if len(S.split()) != len(labels):
1274                            if 'loop_' in S:
1275                                break
1276                            S += file.readline().strip()
1277                        data = S.split()
1278                        if len(data) != len(labels):
1279                            break
1280                        for key in labels:
1281                            if key == 'type_symbol':
1282                                atom[1] = data[labels[key]]
1283                            elif key == 'label':
1284                                atom[0] = data[labels[key]]
1285                            elif key == 'fract_x':
1286                                atom[3] = float(data[labels[key]].split('(')[0])
1287                            elif key == 'fract_y':
1288                                atom[4] = float(data[labels[key]].split('(')[0])
1289                            elif key == 'fract_z':
1290                                atom[5] = float(data[labels[key]].split('(')[0])
1291                            elif key == 'occupancy':
1292                                atom[6] = float(data[labels[key]].split('(')[0])
1293                            elif key == 'thermal_displace_type':
1294                                if data[labels[key]].lower() == 'uiso':
1295                                    atom[9] = 'I'
1296                                    atom[10] = float(data[labels['u_iso_or_equiv']].split('(')[0])
1297                                else:
1298                                    atom[9] = 'A'
1299                                    atom[10] = 0.0
1300                                   
1301                        atom[7],atom[8] = G2spc.SytSym(atom[3:6],SGData)
1302                        atom.append(ran.randint(0,sys.maxint))
1303                        Atoms.append(atom)
1304                        S = file.readline()
1305                else:
1306                    while S:
1307                        S.strip()
1308                        data = S.split()
1309                        if len(data) != len(labels):
1310                            break
1311                        name = data[labels['aniso_label']]
1312                        for atom in Atoms:
1313                            if name == atom[0]:
1314                                for i,uname in enumerate(anisoNames):
1315                                    atom[i+11] = float(data[labels[uname]].split('(')[0])
1316                        S = file.readline()
1317                                                                       
1318        else:           
1319            S = file.readline()
1320    file.close()
1321    if Title:
1322        PhaseName = Title
1323    else:
1324        PhaseName = 'None'
1325    SGlines = G2spc.SGPrint(SGData)
1326    for line in SGlines: print line
1327    cell = [float(a),float(b),float(c),float(alp),float(bet),float(gam)]
1328    Volume = G2lat.calc_V(G2lat.cell2A(cell))
1329    Phase['General'] = {'Name':PhaseName,'Type':'nuclear','SGData':SGData,
1330        'Cell':[False,]+cell+[Volume,]}
1331    Phase['Atoms'] = Atoms
1332    Phase['Drawing'] = {}
1333    Phase['Histograms'] = {}
1334   
1335    return Phase
1336
1337def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
1338    # returns value(esd) string; nTZ=True for no trailing zeros
1339    # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
1340    #get the 2 significant digits in the esd
1341    edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
1342    #get the number of digits to represent them
1343    epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
1344   
1345    mdec = lambda esd: -int(math.log10(abs(esd)))
1346    ndec = lambda esd: int(1.545-math.log10(abs(esd)))
1347    if esd > 0:
1348        fmt = '"%.'+str(ndec(esd))+'f(%d)"'
1349        print fmt,ndec(esd),esd*10**(mdec(esd)+1)
1350        return fmt%(value,int(esd*10**(mdec(esd)+1)))
1351    elif esd < 0:
1352         return str(round(value,mdec(esd)))
1353    else:
1354        text = "%F"%(value)
1355        if nTZ:
1356            return text.rstrip('0')
1357        else:
1358            return text
1359
1360def Fesd(value,esd=0,nTZ=False):
1361#pythonized version of fortran routine in GSAS cifsubs directory - doesn't work correctly
1362    nint = lambda x: int(round(x))
1363    iExp = 0
1364    if value == 0. and esd == 0.:
1365        iDec = 1
1366        iFld = 5
1367    elif value == 0.:
1368        iDec = max(0.,1.545-math.log10(abs(esd)))
1369        iFld = 4+iDec
1370    elif esd == 0.:
1371        iDec = 5
1372        iFld = max(1.,math.log10(abs(value)))+3+iDec
1373    else:
1374        iFld = math.log10(max(abs(esd),abs(value)))
1375        if iFld < -4:
1376            iExp = 1-iFld
1377            iFld -= iExp
1378        elif iFld > 8:
1379            iExp = -iFld
1380            iFld += iExp
1381        if iExp:
1382            value *= 10.**iExp
1383            esd *= 10.**iExp
1384        iDec = min(7,int(max(0.,1.545-math.log10(max(0.000001*abs(value),abs(esd))))))
1385        iFld = max(1,iFld)+3+iDec
1386    if esd <= 0.:
1387        iSigw = 0
1388    else:
1389        iSig = nint(esd*(10.**iDec))
1390        iSigw = 1
1391        if iSig > 0:
1392            iSigw = int(1.+math.log10(1.*iSig))
1393    if iSigw > 2:
1394        xmult = 10.**(iSigw-2)
1395        value = xmult*nint(value/xmult)
1396        iSig = xmult*nint(iSig/xmult)           
1397        iSigw = int(1.+math.log10(1.*iSig))
1398    if iSigw == 0:
1399        fmt = '%.'+str(iDec)+'f'
1400        string = fmt%(value) 
1401    elif iDec > 0:
1402        fmt = '%.'+str(iDec)+'f(%d)'
1403        string = fmt%(value,iSig)
1404    else:
1405        fmt = '%'+str(iFld)+'d(%d)'
1406        string = fmt%(nint(value),iSig)
1407    if iExp:
1408        iDig = 1+math.log10(abs(1.*iExp))
1409        if iExp > 0:
1410            iDig += 1
1411            string += str(-iExp)
1412    return string
1413   
1414
Note: See TracBrowser for help on using the repository browser.