source: trunk/GSASIIIO.py @ 450

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