source: trunk/GSASIIIO.py @ 75

Last change on this file since 75 was 75, checked in by toby, 12 years ago

remove import * from G2comp

File size: 37.3 KB
Line 
1"""GSASIIIO: functions for IO of data
2   Copyright: 2008, Robert B. Von Dreele (Argonne National Laboratory)
3"""
4
5import wx
6import math
7import numpy as np
8import cPickle
9import sys
10import GSASIIpath
11import GSASIIgrid as G2gd
12import os.path as ospath
13
14def sfloat(S):
15    if S.strip():
16        return float(S)
17    else:
18        return 0.0
19
20def sint(S):
21    if S.strip():
22        return int(S)
23    else:
24        return 0
25
26def SelectPowderData(self, filename):
27    """Selects banks of data from a filename of any GSAS powder data format
28    Input - filename: any GSAS powder data formatted file (currently STD, FXYE, FXY & ESD)
29    Returns - a list of banks to be read; each entry in list is a tuple containing:
30    filename: same as input filename
31    Pos: position for start of data; record just after BANK record
32    Bank: the BANK record
33    """
34    File = open(filename,'Ur')
35    Title = '''
36First line of this file:
37'''+File.readline()
38    dlg = wx.MessageDialog(self, Title, 'Is this the file you want?', 
39    wx.YES_NO | wx.ICON_QUESTION)
40    try:
41        result = dlg.ShowModal()
42    finally:
43        dlg.Destroy()
44    if result == wx.ID_NO: return (0,0)
45       
46    self.IparmName = GetInstrumentFile(self,filename)
47    if self.IparmName:
48        Iparm = GetInstrumentData(self.IparmName)
49    else:
50        Iparm = {}                                               #Assume CuKa lab data if no iparm file
51        Iparm['INS   HTYPE '] = 'PXC '
52        Iparm['INS  1 ICONS'] = '  1.540500  1.544300       0.0         0       0.7    0       0.5   '
53        Iparm['INS  1PRCF1 '] = '    3    8      0.01                                                '
54        Iparm['INS  1PRCF11'] = '   2.000000E+00  -2.000000E+00   5.000000E+00   0.000000E+00        '
55        Iparm['INS  1PRCF12'] = '   0.000000E+00   0.000000E+00   0.150000E-01   0.150000E-01        '
56    S = 1
57    Banks = []
58    Pos = []
59    FoundData = []
60    Comments = []
61    wx.BeginBusyCursor()
62    try:
63        while S:
64            S = File.readline()
65            if S[:1] != '#':
66                if S[:4] == 'BANK':
67                    Banks.append(S)
68                    Pos.append(File.tell())
69            else:
70                Comments.append(S[:-1])
71        File.close()
72    finally:
73        wx.EndBusyCursor()
74    if Comments:
75       print 'Comments on file:'
76       for Comment in Comments: print Comment
77    if Banks:
78        result = [0]
79        if len(Banks) >= 2:
80            dlg = wx.MultiChoiceDialog(self, 'Which scans do you want?', 'Select scans', Banks, wx.CHOICEDLG_STYLE)
81            try:
82                if dlg.ShowModal() == wx.ID_OK:
83                    result = dlg.GetSelections()
84                else:
85                    result = []
86            finally:
87                dlg.Destroy()
88        for i in result:
89            FoundData.append((filename,Pos[i],Banks[i]))
90    else:
91        dlg = wx.MessageDialog(self, 'ERROR - this is not a GSAS powder data file', 'No BANK records', wx.OK | wx.ICON_ERROR)
92        try:
93            result = dlg.ShowModal()
94        finally:
95            dlg.Destroy()
96    return FoundData,Iparm,Comments
97
98def GetInstrumentFile(self,filename):
99    import os.path as op
100    dlg = wx.FileDialog(self,'Choose an instrument file','.', '', 'GSAS iparm file (*.prm)|*.prm|All files(*.*)|*.*', wx.OPEN)
101    if self.dirname: 
102        dlg.SetDirectory(self.dirname)
103        Tname = filename[:filename.index('.')]+'.prm'
104        if op.exists(Tname):
105            self.IparmName = Tname       
106    if self.IparmName: dlg.SetFilename(self.IparmName)
107    filename = ''
108    try:
109        if dlg.ShowModal() == wx.ID_OK:
110            filename = dlg.GetPath()
111    finally:
112        dlg.Destroy()
113    return filename
114
115def GetInstrumentData(IparmName):
116    file = open(IparmName, 'Ur')
117    S = 1
118    Iparm = {}
119    while S:
120        S = file.readline()
121        Iparm[S[:12]] = S[12:-1]
122    return Iparm
123   
124def GetPowderPeaks(fileName):
125    sind = lambda x: math.sin(x*math.pi/180.)
126    asind = lambda x: 180.*math.asin(x)/math.pi
127    Cuka = 1.54052
128    File = open(fileName,'Ur')
129    Comments = []
130    peaks = []
131    S = File.readline()
132    while S:
133        if S[:1] == '#':
134            Comments.append(S[:-1])
135        else:
136            item = S.split()
137            if len(item) == 1:
138                peaks.append([float(item[0]),1.0])
139            elif len(item) > 1:
140                peaks.append([float(item[0]),float(item[0])])
141        S = File.readline()
142    File.close()
143    if Comments:
144       print 'Comments on file:'
145       for Comment in Comments: print Comment
146    Peaks = []
147    if peaks[0][0] > peaks[-1][0]:          # d-spacings - assume CuKa
148        for peak in peaks:
149            dsp = peak[0]
150            sth = Cuka/(2.0*dsp)
151            if sth < 1.0:
152                tth = 2.0*asind(sth)
153            else:
154                break
155            Peaks.append([tth,peak[1],True,False,0,0,0,dsp,0.0])
156    else:                                   #2-thetas - assume Cuka (for now)
157        for peak in peaks:
158            tth = peak[0]
159            dsp = Cuka/(2.0*sind(tth/2.0))
160            Peaks.append([tth,peak[1],True,False,0,0,0,dsp,0.0])
161    return Comments,Peaks
162   
163def GetHKLData(filename):
164    print 'Reading: '+filename
165    File = open(filename,'Ur')
166    HKLref = []
167    HKLmin = [1000,1000,1000]
168    HKLmax = [0,0,0]
169    FoMax = 0
170    ifFc = False
171    S = File.readline()
172    while '#' in S[0]:        #get past comments if any
173        S = File.readline()       
174    if '_' in S:         #cif style .hkl file
175        while 'loop_' not in S:         #skip preliminaries if any - can't have 'loop_' in them!
176            S = File.readline()       
177        S = File.readline()             #get past 'loop_' line
178        pos = 0
179        hpos = kpos = lpos = Fosqpos = Fcsqpos = sigpos = -1
180        while S:
181            if '_' in S:
182                if 'index_h' in S:
183                    hpos = pos
184                elif 'index_k' in S:
185                    kpos = pos
186                elif 'index_l' in S:
187                    lpos = pos
188                elif 'F_squared_meas' in S:
189                    Fosqpos = pos
190                elif 'F_squared_calc' in S:
191                    Fcsqpos = pos
192                elif 'F_squared_sigma' in S:
193                    sigpos = pos
194                pos += 1
195            else:
196                data = S.split()
197                if data:                    #avoid blank lines
198                    HKL = np.array([int(data[hpos]),int(data[kpos]),int(data[lpos])])
199                    h,k,l = HKL
200                    Fosq = float(data[Fosqpos])
201                    if sigpos != -1:
202                        sigFosq = float(data[sigpos])
203                    else:
204                        sigFosq = 1.
205                    if Fcsqpos != -1:
206                        Fcsq = float(data[Fcsqpos])
207                        if Fcsq:
208                            ifFc = True
209                    else:
210                        Fcsq = 0.
211                       
212                    HKLmin = [min(h,HKLmin[0]),min(k,HKLmin[1]),min(l,HKLmin[2])]
213                    HKLmax = [max(h,HKLmax[0]),max(k,HKLmax[1]),max(l,HKLmax[2])]
214                    FoMax = max(FoMax,Fosq)
215                    HKLref.append([HKL,Fosq,sigFosq,Fcsq,0,0,0])                 #room for Fc, Fcp, Fcpp & phase
216            S = File.readline()
217    else:                   #dumb h,k,l,Fo,sigFo .hkl file
218        while S:
219            h,k,l,Fo,sigFo = S.split()
220            HKL = np.array([int(h),int(k),int(l)])
221            h,k,l = HKL
222            Fo = float(Fo)
223            sigFo = float(sigFo)
224            HKLmin = [min(h,HKLmin[0]),min(k,HKLmin[1]),min(l,HKLmin[2])]
225            HKLmax = [max(h,HKLmax[0]),max(k,HKLmax[1]),max(l,HKLmax[2])]
226            FoMax = max(FoMax,Fo)
227            HKLref.append([HKL,Fo**2,2.*Fo*sigFo,0,0,0,0])                 #room for Fc, Fcp, Fcpp & phase
228            S = File.readline()
229    File.close()
230    return HKLref,HKLmin,HKLmax,FoMax,ifFc
231
232def GetPowderData(filename,Pos,Bank,DataType):
233    '''Reads one BANK of data from GSAS raw powder data file
234    input:
235    filename: GSAS raw powder file dataname
236    Pos: start of data in file just after BANK record
237    Bank: the BANK record
238    DataType: powder data type, e.g. "PXC" for Powder X-ray CW data
239    returns: list [x,y,e,yc,yb]
240    x: np.array of x-axis values
241    y: np.array of powder pattern intensities
242    w: np.array of w=sig(intensity)^2 values
243    yc: np.array of calc. intensities (zero)
244    yb: np.array of calc. background (zero)
245    yd: np.array of obs-calc profiles
246    '''
247    print 'Reading: '+filename
248    print 'Bank:    '+Bank[:-1]
249    if 'FXYE' in Bank:
250        return GetFXYEdata(filename,Pos,Bank,DataType)
251    elif 'FXY' in Bank:
252        return GetFXYdata(filename,Pos,Bank,DataType)
253    elif 'ESD' in Bank:
254        return GetESDdata(filename,Pos,Bank,DataType)
255    elif 'STD' in Bank:
256        return GetSTDdata(filename,Pos,Bank,DataType)
257    else:
258        return GetSTDdata(filename,Pos,Bank,DataType)
259    return []
260
261def GetFXYEdata(filename,Pos,Bank,DataType):
262    File = open(filename,'Ur')
263    File.seek(Pos)
264    x = []
265    y = []
266    w = []
267    S = File.readline()
268    while S and S[:4] != 'BANK':
269        vals = S.split()
270        if DataType[2] == 'C':
271            x.append(float(vals[0])/100.)               #CW: from centidegrees to degrees
272        elif DataType[2] == 'T':
273            x.append(float(vals[0])/1000.0)             #TOF: from musec to millisec
274        f = float(vals[1])
275        if f <= 0.0:
276            y.append(0.0)
277            w.append(1.0)
278        else:
279            y.append(float(vals[1]))
280            w.append(1.0/float(vals[2])**2)
281        S = File.readline()
282    File.close()
283    N = len(x)
284    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
285   
286def GetFXYdata(filename,Pos,Bank,DataType):
287    File = open(filename,'Ur')
288    File.seek(Pos)
289    x = []
290    y = []
291    w = []
292    S = File.readline()
293    while S and S[:4] != 'BANK':
294        vals = S.split()
295        if DataType[2] == 'C':
296            x.append(float(vals[0])/100.)               #CW: from centidegrees to degrees
297        elif DataType[2] == 'T':
298            x.append(float(vals[0])/1000.0)             #TOF: from musec to millisec
299        f = float(vals[1])
300        if f > 0.0:
301            y.append(f)
302            w.append(1.0/f)
303        else:             
304            y.append(0.0)
305            w.append(1.0)
306        S = File.readline()
307    File.close()
308    N = len(x)
309    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
310   
311def GetESDdata(filename,Pos,Bank,DataType):
312    File = open(filename,'Ur')
313    cons = Bank.split()
314    if DataType[2] == 'C':
315        start = float(cons[5])/100.0               #CW: from centidegrees to degrees
316        step = float(cons[6])/100.0
317    elif DataType[2] == 'T':
318        start = float(cons[5])/1000.0              #TOF: from musec to millisec
319        step = float(cons[6])/1000.0
320    File.seek(Pos)
321    x = []
322    y = []
323    w = []
324    S = File.readline()
325    j = 0
326    while S and S[:4] != 'BANK':
327        for i in range(0,80,16):
328            xi = start+step*j
329            yi = sfloat(S[i:i+8])
330            ei = sfloat(S[i+8:i+16])
331            x.append(xi)
332            if yi > 0.0:
333                y.append(yi)
334                w.append(1.0/ei**2)
335            else:             
336                y.append(0.0)
337                w.append(1.0)
338            j += 1
339        S = File.readline()
340    File.close()
341    N = len(x)
342    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
343
344def GetSTDdata(filename,Pos,Bank,DataType):
345    File = open(filename,'Ur')
346    cons = Bank.split()
347    Nch = cons[2]
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 - not likely!
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,8):
362            xi = start+step*j
363            ni = max(sint(S[i:i+2]),1)
364            yi = max(sfloat(S[i+2:i+8]),0.0)
365            if yi:
366                ei = math.sqrt(yi*ni)
367            else:
368                yi = 0.0
369                ei = 1.0
370            j += 1
371            if j < Nch:
372                x.append(xi)
373                y.append(yi)
374                w.append(1.0/ei**2)
375        S = File.readline()
376    File.close()
377    N = len(x)
378    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
379   
380def GetImageData(imagefile,imageOnly=False):
381    ext = ospath.splitext(imagefile)[1]
382    Comments = []
383    if ext == '.tif':
384        Comments,Data,Size,Image = GetTifData(imagefile)
385    elif ext == '.img':
386        Comments,Data,Size,Image = GetImgData(imagefile)
387        Image[0][0] = 0
388    elif ext == '.mar3450' or ext == '.mar2300':
389        Comments,Data,Size,Image = GetMAR345Data(imagefile)
390    elif ext in ['.sum','.avg']:
391        Comments,Data,Size,Image = GetGEsumData(imagefile)
392    elif ext == '.G2img':
393        return GetG2Image(imagefile)
394    if imageOnly:
395        return Image
396    else:
397        return Comments,Data,Size,Image
398       
399def PutG2Image(filename,image):
400    File = open(filename,'wb')
401    cPickle.dump(image,File,1)
402    File.close()
403    return
404   
405def GetG2Image(filename):
406    File = open(filename,'rb')
407    image = cPickle.load(File)
408    File.close()
409    return image
410   
411def GetGEsumData(filename,imageOnly=False):
412    import array as ar
413    if not imageOnly:
414        print 'Read GE sum file: ',filename   
415    File = open(filename,'rb')
416    size = 2048
417    if '.sum' in filename:
418        head = ['GE detector sum data from APS 1-ID',]
419    if '.avg' in filename:
420        head = ['GE detector avg data from APS 1-ID',]
421    image = np.zeros(shape=(size,size),dtype=np.int32)
422    row = 0
423    pos = 0
424    while row < size:
425        File.seek(pos)
426        if '.sum' in filename:
427            line = ar.array('f',File.read(4*size))
428            pos += 4*size
429        elif '.avg' in filename:
430            line = ar.array('H',File.read(2*size))
431            pos += 2*size
432        image[row] = np.asarray(line)
433        row += 1
434    data = {'pixelSize':(200,200),'wavelength':0.15,'distance':250.0,'center':[204.8,204.8]} 
435    File.close()   
436    if imageOnly:
437        return image
438    else:
439        return head,data,size,image
440       
441def GetImgData(filename,imageOnly=False):
442    import struct as st
443    import array as ar
444    if not imageOnly:
445        print 'Read ADSC img file: ',filename
446    File = open(filename,'rb')
447    head = File.read(511)
448    lines = head.split('\n')
449    head = []
450    center = [0,0]
451    for line in lines[1:-2]:
452        line = line.strip()[:-1]
453        if line:
454            if 'SIZE1' in line:
455                size = int(line.split('=')[1])
456            elif 'WAVELENGTH' in line:
457                wave = float(line.split('=')[1])
458            elif 'BIN' in line:
459                if line.split('=')[1] == '2x2':
460                    pixel=(102,102)
461                else:
462                    pixel = (51,51)
463            elif 'DISTANCE' in line:
464                distance = float(line.split('=')[1])
465            elif 'CENTER_X' in line:
466                center[0] = float(line.split('=')[1])
467            elif 'CENTER_Y' in line:
468                center[1] = float(line.split('=')[1])
469            head.append(line)
470    data = {'pixelSize':pixel,'wavelength':wave,'distance':distance,'center':center}
471    image = []
472    row = 0
473    pos = 512
474    image = np.zeros(shape=(size,size),dtype=np.int32)   
475    while row < size:
476        File.seek(pos)
477        line = ar.array('H',File.read(2*size))
478        image[row] = np.asarray(line)
479        row += 1
480        pos += 2*size
481    File.close()
482    if imageOnly:
483        return image
484    else:
485        return lines[1:-2],data,size,image
486       
487def GetMAR345Data(filename,imageOnly=False):
488    import array as ar
489    import struct as st
490    try:
491        import pack_f as pf
492    except:
493        msg = wx.MessageDialog(None, message="Unable to load the GSAS MAR image decompression, pack_f",
494                               caption="Import Error",
495                               style=wx.ICON_ERROR | wx.OK | wx.STAY_ON_TOP)
496        msg.ShowModal()
497        return None,None,None,None
498
499    if not imageOnly:
500        print 'Read Mar345 file: ',filename
501    File = open(filename,'rb')
502    head = File.read(4095)
503    numbers = st.unpack('<iiiiiiiiii',head[:40])
504    lines = head[128:].split('\n')
505    head = []
506    for line in lines:
507        line = line.strip()
508        if 'PIXEL' in line:
509            values = line.split()
510            pixel = (int(values[2]),int(values[4]))     #in microns
511        elif 'WAVELENGTH' in line:
512            wave = float(line.split()[1])
513        elif 'DISTANCE' in line:
514            distance = float(line.split()[1])           #in mm
515        elif 'CENTER' in line:
516            values = line.split()
517            center = [float(values[2])/10.,float(values[4])/10.]    #make in mm from pixels
518        if line: 
519            head.append(line)
520    data = {'pixelSize':pixel,'wavelength':wave,'distance':distance,'center':center}
521    for line in head:
522        if 'FORMAT' in line[0:6]:
523            items = line.split()
524            size = int(items[1])
525    pos = 4096
526    File.seek(pos)
527    line = File.read(8)
528    while 'CCP4' not in line:       #get past overflow list for now
529        line = File.read(8)
530        pos += 8
531    pos += 37
532    File.seek(pos)
533    raw = File.read()
534    File.close()
535    image = np.zeros(shape=(size,size),dtype=np.int32)
536    image = pf.pack_f(len(raw),raw,size,image)
537    if imageOnly:
538        return image.T              #transpose to get it right way around
539    else:
540        return head,data,size,image.T
541   
542def GetTifData(filename,imageOnly=False):
543    # only works for APS Perkin-Elmer detector data files in "TIFF" format that are readable by Fit2D
544    import struct as st
545    import array as ar
546    if not imageOnly:
547        print 'Read APS PE-detector tiff file: ',filename
548    File = open(filename,'Ur')
549    dataType = 5
550    try:
551        Meta = open(filename+'.metadata','Ur')
552        head = Meta.readlines()
553        for line in head:
554            line = line.strip()
555            if 'dataType' in line:
556                dataType = int(line.split('=')[1])
557        Meta.close()
558    except IOError:
559        print 'no metadata file found - will try to read file anyway'
560        head = 'no metadata file found'
561    tag = File.read(3)
562    if tag != 'II*':
563        lines = ['not a APS PE-detector tiff file',]
564        return lines,0,0
565    size = st.unpack('<i',File.read(4))[0]
566    image = np.zeros(shape=(size,size),dtype=np.int32)
567    row = 0
568    pos = 8
569    while row < size:
570        File.seek(pos)
571        if dataType == 5:
572            line = ar.array('f',File.read(4*size))
573        else:
574            line = ar.array('l',File.read(4*size))
575        image[row] = np.asarray(line)
576        row += 1
577        pos += 4*size
578    data = {'pixelSize':(200,200),'wavelength':0.10,'distance':100.0,'center':[204.8,204.8]}
579    File.close()   
580    if imageOnly:
581        return image
582    else:
583        return head,data,size,image
584   
585def ProjFileOpen(self):
586    file = open(self.GSASprojectfile,'rb')
587    print 'load from file: ',self.GSASprojectfile
588    wx.BeginBusyCursor()
589    try:
590        while True:
591            try:
592                data = cPickle.load(file)
593            except EOFError:
594                break
595            datum = data[0]
596            print 'load: ',datum[0]
597           
598            #temporary fixes to old project files
599            #fix to convert old style list arrays to numpy arrays
600            if 'PWDR' in datum[0] and 'list' in str(type(datum[1][1][0])):     
601                X = datum[1][1]
602                X = [np.array(X[0]),np.array(X[1]),np.array(X[2]),np.array(X[3]),np.array(X[4]),np.array(X[5])]
603                datum[1] = [datum[1][0],X]
604                print 'powder data converted to numpy arrays'
605            #temporary fix to insert 'PWDR' in front of powder names
606            if 'PKS' not in datum[0] and 'IMG' not in datum[0] and 'SNGL' not in datum[0]:
607                if datum[0] not in ['Notebook','Controls','Phases'] and 'PWDR' not in datum[0]:
608                    datum[0] = 'PWDR '+datum[0]
609                    print 'add PWDR to powder names'
610            #end of temporary fixes
611           
612            Id = self.PatternTree.AppendItem(parent=self.root,text=datum[0])
613            self.PatternTree.SetItemPyData(Id,datum[1])
614            for datus in data[1:]:
615                print '    load: ',datus[0]
616               
617                #temporary fix to add azimuthal angle to instrument parameters
618                if 'PWDR' in datum[0] and 'Instrument Parameters' in datus[0]:
619                    if len(datus[1][0]) == 10 or len(datus[1][0]) == 12:
620                        datus[1][0] += (0.0,)                   #add missing azimuthal angle
621                        datus[1][1].append(0.0)
622                        datus[1][2].append(0.0)
623                        datus[1][3].append('Azimuth')
624                        print 'add azimuth to instrument parameters'
625                #end of temporary fix       
626               
627                sub = self.PatternTree.AppendItem(Id,datus[0])
628                self.PatternTree.SetItemPyData(sub,datus[1])
629               
630            #temporary fix to add Comments to powder data sets
631            if 'PWDR' in datum[0] and not G2gd.GetPatternTreeItemId(self,Id, 'Comments'):
632                print 'no comments - add comments'
633                sub = self.PatternTree.AppendItem(Id,'Comments')
634                self.PatternTree.SetItemPyData(sub,['no comments'])
635            #end of temporary fix
636                               
637            if 'IMG' in datum[0]:                   #retreive image default flag & data if set
638                Data = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Image Controls'))
639                if Data['setDefault']:
640                    self.imageDefault = Data
641                #temporary fix to add masks
642                if not G2gd.GetPatternTreeItemId(self,Id, 'Masks'):
643                    Imin,Imax = Data['range'][0]
644                    Masks = {'Points':[],'Rings':[],'Arcs':[],'Polygons':[],'Thresholds':[(Imin,Imax),[Imin,Imax]]}
645                    self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Masks'),Masks)
646                #end of temporary fix
647               
648        file.close()
649       
650        #temporary fix to add Notebook & Controls to project
651        if not G2gd.GetPatternTreeItemId(self,self.root,'Notebook'):
652            sub = self.PatternTree.AppendItem(parent=self.root,text='Notebook')
653            self.PatternTree.SetItemPyData(sub,[''])
654            sub = self.PatternTree.AppendItem(parent=self.root,text='Controls')
655            self.PatternTree.SetItemPyData(sub,[0])
656            print 'add Notebook and Controls to project'
657           
658    finally:
659        wx.EndBusyCursor()
660    print 'project load successful'
661    self.NewPlot = True
662   
663def ProjFileSave(self):
664    if not self.PatternTree.IsEmpty():
665        file = open(self.GSASprojectfile,'wb')
666        print 'save to file: ',self.GSASprojectfile
667        wx.BeginBusyCursor()
668        try:
669            item, cookie = self.PatternTree.GetFirstChild(self.root)
670            while item:
671                data = []
672                name = self.PatternTree.GetItemText(item)
673                print 'save: ',name
674                data.append([name,self.PatternTree.GetItemPyData(item)])
675                item2, cookie2 = self.PatternTree.GetFirstChild(item)
676                while item2:
677                    name = self.PatternTree.GetItemText(item2)
678                    print '    save: ',name
679                    data.append([name,self.PatternTree.GetItemPyData(item2)])
680                    item2, cookie2 = self.PatternTree.GetNextChild(item, cookie2)                           
681                item, cookie = self.PatternTree.GetNextChild(self.root, cookie)                           
682                cPickle.dump(data,file,1)
683            file.close()
684        finally:
685            wx.EndBusyCursor()
686        print 'project save successful'
687       
688def SaveIntegration(self,PickId,data):
689    azms = self.Integrate[1]
690    X = self.Integrate[2].flatten()[:-1]
691    Xminmax = [X[0],X[-1]]
692    N = len(X)
693    Id = self.PatternTree.GetItemParent(PickId)
694    name = self.PatternTree.GetItemText(Id)
695    name = name.replace('IMG ','PWDR ')
696    Comments = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Comments'))
697    names = ['Type','Lam','Zero','Polariz.','U','V','W','X','Y','SH/L','Azimuth'] 
698    codes = [0 for i in range(11)]
699    parms = ['PXC',data['wavelength'],0.0,0.0,1.0,-1.0,0.3,0.0,1.0,0.0,0.0]
700    Azms = [(azms[i+1]+azms[i])/2. for i in range(len(azms)-1)]
701    for i,azm in enumerate(Azms):
702        item, cookie = self.PatternTree.GetFirstChild(self.root)
703        Id = 0
704        while item:
705            Name = self.PatternTree.GetItemText(item)
706            if name == Name:
707                Id = item
708            item, cookie = self.PatternTree.GetNextChild(self.root, cookie)
709        parms[10] = azm
710        Y = self.Integrate[0][i].flatten()
711        W = np.sqrt(Y)
712        if Id:
713            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Comments'),Comments)                   
714            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Limits'),[tuple(Xminmax),Xminmax])
715            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Background'),[['chebyschev',1,3,1.0,0.0,0.0]])
716            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Instrument Parameters'),[tuple(parms),parms,codes,names])
717            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Peak List'),[])
718            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Index Peak List'),[])
719            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Unit Cells List'),[])             
720        else:
721            Id = self.PatternTree.AppendItem(parent=self.root,text=name+" Azm= %.2f"%(azm))
722            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Comments'),Comments)                   
723            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Limits'),[tuple(Xminmax),Xminmax])
724            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Background'),[['chebyschev',1,3,1.0,0.0,0.0]])
725            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Instrument Parameters'),[tuple(parms),parms,codes,names])
726            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Peak List'),[])
727            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Index Peak List'),[])
728            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Unit Cells List'),[])             
729        self.PatternTree.SetItemPyData(Id,[[''],[np.array(X),np.array(Y),np.array(W),np.zeros(N),np.zeros(N),np.zeros(N)]])
730    self.PatternTree.SelectItem(Id)
731    self.PatternTree.Expand(Id)
732    self.PatternId = Id
733           
734def powderFxyeSave(self,powderfile):
735    file = open(powderfile,'w')
736    prm = open(powderfile.strip('fxye')+'prm','w')      #old style GSAS parm file
737    print 'save powder pattern to file: ',powderfile
738    wx.BeginBusyCursor()
739    Inst = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self, \
740                    self.PickId, 'Instrument Parameters'))[1]
741    if len(Inst) == 11:             #single wavelength
742        lam1 = Inst[1]
743        lam2 = 0.0
744        GU,GV,GW = Inst[4:7]
745        LX,LY = Inst[7:9]
746        SL = HL = Inst[9]/2.0   
747    else:                           #Ka1 & Ka2
748        lam1 = Inst[1]
749        lam2 = Inst[2]
750        GU,GV,GW = Inst[6:9]
751        LX,LY = Inst[9:11]
752        SL = HL = Inst[11]/2.0   
753    prm.write( '            123456789012345678901234567890123456789012345678901234567890        '+'\n')
754    prm.write( 'INS   BANK      1                                                               '+'\n')
755    prm.write( 'INS   HTYPE   PXCR                                                              '+'\n')
756    prm.write(('INS  1 ICONS%10.7f%10.7f    0.0000               0.990    0     0.500   '+'\n')%(lam1,lam2))
757    prm.write( 'INS  1 IRAD     0                                                               '+'\n')
758    prm.write( 'INS  1I HEAD                                                                    '+'\n')
759    prm.write( 'INS  1I ITYP    0    0.0000  180.0000         1                                 '+'\n')
760    prm.write( 'INS  1PRCF1     3    8   0.00100                                                '+'\n')
761    prm.write(('INS  1PRCF11     %15.6g%15.6g%15.6g%15.6g   '+'\n')%(GU,GV,GW,0.0))
762    prm.write(('INS  1PRCF12     %15.6g%15.6g%15.6g%15.6g   '+'\n')%(LX,LY,SL,HL))
763    prm.close()
764    try:
765        x,y,w,yc,yb,yd = self.PatternTree.GetItemPyData(self.PickId)[1]
766        file.write(powderfile+'\n')
767        file.write('BANK 1 %d %d CONS %.2f %.2f 0 0 FXYE\n'%(len(x),len(x),\
768            100.*x[0],100.*(x[1]-x[0])))       
769        XYW = zip(x,y,w)
770        for X,Y,W in XYW:
771            file.write("%15.6g %15.6g %15.6g\n" % (100.*X,Y,W))
772        file.close()
773    finally:
774        wx.EndBusyCursor()
775    print 'powder pattern file written'
776       
777def powderXyeSave(self,powderfile):
778    file = open(powderfile,'w')
779    print 'save powder pattern to file: ',powderfile
780    wx.BeginBusyCursor()
781    try:
782        x,y,w,yc,yb,yd = self.PatternTree.GetItemPyData(self.PickId)[1]
783        XYW = zip(x,y,w)
784        for X,Y,W in XYW:
785            file.write("%15.6g %15.6g %15.6g\n" % (X,Y,W))
786        file.close()
787    finally:
788        wx.EndBusyCursor()
789    print 'powder pattern file written'
790   
791def PeakListSave(self,file,peaks):
792    print 'save peak list to file: ',self.peaklistfile
793    if not peaks:
794        dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
795        try:
796            result = dlg.ShowModal()
797        finally:
798            dlg.Destroy()
799        return
800    for peak in peaks:
801        file.write("%10.4f %12.2f %10.3f %10.3f \n" % \
802            (peak[0],peak[2],peak[4],peak[6]))
803    print 'peak list saved'
804             
805def IndexPeakListSave(self,peaks):
806    file = open(self.peaklistfile,'wa')
807    print 'save index peak list to file: ',self.peaklistfile
808    wx.BeginBusyCursor()
809    try:
810        if not peaks:
811            dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
812            try:
813                result = dlg.ShowModal()
814            finally:
815                dlg.Destroy()
816            return
817        for peak in peaks:
818            file.write("%12.6f\n" % (peak[7]))
819        file.close()
820    finally:
821        wx.EndBusyCursor()
822    print 'index peak list saved'
823   
824def ReadEXPPhase(filename):
825    import GSASIIspc as G2spc
826    import GSASIIlattice as G2lat
827    import GSASIIElem as G2el
828    file = open(filename, 'Ur')
829    Phase = {}
830    S = 1
831    Expr = [{},{},{},{},{},{},{},{},{}]
832    while S:
833        S = file.readline()
834        if 'EXPR NPHAS' in S[:12]:
835            Num = S[12:-1].count('0')
836            NPhas = S[12:-1].split()
837        if 'CRS' in S[:3]:
838            N = int(S[3:4])-1
839            Expr[N][S[:12]] = S[12:-1]
840    file.close()
841    PNames = []
842    for n,N in enumerate(NPhas):
843        if N != '0':
844            result = n
845            key = 'CRS'+str(n+1)+'    PNAM'
846            PNames.append(Expr[n][key])
847    if Num < 8:
848        dlg = wx.SingleChoiceDialog(self, 'Which phase to read?', 'Read phase data', PNames, wx.CHOICEDLG_STYLE)
849        try:
850            if dlg.ShowModal() == wx.ID_OK:
851                result = dlg.GetSelection()
852        finally:
853            dlg.Destroy()       
854    EXPphase = Expr[result]
855    keyList = EXPphase.keys()
856    keyList.sort()
857    SGData = {}
858    if NPhas[result] == '1':
859        Ptype = 'nuclear'
860    elif NPhas[result] in ['2','3']:
861        Ptype = 'magnetic'
862    elif NPhas[result] == '4':
863        Ptype = 'macromolecular'
864    elif NPhas[result] == '10':
865        Ptype = 'Pawley'
866    for key in keyList:
867        if 'PNAM' in key:
868           PhaseName = EXPphase[key].strip()
869        elif 'ABC   ' in key:
870            abc = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                       
871        elif 'ANGLES' in key:
872            angles = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                                               
873        elif 'SG SYM' in key:
874            SpGrp = EXPphase[key][:15].strip()
875            E,SGData = G2spc.SpcGroup(SpGrp)
876    Atoms = []
877    if Ptype == 'nuclear':
878        for key in keyList:
879            if 'AT' in key:
880                if key[11:] == 'A':
881                    S = EXPphase[key]
882                elif key[11:] == 'B':
883                    S += EXPphase[key]
884                    Atom = [S[50:58].strip(),S[:10].strip(),'',
885                        float(S[10:20]),float(S[20:30]),float(S[30:40]),
886                        float(S[40:50]),'',int(S[60:62]),S[130:131]]
887                    if Atom[9] == 'I':
888                        Atom += [float(S[68:78]),0.,0.,0.,0.,0.,0.]
889                    elif Atom[9] == 'A':
890                        Atom += [0.0,float(S[68:78]),float(S[78:88]),
891                            float(S[88:98]),float(S[98:108]),
892                            float(S[108:118]),float(S[118:128])]
893                    XYZ = Atom[3:6]
894                    Atom[7],Atom[8] = G2spc.SytSym(XYZ,SGData)
895                    Atoms.append(Atom)
896    elif Ptype == 'macromolecular':
897        for key in keyList:
898            if 'AT' in key[6:8]:
899                S = EXPphase[key]
900                Atom = [int(S[56:60]),S[50:54].strip().upper(),S[54:56],
901                    S[46:51].strip(),S[:8].strip(),'',
902                    float(S[16:24]),float(S[24:32]),float(S[32:40]),
903                    float(S[8:16]),'1',1,'I',float(S[40:46]),0,0,0,0,0,0]
904                XYZ = Atom[6:9]
905                Atom[10],Atom[11] = G2spc.SytSym(XYZ,SGData)
906                Atoms.append(Atom)
907    Volume = G2lat.calc_V(G2lat.cell2A(abc+angles))
908    Phase['General'] = [PhaseName,Ptype,SGData,[False,]+abc+angles+[Volume,],[False,1.0],
909        0,0,0,0,0]
910    Phase['Atoms'] = Atoms
911    return Phase
912       
913def ReadPDBPhase(filename):
914    import GSASIIspc as G2spc
915    import GSASIIlattice as G2lat
916    import GSASIIElem as G2el
917    EightPiSq = 8.*math.pi**2
918    file = open(filename, 'Ur')
919    Phase = {}
920    Title = ''
921    Compnd = ''
922    Atoms = []
923    A = np.zeros(shape=(3,3))
924    S = file.readline()
925    while S:
926        Atom = []
927        if 'TITLE' in S[:5]:
928            Title = S[10:72].strip()
929            S = file.readline()
930        elif 'COMPND    ' in S[:10]:
931            Compnd = S[10:72].strip()
932            S = file.readline()
933        elif 'CRYST' in S[:5]:
934            abc = S[7:34].split()
935            angles = S[34:55].split()
936            cell=[float(abc[0]),float(abc[1]),float(abc[2]),
937                float(angles[0]),float(angles[1]),float(angles[2])]
938            Volume = G2lat.calc_V(G2lat.cell2A(cell))
939            AA,AB = G2lat.cell2AB(cell)
940            SpGrp = S[55:65]
941            E,SGData = G2spc.SpcGroup(SpGrp)
942            S = file.readline()
943        elif 'SCALE' in S[:5]:
944            V = (S[10:41].split())
945            A[int(S[5])-1] = [float(V[0]),float(V[1]),float(V[2])]
946            S = file.readline()
947        elif 'ATOM' in S[:4] or 'HETATM' in S[:6]:
948            XYZ = [float(S[31:39]),float(S[39:47]),float(S[47:55])]
949            XYZ = np.sum(AA*XYZ,axis=1)
950            SytSym,Mult = G2spc.SytSym(XYZ,SGData)
951            Uiso = float(S[61:67])/EightPiSq
952            Type = S[12:14].upper()
953            if Type[0] in '123456789':
954                Type = Type[1:]
955            Atom = [int(S[22:27]),S[17:20].upper(),S[20:22],
956                S[12:17].strip(),Type.strip(),'',XYZ[0],XYZ[1],XYZ[2],
957                float(S[55:61]),SytSym,Mult,'I',Uiso,0,0,0,0,0,0]
958            S = file.readline()
959            if 'ANISOU' in S[:6]:
960                Uij = S[30:72].split()
961                Uij = [float(Uij[0])/10000.,float(Uij[1])/10000.,float(Uij[2])/10000.,
962                    float(Uij[3])/10000.,float(Uij[4])/10000.,float(Uij[5])/10000.]
963                Atom = Atom[:14]+Uij
964                Atom[12] = 'A'
965                S = file.readline()
966            Atoms.append(Atom)
967        else:           
968            S = file.readline()
969    file.close()
970    if Title:
971        PhaseName = Title
972    elif Compnd:
973        PhaseName = Compnd
974    else:
975        PhaseName = 'None'
976    Phase['General'] = [PhaseName,'macromolecular',SGData,[False,]+cell+[Volume,],[False,1.0],
977        0,0,0,0,0]
978    Phase['Atoms'] = Atoms
979   
980    return Phase
981   
982def ReadCIFAtoms(self,data):
983    print data
Note: See TracBrowser for help on using the repository browser.