source: trunk/GSASIIIO.py @ 391

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

set default atom frac=1 for cif files
store covMatrix not normalized by diagonal
make esds on cell parameters

  • Property svn:keywords set to Date Author Revision URL Id
File size: 55.6 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-13 15:48:51 +0000 (Thu, 13 Oct 2011) $
6# $Author: vondreele $
7# $Revision: 391 $
8# $URL: trunk/GSASIIIO.py $
9# $Id: GSASIIIO.py 391 2011-10-13 15:48:51Z 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 GSASIIElem as G2el
22import os.path as ospath
23
24def sfloat(S):
25    if S.strip():
26        return float(S)
27    else:
28        return 0.0
29
30def sint(S):
31    if S.strip():
32        return int(S)
33    else:
34        return 0
35
36def SelectPowderData(self, filename):
37    """Selects banks of data from a filename of any GSAS powder data format
38    Input - filename: any GSAS powder data formatted file (currently STD, FXYE, FXY & ESD)
39    Returns - a list of banks to be read; each entry in list is a tuple containing:
40    filename: same as input filename
41    Pos: position for start of data; record just after BANK record
42    Bank: the BANK record
43    """
44    File = open(filename,'Ur')
45    Title = '''
46First line of this file:
47'''+File.readline()
48    dlg = wx.MessageDialog(self, Title, 'Is this the file you want?', 
49        wx.YES_NO | wx.ICON_QUESTION)
50    try:
51        result = dlg.ShowModal()
52    finally:
53        dlg.Destroy()
54    if result == wx.ID_NO: return (0,0)
55    Temperature = 300
56   
57    if '.xye' in filename:      #Topas style xye file (e.g. 2-th, I, sig) - no iparm file/no BANK record
58        dlg = wx.MessageDialog(self,'''Is this laboratory Cu Ka1/Ka2 data?
59(No = 0.6A wavelength synchrotron data)
60Change wavelength in Instrument Parameters if needed''','Data type?',
61            wx.YES_NO | wx.ICON_QUESTION)
62        try:
63            result = dlg.ShowModal()
64        finally:
65            dlg.Destroy()
66        print result
67        if result == wx.ID_YES:
68            Iparm = {}                                               #Assume CuKa lab data
69            Iparm['INS   HTYPE '] = 'PXC '
70            Iparm['INS  1 ICONS'] = '  1.540500  1.544300       0.0         0       0.7    0       0.5   '
71            Iparm['INS  1PRCF1 '] = '    3    8      0.01                                                '
72            Iparm['INS  1PRCF11'] = '   2.000000E+00  -2.000000E+00   5.000000E+00   0.000000E+00        '
73            Iparm['INS  1PRCF12'] = '   0.000000E+00   0.000000E+00   0.150000E-01   0.150000E-01        '
74        else:
75            Iparm = {}                                               #Assume 0.6A synchrotron data
76            Iparm['INS   HTYPE '] = 'PXC '
77            Iparm['INS  1 ICONS'] = '  0.600000  0.000000       0.0         0      0.99    0       0.5   '
78            Iparm['INS  1PRCF1 '] = '    3    8      0.01                                                '
79            Iparm['INS  1PRCF11'] = '   1.000000E+00  -1.000000E+00   0.300000E+00   0.000000E+00        '
80            Iparm['INS  1PRCF12'] = '   0.000000E+00   0.000000E+00   0.100000E-01   0.100000E-01        '
81                       
82       
83    else:                       #GSAS style fxye or fxy file (e.g. 100*2-th, I, sig)
84        self.IparmName = GetInstrumentFile(self,filename)
85        if self.IparmName:
86            Iparm = GetInstrumentData(self.IparmName)
87        else:
88            Iparm = {}                                               #Assume CuKa lab data if no iparm file
89            Iparm['INS   HTYPE '] = 'PXC '
90            Iparm['INS  1 ICONS'] = '  1.540500  1.544300       0.0         0       0.7    0       0.5   '
91            Iparm['INS  1PRCF1 '] = '    3    8      0.01                                                '
92            Iparm['INS  1PRCF11'] = '   2.000000E+00  -2.000000E+00   5.000000E+00   0.000000E+00        '
93            Iparm['INS  1PRCF12'] = '   0.000000E+00   0.000000E+00   0.150000E-01   0.150000E-01        '
94    S = 1
95    Banks = []
96    Pos = []
97    FoundData = []
98    Comments = []
99    wx.BeginBusyCursor()
100    try:
101        while S:
102            S = File.readline()
103            if S[:1] != '#':
104                if S[:4] == 'BANK':
105                    Banks.append(S)
106                    Pos.append(File.tell())
107                elif '.xye' in filename:    #No BANK in a xye file
108                    Banks.append('BANK 1 XYE')
109                    Pos.append(File.tell())
110                    break
111            else:
112                Comments.append(S[:-1])
113                if 'Temp' in S.split('=')[0]:
114                    Temperature = float(S.split('=')[1])
115        File.close()
116    finally:
117        wx.EndBusyCursor()
118    if Comments:
119       print 'Comments on file:'
120       for Comment in Comments: print Comment
121    if Banks:
122        result = [0]
123        if len(Banks) >= 2:
124            dlg = wx.MultiChoiceDialog(self, 'Which scans do you want?', 'Select scans', Banks, wx.CHOICEDLG_STYLE)
125            try:
126                if dlg.ShowModal() == wx.ID_OK:
127                    result = dlg.GetSelections()
128                else:
129                    result = []
130            finally:
131                dlg.Destroy()
132        for i in result:
133            FoundData.append((filename,Pos[i],Banks[i]))
134    else:
135        dlg = wx.MessageDialog(self, 'ERROR - this is not a GSAS powder data file', 'No BANK records', wx.OK | wx.ICON_ERROR)
136        try:
137            result = dlg.ShowModal()
138        finally:
139            dlg.Destroy()
140    return FoundData,Iparm,Comments,Temperature
141
142def GetInstrumentFile(self,filename):
143    import os.path as op
144    dlg = wx.FileDialog(self,'Choose an instrument file','.', '', 'GSAS iparm file (*.prm)|*.prm|All files(*.*)|*.*', wx.OPEN)
145    if self.dirname: 
146        dlg.SetDirectory(self.dirname)
147        Tname = filename[:filename.index('.')]+'.prm'
148        if op.exists(Tname):
149            self.IparmName = Tname       
150    if self.IparmName: dlg.SetFilename(self.IparmName)
151    filename = ''
152    try:
153        if dlg.ShowModal() == wx.ID_OK:
154            filename = dlg.GetPath()
155    finally:
156        dlg.Destroy()
157    return filename
158
159def GetInstrumentData(IparmName):
160    file = open(IparmName, 'Ur')
161    S = 1
162    Iparm = {}
163    while S:
164        S = file.readline()
165        Iparm[S[:12]] = S[12:-1]
166    return Iparm
167   
168def GetPowderPeaks(fileName):
169    sind = lambda x: math.sin(x*math.pi/180.)
170    asind = lambda x: 180.*math.asin(x)/math.pi
171    Cuka = 1.54052
172    File = open(fileName,'Ur')
173    Comments = []
174    peaks = []
175    S = File.readline()
176    while S:
177        if S[:1] == '#':
178            Comments.append(S[:-1])
179        else:
180            item = S.split()
181            if len(item) == 1:
182                peaks.append([float(item[0]),1.0])
183            elif len(item) > 1:
184                peaks.append([float(item[0]),float(item[0])])
185        S = File.readline()
186    File.close()
187    if Comments:
188       print 'Comments on file:'
189       for Comment in Comments: print Comment
190    Peaks = []
191    if peaks[0][0] > peaks[-1][0]:          # d-spacings - assume CuKa
192        for peak in peaks:
193            dsp = peak[0]
194            sth = Cuka/(2.0*dsp)
195            if sth < 1.0:
196                tth = 2.0*asind(sth)
197            else:
198                break
199            Peaks.append([tth,peak[1],True,False,0,0,0,dsp,0.0])
200    else:                                   #2-thetas - assume Cuka (for now)
201        for peak in peaks:
202            tth = peak[0]
203            dsp = Cuka/(2.0*sind(tth/2.0))
204            Peaks.append([tth,peak[1],True,False,0,0,0,dsp,0.0])
205    return Comments,Peaks
206
207def GetPawleyPeaks(filename):
208    rt2ln2x2 = 2.35482
209    File = open(filename,'Ur')
210    PawleyPeaks = []
211    S = File.readline()         #skip header
212    S = File.readline()
213    item = S.split()
214    while S:
215        h,k,l = int(item[0]),int(item[1]),int(item[2])
216        mult = int(item[3])
217        tth = float(item[5])
218        sig = float(item[6])/rt2ln2x2
219        Iobs = float(item[7])*mult
220        PawleyPeaks.append([h,k,l,mult,tth,False,Iobs,0.0])
221        S = File.readline()
222        item = S.split()
223        if item[3] == '-100.0000':       #find trailer
224            break
225    File.close()
226    return PawleyPeaks
227   
228def GetHKLData(filename):
229    print 'Reading: '+filename
230    File = open(filename,'Ur')
231    HKLref = []
232    HKLmin = [1000,1000,1000]
233    HKLmax = [0,0,0]
234    FoMax = 0
235    ifFc = False
236    S = File.readline()
237    while '#' in S[0]:        #get past comments if any
238        S = File.readline()       
239    if '_' in S:         #cif style .hkl file
240        while 'loop_' not in S:         #skip preliminaries if any - can't have 'loop_' in them!
241            S = File.readline()       
242        S = File.readline()             #get past 'loop_' line
243        pos = 0
244        hpos = kpos = lpos = Fosqpos = Fcsqpos = sigpos = -1
245        while S:
246            if '_' in S:
247                if 'index_h' in S:
248                    hpos = pos
249                elif 'index_k' in S:
250                    kpos = pos
251                elif 'index_l' in S:
252                    lpos = pos
253                elif 'F_squared_meas' in S:
254                    Fosqpos = pos
255                elif 'F_squared_calc' in S:
256                    Fcsqpos = pos
257                elif 'F_squared_sigma' in S:
258                    sigpos = pos
259                pos += 1
260            else:
261                data = S.split()
262                if data:                    #avoid blank lines
263                    HKL = np.array([int(data[hpos]),int(data[kpos]),int(data[lpos])])
264                    h,k,l = HKL
265                    Fosq = float(data[Fosqpos])
266                    if sigpos != -1:
267                        sigFosq = float(data[sigpos])
268                    else:
269                        sigFosq = 1.
270                    if Fcsqpos != -1:
271                        Fcsq = float(data[Fcsqpos])
272                        if Fcsq:
273                            ifFc = True
274                    else:
275                        Fcsq = 0.
276                       
277                    HKLmin = [min(h,HKLmin[0]),min(k,HKLmin[1]),min(l,HKLmin[2])]
278                    HKLmax = [max(h,HKLmax[0]),max(k,HKLmax[1]),max(l,HKLmax[2])]
279                    FoMax = max(FoMax,Fosq)
280                    HKLref.append([HKL,Fosq,sigFosq,Fcsq,0,0,0])                 #room for Fcp, Fcpp & phase
281            S = File.readline()
282    else:                   #dumb h,k,l,Fo,sigFo .hkl file
283        while S:
284            h,k,l,Fo,sigFo = S.split()
285            HKL = np.array([int(h),int(k),int(l)])
286            h,k,l = HKL
287            Fo = float(Fo)
288            sigFo = float(sigFo)
289            HKLmin = [min(h,HKLmin[0]),min(k,HKLmin[1]),min(l,HKLmin[2])]
290            HKLmax = [max(h,HKLmax[0]),max(k,HKLmax[1]),max(l,HKLmax[2])]
291            FoMax = max(FoMax,Fo)
292            HKLref.append([HKL,Fo**2,2.*Fo*sigFo,0,0,0,0])                 #room for Fc, Fcp, Fcpp & phase
293            S = File.readline()
294    File.close()
295    return HKLref,HKLmin,HKLmax,FoMax,ifFc
296
297def GetPowderData(filename,Pos,Bank,DataType):
298    '''Reads one BANK of data from GSAS raw powder data file
299    input:
300    filename: GSAS raw powder file dataname
301    Pos: start of data in file just after BANK record
302    Bank: the BANK record
303    DataType: powder data type, e.g. "PXC" for Powder X-ray CW data
304    returns: list [x,y,e,yc,yb]
305    x: np.array of x-axis values
306    y: np.array of powder pattern intensities
307    w: np.array of w=sig(intensity)^2 values
308    yc: np.array of calc. intensities (zero)
309    yb: np.array of calc. background (zero)
310    yd: np.array of obs-calc profiles
311    '''
312    print 'Reading: '+filename
313    print 'Bank:    '+Bank[:-1]
314    if 'FXYE' in Bank:
315        return GetFXYEdata(filename,Pos,Bank,DataType)
316    elif ' XYE' in Bank:
317        return GetXYEdata(filename,Pos,Bank,DataType)
318    elif 'FXY' in Bank:
319        return GetFXYdata(filename,Pos,Bank,DataType)
320    elif 'ESD' in Bank:
321        return GetESDdata(filename,Pos,Bank,DataType)
322    elif 'STD' in Bank:
323        return GetSTDdata(filename,Pos,Bank,DataType)
324    else:
325        return GetSTDdata(filename,Pos,Bank,DataType)
326    return []
327
328def GetFXYEdata(filename,Pos,Bank,DataType):
329    File = open(filename,'Ur')
330    File.seek(Pos)
331    x = []
332    y = []
333    w = []
334    S = File.readline()
335    while S and S[:4] != 'BANK':
336        vals = S.split()
337        if DataType[2] == 'C':
338            x.append(float(vals[0])/100.)               #CW: from centidegrees to degrees
339        elif DataType[2] == 'T':
340            x.append(float(vals[0])/1000.0)             #TOF: from musec to millisec
341        f = float(vals[1])
342        if f <= 0.0:
343            y.append(0.0)
344            w.append(1.0)
345        else:
346            y.append(float(vals[1]))
347            w.append(1.0/float(vals[2])**2)
348        S = File.readline()
349    File.close()
350    N = len(x)
351    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
352   
353def GetXYEdata(filename,Pos,Bank,DataType):
354    File = open(filename,'Ur')
355    File.seek(Pos)
356    x = []
357    y = []
358    w = []
359    S = File.readline()
360    while S:
361        vals = S.split()
362        try:
363            x.append(float(vals[0]))
364            f = float(vals[1])
365            if f <= 0.0:
366                y.append(0.0)
367                w.append(1.0)
368            else:
369                y.append(float(vals[1]))
370                w.append(1.0/float(vals[2])**2)
371            S = File.readline()
372        except ValueError:
373            break
374    File.close()
375    N = len(x)
376    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
377   
378   
379def GetFXYdata(filename,Pos,Bank,DataType):
380    File = open(filename,'Ur')
381    File.seek(Pos)
382    x = []
383    y = []
384    w = []
385    S = File.readline()
386    while S and S[:4] != 'BANK':
387        vals = S.split()
388        if DataType[2] == 'C':
389            x.append(float(vals[0])/100.)               #CW: from centidegrees to degrees
390        elif DataType[2] == 'T':
391            x.append(float(vals[0])/1000.0)             #TOF: from musec to millisec
392        f = float(vals[1])
393        if f > 0.0:
394            y.append(f)
395            w.append(1.0/f)
396        else:             
397            y.append(0.0)
398            w.append(1.0)
399        S = File.readline()
400    File.close()
401    N = len(x)
402    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
403   
404def GetESDdata(filename,Pos,Bank,DataType):
405    File = open(filename,'Ur')
406    cons = Bank.split()
407    if DataType[2] == 'C':
408        start = float(cons[5])/100.0               #CW: from centidegrees to degrees
409        step = float(cons[6])/100.0
410    elif DataType[2] == 'T':
411        start = float(cons[5])/1000.0              #TOF: from musec to millisec
412        step = float(cons[6])/1000.0
413    File.seek(Pos)
414    x = []
415    y = []
416    w = []
417    S = File.readline()
418    j = 0
419    while S and S[:4] != 'BANK':
420        for i in range(0,80,16):
421            xi = start+step*j
422            yi = sfloat(S[i:i+8])
423            ei = sfloat(S[i+8:i+16])
424            x.append(xi)
425            if yi > 0.0:
426                y.append(yi)
427                w.append(1.0/ei**2)
428            else:             
429                y.append(0.0)
430                w.append(1.0)
431            j += 1
432        S = File.readline()
433    File.close()
434    N = len(x)
435    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
436
437def GetSTDdata(filename,Pos,Bank,DataType):
438    File = open(filename,'Ur')
439    cons = Bank.split()
440    Nch = cons[2]
441    if DataType[2] == 'C':
442        start = float(cons[5])/100.0               #CW: from centidegrees to degrees
443        step = float(cons[6])/100.0
444    elif DataType[2] == 'T':
445        start = float(cons[5])/1000.0              #TOF: from musec to millisec - not likely!
446        step = float(cons[6])/1000.0
447    File.seek(Pos)
448    x = []
449    y = []
450    w = []
451    S = File.readline()
452    j = 0
453    while S and S[:4] != 'BANK':
454        for i in range(0,80,8):
455            xi = start+step*j
456            ni = max(sint(S[i:i+2]),1)
457            yi = max(sfloat(S[i+2:i+8]),0.0)
458            if yi:
459                ei = math.sqrt(yi*ni)
460            else:
461                yi = 0.0
462                ei = 1.0
463            j += 1
464            if j < Nch:
465                x.append(xi)
466                y.append(yi)
467                w.append(1.0/ei**2)
468        S = File.readline()
469    File.close()
470    N = len(x)
471    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
472   
473def CheckImageFile(self,imagefile):
474    if not ospath.exists(imagefile):
475        dlg = wx.FileDialog(self, 'Bad image file name; choose name', '.', '',\
476        'Any image file (*.tif;*.tiff;*.mar*;*.avg;*.sum;*.img)\
477        |*.tif;*.tiff;*.mar*;*.avg;*.sum;*.img|\
478        Any detector tif (*.tif;*.tiff)|*.tif;*.tiff|\
479        MAR file (*.mar*)|*.mar*|\
480        GE Image (*.avg;*.sum)|*.avg;*.sum|\
481        ADSC Image (*.img)|*.img|\
482        All files (*.*)|*.*',wx.OPEN)
483        if self.dirname:
484            dlg.SetDirectory(self.dirname)
485        try:
486            dlg.SetFilename(ospath.split(imagefile)[1])
487            if dlg.ShowModal() == wx.ID_OK:
488                self.dirname = dlg.GetDirectory()
489                imagefile = dlg.GetPath()
490            else:
491                imagefile = False
492        finally:
493            dlg.Destroy()
494    return imagefile
495       
496def GetImageData(self,imagefile,imageOnly=False):       
497    ext = ospath.splitext(imagefile)[1]
498    Comments = []
499    if ext == '.tif' or ext == '.tiff':
500        Comments,Data,Npix,Image = GetTifData(imagefile)
501    elif ext == '.img':
502        Comments,Data,Npix,Image = GetImgData(imagefile)
503        Image[0][0] = 0
504    elif ext == '.mar3450' or ext == '.mar2300':
505        Comments,Data,Npix,Image = GetMAR345Data(imagefile)
506    elif ext in ['.sum','.avg','']:
507        Comments,Data,Npix,Image = GetGEsumData(imagefile)
508    elif ext == '.G2img':
509        Comments,Data,Npix,Image = GetG2Image(imagefile)
510    if imageOnly:
511        return Image
512    else:
513        return Comments,Data,Npix,Image
514       
515def PutG2Image(filename,Comments,Data,Npix,image):
516    File = open(filename,'wb')
517    cPickle.dump([Comments,Data,Npix,image],File,1)
518    File.close()
519    return
520   
521def GetG2Image(filename):
522    File = open(filename,'rb')
523    Comments,Data,Npix,image = cPickle.load(File)
524    File.close()
525    return Comments,Data,Npix,image
526   
527def GetGEsumData(filename,imageOnly=False):
528    import struct as st
529    import array as ar
530    if not imageOnly:
531        print 'Read GE sum file: ',filename   
532    File = open(filename,'rb')
533    if '.sum' in filename:
534        head = ['GE detector sum data from APS 1-ID',]
535        sizexy = [2048,2048]
536    elif '.avg' in filename:
537        head = ['GE detector avg data from APS 1-ID',]
538        sizexy = [2048,2048]
539    else:
540        head = ['GE detector raw data from APS 1-ID',]
541        File.seek(18)
542        size,nframes = st.unpack('<ih',File.read(6))
543        sizexy = [2048,2048]
544        pos = 8192
545        File.seek(pos)
546    Npix = sizexy[0]*sizexy[1]
547    if '.sum' in filename:
548        image = np.array(ar.array('f',File.read(4*Npix)),dtype=np.int32)
549    elif '.avg' in filename:
550        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
551    else:
552        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
553        while nframes > 1:
554            image += np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
555            nframes -= 1
556    image = np.reshape(image,(sizexy[1],sizexy[0]))
557    data = {'pixelSize':(200,200),'wavelength':0.15,'distance':250.0,'center':[204.8,204.8],'size':sizexy} 
558    File.close()   
559    if imageOnly:
560        return image
561    else:
562        return head,data,Npix,image
563       
564def GetImgData(filename,imageOnly=False):
565    import struct as st
566    import array as ar
567    if not imageOnly:
568        print 'Read ADSC img file: ',filename
569    File = open(filename,'rb')
570    head = File.read(511)
571    lines = head.split('\n')
572    head = []
573    center = [0,0]
574    for line in lines[1:-2]:
575        line = line.strip()[:-1]
576        if line:
577            if 'SIZE1' in line:
578                size = int(line.split('=')[1])
579                Npix = size*size
580            elif 'WAVELENGTH' in line:
581                wave = float(line.split('=')[1])
582            elif 'BIN' in line:
583                if line.split('=')[1] == '2x2':
584                    pixel=(102,102)
585                else:
586                    pixel = (51,51)
587            elif 'DISTANCE' in line:
588                distance = float(line.split('=')[1])
589            elif 'CENTER_X' in line:
590                center[0] = float(line.split('=')[1])
591            elif 'CENTER_Y' in line:
592                center[1] = float(line.split('=')[1])
593            head.append(line)
594    data = {'pixelSize':pixel,'wavelength':wave,'distance':distance,'center':center,'size':[size,size]}
595    image = []
596    row = 0
597    pos = 512
598    File.seek(pos)
599    image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
600    image = np.reshape(image,(sizexy[1],sizexy[0]))
601#    image = np.zeros(shape=(size,size),dtype=np.int32)   
602#    while row < size:
603#        File.seek(pos)
604#        line = ar.array('H',File.read(2*size))
605#        image[row] = np.asarray(line)
606#        row += 1
607#        pos += 2*size
608    File.close()
609    if imageOnly:
610        return image
611    else:
612        return lines[1:-2],data,Npix,image
613       
614def GetMAR345Data(filename,imageOnly=False):
615    import array as ar
616    import struct as st
617    try:
618        import pack_f as pf
619    except:
620        msg = wx.MessageDialog(None, message="Unable to load the GSAS MAR image decompression, pack_f",
621                               caption="Import Error",
622                               style=wx.ICON_ERROR | wx.OK | wx.STAY_ON_TOP)
623        msg.ShowModal()
624        return None,None,None,None
625
626    if not imageOnly:
627        print 'Read Mar345 file: ',filename
628    File = open(filename,'rb')
629    head = File.read(4095)
630    numbers = st.unpack('<iiiiiiiiii',head[:40])
631    lines = head[128:].split('\n')
632    head = []
633    for line in lines:
634        line = line.strip()
635        if 'PIXEL' in line:
636            values = line.split()
637            pixel = (int(values[2]),int(values[4]))     #in microns
638        elif 'WAVELENGTH' in line:
639            wave = float(line.split()[1])
640        elif 'DISTANCE' in line:
641            distance = float(line.split()[1])           #in mm
642        elif 'CENTER' in line:
643            values = line.split()
644            center = [float(values[2])/10.,float(values[4])/10.]    #make in mm from pixels
645        if line: 
646            head.append(line)
647    data = {'pixelSize':pixel,'wavelength':wave,'distance':distance,'center':center}
648    for line in head:
649        if 'FORMAT' in line[0:6]:
650            items = line.split()
651            size = int(items[1])
652            Npix = size*size
653    pos = 4096
654    data['size'] = [size,size]
655    File.seek(pos)
656    line = File.read(8)
657    while 'CCP4' not in line:       #get past overflow list for now
658        line = File.read(8)
659        pos += 8
660    pos += 37
661    File.seek(pos)
662    raw = File.read()
663    File.close()
664    image = np.zeros(shape=(size,size),dtype=np.int32)
665    image = pf.pack_f(len(raw),raw,size,image)
666    if imageOnly:
667        return image.T              #transpose to get it right way around
668    else:
669        return head,data,Npix,image.T
670       
671def GetTifData(filename,imageOnly=False):
672    import struct as st
673    import array as ar
674    File = open(filename,'rb')
675    dataType = 5
676    try:
677        Meta = open(filename+'.metadata','Ur')
678        head = Meta.readlines()
679        for line in head:
680            line = line.strip()
681            if 'dataType' in line:
682                dataType = int(line.split('=')[1])
683        Meta.close()
684    except IOError:
685        print 'no metadata file found - will try to read file anyway'
686        head = ['no metadata file found',]
687       
688    tag = File.read(2)
689    byteOrd = '<'
690    if tag == 'II' and int(st.unpack('<h',File.read(2))[0]) == 42:     #little endian
691        IFD = int(st.unpack(byteOrd+'i',File.read(4))[0])
692    elif tag == 'MM' and int(st.unpack('>h',File.read(2))[0]) == 42:   #big endian
693        byteOrd = '>'
694        IFD = int(st.unpack(byteOrd+'i',File.read(4))[0])       
695    else:
696        lines = ['not a detector tiff file',]
697        return lines,0,0,0
698    File.seek(IFD)                                                  #get number of directory entries
699    NED = int(st.unpack(byteOrd+'h',File.read(2))[0])
700    print IFD,NED
701    IFD = {}
702    for ied in range(NED):
703        Tag,Type = st.unpack(byteOrd+'Hh',File.read(4))
704        nVal = st.unpack(byteOrd+'i',File.read(4))[0]
705        if Type == 1:
706            Value = st.unpack(byteOrd+nVal*'b',File.read(nVal))
707        elif Type == 2:
708            Value = st.unpack(byteOrd+'i',File.read(4))
709        elif Type == 3:
710            Value = st.unpack(byteOrd+nVal*'h',File.read(nVal*2))
711            x = st.unpack(byteOrd+nVal*'h',File.read(nVal*2))
712        elif Type == 4:
713            Value = st.unpack(byteOrd+nVal*'i',File.read(nVal*4))
714        elif Type == 5:
715            Value = st.unpack(byteOrd+nVal*'i',File.read(nVal*4))
716        elif Type == 11:
717            Value = st.unpack(byteOrd+nVal*'f',File.read(nVal*4))
718        IFD[Tag] = [Type,nVal,Value]
719#    for key in IFD:
720#        print key,IFD[key]
721    sizexy = [IFD[256][2][0],IFD[257][2][0]]
722    [nx,ny] = sizexy
723    Npix = nx*ny
724    if 272 in IFD:
725        ifd = IFD[272]
726        File.seek(ifd[2][0])
727        S = File.read(ifd[1])
728        if 'PILATUS' in S:
729            tifType = 'Pilatus'
730            dataType = 0
731            pixy = (172,172)
732            File.seek(4096)
733            if not imageOnly:
734                print 'Read Pilatus tiff file: ',filename
735            image = ar.array('L',File.read(4*Npix))
736            image = np.array(np.asarray(image),dtype=np.int32)
737    elif 262 in IFD and IFD[262][2][0] > 4:
738        tifType = 'DND'
739        pixy = (158,158)
740        File.seek(512)
741        if not imageOnly:
742            print 'Read DND SAX/WAX-detector tiff file: ',filename
743        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
744    elif sizexy == [1536,1536]:
745        tifType = 'APS Gold'
746        pixy = (150,150)
747        File.seek(64)
748        if not imageOnly:
749            print 'Read Gold tiff file:',filename
750        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
751    elif sizexy == [2048,2048] or sizexy == [1024,1024]:
752        if IFD[273][2][0] == 8:
753            if IFD[258][2][0] == 32:
754                tifType = 'PE'
755                pixy = (200,200)
756                File.seek(8)
757                if not imageOnly:
758                    print 'Read APS PE-detector tiff file: ',filename
759                if dataType == 5:
760                    image = np.array(ar.array('f',File.read(4*Npix)),dtype=np.float32)
761                else:
762                    image = np.array(ar.array('I',File.read(4*Npix)),dtype=np.int32)
763        elif IFD[273][2][0] == 4096:
764            tifType = 'MAR'
765            pixy = (158,158)
766            File.seek(4096)
767            if not imageOnly:
768                print 'Read MAR CCD tiff file: ',filename
769            image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
770#    elif sizexy == [960,960]:
771#        tiftype = 'PE-BE'
772#        pixy = (200,200)
773#        File.seek(8)
774#        if not imageOnly:
775#            print 'Read Gold tiff file:',filename
776#        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
777           
778    else:
779        lines = ['not a known detector tiff file',]
780        return lines,0,0,0
781       
782    image = np.reshape(image,(sizexy[1],sizexy[0]))
783    center = [pixy[0]*sizexy[0]/2000,pixy[1]*sizexy[1]/2000]
784    data = {'pixelSize':pixy,'wavelength':0.10,'distance':100.0,'center':center,'size':sizexy}
785    File.close()   
786    if imageOnly:
787        return image
788    else:
789        return head,data,Npix,image
790   
791def ProjFileOpen(self):
792    file = open(self.GSASprojectfile,'rb')
793    print 'load from file: ',self.GSASprojectfile
794    wx.BeginBusyCursor()
795    try:
796        while True:
797            try:
798                data = cPickle.load(file)
799            except EOFError:
800                break
801            datum = data[0]
802            print 'load: ',datum[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                print '    load: ',datus[0]
811                sub = self.PatternTree.AppendItem(Id,datus[0])
812                self.PatternTree.SetItemPyData(sub,datus[1])
813            if 'IMG' in datum[0]:                   #retreive image default flag & data if set
814                Data = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Image Controls'))
815                if Data['setDefault']:
816                    self.imageDefault = Data               
817        file.close()
818       
819    finally:
820        wx.EndBusyCursor()
821    print 'project load successful'
822    self.NewPlot = True
823   
824def ProjFileSave(self):
825    if not self.PatternTree.IsEmpty():
826        file = open(self.GSASprojectfile,'wb')
827        print 'save to file: ',self.GSASprojectfile
828        wx.BeginBusyCursor()
829        try:
830            item, cookie = self.PatternTree.GetFirstChild(self.root)
831            while item:
832                data = []
833                name = self.PatternTree.GetItemText(item)
834                print 'save: ',name
835                data.append([name,self.PatternTree.GetItemPyData(item)])
836                item2, cookie2 = self.PatternTree.GetFirstChild(item)
837                while item2:
838                    name = self.PatternTree.GetItemText(item2)
839                    print '    save: ',name
840                    data.append([name,self.PatternTree.GetItemPyData(item2)])
841                    item2, cookie2 = self.PatternTree.GetNextChild(item, cookie2)                           
842                item, cookie = self.PatternTree.GetNextChild(self.root, cookie)                           
843                cPickle.dump(data,file,1)
844            file.close()
845        finally:
846            wx.EndBusyCursor()
847        print 'project save successful'
848       
849def SaveIntegration(self,PickId,data):
850    azms = self.Integrate[1][:-1]
851    X = self.Integrate[2][:-1]
852    Xminmax = [X[0],X[-1]]
853    N = len(X)
854    Id = self.PatternTree.GetItemParent(PickId)
855    name = self.PatternTree.GetItemText(Id)
856    name = name.replace('IMG ','PWDR ')
857    Comments = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Comments'))
858    names = ['Type','Lam','Zero','Polariz.','U','V','W','X','Y','SH/L','Azimuth'] 
859    codes = [0 for i in range(11)]
860    LRazm = data['LRazimuth']
861    if data['fullIntegrate'] and data['outAzimuths'] == 1:
862        Azms = [45.0,]                              #a poor man's average?
863    for i,azm in enumerate(azms):
864        item, cookie = self.PatternTree.GetFirstChild(self.root)
865        Id = 0
866        while item:
867            Name = self.PatternTree.GetItemText(item)
868            if name == Name:
869                Id = item
870            item, cookie = self.PatternTree.GetNextChild(self.root, cookie)
871        parms = ['PXC',data['wavelength'],0.0,0.99,1.0,-0.10,0.4,0.30,1.0,0.0001,azm]    #set polarization for synchrotron radiation!
872        Y = self.Integrate[0][i]
873        W = 1./Y                    #probably not true
874        Sample = {'Scale':[1.0,True],'Type':'Debye-Scherrer','Absorption':[0.0,False],'DisplaceX':[0.0,False],
875            'DisplaceY':[0.0,False],'Diffuse':[],'Temperature':300.,'Pressure':1.0,'Humidity':0.0,'Voltage':0.0,'Force':0.0}
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"%(azm))
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.