source: trunk/GSASIIIO.py @ 230

Last change on this file since 230 was 230, checked in by vondreele, 13 years ago

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