source: trunk/GSASIIIO.py @ 251

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

Add NaCl? even hkl to ImageCalibrants?.py
Fix bug in polygon mask in GSASIIimage.py
Fix copy bugs in GSASIIimgGUI.py
Make float image for float tiff file in GSASIIIO.py

  • Property svn:keywords set to Date Author Revision URL Id
File size: 47.0 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-03-01 19:16:27 +0000 (Tue, 01 Mar 2011) $
6# $Author: vondreele $
7# $Revision: 251 $
8# $URL: trunk/GSASIIIO.py $
9# $Id: GSASIIIO.py 251 2011-03-01 19:16:27Z 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':
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*'i',File.read(nVal*4))
651        elif Type == 4:
652            Value = st.unpack(byteOrd+nVal*'i',File.read(nVal*4))
653        elif Type == 5:
654            Value = st.unpack(byteOrd+nVal*'i',File.read(nVal*4))
655        elif Type == 11:
656            Value = st.unpack(byteOrd+nVal*'f',File.read(nVal*4))
657        IFD[Tag] = [Type,nVal,Value]
658    sizexy = [IFD[256][2][0],IFD[257][2][0]]
659    [nx,ny] = sizexy
660    Npix = nx*ny
661    if 272 in IFD:
662        ifd = IFD[272]
663        File.seek(ifd[2][0])
664        S = File.read(ifd[1])
665        if 'PILATUS' in S:
666            tifType = 'Pilatus'
667            dataType = 0
668            pixy = (172,172)
669            File.seek(4096)
670            if not imageOnly:
671                print 'Read Pilatus tiff file: ',filename
672            image = ar.array('L',File.read(4*Npix))
673            image = np.array(np.asarray(image),dtype=np.int32)
674    elif 262 in IFD and IFD[262][2][0] > 4:
675        tifType = 'DND'
676        pixy = (158,158)
677        File.seek(512)
678        if not imageOnly:
679            print 'Read DND SAX/WAX-detector tiff file: ',filename
680        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
681    elif sizexy == [1536,1536]:
682        tifType = 'APS Gold'
683        pixy = (150,150)
684        File.seek(64)
685        if not imageOnly:
686            print 'Read Gold tiff file:',filename
687        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
688    elif sizexy == [2048,2048] or sizexy == [1024,1024]:
689        if IFD[273][2][0] == 8:
690            if IFD[258][2][0] == 32:
691                tifType = 'PE'
692                pixy = (200,200)
693                File.seek(8)
694                if not imageOnly:
695                    print 'Read APS PE-detector tiff file: ',filename
696                if dataType == 5:
697                    image = np.array(ar.array('f',File.read(4*Npix)),dtype=np.float32)
698                else:
699                    image = np.array(ar.array('I',File.read(4*Npix)),dtype=np.int32)
700        elif IFD[273][2][0] == 4096:
701            tifType = 'MAR'
702            pixy = (158,158)
703            File.seek(4096)
704            if not imageOnly:
705                print 'Read MAR CCD tiff file: ',filename
706            image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
707           
708    else:
709        lines = ['not a known detector tiff file',]
710        return lines,0,0,0
711       
712    image = np.reshape(image,(sizexy[1],sizexy[0]))
713    center = [pixy[0]*sizexy[0]/2000,pixy[1]*sizexy[1]/2000]
714    data = {'pixelSize':pixy,'wavelength':0.10,'distance':100.0,'center':center,'size':sizexy}
715    File.close()   
716    if imageOnly:
717        return image
718    else:
719        return head,data,Npix,image
720   
721def ProjFileOpen(self):
722    file = open(self.GSASprojectfile,'rb')
723    print 'load from file: ',self.GSASprojectfile
724    wx.BeginBusyCursor()
725    try:
726        while True:
727            try:
728                data = cPickle.load(file)
729            except EOFError:
730                break
731            datum = data[0]
732            print 'load: ',datum[0]
733           
734            Id = self.PatternTree.AppendItem(parent=self.root,text=datum[0])
735            if 'PWDR' in datum[0]:               
736                self.PatternTree.SetItemPyData(Id,datum[1][:3])     #temp. trim off junk
737            else:
738                self.PatternTree.SetItemPyData(Id,datum[1])
739            for datus in data[1:]:
740                print '    load: ',datus[0]
741                               
742                sub = self.PatternTree.AppendItem(Id,datus[0])
743                self.PatternTree.SetItemPyData(sub,datus[1])
744                                               
745            if 'IMG' in datum[0]:                   #retreive image default flag & data if set
746                Data = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Image Controls'))
747                if Data['setDefault']:
748                    self.imageDefault = Data               
749        file.close()
750       
751    finally:
752        wx.EndBusyCursor()
753    print 'project load successful'
754    self.NewPlot = True
755   
756def ProjFileSave(self):
757    if not self.PatternTree.IsEmpty():
758        file = open(self.GSASprojectfile,'wb')
759        print 'save to file: ',self.GSASprojectfile
760        wx.BeginBusyCursor()
761        try:
762            item, cookie = self.PatternTree.GetFirstChild(self.root)
763            while item:
764                data = []
765                name = self.PatternTree.GetItemText(item)
766                print 'save: ',name
767                data.append([name,self.PatternTree.GetItemPyData(item)])
768                item2, cookie2 = self.PatternTree.GetFirstChild(item)
769                while item2:
770                    name = self.PatternTree.GetItemText(item2)
771                    print '    save: ',name
772                    data.append([name,self.PatternTree.GetItemPyData(item2)])
773                    item2, cookie2 = self.PatternTree.GetNextChild(item, cookie2)                           
774                item, cookie = self.PatternTree.GetNextChild(self.root, cookie)                           
775                cPickle.dump(data,file,1)
776            file.close()
777        finally:
778            wx.EndBusyCursor()
779        print 'project save successful'
780       
781def SaveIntegration(self,PickId,data):
782    azms = self.Integrate[1]
783    X = self.Integrate[2][:-1]
784    Xminmax = [X[0],X[-1]]
785    N = len(X)
786    Id = self.PatternTree.GetItemParent(PickId)
787    name = self.PatternTree.GetItemText(Id)
788    name = name.replace('IMG ','PWDR ')
789    Comments = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Comments'))
790    names = ['Type','Lam','Zero','Polariz.','U','V','W','X','Y','SH/L','Azimuth'] 
791    codes = [0 for i in range(11)]
792    parms = ['PXC',data['wavelength'],0.0,0.0,1.0,-1.0,0.3,0.0,1.0,0.0,0.0]
793    Azms = [(azms[i+1]+azms[i])/2. for i in range(len(azms)-1)]
794    for i,azm in enumerate(Azms):
795        item, cookie = self.PatternTree.GetFirstChild(self.root)
796        Id = 0
797        while item:
798            Name = self.PatternTree.GetItemText(item)
799            if name == Name:
800                Id = item
801            item, cookie = self.PatternTree.GetNextChild(self.root, cookie)
802        parms[10] = azm
803        Y = self.Integrate[0][i]
804        W = 1./Y                    #probably not true
805        Sample = {'Scale':[1.0,True],'Type':'Debye-Scherrer','Absorption':[0.0,False],'DisplaceX':[0.0,False],
806            'DisplaceY':[0.0,False],'Diffuse':[],'Temperature':300.,'Pressure':1.0,'Humidity':0.0,'Voltage':0.0,'Force':0.0}
807        if Id:
808            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Comments'),Comments)                   
809            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Limits'),[tuple(Xminmax),Xminmax])
810            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Background'),[['chebyschev',1,3,1.0,0.0,0.0]])
811            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Instrument Parameters'),[tuple(parms),parms,codes,names])
812            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Peak List'),[])
813            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Index Peak List'),[])
814            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Unit Cells List'),[])             
815        else:
816            Id = self.PatternTree.AppendItem(parent=self.root,text=name+" Azm= %.2f"%(azm))
817            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Comments'),Comments)                   
818            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Limits'),[tuple(Xminmax),Xminmax])
819            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Background'),[['chebyschev',1,3,1.0,0.0,0.0]])
820            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Instrument Parameters'),[tuple(parms),parms,codes,names])
821            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Sample Parameters'),Sample)
822            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Peak List'),[])
823            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Index Peak List'),[])
824            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Unit Cells List'),[])             
825        self.PatternTree.SetItemPyData(Id,[[''],[np.array(X),np.array(Y),np.array(W),np.zeros(N),np.zeros(N),np.zeros(N)]])
826    self.PatternTree.SelectItem(Id)
827    self.PatternTree.Expand(Id)
828    self.PatternId = Id
829           
830def powderFxyeSave(self,exports,powderfile):
831    head,tail = ospath.split(powderfile)
832    name,ext = tail.split('.')
833    wx.BeginBusyCursor()
834    for i,export in enumerate(exports):
835        filename = ospath.join(head,name+'-%03d.'%(i)+ext)
836        prmname = filename.strip(ext)+'prm'
837        prm = open(prmname,'w')      #old style GSAS parm file
838        PickId = G2gd.GetPatternTreeItemId(self, self.root, export)
839        Values,Names = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self, \
840            PickId, 'Instrument Parameters'))[1::2]     #get values & names
841        Inst = dict(zip(Names,Values))
842        print Inst['Type']
843        prm.write( '            123456789012345678901234567890123456789012345678901234567890        '+'\n')
844        prm.write( 'INS   BANK      1                                                               '+'\n')
845        prm.write(('INS   HTYPE   %sR                                                              '+'\n')%(Inst['Type']))
846        if 'Lam1' in Inst:              #Ka1 & Ka2
847            prm.write(('INS  1 ICONS%10.7f%10.7f    0.0000               0.990    0     0.500   '+'\n')%(Inst['Lam1'],Inst['Lam2']))
848        elif 'Lam' in Inst:             #single wavelength
849            prm.write(('INS  1 ICONS%10.7f%10.7f    0.0000               0.990    0     0.500   '+'\n')%(Inst['Lam'],0.0))
850        prm.write( 'INS  1 IRAD     0                                                               '+'\n')
851        prm.write( 'INS  1I HEAD                                                                    '+'\n')
852        prm.write( 'INS  1I ITYP    0    0.0000  180.0000         1                                 '+'\n')
853        prm.write(('INS  1DETAZM%10.3f                                                          '+'\n')%(Inst['Azimuth']))
854        prm.write( 'INS  1PRCF1     3    8   0.00100                                                '+'\n')
855        prm.write(('INS  1PRCF11     %15.6g%15.6g%15.6g%15.6g   '+'\n')%(Inst['U'],Inst['V'],Inst['W'],0.0))
856        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.))
857        prm.close()
858        file = open(filename,'w')
859        print 'save powder pattern to file: ',filename
860        try:
861            x,y,w,yc,yb,yd = self.PatternTree.GetItemPyData(PickId)[1]
862            file.write(powderfile+'\n')
863            file.write('BANK 1 %d %d CONS %.2f %.2f 0 0 FXYE\n'%(len(x),len(x),\
864                100.*x[0],100.*(x[1]-x[0])))
865            s = list(np.sqrt(1./np.array(w)))       
866            XYW = zip(x,y,s)
867            for X,Y,S in XYW:
868                file.write("%15.6g %15.6g %15.6g\n" % (100.*X,Y,max(S,1.0)))
869            file.close()
870        finally:
871            wx.EndBusyCursor()
872        print 'powder pattern file written'
873       
874def powderXyeSave(self,exports,powderfile):
875    head,tail = ospath.split(powderfile)
876    name,ext = tail.split('.')
877    for i,export in enumerate(exports):
878        filename = ospath.join(head,name+'-%03d.'%(i)+ext)
879        PickId = G2gd.GetPatternTreeItemId(self, self.root, export)
880        file = open(filename,'w')
881        file.write('#%s\n'%(export))
882        print 'save powder pattern to file: ',filename
883        wx.BeginBusyCursor()
884        try:
885            x,y,w,yc,yb,yd = self.PatternTree.GetItemPyData(PickId)[1]
886            XYW = zip(x,y,w)
887            for X,Y,W in XYW:
888                file.write("%15.6g %15.6g %15.6g\n" % (X,Y,W))
889            file.close()
890        finally:
891            wx.EndBusyCursor()
892        print 'powder pattern file written'
893   
894def PeakListSave(self,file,peaks):
895    print 'save peak list to file: ',self.peaklistfile
896    if not peaks:
897        dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
898        try:
899            result = dlg.ShowModal()
900        finally:
901            dlg.Destroy()
902        return
903    for peak in peaks:
904        file.write("%10.4f %12.2f %10.3f %10.3f \n" % \
905            (peak[0],peak[2],peak[4],peak[6]))
906    print 'peak list saved'
907             
908def IndexPeakListSave(self,peaks):
909    file = open(self.peaklistfile,'wa')
910    print 'save index peak list to file: ',self.peaklistfile
911    wx.BeginBusyCursor()
912    try:
913        if not peaks:
914            dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
915            try:
916                result = dlg.ShowModal()
917            finally:
918                dlg.Destroy()
919            return
920        for peak in peaks:
921            file.write("%12.6f\n" % (peak[7]))
922        file.close()
923    finally:
924        wx.EndBusyCursor()
925    print 'index peak list saved'
926   
927def ReadEXPPhase(self,filename):
928    file = open(filename, 'Ur')
929    Phase = {}
930    S = 1
931    Expr = [{},{},{},{},{},{},{},{},{}]
932    while S:
933        S = file.readline()
934        if 'EXPR NPHAS' in S[:12]:
935            Num = S[12:-1].count('0')
936            NPhas = S[12:-1].split()
937        if 'CRS' in S[:3]:
938            N = int(S[3:4])-1
939            Expr[N][S[:12]] = S[12:-1]
940    file.close()
941    PNames = []
942    for n,N in enumerate(NPhas):
943        if N != '0':
944            result = n
945            key = 'CRS'+str(n+1)+'    PNAM'
946            PNames.append(Expr[n][key])
947    if Num < 8:
948        dlg = wx.SingleChoiceDialog(self, 'Which phase to read?', 'Read phase data', PNames, wx.CHOICEDLG_STYLE)
949        try:
950            if dlg.ShowModal() == wx.ID_OK:
951                result = dlg.GetSelection()
952        finally:
953            dlg.Destroy()       
954    EXPphase = Expr[result]
955    keyList = EXPphase.keys()
956    keyList.sort()
957    SGData = {}
958    if NPhas[result] == '1':
959        Ptype = 'nuclear'
960    elif NPhas[result] in ['2','3']:
961        Ptype = 'magnetic'
962    elif NPhas[result] == '4':
963        Ptype = 'macromolecular'
964    elif NPhas[result] == '10':
965        Ptype = 'Pawley'
966    for key in keyList:
967        if 'PNAM' in key:
968           PhaseName = EXPphase[key].strip()
969        elif 'ABC   ' in key:
970            abc = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                       
971        elif 'ANGLES' in key:
972            angles = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                                               
973        elif 'SG SYM' in key:
974            SpGrp = EXPphase[key][:15].strip()
975            E,SGData = G2spc.SpcGroup(SpGrp)
976    Atoms = []
977    if Ptype == 'nuclear':
978        for key in keyList:
979            if 'AT' in key:
980                if key[11:] == 'A':
981                    S = EXPphase[key]
982                elif key[11:] == 'B':
983                    S += EXPphase[key]
984                    Atom = [S[50:58].strip(),S[:10].strip(),'',
985                        float(S[10:20]),float(S[20:30]),float(S[30:40]),
986                        float(S[40:50]),'',int(S[60:62]),S[130:131]]
987                    if Atom[9] == 'I':
988                        Atom += [float(S[68:78]),0.,0.,0.,0.,0.,0.]
989                    elif Atom[9] == 'A':
990                        Atom += [0.0,float(S[68:78]),float(S[78:88]),
991                            float(S[88:98]),float(S[98:108]),
992                            float(S[108:118]),float(S[118:128])]
993                    XYZ = Atom[3:6]
994                    Atom[7],Atom[8] = G2spc.SytSym(XYZ,SGData)
995                    Atom.append(ran.randint(0,sys.maxint))
996                    Atoms.append(Atom)
997    elif Ptype == 'macromolecular':
998        for key in keyList:
999            if 'AT' in key[6:8]:
1000                S = EXPphase[key]
1001                Atom = [int(S[56:60]),S[50:54].strip().upper(),S[54:56],
1002                    S[46:51].strip(),S[:8].strip(),'',
1003                    float(S[16:24]),float(S[24:32]),float(S[32:40]),
1004                    float(S[8:16]),'1',1,'I',float(S[40:46]),0,0,0,0,0,0]
1005                XYZ = Atom[6:9]
1006                Atom[10],Atom[11] = G2spc.SytSym(XYZ,SGData)
1007                Atom.append(ran.randint(0,sys.maxint))
1008                Atoms.append(Atom)
1009    Volume = G2lat.calc_V(G2lat.cell2A(abc+angles))
1010    Phase['General'] = {'Name':PhaseName,'Type':Ptype,'SGData':SGData,
1011        'Cell':[False,]+abc+angles+[Volume,]}
1012    Phase['Atoms'] = Atoms
1013    Phase['Drawing'] = {}
1014    Phase['Histograms'] = {}
1015    return Phase
1016       
1017def ReadPDBPhase(filename):
1018    EightPiSq = 8.*math.pi**2
1019    file = open(filename, 'Ur')
1020    Phase = {}
1021    Title = ''
1022    Compnd = ''
1023    Atoms = []
1024    A = np.zeros(shape=(3,3))
1025    S = file.readline()
1026    while S:
1027        Atom = []
1028        if 'TITLE' in S[:5]:
1029            Title = S[10:72].strip()
1030            S = file.readline()
1031        elif 'COMPND    ' in S[:10]:
1032            Compnd = S[10:72].strip()
1033            S = file.readline()
1034        elif 'CRYST' in S[:5]:
1035            abc = S[7:34].split()
1036            angles = S[34:55].split()
1037            cell=[float(abc[0]),float(abc[1]),float(abc[2]),
1038                float(angles[0]),float(angles[1]),float(angles[2])]
1039            Volume = G2lat.calc_V(G2lat.cell2A(cell))
1040            AA,AB = G2lat.cell2AB(cell)
1041            SpGrp = S[55:65]
1042            E,SGData = G2spc.SpcGroup(SpGrp)
1043            if E: 
1044                print ' ERROR in space group symbol ',SpGrp,' in file ',filename
1045                print ' N.B.: make sure spaces separate axial fields in symbol' 
1046                print G2spc.SGErrors(E)
1047                return None
1048            SGlines = G2spc.SGPrint(SGData)
1049            for line in SGlines: print line
1050            S = file.readline()
1051        elif 'SCALE' in S[:5]:
1052            V = (S[10:41].split())
1053            A[int(S[5])-1] = [float(V[0]),float(V[1]),float(V[2])]
1054            S = file.readline()
1055        elif 'ATOM' in S[:4] or 'HETATM' in S[:6]:
1056            XYZ = [float(S[31:39]),float(S[39:47]),float(S[47:55])]
1057            XYZ = np.inner(AB,XYZ)
1058            SytSym,Mult = G2spc.SytSym(XYZ,SGData)
1059            Uiso = float(S[61:67])/EightPiSq
1060            Type = S[12:14].upper()
1061            if Type[0] in '123456789':
1062                Type = Type[1:]
1063            Atom = [S[22:27].strip(),S[17:20].upper(),S[20:22],
1064                S[12:17].strip(),Type.strip(),'',XYZ[0],XYZ[1],XYZ[2],
1065                float(S[55:61]),SytSym,Mult,'I',Uiso,0,0,0,0,0,0]
1066            S = file.readline()
1067            if 'ANISOU' in S[:6]:
1068                Uij = S[30:72].split()
1069                Uij = [float(Uij[0])/10000.,float(Uij[1])/10000.,float(Uij[2])/10000.,
1070                    float(Uij[3])/10000.,float(Uij[4])/10000.,float(Uij[5])/10000.]
1071                Atom = Atom[:14]+Uij
1072                Atom[12] = 'A'
1073                S = file.readline()
1074            Atom.append(ran.randint(0,sys.maxint))
1075            Atoms.append(Atom)
1076        else:           
1077            S = file.readline()
1078    file.close()
1079    if Title:
1080        PhaseName = Title
1081    elif Compnd:
1082        PhaseName = Compnd
1083    else:
1084        PhaseName = 'None'
1085    Phase['General'] = {'Name':PhaseName,'Type':'macromolecular','SGData':SGData,
1086        'Cell':[False,]+cell+[Volume,]}
1087    Phase['Atoms'] = Atoms
1088    Phase['Drawing'] = {}
1089    Phase['Histograms'] = {}
1090   
1091    return Phase
1092   
1093def ReadCIFPhase(filename):
1094    anisoNames = ['aniso_u_11','aniso_u_22','aniso_u_33','aniso_u_12','aniso_u_13','aniso_u_23']
1095    file = open(filename, 'Ur')
1096    Phase = {}
1097    Title = ospath.split(filename)[-1]
1098    print '\n Reading cif file: ',Title
1099    Compnd = ''
1100    Atoms = []
1101    A = np.zeros(shape=(3,3))
1102    S = file.readline()
1103    while S:
1104        if '_symmetry_space_group_name_H-M' in S:
1105            SpGrp = S.split("_symmetry_space_group_name_H-M")[1].strip().strip('"').strip("'")
1106            E,SGData = G2spc.SpcGroup(SpGrp)
1107            if E:
1108                print ' ERROR in space group symbol ',SpGrp,' in file ',filename
1109                print ' N.B.: make sure spaces separate axial fields in symbol' 
1110                print G2spc.SGErrors(E)
1111                return None
1112            S = file.readline()
1113        elif '_cell' in S:
1114            if '_cell_length_a' in S:
1115                a = S.split('_cell_length_a')[1].strip().strip('"').strip("'").split('(')[0]
1116            elif '_cell_length_b' in S:
1117                b = S.split('_cell_length_b')[1].strip().strip('"').strip("'").split('(')[0]
1118            elif '_cell_length_c' in S:
1119                c = S.split('_cell_length_c')[1].strip().strip('"').strip("'").split('(')[0]
1120            elif '_cell_angle_alpha' in S:
1121                alp = S.split('_cell_angle_alpha')[1].strip().strip('"').strip("'").split('(')[0]
1122            elif '_cell_angle_beta' in S:
1123                bet = S.split('_cell_angle_beta')[1].strip().strip('"').strip("'").split('(')[0]
1124            elif '_cell_angle_gamma' in S:
1125                gam = S.split('_cell_angle_gamma')[1].strip().strip('"').strip("'").split('(')[0]
1126            S = file.readline()
1127        elif 'loop_' in S:
1128            labels = {}
1129            i = 0
1130            while S:
1131                S = file.readline()
1132                if '_atom_site' in S.strip()[:10]:
1133                    labels[S.strip().split('_atom_site_')[1].lower()] = i
1134                    i += 1
1135                else:
1136                    break
1137            if labels:
1138                if 'aniso_label' not in labels:
1139                    while S:
1140                        atom = ['','','',0,0,0,0,'','','I',0.01,0,0,0,0,0,0]
1141                        S.strip()
1142                        if len(S.split()) != len(labels):
1143                            if 'loop_' in S:
1144                                break
1145                            S += file.readline().strip()
1146                        data = S.split()
1147                        if len(data) != len(labels):
1148                            break
1149                        for key in labels:
1150                            if key == 'type_symbol':
1151                                atom[1] = data[labels[key]]
1152                            elif key == 'label':
1153                                atom[0] = data[labels[key]]
1154                            elif key == 'fract_x':
1155                                atom[3] = float(data[labels[key]].split('(')[0])
1156                            elif key == 'fract_y':
1157                                atom[4] = float(data[labels[key]].split('(')[0])
1158                            elif key == 'fract_z':
1159                                atom[5] = float(data[labels[key]].split('(')[0])
1160                            elif key == 'occupancy':
1161                                atom[6] = float(data[labels[key]].split('(')[0])
1162                            elif key == 'thermal_displace_type':
1163                                if data[labels[key]].lower() == 'uiso':
1164                                    atom[9] = 'I'
1165                                    atom[10] = float(data[labels['u_iso_or_equiv']].split('(')[0])
1166                                else:
1167                                    atom[9] = 'A'
1168                                    atom[10] = 0.0
1169                                   
1170                        atom[7],atom[8] = G2spc.SytSym(atom[3:6],SGData)
1171                        atom.append(ran.randint(0,sys.maxint))
1172                        Atoms.append(atom)
1173                        S = file.readline()
1174                else:
1175                    while S:
1176                        S.strip()
1177                        data = S.split()
1178                        if len(data) != len(labels):
1179                            break
1180                        name = data[labels['aniso_label']]
1181                        for atom in Atoms:
1182                            if name == atom[0]:
1183                                for i,uname in enumerate(anisoNames):
1184                                    atom[i+11] = float(data[labels[uname]].split('(')[0])
1185                        S = file.readline()
1186                                                                       
1187        else:           
1188            S = file.readline()
1189    file.close()
1190    if Title:
1191        PhaseName = Title
1192    else:
1193        PhaseName = 'None'
1194    SGlines = G2spc.SGPrint(SGData)
1195    for line in SGlines: print line
1196    cell = [float(a),float(b),float(c),float(alp),float(bet),float(gam)]
1197    Volume = G2lat.calc_V(G2lat.cell2A(cell))
1198    Phase['General'] = {'Name':PhaseName,'Type':'nuclear','SGData':SGData,
1199        'Cell':[False,]+cell+[Volume,]}
1200    Phase['Atoms'] = Atoms
1201    Phase['Drawing'] = {}
1202    Phase['Histograms'] = {}
1203   
1204    return Phase
Note: See TracBrowser for help on using the repository browser.