source: trunk/GSASIIIO.py @ 273

Last change on this file since 273 was 273, checked in by vondreele, 13 years ago

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