source: trunk/GSASIIIO.py @ 303

Last change on this file since 303 was 303, checked in by vondreele, 12 years ago

Allow reading of xye (Topas style) files
Begin implementation of spherical harmonics texture
Refactor indexing; replace cell refinement from peak positions

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