source: trunk/GSASIIIO.py @ 430

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

NB: There are a few stray print statements in this version
Some work on help stuff in GSASIIgrid.py
change default sum for constraints to 1.0
get data plots to include size & preferred orientation
modify plot toolbar - eliminate subplot button & add a Help button; doesn't do anything except recognize what plot it is from
Edited gsasII.html with Word - probably messed up all the bookmarks - but added a lot of text

  • 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-11-30 23:02:08 +0000 (Wed, 30 Nov 2011) $
6# $Author: vondreele $
7# $Revision: 430 $
8# $URL: trunk/GSASIIIO.py $
9# $Id: GSASIIIO.py 430 2011-11-30 23:02:08Z 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    print cons
442    Nch = int(cons[2])
443    if DataType[2] == 'C':
444        start = float(cons[5])/100.0               #CW: from centidegrees to degrees
445        step = float(cons[6])/100.0
446    elif DataType[2] == 'T':
447        start = float(cons[5])/1000.0              #TOF: from musec to millisec - not likely!
448        step = float(cons[6])/1000.0
449    File.seek(Pos)
450    x = []
451    y = []
452    w = []
453    S = File.readline()
454    j = 0
455    while S and S[:4] != 'BANK':
456        for i in range(0,80,8):
457            xi = start+step*j
458            ni = max(sint(S[i:i+2]),1)
459            yi = max(sfloat(S[i+2:i+8]),0.0)
460            if yi:
461                vi = yi/ni
462            else:
463                yi = 0.0
464                vi = 1.0
465            j += 1
466            if j < Nch:
467                x.append(xi)
468                y.append(yi)
469                w.append(1.0/vi)
470            print '%8.3f %8.0f %d %8g'%(xi,yi,ni,1./vi)
471        S = File.readline()
472    File.close()
473    N = len(x)
474    print N
475    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
476   
477def CheckImageFile(self,imagefile):
478    if not ospath.exists(imagefile):
479        dlg = wx.FileDialog(self, 'Bad image file name; choose name', '.', '',\
480        'Any image file (*.tif;*.tiff;*.mar*;*.avg;*.sum;*.img)\
481        |*.tif;*.tiff;*.mar*;*.avg;*.sum;*.img|\
482        Any detector tif (*.tif;*.tiff)|*.tif;*.tiff|\
483        MAR file (*.mar*)|*.mar*|\
484        GE Image (*.avg;*.sum)|*.avg;*.sum|\
485        ADSC Image (*.img)|*.img|\
486        All files (*.*)|*.*',wx.OPEN)
487        if self.dirname:
488            dlg.SetDirectory(self.dirname)
489        try:
490            dlg.SetFilename(ospath.split(imagefile)[1])
491            if dlg.ShowModal() == wx.ID_OK:
492                self.dirname = dlg.GetDirectory()
493                imagefile = dlg.GetPath()
494            else:
495                imagefile = False
496        finally:
497            dlg.Destroy()
498    return imagefile
499       
500def GetImageData(self,imagefile,imageOnly=False):       
501    ext = ospath.splitext(imagefile)[1]
502    Comments = []
503    if ext == '.tif' or ext == '.tiff':
504        Comments,Data,Npix,Image = GetTifData(imagefile)
505    elif ext == '.img':
506        Comments,Data,Npix,Image = GetImgData(imagefile)
507        Image[0][0] = 0
508    elif ext == '.mar3450' or ext == '.mar2300':
509        Comments,Data,Npix,Image = GetMAR345Data(imagefile)
510    elif ext in ['.sum','.avg','']:
511        Comments,Data,Npix,Image = GetGEsumData(imagefile)
512    elif ext == '.G2img':
513        Comments,Data,Npix,Image = GetG2Image(imagefile)
514    if imageOnly:
515        return Image
516    else:
517        return Comments,Data,Npix,Image
518       
519def PutG2Image(filename,Comments,Data,Npix,image):
520    File = open(filename,'wb')
521    cPickle.dump([Comments,Data,Npix,image],File,1)
522    File.close()
523    return
524   
525def GetG2Image(filename):
526    File = open(filename,'rb')
527    Comments,Data,Npix,image = cPickle.load(File)
528    File.close()
529    return Comments,Data,Npix,image
530   
531def GetGEsumData(filename,imageOnly=False):
532    import struct as st
533    import array as ar
534    if not imageOnly:
535        print 'Read GE sum file: ',filename   
536    File = open(filename,'rb')
537    if '.sum' in filename:
538        head = ['GE detector sum data from APS 1-ID',]
539        sizexy = [2048,2048]
540    elif '.avg' in filename:
541        head = ['GE detector avg data from APS 1-ID',]
542        sizexy = [2048,2048]
543    else:
544        head = ['GE detector raw data from APS 1-ID',]
545        File.seek(18)
546        size,nframes = st.unpack('<ih',File.read(6))
547        sizexy = [2048,2048]
548        pos = 8192
549        File.seek(pos)
550    Npix = sizexy[0]*sizexy[1]
551    if '.sum' in filename:
552        image = np.array(ar.array('f',File.read(4*Npix)),dtype=np.int32)
553    elif '.avg' in filename:
554        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
555    else:
556        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
557        while nframes > 1:
558            image += np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
559            nframes -= 1
560    image = np.reshape(image,(sizexy[1],sizexy[0]))
561    data = {'pixelSize':(200,200),'wavelength':0.15,'distance':250.0,'center':[204.8,204.8],'size':sizexy} 
562    File.close()   
563    if imageOnly:
564        return image
565    else:
566        return head,data,Npix,image
567       
568def GetImgData(filename,imageOnly=False):
569    import struct as st
570    import array as ar
571    if not imageOnly:
572        print 'Read ADSC img file: ',filename
573    File = open(filename,'rb')
574    head = File.read(511)
575    lines = head.split('\n')
576    head = []
577    center = [0,0]
578    for line in lines[1:-2]:
579        line = line.strip()[:-1]
580        if line:
581            if 'SIZE1' in line:
582                size = int(line.split('=')[1])
583                Npix = size*size
584            elif 'WAVELENGTH' in line:
585                wave = float(line.split('=')[1])
586            elif 'BIN' in line:
587                if line.split('=')[1] == '2x2':
588                    pixel=(102,102)
589                else:
590                    pixel = (51,51)
591            elif 'DISTANCE' in line:
592                distance = float(line.split('=')[1])
593            elif 'CENTER_X' in line:
594                center[0] = float(line.split('=')[1])
595            elif 'CENTER_Y' in line:
596                center[1] = float(line.split('=')[1])
597            head.append(line)
598    data = {'pixelSize':pixel,'wavelength':wave,'distance':distance,'center':center,'size':[size,size]}
599    image = []
600    row = 0
601    pos = 512
602    File.seek(pos)
603    image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
604    image = np.reshape(image,(sizexy[1],sizexy[0]))
605#    image = np.zeros(shape=(size,size),dtype=np.int32)   
606#    while row < size:
607#        File.seek(pos)
608#        line = ar.array('H',File.read(2*size))
609#        image[row] = np.asarray(line)
610#        row += 1
611#        pos += 2*size
612    File.close()
613    if imageOnly:
614        return image
615    else:
616        return lines[1:-2],data,Npix,image
617       
618def GetMAR345Data(filename,imageOnly=False):
619    import array as ar
620    import struct as st
621    try:
622        import pack_f as pf
623    except:
624        msg = wx.MessageDialog(None, message="Unable to load the GSAS MAR image decompression, pack_f",
625                               caption="Import Error",
626                               style=wx.ICON_ERROR | wx.OK | wx.STAY_ON_TOP)
627        msg.ShowModal()
628        return None,None,None,None
629
630    if not imageOnly:
631        print 'Read Mar345 file: ',filename
632    File = open(filename,'rb')
633    head = File.read(4095)
634    numbers = st.unpack('<iiiiiiiiii',head[:40])
635    lines = head[128:].split('\n')
636    head = []
637    for line in lines:
638        line = line.strip()
639        if 'PIXEL' in line:
640            values = line.split()
641            pixel = (int(values[2]),int(values[4]))     #in microns
642        elif 'WAVELENGTH' in line:
643            wave = float(line.split()[1])
644        elif 'DISTANCE' in line:
645            distance = float(line.split()[1])           #in mm
646        elif 'CENTER' in line:
647            values = line.split()
648            center = [float(values[2])/10.,float(values[4])/10.]    #make in mm from pixels
649        if line: 
650            head.append(line)
651    data = {'pixelSize':pixel,'wavelength':wave,'distance':distance,'center':center}
652    for line in head:
653        if 'FORMAT' in line[0:6]:
654            items = line.split()
655            size = int(items[1])
656            Npix = size*size
657    pos = 4096
658    data['size'] = [size,size]
659    File.seek(pos)
660    line = File.read(8)
661    while 'CCP4' not in line:       #get past overflow list for now
662        line = File.read(8)
663        pos += 8
664    pos += 37
665    File.seek(pos)
666    raw = File.read()
667    File.close()
668    image = np.zeros(shape=(size,size),dtype=np.int32)
669    image = pf.pack_f(len(raw),raw,size,image)
670    if imageOnly:
671        return image.T              #transpose to get it right way around
672    else:
673        return head,data,Npix,image.T
674       
675def GetTifData(filename,imageOnly=False):
676    import struct as st
677    import array as ar
678    File = open(filename,'rb')
679    dataType = 5
680    try:
681        Meta = open(filename+'.metadata','Ur')
682        head = Meta.readlines()
683        for line in head:
684            line = line.strip()
685            if 'dataType' in line:
686                dataType = int(line.split('=')[1])
687        Meta.close()
688    except IOError:
689        print 'no metadata file found - will try to read file anyway'
690        head = ['no metadata file found',]
691       
692    tag = File.read(2)
693    byteOrd = '<'
694    if tag == 'II' and int(st.unpack('<h',File.read(2))[0]) == 42:     #little endian
695        IFD = int(st.unpack(byteOrd+'i',File.read(4))[0])
696    elif tag == 'MM' and int(st.unpack('>h',File.read(2))[0]) == 42:   #big endian
697        byteOrd = '>'
698        IFD = int(st.unpack(byteOrd+'i',File.read(4))[0])       
699    else:
700        lines = ['not a detector tiff file',]
701        return lines,0,0,0
702    File.seek(IFD)                                                  #get number of directory entries
703    NED = int(st.unpack(byteOrd+'h',File.read(2))[0])
704    IFD = {}
705    for ied in range(NED):
706        Tag,Type = st.unpack(byteOrd+'Hh',File.read(4))
707        nVal = st.unpack(byteOrd+'i',File.read(4))[0]
708        if Type == 1:
709            Value = st.unpack(byteOrd+nVal*'b',File.read(nVal))
710        elif Type == 2:
711            Value = st.unpack(byteOrd+'i',File.read(4))
712        elif Type == 3:
713            Value = st.unpack(byteOrd+nVal*'h',File.read(nVal*2))
714            x = st.unpack(byteOrd+nVal*'h',File.read(nVal*2))
715        elif Type == 4:
716            Value = st.unpack(byteOrd+nVal*'i',File.read(nVal*4))
717        elif Type == 5:
718            Value = st.unpack(byteOrd+nVal*'i',File.read(nVal*4))
719        elif Type == 11:
720            Value = st.unpack(byteOrd+nVal*'f',File.read(nVal*4))
721        IFD[Tag] = [Type,nVal,Value]
722#    for key in IFD:
723#        print key,IFD[key]
724    sizexy = [IFD[256][2][0],IFD[257][2][0]]
725    [nx,ny] = sizexy
726    Npix = nx*ny
727    if 272 in IFD:
728        ifd = IFD[272]
729        File.seek(ifd[2][0])
730        S = File.read(ifd[1])
731        if 'PILATUS' in S:
732            tifType = 'Pilatus'
733            dataType = 0
734            pixy = (172,172)
735            File.seek(4096)
736            if not imageOnly:
737                print 'Read Pilatus tiff file: ',filename
738            image = ar.array('L',File.read(4*Npix))
739            image = np.array(np.asarray(image),dtype=np.int32)
740    elif 262 in IFD and IFD[262][2][0] > 4:
741        tifType = 'DND'
742        pixy = (158,158)
743        File.seek(512)
744        if not imageOnly:
745            print 'Read DND SAX/WAX-detector tiff file: ',filename
746        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
747    elif sizexy == [1536,1536]:
748        tifType = 'APS Gold'
749        pixy = (150,150)
750        File.seek(64)
751        if not imageOnly:
752            print 'Read Gold tiff file:',filename
753        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
754    elif sizexy == [2048,2048] or sizexy == [1024,1024]:
755        if IFD[273][2][0] == 8:
756            if IFD[258][2][0] == 32:
757                tifType = 'PE'
758                pixy = (200,200)
759                File.seek(8)
760                if not imageOnly:
761                    print 'Read APS PE-detector tiff file: ',filename
762                if dataType == 5:
763                    image = np.array(ar.array('f',File.read(4*Npix)),dtype=np.float32)
764                else:
765                    image = np.array(ar.array('I',File.read(4*Npix)),dtype=np.int32)
766        elif IFD[273][2][0] == 4096:
767            tifType = 'MAR'
768            pixy = (158,158)
769            File.seek(4096)
770            if not imageOnly:
771                print 'Read MAR CCD tiff file: ',filename
772            image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
773#    elif sizexy == [960,960]:
774#        tiftype = 'PE-BE'
775#        pixy = (200,200)
776#        File.seek(8)
777#        if not imageOnly:
778#            print 'Read Gold tiff file:',filename
779#        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
780           
781    else:
782        lines = ['not a known detector tiff file',]
783        return lines,0,0,0
784       
785    image = np.reshape(image,(sizexy[1],sizexy[0]))
786    center = [pixy[0]*sizexy[0]/2000,pixy[1]*sizexy[1]/2000]
787    data = {'pixelSize':pixy,'wavelength':0.10,'distance':100.0,'center':center,'size':sizexy}
788    File.close()   
789    if imageOnly:
790        return image
791    else:
792        return head,data,Npix,image
793   
794def ProjFileOpen(self):
795    file = open(self.GSASprojectfile,'rb')
796    print 'load from file: ',self.GSASprojectfile
797    wx.BeginBusyCursor()
798    try:
799        while True:
800            try:
801                data = cPickle.load(file)
802            except EOFError:
803                break
804            datum = data[0]
805           
806            Id = self.PatternTree.AppendItem(parent=self.root,text=datum[0])
807            if 'PWDR' in datum[0]:               
808                self.PatternTree.SetItemPyData(Id,datum[1][:3])     #temp. trim off junk
809            else:
810                self.PatternTree.SetItemPyData(Id,datum[1])
811            for datus in data[1:]:
812                sub = self.PatternTree.AppendItem(Id,datus[0])
813                self.PatternTree.SetItemPyData(sub,datus[1])
814            if 'IMG' in datum[0]:                   #retreive image default flag & data if set
815                Data = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Image Controls'))
816                if Data['setDefault']:
817                    self.imageDefault = Data               
818        file.close()
819       
820    finally:
821        wx.EndBusyCursor()
822    print 'project load successful'
823    self.NewPlot = True
824   
825def ProjFileSave(self):
826    if not self.PatternTree.IsEmpty():
827        file = open(self.GSASprojectfile,'wb')
828        print 'save to file: ',self.GSASprojectfile
829        wx.BeginBusyCursor()
830        try:
831            item, cookie = self.PatternTree.GetFirstChild(self.root)
832            while item:
833                data = []
834                name = self.PatternTree.GetItemText(item)
835                data.append([name,self.PatternTree.GetItemPyData(item)])
836                item2, cookie2 = self.PatternTree.GetFirstChild(item)
837                while item2:
838                    name = self.PatternTree.GetItemText(item2)
839                    data.append([name,self.PatternTree.GetItemPyData(item2)])
840                    item2, cookie2 = self.PatternTree.GetNextChild(item, cookie2)                           
841                item, cookie = self.PatternTree.GetNextChild(self.root, cookie)                           
842                cPickle.dump(data,file,1)
843            file.close()
844        finally:
845            wx.EndBusyCursor()
846        print 'project save successful'
847       
848def SaveIntegration(self,PickId,data):
849    azms = self.Integrate[1]
850    X = self.Integrate[2][:-1]
851    Xminmax = [X[0],X[-1]]
852    N = len(X)
853    Id = self.PatternTree.GetItemParent(PickId)
854    name = self.PatternTree.GetItemText(Id)
855    name = name.replace('IMG ','PWDR ')
856    Comments = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Comments'))
857    names = ['Type','Lam','Zero','Polariz.','U','V','W','X','Y','SH/L','Azimuth'] 
858    codes = [0 for i in range(11)]
859    LRazm = data['LRazimuth']
860    Azms = []
861    if data['fullIntegrate'] and data['outAzimuths'] == 1:
862        Azms = [45.0,]                              #a poor man's average?
863    else:
864        for i,azm in enumerate(azms[:-1]):
865            Azms.append((azms[i+1]+azm)/2.)
866    for i,azm in enumerate(azms[:-1]):
867        item, cookie = self.PatternTree.GetFirstChild(self.root)
868        Id = 0
869        while item:
870            Name = self.PatternTree.GetItemText(item)
871            if name == Name:
872                Id = item
873            item, cookie = self.PatternTree.GetNextChild(self.root, cookie)
874        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!
875        Y = self.Integrate[0][i]
876        W = 1./Y                    #probably not true
877        Sample = G2pdG.SetDefaultSample()
878        if Id:
879            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Comments'),Comments)                   
880            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Limits'),[tuple(Xminmax),Xminmax])
881            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Background'),[['chebyschev',1,3,1.0,0.0,0.0]])
882            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Instrument Parameters'),[tuple(parms),parms[:],codes,names])
883            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Peak List'),[])
884            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Index Peak List'),[])
885            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Unit Cells List'),[])             
886            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Reflection Lists'),{})             
887        else:
888            Id = self.PatternTree.AppendItem(parent=self.root,text=name+" Azm= %.2f"%(Azms[i]))
889            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Comments'),Comments)                   
890            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Limits'),[tuple(Xminmax),Xminmax])
891            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Background'),[['chebyschev',1,3,1.0,0.0,0.0]])
892            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Instrument Parameters'),[tuple(parms),parms[:],codes,names])
893            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Sample Parameters'),Sample)
894            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Peak List'),[])
895            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Index Peak List'),[])
896            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Unit Cells List'),[])
897            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Reflection Lists'),{})             
898        self.PatternTree.SetItemPyData(Id,[[''],[np.array(X),np.array(Y),np.array(W),np.zeros(N),np.zeros(N),np.zeros(N)]])
899    self.PatternTree.SelectItem(Id)
900    self.PatternTree.Expand(Id)
901    self.PatternId = Id
902           
903def powderFxyeSave(self,exports,powderfile):
904    head,tail = ospath.split(powderfile)
905    name,ext = tail.split('.')
906    wx.BeginBusyCursor()
907    for i,export in enumerate(exports):
908        filename = ospath.join(head,name+'-%03d.'%(i)+ext)
909        prmname = filename.strip(ext)+'prm'
910        prm = open(prmname,'w')      #old style GSAS parm file
911        PickId = G2gd.GetPatternTreeItemId(self, self.root, export)
912        Values,Names = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self, \
913            PickId, 'Instrument Parameters'))[1::2]     #get values & names
914        Inst = dict(zip(Names,Values))
915        print Inst['Type']
916        prm.write( '            123456789012345678901234567890123456789012345678901234567890        '+'\n')
917        prm.write( 'INS   BANK      1                                                               '+'\n')
918        prm.write(('INS   HTYPE   %sR                                                              '+'\n')%(Inst['Type']))
919        if 'Lam1' in Inst:              #Ka1 & Ka2
920            prm.write(('INS  1 ICONS%10.7f%10.7f    0.0000               0.990    0     0.500   '+'\n')%(Inst['Lam1'],Inst['Lam2']))
921        elif 'Lam' in Inst:             #single wavelength
922            prm.write(('INS  1 ICONS%10.7f%10.7f    0.0000               0.990    0     0.500   '+'\n')%(Inst['Lam'],0.0))
923        prm.write( 'INS  1 IRAD     0                                                               '+'\n')
924        prm.write( 'INS  1I HEAD                                                                    '+'\n')
925        prm.write( 'INS  1I ITYP    0    0.0000  180.0000         1                                 '+'\n')
926        prm.write(('INS  1DETAZM%10.3f                                                          '+'\n')%(Inst['Azimuth']))
927        prm.write( 'INS  1PRCF1     3    8   0.00100                                                '+'\n')
928        prm.write(('INS  1PRCF11     %15.6g%15.6g%15.6g%15.6g   '+'\n')%(Inst['U'],Inst['V'],Inst['W'],0.0))
929        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.))
930        prm.close()
931        file = open(filename,'w')
932        print 'save powder pattern to file: ',filename
933        try:
934            x,y,w,yc,yb,yd = self.PatternTree.GetItemPyData(PickId)[1]
935            file.write(powderfile+'\n')
936            file.write('BANK 1 %d %d CONS %.2f %.2f 0 0 FXYE\n'%(len(x),len(x),\
937                100.*x[0],100.*(x[1]-x[0])))
938            s = list(np.sqrt(1./np.array(w)))       
939            XYW = zip(x,y,s)
940            for X,Y,S in XYW:
941                file.write("%15.6g %15.6g %15.6g\n" % (100.*X,Y,max(S,1.0)))
942            file.close()
943        finally:
944            wx.EndBusyCursor()
945        print 'powder pattern file written'
946       
947def powderXyeSave(self,exports,powderfile):
948    head,tail = ospath.split(powderfile)
949    name,ext = tail.split('.')
950    for i,export in enumerate(exports):
951        filename = ospath.join(head,name+'-%03d.'%(i)+ext)
952        PickId = G2gd.GetPatternTreeItemId(self, self.root, export)
953        file = open(filename,'w')
954        file.write('#%s\n'%(export))
955        print 'save powder pattern to file: ',filename
956        wx.BeginBusyCursor()
957        try:
958            x,y,w,yc,yb,yd = self.PatternTree.GetItemPyData(PickId)[1]
959            s = list(np.sqrt(1./np.array(w)))       
960            XYW = zip(x,y,s)
961            for X,Y,W in XYW:
962                file.write("%15.6g %15.6g %15.6g\n" % (X,Y,W))
963            file.close()
964        finally:
965            wx.EndBusyCursor()
966        print 'powder pattern file written'
967       
968def PDFSave(self,exports):   
969    for export in exports:
970        PickId = G2gd.GetPatternTreeItemId(self, self.root, export)
971        SQname = 'S(Q)'+export[4:]
972        GRname = 'G(R)'+export[4:]
973        sqfilename = ospath.join(self.dirname,export.replace(' ','_')[5:]+'.sq')
974        grfilename = ospath.join(self.dirname,export.replace(' ','_')[5:]+'.gr')
975        sqId = G2gd.GetPatternTreeItemId(self, PickId, SQname)
976        grId = G2gd.GetPatternTreeItemId(self, PickId, GRname)
977        sqdata = np.array(self.PatternTree.GetItemPyData(sqId)[1][:2]).T
978        grdata = np.array(self.PatternTree.GetItemPyData(grId)[1][:2]).T
979        sqfile = open(sqfilename,'w')
980        grfile = open(grfilename,'w')
981        sqfile.write('#T S(Q) %s\n'%(export))
982        grfile.write('#T G(R) %s\n'%(export))
983        sqfile.write('#L Q     S(Q)\n')
984        grfile.write('#L R     G(R)\n')
985        for q,sq in sqdata:
986            sqfile.write("%15.6g %15.6g\n" % (q,sq))
987        sqfile.close()
988        for r,gr in grdata:
989            grfile.write("%15.6g %15.6g\n" % (r,gr))
990        grfile.close()
991   
992def PeakListSave(self,file,peaks):
993    print 'save peak list to file: ',self.peaklistfile
994    if not peaks:
995        dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
996        try:
997            result = dlg.ShowModal()
998        finally:
999            dlg.Destroy()
1000        return
1001    for peak in peaks:
1002        file.write("%10.4f %12.2f %10.3f %10.3f \n" % \
1003            (peak[0],peak[2],peak[4],peak[6]))
1004    print 'peak list saved'
1005             
1006def IndexPeakListSave(self,peaks):
1007    file = open(self.peaklistfile,'wa')
1008    print 'save index peak list to file: ',self.peaklistfile
1009    wx.BeginBusyCursor()
1010    try:
1011        if not peaks:
1012            dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
1013            try:
1014                result = dlg.ShowModal()
1015            finally:
1016                dlg.Destroy()
1017            return
1018        for peak in peaks:
1019            file.write("%12.6f\n" % (peak[7]))
1020        file.close()
1021    finally:
1022        wx.EndBusyCursor()
1023    print 'index peak list saved'
1024   
1025def ReadEXPPhase(self,filename):
1026    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1027    textureData = {'Order':0,'Model':'cylindrical','Sample omega':[False,0.0],
1028        'Sample chi':[False,0.0],'Sample phi':[False,0.0],'SH Coeff':[False,{}],
1029        'SHShow':False,'PFhkl':[0,0,1],'PFxyz':[0,0,1],'PlotType':'Pole figure'}
1030    shNcof = 0
1031    file = open(filename, 'Ur')
1032    S = 1
1033    Expr = [{},{},{},{},{},{},{},{},{}]
1034    while S:
1035        S = file.readline()
1036        if 'EXPR NPHAS' in S[:12]:
1037            Num = S[12:-1].count('0')
1038            NPhas = S[12:-1].split()
1039        if 'CRS' in S[:3]:
1040            N = int(S[3:4])-1
1041            Expr[N][S[:12]] = S[12:-1]
1042    file.close()
1043    PNames = []
1044    for n,N in enumerate(NPhas):
1045        if N != '0':
1046            result = n
1047            key = 'CRS'+str(n+1)+'    PNAM'
1048            PNames.append(Expr[n][key])
1049    if Num < 8:
1050        dlg = wx.SingleChoiceDialog(self, 'Which phase to read?', 'Read phase data', PNames, wx.CHOICEDLG_STYLE)
1051        try:
1052            if dlg.ShowModal() == wx.ID_OK:
1053                result = dlg.GetSelection()
1054        finally:
1055            dlg.Destroy()       
1056    EXPphase = Expr[result]
1057    keyList = EXPphase.keys()
1058    keyList.sort()
1059    SGData = {}
1060    if NPhas[result] == '1':
1061        Ptype = 'nuclear'
1062    elif NPhas[result] in ['2','3']:
1063        Ptype = 'magnetic'
1064    elif NPhas[result] == '4':
1065        Ptype = 'macromolecular'
1066    elif NPhas[result] == '10':
1067        Ptype = 'Pawley'
1068    for key in keyList:
1069        if 'PNAM' in key:
1070           PhaseName = EXPphase[key].strip()
1071        elif 'ABC   ' in key:
1072            abc = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                       
1073        elif 'ANGLES' in key:
1074            angles = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                                               
1075        elif 'SG SYM' in key:
1076            SpGrp = EXPphase[key][:15].strip()
1077            E,SGData = G2spc.SpcGroup(SpGrp)
1078        elif 'OD    ' in key:
1079            SHdata = EXPphase[key].split()
1080            textureData['Order'] = int(SHdata[0])
1081            textureData['Model'] = shModels[int(SHdata[2])]
1082            textureData['Sample omega'] = [False,float(SHdata[6])]
1083            textureData['Sample chi'] = [False,float(SHdata[7])]
1084            textureData['Sample phi'] = [False,float(SHdata[8])]
1085            shNcof = int(SHdata[1])
1086    Atoms = []
1087    if Ptype == 'nuclear':
1088        for key in keyList:
1089            if 'AT' in key:
1090                if key[11:] == 'A':
1091                    S = EXPphase[key]
1092                elif key[11:] == 'B':
1093                    S += EXPphase[key]
1094                    Atom = [S[50:58].strip(),S[:10].strip(),'',
1095                        float(S[10:20]),float(S[20:30]),float(S[30:40]),
1096                        float(S[40:50]),'',int(S[60:62]),S[130:131]]
1097                    if Atom[9] == 'I':
1098                        Atom += [float(S[68:78]),0.,0.,0.,0.,0.,0.]
1099                    elif Atom[9] == 'A':
1100                        Atom += [0.0,float(S[68:78]),float(S[78:88]),
1101                            float(S[88:98]),float(S[98:108]),
1102                            float(S[108:118]),float(S[118:128])]
1103                    XYZ = Atom[3:6]
1104                    Atom[7],Atom[8] = G2spc.SytSym(XYZ,SGData)
1105                    Atom.append(ran.randint(0,sys.maxint))
1106                    Atoms.append(Atom)
1107    elif Ptype == 'macromolecular':
1108        for key in keyList:
1109            if 'AT' in key[6:8]:
1110                S = EXPphase[key]
1111                Atom = [int(S[56:60]),S[50:54].strip().upper(),S[54:56],
1112                    S[46:51].strip(),S[:8].strip(),'',
1113                    float(S[16:24]),float(S[24:32]),float(S[32:40]),
1114                    float(S[8:16]),'1',1,'I',float(S[40:46]),0,0,0,0,0,0]
1115                XYZ = Atom[6:9]
1116                Atom[10],Atom[11] = G2spc.SytSym(XYZ,SGData)
1117                Atom.append(ran.randint(0,sys.maxint))
1118                Atoms.append(Atom)
1119    Volume = G2lat.calc_V(G2lat.cell2A(abc+angles))
1120    if shNcof:
1121        shCoef = {}
1122        nRec = [i+1 for i in range((shNcof-1)/6+1)]
1123        for irec in nRec:
1124            ODkey = keyList[0][:6]+'OD'+'%3dA'%(irec)
1125            indx = EXPphase[ODkey].split()
1126            ODkey = ODkey[:-1]+'B'
1127            vals = EXPphase[ODkey].split()
1128            for i,val in enumerate(vals):
1129                key = 'C(%s,%s,%s)'%(indx[3*i],indx[3*i+1],indx[3*i+2])
1130                shCoef[key] = float(val)
1131        textureData['SH Coeff'] = [False,shCoef]
1132       
1133    Phase = {
1134            'General':{
1135                'Name':PhaseName,
1136                'Type':Ptype,
1137                'SGData':SGData,
1138                'Cell':[False,]+abc+angles+[Volume,],
1139                'Pawley dmin':1.0},
1140            'Atoms':Atoms,
1141            'Drawing':{},
1142            'Histograms':{},
1143            'Pawley ref':[],
1144            'Models':{},
1145            'SH Texture':textureData
1146            }
1147           
1148    return Phase
1149       
1150def ReadPDBPhase(filename):
1151    EightPiSq = 8.*math.pi**2
1152    file = open(filename, 'Ur')
1153    Phase = {}
1154    Title = ''
1155    Compnd = ''
1156    Atoms = []
1157    A = np.zeros(shape=(3,3))
1158    S = file.readline()
1159    while S:
1160        Atom = []
1161        if 'TITLE' in S[:5]:
1162            Title = S[10:72].strip()
1163            S = file.readline()
1164        elif 'COMPND    ' in S[:10]:
1165            Compnd = S[10:72].strip()
1166            S = file.readline()
1167        elif 'CRYST' in S[:5]:
1168            abc = S[7:34].split()
1169            angles = S[34:55].split()
1170            cell=[float(abc[0]),float(abc[1]),float(abc[2]),
1171                float(angles[0]),float(angles[1]),float(angles[2])]
1172            Volume = G2lat.calc_V(G2lat.cell2A(cell))
1173            AA,AB = G2lat.cell2AB(cell)
1174            SpGrp = S[55:65]
1175            E,SGData = G2spc.SpcGroup(SpGrp)
1176            if E: 
1177                print ' ERROR in space group symbol ',SpGrp,' in file ',filename
1178                print ' N.B.: make sure spaces separate axial fields in symbol' 
1179                print G2spc.SGErrors(E)
1180                return None
1181            SGlines = G2spc.SGPrint(SGData)
1182            for line in SGlines: print line
1183            S = file.readline()
1184        elif 'SCALE' in S[:5]:
1185            V = (S[10:41].split())
1186            A[int(S[5])-1] = [float(V[0]),float(V[1]),float(V[2])]
1187            S = file.readline()
1188        elif 'ATOM' in S[:4] or 'HETATM' in S[:6]:
1189            XYZ = [float(S[31:39]),float(S[39:47]),float(S[47:55])]
1190            XYZ = np.inner(AB,XYZ)
1191            SytSym,Mult = G2spc.SytSym(XYZ,SGData)
1192            Uiso = float(S[61:67])/EightPiSq
1193            Type = S[12:14].upper()
1194            if Type[0] in '123456789':
1195                Type = Type[1:]
1196            Atom = [S[22:27].strip(),S[17:20].upper(),S[20:22],
1197                S[12:17].strip(),Type.strip(),'',XYZ[0],XYZ[1],XYZ[2],
1198                float(S[55:61]),SytSym,Mult,'I',Uiso,0,0,0,0,0,0]
1199            S = file.readline()
1200            if 'ANISOU' in S[:6]:
1201                Uij = S[30:72].split()
1202                Uij = [float(Uij[0])/10000.,float(Uij[1])/10000.,float(Uij[2])/10000.,
1203                    float(Uij[3])/10000.,float(Uij[4])/10000.,float(Uij[5])/10000.]
1204                Atom = Atom[:14]+Uij
1205                Atom[12] = 'A'
1206                S = file.readline()
1207            Atom.append(ran.randint(0,sys.maxint))
1208            Atoms.append(Atom)
1209        else:           
1210            S = file.readline()
1211    file.close()
1212    if Title:
1213        PhaseName = Title
1214    elif Compnd:
1215        PhaseName = Compnd
1216    else:
1217        PhaseName = 'None'
1218    Phase['General'] = {'Name':PhaseName,'Type':'macromolecular','SGData':SGData,
1219        'Cell':[False,]+cell+[Volume,]}
1220    Phase['Atoms'] = Atoms
1221    Phase['Drawing'] = {}
1222    Phase['Histograms'] = {}
1223   
1224    return Phase
1225   
1226def ReadCIFPhase(filename):
1227    anisoNames = ['aniso_u_11','aniso_u_22','aniso_u_33','aniso_u_12','aniso_u_13','aniso_u_23']
1228    file = open(filename, 'Ur')
1229    Phase = {}
1230    Title = ospath.split(filename)[-1]
1231    print '\n Reading cif file: ',Title
1232    Compnd = ''
1233    Atoms = []
1234    A = np.zeros(shape=(3,3))
1235    S = file.readline()
1236    while S:
1237        if '_symmetry_space_group_name_H-M' in S:
1238            SpGrp = S.split("_symmetry_space_group_name_H-M")[1].strip().strip('"').strip("'")
1239            E,SGData = G2spc.SpcGroup(SpGrp)
1240            if E:
1241                print ' ERROR in space group symbol ',SpGrp,' in file ',filename
1242                print ' N.B.: make sure spaces separate axial fields in symbol' 
1243                print G2spc.SGErrors(E)
1244                return None
1245            S = file.readline()
1246        elif '_cell' in S:
1247            if '_cell_length_a' in S:
1248                a = S.split('_cell_length_a')[1].strip().strip('"').strip("'").split('(')[0]
1249            elif '_cell_length_b' in S:
1250                b = S.split('_cell_length_b')[1].strip().strip('"').strip("'").split('(')[0]
1251            elif '_cell_length_c' in S:
1252                c = S.split('_cell_length_c')[1].strip().strip('"').strip("'").split('(')[0]
1253            elif '_cell_angle_alpha' in S:
1254                alp = S.split('_cell_angle_alpha')[1].strip().strip('"').strip("'").split('(')[0]
1255            elif '_cell_angle_beta' in S:
1256                bet = S.split('_cell_angle_beta')[1].strip().strip('"').strip("'").split('(')[0]
1257            elif '_cell_angle_gamma' in S:
1258                gam = S.split('_cell_angle_gamma')[1].strip().strip('"').strip("'").split('(')[0]
1259            S = file.readline()
1260        elif 'loop_' in S:
1261            labels = {}
1262            i = 0
1263            while S:
1264                S = file.readline()
1265                if '_atom_site' in S.strip()[:10]:
1266                    labels[S.strip().split('_atom_site_')[1].lower()] = i
1267                    i += 1
1268                else:
1269                    break
1270            if labels:
1271                if 'aniso_label' not in labels:
1272                    while S:
1273                        atom = ['','','',0,0,0,1.0,'','','I',0.01,0,0,0,0,0,0]
1274                        S.strip()
1275                        if len(S.split()) != len(labels):
1276                            if 'loop_' in S:
1277                                break
1278                            S += file.readline().strip()
1279                        data = S.split()
1280                        if len(data) != len(labels):
1281                            break
1282                        for key in labels:
1283                            if key == 'type_symbol':
1284                                atom[1] = data[labels[key]]
1285                            elif key == 'label':
1286                                atom[0] = data[labels[key]]
1287                            elif key == 'fract_x':
1288                                atom[3] = float(data[labels[key]].split('(')[0])
1289                            elif key == 'fract_y':
1290                                atom[4] = float(data[labels[key]].split('(')[0])
1291                            elif key == 'fract_z':
1292                                atom[5] = float(data[labels[key]].split('(')[0])
1293                            elif key == 'occupancy':
1294                                atom[6] = float(data[labels[key]].split('(')[0])
1295                            elif key == 'thermal_displace_type':
1296                                if data[labels[key]].lower() == 'uiso':
1297                                    atom[9] = 'I'
1298                                    atom[10] = float(data[labels['u_iso_or_equiv']].split('(')[0])
1299                                else:
1300                                    atom[9] = 'A'
1301                                    atom[10] = 0.0
1302                                   
1303                        atom[7],atom[8] = G2spc.SytSym(atom[3:6],SGData)
1304                        atom.append(ran.randint(0,sys.maxint))
1305                        Atoms.append(atom)
1306                        S = file.readline()
1307                else:
1308                    while S:
1309                        S.strip()
1310                        data = S.split()
1311                        if len(data) != len(labels):
1312                            break
1313                        name = data[labels['aniso_label']]
1314                        for atom in Atoms:
1315                            if name == atom[0]:
1316                                for i,uname in enumerate(anisoNames):
1317                                    atom[i+11] = float(data[labels[uname]].split('(')[0])
1318                        S = file.readline()
1319                                                                       
1320        else:           
1321            S = file.readline()
1322    file.close()
1323    if Title:
1324        PhaseName = Title
1325    else:
1326        PhaseName = 'None'
1327    SGlines = G2spc.SGPrint(SGData)
1328    for line in SGlines: print line
1329    cell = [float(a),float(b),float(c),float(alp),float(bet),float(gam)]
1330    Volume = G2lat.calc_V(G2lat.cell2A(cell))
1331    Phase['General'] = {'Name':PhaseName,'Type':'nuclear','SGData':SGData,
1332        'Cell':[False,]+cell+[Volume,]}
1333    Phase['Atoms'] = Atoms
1334    Phase['Drawing'] = {}
1335    Phase['Histograms'] = {}
1336   
1337    return Phase
1338
1339def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
1340    # returns value(esd) string; nTZ=True for no trailing zeros
1341    # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
1342    #get the 2 significant digits in the esd
1343    edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
1344    #get the number of digits to represent them
1345    epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
1346   
1347    mdec = lambda esd: -int(math.log10(abs(esd)))
1348    ndec = lambda esd: int(1.545-math.log10(abs(esd)))
1349    if esd > 0:
1350        fmt = '"%.'+str(ndec(esd))+'f(%d)"'
1351        print fmt,ndec(esd),esd*10**(mdec(esd)+1)
1352        return fmt%(value,int(esd*10**(mdec(esd)+1)))
1353    elif esd < 0:
1354         return str(round(value,mdec(esd)))
1355    else:
1356        text = "%F"%(value)
1357        if nTZ:
1358            return text.rstrip('0')
1359        else:
1360            return text
1361
1362def Fesd(value,esd=0,nTZ=False):
1363#pythonized version of fortran routine in GSAS cifsubs directory - doesn't work correctly
1364    nint = lambda x: int(round(x))
1365    iExp = 0
1366    if value == 0. and esd == 0.:
1367        iDec = 1
1368        iFld = 5
1369    elif value == 0.:
1370        iDec = max(0.,1.545-math.log10(abs(esd)))
1371        iFld = 4+iDec
1372    elif esd == 0.:
1373        iDec = 5
1374        iFld = max(1.,math.log10(abs(value)))+3+iDec
1375    else:
1376        iFld = math.log10(max(abs(esd),abs(value)))
1377        if iFld < -4:
1378            iExp = 1-iFld
1379            iFld -= iExp
1380        elif iFld > 8:
1381            iExp = -iFld
1382            iFld += iExp
1383        if iExp:
1384            value *= 10.**iExp
1385            esd *= 10.**iExp
1386        iDec = min(7,int(max(0.,1.545-math.log10(max(0.000001*abs(value),abs(esd))))))
1387        iFld = max(1,iFld)+3+iDec
1388    if esd <= 0.:
1389        iSigw = 0
1390    else:
1391        iSig = nint(esd*(10.**iDec))
1392        iSigw = 1
1393        if iSig > 0:
1394            iSigw = int(1.+math.log10(1.*iSig))
1395    if iSigw > 2:
1396        xmult = 10.**(iSigw-2)
1397        value = xmult*nint(value/xmult)
1398        iSig = xmult*nint(iSig/xmult)           
1399        iSigw = int(1.+math.log10(1.*iSig))
1400    if iSigw == 0:
1401        fmt = '%.'+str(iDec)+'f'
1402        string = fmt%(value) 
1403    elif iDec > 0:
1404        fmt = '%.'+str(iDec)+'f(%d)'
1405        string = fmt%(value,iSig)
1406    else:
1407        fmt = '%'+str(iFld)+'d(%d)'
1408        string = fmt%(nint(value),iSig)
1409    if iExp:
1410        iDig = 1+math.log10(abs(1.*iExp))
1411        if iExp > 0:
1412            iDig += 1
1413            string += str(-iExp)
1414    return string
1415   
1416
Note: See TracBrowser for help on using the repository browser.