source: trunk/GSASIIIO.py @ 400

Last change on this file since 400 was 400, checked in by vondreele, 11 years ago

new default for Vcov contour plot - RdYlGn?
faster cleanup on changing/reloading projects
cleanup data delete
implement sample parameter copy
improve Vcov plotting routine
implement plot of vcov from seq refinements

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