source: trunk/GSASIIIO.py @ 267

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

removed azm rotation - didn't work
fix azimuthal polarization

  • Property svn:keywords set to Date Author Revision URL Id
File size: 47.2 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-04-20 18:09:53 +0000 (Wed, 20 Apr 2011) $
6# $Author: vondreele $
7# $Revision: 267 $
8# $URL: trunk/GSASIIIO.py $
9# $Id: GSASIIIO.py 267 2011-04-20 18:09:53Z 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 PeakListSave(self,file,peaks):
899    print 'save peak list to file: ',self.peaklistfile
900    if not peaks:
901        dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
902        try:
903            result = dlg.ShowModal()
904        finally:
905            dlg.Destroy()
906        return
907    for peak in peaks:
908        file.write("%10.4f %12.2f %10.3f %10.3f \n" % \
909            (peak[0],peak[2],peak[4],peak[6]))
910    print 'peak list saved'
911             
912def IndexPeakListSave(self,peaks):
913    file = open(self.peaklistfile,'wa')
914    print 'save index peak list to file: ',self.peaklistfile
915    wx.BeginBusyCursor()
916    try:
917        if not peaks:
918            dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
919            try:
920                result = dlg.ShowModal()
921            finally:
922                dlg.Destroy()
923            return
924        for peak in peaks:
925            file.write("%12.6f\n" % (peak[7]))
926        file.close()
927    finally:
928        wx.EndBusyCursor()
929    print 'index peak list saved'
930   
931def ReadEXPPhase(self,filename):
932    file = open(filename, 'Ur')
933    Phase = {}
934    S = 1
935    Expr = [{},{},{},{},{},{},{},{},{}]
936    while S:
937        S = file.readline()
938        if 'EXPR NPHAS' in S[:12]:
939            Num = S[12:-1].count('0')
940            NPhas = S[12:-1].split()
941        if 'CRS' in S[:3]:
942            N = int(S[3:4])-1
943            Expr[N][S[:12]] = S[12:-1]
944    file.close()
945    PNames = []
946    for n,N in enumerate(NPhas):
947        if N != '0':
948            result = n
949            key = 'CRS'+str(n+1)+'    PNAM'
950            PNames.append(Expr[n][key])
951    if Num < 8:
952        dlg = wx.SingleChoiceDialog(self, 'Which phase to read?', 'Read phase data', PNames, wx.CHOICEDLG_STYLE)
953        try:
954            if dlg.ShowModal() == wx.ID_OK:
955                result = dlg.GetSelection()
956        finally:
957            dlg.Destroy()       
958    EXPphase = Expr[result]
959    keyList = EXPphase.keys()
960    keyList.sort()
961    SGData = {}
962    if NPhas[result] == '1':
963        Ptype = 'nuclear'
964    elif NPhas[result] in ['2','3']:
965        Ptype = 'magnetic'
966    elif NPhas[result] == '4':
967        Ptype = 'macromolecular'
968    elif NPhas[result] == '10':
969        Ptype = 'Pawley'
970    for key in keyList:
971        if 'PNAM' in key:
972           PhaseName = EXPphase[key].strip()
973        elif 'ABC   ' in key:
974            abc = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                       
975        elif 'ANGLES' in key:
976            angles = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                                               
977        elif 'SG SYM' in key:
978            SpGrp = EXPphase[key][:15].strip()
979            E,SGData = G2spc.SpcGroup(SpGrp)
980    Atoms = []
981    if Ptype == 'nuclear':
982        for key in keyList:
983            if 'AT' in key:
984                if key[11:] == 'A':
985                    S = EXPphase[key]
986                elif key[11:] == 'B':
987                    S += EXPphase[key]
988                    Atom = [S[50:58].strip(),S[:10].strip(),'',
989                        float(S[10:20]),float(S[20:30]),float(S[30:40]),
990                        float(S[40:50]),'',int(S[60:62]),S[130:131]]
991                    if Atom[9] == 'I':
992                        Atom += [float(S[68:78]),0.,0.,0.,0.,0.,0.]
993                    elif Atom[9] == 'A':
994                        Atom += [0.0,float(S[68:78]),float(S[78:88]),
995                            float(S[88:98]),float(S[98:108]),
996                            float(S[108:118]),float(S[118:128])]
997                    XYZ = Atom[3:6]
998                    Atom[7],Atom[8] = G2spc.SytSym(XYZ,SGData)
999                    Atom.append(ran.randint(0,sys.maxint))
1000                    Atoms.append(Atom)
1001    elif Ptype == 'macromolecular':
1002        for key in keyList:
1003            if 'AT' in key[6:8]:
1004                S = EXPphase[key]
1005                Atom = [int(S[56:60]),S[50:54].strip().upper(),S[54:56],
1006                    S[46:51].strip(),S[:8].strip(),'',
1007                    float(S[16:24]),float(S[24:32]),float(S[32:40]),
1008                    float(S[8:16]),'1',1,'I',float(S[40:46]),0,0,0,0,0,0]
1009                XYZ = Atom[6:9]
1010                Atom[10],Atom[11] = G2spc.SytSym(XYZ,SGData)
1011                Atom.append(ran.randint(0,sys.maxint))
1012                Atoms.append(Atom)
1013    Volume = G2lat.calc_V(G2lat.cell2A(abc+angles))
1014    Phase['General'] = {'Name':PhaseName,'Type':Ptype,'SGData':SGData,
1015        'Cell':[False,]+abc+angles+[Volume,]}
1016    Phase['Atoms'] = Atoms
1017    Phase['Drawing'] = {}
1018    Phase['Histograms'] = {}
1019    return Phase
1020       
1021def ReadPDBPhase(filename):
1022    EightPiSq = 8.*math.pi**2
1023    file = open(filename, 'Ur')
1024    Phase = {}
1025    Title = ''
1026    Compnd = ''
1027    Atoms = []
1028    A = np.zeros(shape=(3,3))
1029    S = file.readline()
1030    while S:
1031        Atom = []
1032        if 'TITLE' in S[:5]:
1033            Title = S[10:72].strip()
1034            S = file.readline()
1035        elif 'COMPND    ' in S[:10]:
1036            Compnd = S[10:72].strip()
1037            S = file.readline()
1038        elif 'CRYST' in S[:5]:
1039            abc = S[7:34].split()
1040            angles = S[34:55].split()
1041            cell=[float(abc[0]),float(abc[1]),float(abc[2]),
1042                float(angles[0]),float(angles[1]),float(angles[2])]
1043            Volume = G2lat.calc_V(G2lat.cell2A(cell))
1044            AA,AB = G2lat.cell2AB(cell)
1045            SpGrp = S[55:65]
1046            E,SGData = G2spc.SpcGroup(SpGrp)
1047            if E: 
1048                print ' ERROR in space group symbol ',SpGrp,' in file ',filename
1049                print ' N.B.: make sure spaces separate axial fields in symbol' 
1050                print G2spc.SGErrors(E)
1051                return None
1052            SGlines = G2spc.SGPrint(SGData)
1053            for line in SGlines: print line
1054            S = file.readline()
1055        elif 'SCALE' in S[:5]:
1056            V = (S[10:41].split())
1057            A[int(S[5])-1] = [float(V[0]),float(V[1]),float(V[2])]
1058            S = file.readline()
1059        elif 'ATOM' in S[:4] or 'HETATM' in S[:6]:
1060            XYZ = [float(S[31:39]),float(S[39:47]),float(S[47:55])]
1061            XYZ = np.inner(AB,XYZ)
1062            SytSym,Mult = G2spc.SytSym(XYZ,SGData)
1063            Uiso = float(S[61:67])/EightPiSq
1064            Type = S[12:14].upper()
1065            if Type[0] in '123456789':
1066                Type = Type[1:]
1067            Atom = [S[22:27].strip(),S[17:20].upper(),S[20:22],
1068                S[12:17].strip(),Type.strip(),'',XYZ[0],XYZ[1],XYZ[2],
1069                float(S[55:61]),SytSym,Mult,'I',Uiso,0,0,0,0,0,0]
1070            S = file.readline()
1071            if 'ANISOU' in S[:6]:
1072                Uij = S[30:72].split()
1073                Uij = [float(Uij[0])/10000.,float(Uij[1])/10000.,float(Uij[2])/10000.,
1074                    float(Uij[3])/10000.,float(Uij[4])/10000.,float(Uij[5])/10000.]
1075                Atom = Atom[:14]+Uij
1076                Atom[12] = 'A'
1077                S = file.readline()
1078            Atom.append(ran.randint(0,sys.maxint))
1079            Atoms.append(Atom)
1080        else:           
1081            S = file.readline()
1082    file.close()
1083    if Title:
1084        PhaseName = Title
1085    elif Compnd:
1086        PhaseName = Compnd
1087    else:
1088        PhaseName = 'None'
1089    Phase['General'] = {'Name':PhaseName,'Type':'macromolecular','SGData':SGData,
1090        'Cell':[False,]+cell+[Volume,]}
1091    Phase['Atoms'] = Atoms
1092    Phase['Drawing'] = {}
1093    Phase['Histograms'] = {}
1094   
1095    return Phase
1096   
1097def ReadCIFPhase(filename):
1098    anisoNames = ['aniso_u_11','aniso_u_22','aniso_u_33','aniso_u_12','aniso_u_13','aniso_u_23']
1099    file = open(filename, 'Ur')
1100    Phase = {}
1101    Title = ospath.split(filename)[-1]
1102    print '\n Reading cif file: ',Title
1103    Compnd = ''
1104    Atoms = []
1105    A = np.zeros(shape=(3,3))
1106    S = file.readline()
1107    while S:
1108        if '_symmetry_space_group_name_H-M' in S:
1109            SpGrp = S.split("_symmetry_space_group_name_H-M")[1].strip().strip('"').strip("'")
1110            E,SGData = G2spc.SpcGroup(SpGrp)
1111            if E:
1112                print ' ERROR in space group symbol ',SpGrp,' in file ',filename
1113                print ' N.B.: make sure spaces separate axial fields in symbol' 
1114                print G2spc.SGErrors(E)
1115                return None
1116            S = file.readline()
1117        elif '_cell' in S:
1118            if '_cell_length_a' in S:
1119                a = S.split('_cell_length_a')[1].strip().strip('"').strip("'").split('(')[0]
1120            elif '_cell_length_b' in S:
1121                b = S.split('_cell_length_b')[1].strip().strip('"').strip("'").split('(')[0]
1122            elif '_cell_length_c' in S:
1123                c = S.split('_cell_length_c')[1].strip().strip('"').strip("'").split('(')[0]
1124            elif '_cell_angle_alpha' in S:
1125                alp = S.split('_cell_angle_alpha')[1].strip().strip('"').strip("'").split('(')[0]
1126            elif '_cell_angle_beta' in S:
1127                bet = S.split('_cell_angle_beta')[1].strip().strip('"').strip("'").split('(')[0]
1128            elif '_cell_angle_gamma' in S:
1129                gam = S.split('_cell_angle_gamma')[1].strip().strip('"').strip("'").split('(')[0]
1130            S = file.readline()
1131        elif 'loop_' in S:
1132            labels = {}
1133            i = 0
1134            while S:
1135                S = file.readline()
1136                if '_atom_site' in S.strip()[:10]:
1137                    labels[S.strip().split('_atom_site_')[1].lower()] = i
1138                    i += 1
1139                else:
1140                    break
1141            if labels:
1142                if 'aniso_label' not in labels:
1143                    while S:
1144                        atom = ['','','',0,0,0,0,'','','I',0.01,0,0,0,0,0,0]
1145                        S.strip()
1146                        if len(S.split()) != len(labels):
1147                            if 'loop_' in S:
1148                                break
1149                            S += file.readline().strip()
1150                        data = S.split()
1151                        if len(data) != len(labels):
1152                            break
1153                        for key in labels:
1154                            if key == 'type_symbol':
1155                                atom[1] = data[labels[key]]
1156                            elif key == 'label':
1157                                atom[0] = data[labels[key]]
1158                            elif key == 'fract_x':
1159                                atom[3] = float(data[labels[key]].split('(')[0])
1160                            elif key == 'fract_y':
1161                                atom[4] = float(data[labels[key]].split('(')[0])
1162                            elif key == 'fract_z':
1163                                atom[5] = float(data[labels[key]].split('(')[0])
1164                            elif key == 'occupancy':
1165                                atom[6] = float(data[labels[key]].split('(')[0])
1166                            elif key == 'thermal_displace_type':
1167                                if data[labels[key]].lower() == 'uiso':
1168                                    atom[9] = 'I'
1169                                    atom[10] = float(data[labels['u_iso_or_equiv']].split('(')[0])
1170                                else:
1171                                    atom[9] = 'A'
1172                                    atom[10] = 0.0
1173                                   
1174                        atom[7],atom[8] = G2spc.SytSym(atom[3:6],SGData)
1175                        atom.append(ran.randint(0,sys.maxint))
1176                        Atoms.append(atom)
1177                        S = file.readline()
1178                else:
1179                    while S:
1180                        S.strip()
1181                        data = S.split()
1182                        if len(data) != len(labels):
1183                            break
1184                        name = data[labels['aniso_label']]
1185                        for atom in Atoms:
1186                            if name == atom[0]:
1187                                for i,uname in enumerate(anisoNames):
1188                                    atom[i+11] = float(data[labels[uname]].split('(')[0])
1189                        S = file.readline()
1190                                                                       
1191        else:           
1192            S = file.readline()
1193    file.close()
1194    if Title:
1195        PhaseName = Title
1196    else:
1197        PhaseName = 'None'
1198    SGlines = G2spc.SGPrint(SGData)
1199    for line in SGlines: print line
1200    cell = [float(a),float(b),float(c),float(alp),float(bet),float(gam)]
1201    Volume = G2lat.calc_V(G2lat.cell2A(cell))
1202    Phase['General'] = {'Name':PhaseName,'Type':'nuclear','SGData':SGData,
1203        'Cell':[False,]+cell+[Volume,]}
1204    Phase['Atoms'] = Atoms
1205    Phase['Drawing'] = {}
1206    Phase['Histograms'] = {}
1207   
1208    return Phase
Note: See TracBrowser for help on using the repository browser.