source: trunk/GSASIIIO.py @ 193

Last change on this file since 193 was 193, checked in by vondreele, 12 years ago

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