source: trunk/GSASIIIO.py @ 85

Last change on this file since 85 was 85, checked in by vondreel, 12 years ago

threshold masks operational
add index peak list load/reload menu & remove delete key method
plotting of threshold masks
remove "U" in image open command; implement reading of MAR CCD tif images

File size: 37.9 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    File = open(filename,'rb')
547    dataType = 5
548    try:
549        Meta = open(filename+'.metadata','Ur')
550        head = Meta.readlines()
551        for line in head:
552            line = line.strip()
553            if 'dataType' in line:
554                dataType = int(line.split('=')[1])
555        Meta.close()
556    except IOError:
557        print 'no metadata file found - will try to read file anyway'
558        head = 'no metadata file found'
559    tag = File.read(3)
560    if tag != 'II*':
561        lines = ['not a detector tiff file',]
562        return lines,0,0
563    size,Ityp = st.unpack('<ii',File.read(8))
564    if Ityp == 0:
565        tifType = 'Pilatus'
566        pixy = (172,172)
567        pos = 4096
568        if not imageOnly:
569            print 'Read Pilatus tiff file: ',filename
570    elif Ityp == 1:
571        tifType = 'PE'
572        pixy = (200,200)
573        pos = 8
574        if not imageOnly:
575            print 'Read APS PE-detector tiff file: ',filename
576    elif Ityp == 3328:
577        tifType = 'MAR'
578        pixy = (79,79)
579        pos = 4096
580        if not imageOnly:
581            print 'Read MAR CCD tiff file: ',filename
582    else:
583        lines = 'unknown tif type'
584        return lines,0,0
585    image = np.zeros(shape=(size,size),dtype=np.int32)
586    row = 0
587    while row < size:
588        File.seek(pos)
589        if 'PE' in tifType: 
590            if dataType == 5:
591                line = ar.array('f',File.read(4*size))
592            else:
593                line = ar.array('l',File.read(4*size))
594            pos += 4*size
595        elif 'MAR' in tifType:
596            line = ar.array('H',File.read(2*size))
597            pos += 2*size
598        image[row] = np.asarray(line)
599        row += 1
600    data = {'pixelSize':pixy,'wavelength':0.10,'distance':100.0,'center':[204.8,204.8]}
601    File.close()   
602    if imageOnly:
603        return image
604    else:
605        return head,data,size,image
606   
607def ProjFileOpen(self):
608    file = open(self.GSASprojectfile,'rb')
609    print 'load from file: ',self.GSASprojectfile
610    wx.BeginBusyCursor()
611    try:
612        while True:
613            try:
614                data = cPickle.load(file)
615            except EOFError:
616                break
617            datum = data[0]
618            print 'load: ',datum[0]
619           
620            #temporary fixes to old project files
621            #fix to convert old style list arrays to numpy arrays
622            if 'PWDR' in datum[0] and 'list' in str(type(datum[1][1][0])):     
623                X = datum[1][1]
624                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])]
625                datum[1] = [datum[1][0],X]
626                print 'powder data converted to numpy arrays'
627            #temporary fix to insert 'PWDR' in front of powder names
628            if 'PKS' not in datum[0] and 'IMG' not in datum[0] and 'SNGL' not in datum[0]:
629                if datum[0] not in ['Notebook','Controls','Phases'] and 'PWDR' not in datum[0]:
630                    datum[0] = 'PWDR '+datum[0]
631                    print 'add PWDR to powder names'
632            #end of temporary fixes
633           
634            Id = self.PatternTree.AppendItem(parent=self.root,text=datum[0])
635            self.PatternTree.SetItemPyData(Id,datum[1])
636            for datus in data[1:]:
637                print '    load: ',datus[0]
638               
639                #temporary fix to add azimuthal angle to instrument parameters
640                if 'PWDR' in datum[0] and 'Instrument Parameters' in datus[0]:
641                    if len(datus[1][0]) == 10 or len(datus[1][0]) == 12:
642                        datus[1][0] += (0.0,)                   #add missing azimuthal angle
643                        datus[1][1].append(0.0)
644                        datus[1][2].append(0.0)
645                        datus[1][3].append('Azimuth')
646                        print 'add azimuth to instrument parameters'
647                #end of temporary fix       
648               
649                sub = self.PatternTree.AppendItem(Id,datus[0])
650                self.PatternTree.SetItemPyData(sub,datus[1])
651               
652            #temporary fix to add Comments to powder data sets
653            if 'PWDR' in datum[0] and not G2gd.GetPatternTreeItemId(self,Id, 'Comments'):
654                print 'no comments - add comments'
655                sub = self.PatternTree.AppendItem(Id,'Comments')
656                self.PatternTree.SetItemPyData(sub,['no comments'])
657            #end of temporary fix
658                               
659            if 'IMG' in datum[0]:                   #retreive image default flag & data if set
660                Data = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Image Controls'))
661                if Data['setDefault']:
662                    self.imageDefault = Data
663                #temporary fix to add masks
664                if not G2gd.GetPatternTreeItemId(self,Id, 'Masks'):
665                    Imin,Imax = Data['range'][0]
666                    Masks = {'Points':[],'Rings':[],'Arcs':[],'Polygons':[],'Thresholds':[(Imin,Imax),[Imin,Imax]]}
667                    self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Masks'),Masks)
668                #end of temporary fix
669               
670        file.close()
671       
672        #temporary fix to add Notebook & Controls to project
673        if not G2gd.GetPatternTreeItemId(self,self.root,'Notebook'):
674            sub = self.PatternTree.AppendItem(parent=self.root,text='Notebook')
675            self.PatternTree.SetItemPyData(sub,[''])
676            sub = self.PatternTree.AppendItem(parent=self.root,text='Controls')
677            self.PatternTree.SetItemPyData(sub,[0])
678            print 'add Notebook and Controls to project'
679           
680    finally:
681        wx.EndBusyCursor()
682    print 'project load successful'
683    self.NewPlot = True
684   
685def ProjFileSave(self):
686    if not self.PatternTree.IsEmpty():
687        file = open(self.GSASprojectfile,'wb')
688        print 'save to file: ',self.GSASprojectfile
689        wx.BeginBusyCursor()
690        try:
691            item, cookie = self.PatternTree.GetFirstChild(self.root)
692            while item:
693                data = []
694                name = self.PatternTree.GetItemText(item)
695                print 'save: ',name
696                data.append([name,self.PatternTree.GetItemPyData(item)])
697                item2, cookie2 = self.PatternTree.GetFirstChild(item)
698                while item2:
699                    name = self.PatternTree.GetItemText(item2)
700                    print '    save: ',name
701                    data.append([name,self.PatternTree.GetItemPyData(item2)])
702                    item2, cookie2 = self.PatternTree.GetNextChild(item, cookie2)                           
703                item, cookie = self.PatternTree.GetNextChild(self.root, cookie)                           
704                cPickle.dump(data,file,1)
705            file.close()
706        finally:
707            wx.EndBusyCursor()
708        print 'project save successful'
709       
710def SaveIntegration(self,PickId,data):
711    azms = self.Integrate[1]
712    X = self.Integrate[2].flatten()[:-1]
713    Xminmax = [X[0],X[-1]]
714    N = len(X)
715    Id = self.PatternTree.GetItemParent(PickId)
716    name = self.PatternTree.GetItemText(Id)
717    name = name.replace('IMG ','PWDR ')
718    Comments = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Comments'))
719    names = ['Type','Lam','Zero','Polariz.','U','V','W','X','Y','SH/L','Azimuth'] 
720    codes = [0 for i in range(11)]
721    parms = ['PXC',data['wavelength'],0.0,0.0,1.0,-1.0,0.3,0.0,1.0,0.0,0.0]
722    Azms = [(azms[i+1]+azms[i])/2. for i in range(len(azms)-1)]
723    for i,azm in enumerate(Azms):
724        item, cookie = self.PatternTree.GetFirstChild(self.root)
725        Id = 0
726        while item:
727            Name = self.PatternTree.GetItemText(item)
728            if name == Name:
729                Id = item
730            item, cookie = self.PatternTree.GetNextChild(self.root, cookie)
731        parms[10] = azm
732        Y = self.Integrate[0][i].flatten()
733        W = np.sqrt(Y)
734        if Id:
735            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Comments'),Comments)                   
736            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Limits'),[tuple(Xminmax),Xminmax])
737            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Background'),[['chebyschev',1,3,1.0,0.0,0.0]])
738            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Instrument Parameters'),[tuple(parms),parms,codes,names])
739            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Peak List'),[])
740            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Index Peak List'),[])
741            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Unit Cells List'),[])             
742        else:
743            Id = self.PatternTree.AppendItem(parent=self.root,text=name+" Azm= %.2f"%(azm))
744            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Comments'),Comments)                   
745            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Limits'),[tuple(Xminmax),Xminmax])
746            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Background'),[['chebyschev',1,3,1.0,0.0,0.0]])
747            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Instrument Parameters'),[tuple(parms),parms,codes,names])
748            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Peak List'),[])
749            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Index Peak List'),[])
750            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Unit Cells List'),[])             
751        self.PatternTree.SetItemPyData(Id,[[''],[np.array(X),np.array(Y),np.array(W),np.zeros(N),np.zeros(N),np.zeros(N)]])
752    self.PatternTree.SelectItem(Id)
753    self.PatternTree.Expand(Id)
754    self.PatternId = Id
755           
756def powderFxyeSave(self,powderfile):
757    file = open(powderfile,'w')
758    prm = open(powderfile.strip('fxye')+'prm','w')      #old style GSAS parm file
759    print 'save powder pattern to file: ',powderfile
760    wx.BeginBusyCursor()
761    Inst = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self, \
762                    self.PickId, 'Instrument Parameters'))[1]
763    if len(Inst) == 11:             #single wavelength
764        lam1 = Inst[1]
765        lam2 = 0.0
766        GU,GV,GW = Inst[4:7]
767        LX,LY = Inst[7:9]
768        SL = HL = Inst[9]/2.0   
769    else:                           #Ka1 & Ka2
770        lam1 = Inst[1]
771        lam2 = Inst[2]
772        GU,GV,GW = Inst[6:9]
773        LX,LY = Inst[9:11]
774        SL = HL = Inst[11]/2.0   
775    prm.write( '            123456789012345678901234567890123456789012345678901234567890        '+'\n')
776    prm.write( 'INS   BANK      1                                                               '+'\n')
777    prm.write( 'INS   HTYPE   PXCR                                                              '+'\n')
778    prm.write(('INS  1 ICONS%10.7f%10.7f    0.0000               0.990    0     0.500   '+'\n')%(lam1,lam2))
779    prm.write( 'INS  1 IRAD     0                                                               '+'\n')
780    prm.write( 'INS  1I HEAD                                                                    '+'\n')
781    prm.write( 'INS  1I ITYP    0    0.0000  180.0000         1                                 '+'\n')
782    prm.write( 'INS  1PRCF1     3    8   0.00100                                                '+'\n')
783    prm.write(('INS  1PRCF11     %15.6g%15.6g%15.6g%15.6g   '+'\n')%(GU,GV,GW,0.0))
784    prm.write(('INS  1PRCF12     %15.6g%15.6g%15.6g%15.6g   '+'\n')%(LX,LY,SL,HL))
785    prm.close()
786    try:
787        x,y,w,yc,yb,yd = self.PatternTree.GetItemPyData(self.PickId)[1]
788        file.write(powderfile+'\n')
789        file.write('BANK 1 %d %d CONS %.2f %.2f 0 0 FXYE\n'%(len(x),len(x),\
790            100.*x[0],100.*(x[1]-x[0])))       
791        XYW = zip(x,y,w)
792        for X,Y,W in XYW:
793            file.write("%15.6g %15.6g %15.6g\n" % (100.*X,Y,W))
794        file.close()
795    finally:
796        wx.EndBusyCursor()
797    print 'powder pattern file written'
798       
799def powderXyeSave(self,powderfile):
800    file = open(powderfile,'w')
801    print 'save powder pattern to file: ',powderfile
802    wx.BeginBusyCursor()
803    try:
804        x,y,w,yc,yb,yd = self.PatternTree.GetItemPyData(self.PickId)[1]
805        XYW = zip(x,y,w)
806        for X,Y,W in XYW:
807            file.write("%15.6g %15.6g %15.6g\n" % (X,Y,W))
808        file.close()
809    finally:
810        wx.EndBusyCursor()
811    print 'powder pattern file written'
812   
813def PeakListSave(self,file,peaks):
814    print 'save peak list to file: ',self.peaklistfile
815    if not peaks:
816        dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
817        try:
818            result = dlg.ShowModal()
819        finally:
820            dlg.Destroy()
821        return
822    for peak in peaks:
823        file.write("%10.4f %12.2f %10.3f %10.3f \n" % \
824            (peak[0],peak[2],peak[4],peak[6]))
825    print 'peak list saved'
826             
827def IndexPeakListSave(self,peaks):
828    file = open(self.peaklistfile,'wa')
829    print 'save index peak list to file: ',self.peaklistfile
830    wx.BeginBusyCursor()
831    try:
832        if not peaks:
833            dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
834            try:
835                result = dlg.ShowModal()
836            finally:
837                dlg.Destroy()
838            return
839        for peak in peaks:
840            file.write("%12.6f\n" % (peak[7]))
841        file.close()
842    finally:
843        wx.EndBusyCursor()
844    print 'index peak list saved'
845   
846def ReadEXPPhase(filename):
847    import GSASIIspc as G2spc
848    import GSASIIlattice as G2lat
849    import GSASIIElem as G2el
850    file = open(filename, 'Ur')
851    Phase = {}
852    S = 1
853    Expr = [{},{},{},{},{},{},{},{},{}]
854    while S:
855        S = file.readline()
856        if 'EXPR NPHAS' in S[:12]:
857            Num = S[12:-1].count('0')
858            NPhas = S[12:-1].split()
859        if 'CRS' in S[:3]:
860            N = int(S[3:4])-1
861            Expr[N][S[:12]] = S[12:-1]
862    file.close()
863    PNames = []
864    for n,N in enumerate(NPhas):
865        if N != '0':
866            result = n
867            key = 'CRS'+str(n+1)+'    PNAM'
868            PNames.append(Expr[n][key])
869    if Num < 8:
870        dlg = wx.SingleChoiceDialog(self, 'Which phase to read?', 'Read phase data', PNames, wx.CHOICEDLG_STYLE)
871        try:
872            if dlg.ShowModal() == wx.ID_OK:
873                result = dlg.GetSelection()
874        finally:
875            dlg.Destroy()       
876    EXPphase = Expr[result]
877    keyList = EXPphase.keys()
878    keyList.sort()
879    SGData = {}
880    if NPhas[result] == '1':
881        Ptype = 'nuclear'
882    elif NPhas[result] in ['2','3']:
883        Ptype = 'magnetic'
884    elif NPhas[result] == '4':
885        Ptype = 'macromolecular'
886    elif NPhas[result] == '10':
887        Ptype = 'Pawley'
888    for key in keyList:
889        if 'PNAM' in key:
890           PhaseName = EXPphase[key].strip()
891        elif 'ABC   ' in key:
892            abc = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                       
893        elif 'ANGLES' in key:
894            angles = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                                               
895        elif 'SG SYM' in key:
896            SpGrp = EXPphase[key][:15].strip()
897            E,SGData = G2spc.SpcGroup(SpGrp)
898    Atoms = []
899    if Ptype == 'nuclear':
900        for key in keyList:
901            if 'AT' in key:
902                if key[11:] == 'A':
903                    S = EXPphase[key]
904                elif key[11:] == 'B':
905                    S += EXPphase[key]
906                    Atom = [S[50:58].strip(),S[:10].strip(),'',
907                        float(S[10:20]),float(S[20:30]),float(S[30:40]),
908                        float(S[40:50]),'',int(S[60:62]),S[130:131]]
909                    if Atom[9] == 'I':
910                        Atom += [float(S[68:78]),0.,0.,0.,0.,0.,0.]
911                    elif Atom[9] == 'A':
912                        Atom += [0.0,float(S[68:78]),float(S[78:88]),
913                            float(S[88:98]),float(S[98:108]),
914                            float(S[108:118]),float(S[118:128])]
915                    XYZ = Atom[3:6]
916                    Atom[7],Atom[8] = G2spc.SytSym(XYZ,SGData)
917                    Atoms.append(Atom)
918    elif Ptype == 'macromolecular':
919        for key in keyList:
920            if 'AT' in key[6:8]:
921                S = EXPphase[key]
922                Atom = [int(S[56:60]),S[50:54].strip().upper(),S[54:56],
923                    S[46:51].strip(),S[:8].strip(),'',
924                    float(S[16:24]),float(S[24:32]),float(S[32:40]),
925                    float(S[8:16]),'1',1,'I',float(S[40:46]),0,0,0,0,0,0]
926                XYZ = Atom[6:9]
927                Atom[10],Atom[11] = G2spc.SytSym(XYZ,SGData)
928                Atoms.append(Atom)
929    Volume = G2lat.calc_V(G2lat.cell2A(abc+angles))
930    Phase['General'] = [PhaseName,Ptype,SGData,[False,]+abc+angles+[Volume,],[False,1.0],
931        0,0,0,0,0]
932    Phase['Atoms'] = Atoms
933    return Phase
934       
935def ReadPDBPhase(filename):
936    import GSASIIspc as G2spc
937    import GSASIIlattice as G2lat
938    import GSASIIElem as G2el
939    EightPiSq = 8.*math.pi**2
940    file = open(filename, 'Ur')
941    Phase = {}
942    Title = ''
943    Compnd = ''
944    Atoms = []
945    A = np.zeros(shape=(3,3))
946    S = file.readline()
947    while S:
948        Atom = []
949        if 'TITLE' in S[:5]:
950            Title = S[10:72].strip()
951            S = file.readline()
952        elif 'COMPND    ' in S[:10]:
953            Compnd = S[10:72].strip()
954            S = file.readline()
955        elif 'CRYST' in S[:5]:
956            abc = S[7:34].split()
957            angles = S[34:55].split()
958            cell=[float(abc[0]),float(abc[1]),float(abc[2]),
959                float(angles[0]),float(angles[1]),float(angles[2])]
960            Volume = G2lat.calc_V(G2lat.cell2A(cell))
961            AA,AB = G2lat.cell2AB(cell)
962            SpGrp = S[55:65]
963            E,SGData = G2spc.SpcGroup(SpGrp)
964            S = file.readline()
965        elif 'SCALE' in S[:5]:
966            V = (S[10:41].split())
967            A[int(S[5])-1] = [float(V[0]),float(V[1]),float(V[2])]
968            S = file.readline()
969        elif 'ATOM' in S[:4] or 'HETATM' in S[:6]:
970            XYZ = [float(S[31:39]),float(S[39:47]),float(S[47:55])]
971            XYZ = np.sum(AA*XYZ,axis=1)
972            SytSym,Mult = G2spc.SytSym(XYZ,SGData)
973            Uiso = float(S[61:67])/EightPiSq
974            Type = S[12:14].upper()
975            if Type[0] in '123456789':
976                Type = Type[1:]
977            Atom = [int(S[22:27]),S[17:20].upper(),S[20:22],
978                S[12:17].strip(),Type.strip(),'',XYZ[0],XYZ[1],XYZ[2],
979                float(S[55:61]),SytSym,Mult,'I',Uiso,0,0,0,0,0,0]
980            S = file.readline()
981            if 'ANISOU' in S[:6]:
982                Uij = S[30:72].split()
983                Uij = [float(Uij[0])/10000.,float(Uij[1])/10000.,float(Uij[2])/10000.,
984                    float(Uij[3])/10000.,float(Uij[4])/10000.,float(Uij[5])/10000.]
985                Atom = Atom[:14]+Uij
986                Atom[12] = 'A'
987                S = file.readline()
988            Atoms.append(Atom)
989        else:           
990            S = file.readline()
991    file.close()
992    if Title:
993        PhaseName = Title
994    elif Compnd:
995        PhaseName = Compnd
996    else:
997        PhaseName = 'None'
998    Phase['General'] = [PhaseName,'macromolecular',SGData,[False,]+cell+[Volume,],[False,1.0],
999        0,0,0,0,0]
1000    Phase['Atoms'] = Atoms
1001   
1002    return Phase
1003   
1004def ReadCIFAtoms(self,data):
1005    print data
Note: See TracBrowser for help on using the repository browser.