source: trunk/GSASIIIO.py @ 342

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

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