source: trunk/GSASIIIO.py @ 37

Last change on this file since 37 was 37, checked in by vondreel, 13 years ago

import sys missing

File size: 29.0 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 GSASIIgrid as G2gd
11# needed for MAR345
12import os.path as ospath
13# determine a binary path pased on the host OS and the python version, path is relative to
14# location of this file
15if sys.platform == "win32":
16    bindir = 'binwin%d.%d' % sys.version_info[0:2]
17elif sys.platform == "darwin":
18    bindir = 'binmac%d.%d' % sys.version_info[0:2]
19else:
20    bindir = 'bin'
21if ospath.exists(ospath.join(sys.path[0],bindir)) and ospath.join(sys.path[0],bindir) not in sys.path: 
22    sys.path.insert(0,ospath.join(sys.path[0],bindir))
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       
56    self.IparmName = GetInstrumentFile(self,filename)
57    if self.IparmName:
58        Iparm = GetInstrumentData(self.IparmName)
59    else:
60        Iparm = {}                                               #Assume CuKa lab data if no iparm file
61        Iparm['INS   HTYPE '] = 'PXC '
62        Iparm['INS  1 ICONS'] = '  1.540500  1.544300       0.0         0       0.7    0       0.5   '
63        Iparm['INS  1PRCF1 '] = '    3    8      0.01                                                '
64        Iparm['INS  1PRCF11'] = '   2.000000E+00  -2.000000E+00   5.000000E+00   0.000000E+00        '
65        Iparm['INS  1PRCF12'] = '   0.000000E+00   0.000000E+00   0.150000E-01   0.150000E-01        '
66    S = 1
67    Banks = []
68    Pos = []
69    FoundData = []
70    Comments = []
71    wx.BeginBusyCursor()
72    try:
73        while S:
74            S = File.readline()
75            if S[:1] != '#':
76                if S[:4] == 'BANK':
77                    Banks.append(S)
78                    Pos.append(File.tell())
79            else:
80                Comments.append(S[:-1])
81        File.close()
82    finally:
83        wx.EndBusyCursor()
84    if Comments:
85       print 'Comments on file:'
86       for Comment in Comments: print Comment
87    if Banks:
88        result = [0]
89        if len(Banks) >= 2:
90            dlg = wx.MultiChoiceDialog(self, 'Which scans do you want?', 'Select scans', Banks, wx.CHOICEDLG_STYLE)
91            try:
92                if dlg.ShowModal() == wx.ID_OK:
93                    result = dlg.GetSelections()
94                else:
95                    result = []
96            finally:
97                dlg.Destroy()
98        for i in result:
99            FoundData.append((filename,Pos[i],Banks[i]))
100    else:
101        dlg = wx.MessageDialog(self, 'ERROR - this is not a GSAS powder data file', 'No BANK records', wx.OK | wx.ICON_ERROR)
102        try:
103            result = dlg.ShowModal()
104        finally:
105            dlg.Destroy()
106    return FoundData,Iparm,Comments
107
108def GetInstrumentFile(self,filename):
109    import os.path as op
110    dlg = wx.FileDialog(self,'Choose an instrument file','.', '', 'GSAS iparm file (*.prm)|*.prm|All files(*.*)|*.*', wx.OPEN)
111    if self.dirname: 
112        dlg.SetDirectory(self.dirname)
113        Tname = filename[:filename.index('.')]+'.prm'
114        if op.exists(Tname):
115            self.IparmName = Tname       
116    if self.IparmName: dlg.SetFilename(self.IparmName)
117    filename = ''
118    try:
119        if dlg.ShowModal() == wx.ID_OK:
120            filename = dlg.GetPath()
121    finally:
122        dlg.Destroy()
123    return filename
124
125def GetInstrumentData(IparmName):
126    file = open(IparmName, 'Ur')
127    S = 1
128    Iparm = {}
129    while S:
130        S = file.readline()
131        Iparm[S[:12]] = S[12:-1]
132    return Iparm
133   
134def GetPowderPeaks(fileName):
135    sind = lambda x: math.sin(x*math.pi/180.)
136    asind = lambda x: 180.*math.asin(x)/math.pi
137    Cuka = 1.54052
138    File = open(fileName,'Ur')
139    Comments = []
140    peaks = []
141    S = File.readline()
142    while S:
143        if S[:1] == '#':
144            Comments.append(S[:-1])
145        else:
146            item = S.split()
147            if len(item) == 1:
148                peaks.append([float(item[0]),1.0])
149            elif len(item) > 1:
150                peaks.append([float(item[0]),float(item[0])])
151        S = File.readline()
152    File.close()
153    if Comments:
154       print 'Comments on file:'
155       for Comment in Comments: print Comment
156    Peaks = []
157    if peaks[0][0] > peaks[-1][0]:          # d-spacings - assume CuKa
158        for peak in peaks:
159            dsp = peak[0]
160            sth = Cuka/(2.0*dsp)
161            if sth < 1.0:
162                tth = 2.0*asind(sth)
163            else:
164                break
165            Peaks.append([tth,peak[1],True,False,0,0,0,dsp,0.0])
166    else:                                   #2-thetas - assume Cuka (for now)
167        for peak in peaks:
168            tth = peak[0]
169            dsp = Cuka/(2.0*sind(tth/2.0))
170            Peaks.append([tth,peak[1],True,False,0,0,0,dsp,0.0])
171    return Comments,Peaks
172   
173def GetHKLData(filename):
174    print 'Reading: '+filename
175    File = open(filename,'Ur')
176    HKLref = []
177    HKLmin = [1000,1000,1000]
178    HKLmax = [0,0,0]
179    FoMax = 0
180    ifFc = False
181    S = File.readline()
182    while '#' in S[0]:        #get past comments if any
183        S = File.readline()       
184    if '_' in S:         #cif style .hkl file
185        while 'loop_' not in S:         #skip preliminaries if any - can't have 'loop_' in them!
186            S = File.readline()       
187        S = File.readline()             #get past 'loop_' line
188        pos = 0
189        hpos = kpos = lpos = Fosqpos = Fcsqpos = sigpos = -1
190        while S:
191            if '_' in S:
192                if 'index_h' in S:
193                    hpos = pos
194                elif 'index_k' in S:
195                    kpos = pos
196                elif 'index_l' in S:
197                    lpos = pos
198                elif 'F_squared_meas' in S:
199                    Fosqpos = pos
200                elif 'F_squared_calc' in S:
201                    Fcsqpos = pos
202                elif 'F_squared_sigma' in S:
203                    sigpos = pos
204                pos += 1
205            else:
206                data = S.split()
207                if data:                    #avoid blank lines
208                    h = int(data[hpos])
209                    k = int(data[kpos])
210                    l = int(data[lpos])
211                    Fosq = float(data[Fosqpos])
212                    if sigpos != -1:
213                        sigFosq = float(data[sigpos])
214                    else:
215                        sigFosq = 1.
216                    if Fcsqpos != -1:
217                        Fcsq = float(data[Fcsqpos])
218                        if Fcsq:
219                            ifFc = True
220                    else:
221                        Fcsq = 0.
222                       
223                    HKLmin = [min(h,HKLmin[0]),min(k,HKLmin[1]),min(l,HKLmin[2])]
224                    HKLmax = [max(h,HKLmax[0]),max(k,HKLmax[1]),max(l,HKLmax[2])]
225                    FoMax = max(FoMax,Fosq)
226                    HKLref.append([h,k,l,Fosq,sigFosq,Fcsq,0,0,0])                 #room for Fc, Fcp, Fcpp & phase
227            S = File.readline()
228    else:                   #dumb h,k,l,Fo,sigFo .hkl file
229        while S:
230            h,k,l,Fo,sigFo = S.split()
231            h = int(h)
232            k = int(k)
233            l = int(l)
234            Fo = float(Fo)
235            sigFo = float(sigFo)
236            HKLmin = [min(h,HKLmin[0]),min(k,HKLmin[1]),min(l,HKLmin[2])]
237            HKLmax = [max(h,HKLmax[0]),max(k,HKLmax[1]),max(l,HKLmax[2])]
238            FoMax = max(FoMax,Fo)
239            HKLref.append([h,k,l,Fo**2,2.*Fo*sigFo,0,0,0,0])                 #room for Fc, Fcp, Fcpp & phase
240            S = File.readline()
241    File.close()
242    return HKLref,HKLmin,HKLmax,FoMax,ifFc
243
244def GetPowderData(filename,Pos,Bank,DataType):
245    '''Reads one BANK of data from GSAS raw powder data file
246    input:
247    filename: GSAS raw powder file dataname
248    Pos: start of data in file just after BANK record
249    Bank: the BANK record
250    DataType: powder data type, e.g. "PXC" for Powder X-ray CW data
251    returns: list [x,y,e,yc,yb]
252    x: array of x-axis values
253    y: array of powder pattern intensities
254    w: array of w=sig(intensity)^2 values
255    yc: array of calc. intensities (zero)
256    yb: array of calc. background (zero)
257    yd: array of obs-calc profiles
258    '''
259    print 'Reading: '+filename
260    print 'Bank:    '+Bank[:-1]
261    if 'FXYE' in Bank:
262        return GetFXYEdata(filename,Pos,Bank,DataType)
263    elif 'FXY' in Bank:
264        return GetFXYdata(filename,Pos,Bank,DataType)
265    elif 'ESD' in Bank:
266        return GetESDdata(filename,Pos,Bank,DataType)
267    elif 'STD' in Bank:
268        return GetSTDdata(filename,Pos,Bank,DataType)
269    else:
270        return GetSTDdata(filename,Pos,Bank,DataType)
271    return []
272
273def GetFXYEdata(filename,Pos,Bank,DataType):
274    File = open(filename,'Ur')
275    File.seek(Pos)
276    x = []
277    y = []
278    w = []
279    yc = []
280    yb = []
281    yd = []
282    S = File.readline()
283    while S and S[:4] != 'BANK':
284        vals = S.split()
285        if DataType[2] == 'C':
286            x.append(float(vals[0])/100.)               #CW: from centidegrees to degrees
287        elif DataType[2] == 'T':
288            x.append(float(vals[0])/1000.0)             #TOF: from musec to millisec
289        f = float(vals[1])
290        if f <= 0.0:
291            y.append(0.0)
292            w.append(1.0)
293        else:
294            y.append(float(vals[1]))
295            w.append(1.0/float(vals[2])**2)
296        yc.append(0.0)
297        yb.append(0.0)
298        yd.append(0.0)
299        S = File.readline()
300    File.close()
301    return [x,y,w,yc,yb,yd]
302   
303def GetFXYdata(filename,Pos,Bank,DataType):
304    File = open(filename,'Ur')
305    File.seek(Pos)
306    x = []
307    y = []
308    w = []
309    yc = []
310    yb = []
311    yd = []
312    S = File.readline()
313    while S and S[:4] != 'BANK':
314        vals = S.split()
315        if DataType[2] == 'C':
316            x.append(float(vals[0])/100.)               #CW: from centidegrees to degrees
317        elif DataType[2] == 'T':
318            x.append(float(vals[0])/1000.0)             #TOF: from musec to millisec
319        f = float(vals[1])
320        if f > 0.0:
321            y.append(f)
322            w.append(1.0/f)
323        else:             
324            y.append(0.0)
325            w.append(1.0)
326        S = File.readline()
327        yc.append(0.0)
328        yb.append(0.0)
329        yd.append(0.0)
330    File.close()
331    return [x,y,w,yc,yb,yd]
332   
333def GetESDdata(filename,Pos,Bank,DataType):
334    File = open(filename,'Ur')
335    cons = Bank.split()
336    if DataType[2] == 'C':
337        start = float(cons[5])/100.0               #CW: from centidegrees to degrees
338        step = float(cons[6])/100.0
339    elif DataType[2] == 'T':
340        start = float(cons[5])/1000.0              #TOF: from musec to millisec
341        step = float(cons[6])/1000.0
342    File.seek(Pos)
343    x = []
344    y = []
345    w = []
346    yc = []
347    yb = []
348    yd = []
349    S = File.readline()
350    j = 0
351    while S and S[:4] != 'BANK':
352        for i in range(0,80,16):
353            xi = start+step*j
354            yi = sfloat(S[i:i+8])
355            ei = sfloat(S[i+8:i+16])
356            x.append(xi)
357            if yi > 0.0:
358                y.append(yi)
359                w.append(1.0/ei**2)
360            else:             
361                y.append(0.0)
362                w.append(1.0)
363            yc.append(0.0)
364            yb.append(0.0)
365            yd.append(0.0)
366            j += 1
367        S = File.readline()
368    File.close()
369    return [x,y,w,yc,yb,yd]
370
371def GetSTDdata(filename,Pos,Bank,DataType):
372    File = open(filename,'Ur')
373    cons = Bank.split()
374    Nch = cons[2]
375    if DataType[2] == 'C':
376        start = float(cons[5])/100.0               #CW: from centidegrees to degrees
377        step = float(cons[6])/100.0
378    elif DataType[2] == 'T':
379        start = float(cons[5])/1000.0              #TOF: from musec to millisec - not likely!
380        step = float(cons[6])/1000.0
381    File.seek(Pos)
382    x = []
383    y = []
384    w = []
385    yc = []
386    yb = []
387    yd = []
388    S = File.readline()
389    j = 0
390    while S and S[:4] != 'BANK':
391        for i in range(0,80,8):
392            xi = start+step*j
393            ni = max(sint(S[i:i+2]),1)
394            yi = max(sfloat(S[i+2:i+8]),0.0)
395            if yi:
396                ei = math.sqrt(yi*ni)
397            else:
398                yi = 0.0
399                ei = 1.0
400            j += 1
401            if j < Nch:
402                x.append(xi)
403                y.append(yi)
404                w.append(1.0/ei**2)
405                yc.append(0.0)
406                yb.append(0.0)
407                yd.append(0.0)
408        S = File.readline()
409    File.close()
410    return [x,y,w,yc,yb,yd]
411   
412def GetGEsumData(filename):
413    import array as ar
414    print 'Read GE sum file: ',filename   
415    File = open(filename,'rb')
416    size = 2048
417    if '.sum' in filename:
418        head = ['GE detector sum data from APS 1-ID',]
419    if '.avg' in filename:
420        head = ['GE detector avg data from APS 1-ID',]
421    image = np.zeros(shape=(size,size),dtype=np.int32)
422    row = 0
423    pos = 0
424    while row < size:
425        File.seek(pos)
426        if '.sum' in filename:
427            line = ar.array('f',File.read(4*size))
428            pos += 4*size
429        elif '.avg' in filename:
430            line = ar.array('H',File.read(2*size))
431            pos += 2*size
432        image[row] = np.asarray(line)
433        row += 1
434    data = {'pixelSize':(200,200),'wavelength':0.15,'distance':250.0,'center':[204.8,204.8]} 
435    return head,data,size,image
436    File.close()   
437       
438def GetImgData(filename):
439    import struct as st
440    import array as ar
441    print 'Read ADSC img file: ',filename
442    File = open(filename,'rb')
443    head = File.read(511)
444    lines = head.split('\n')
445    head = []
446    center = [0,0]
447    for line in lines[1:-2]:
448        line = line.strip()[:-1]
449        if line:
450            if 'SIZE1' in line:
451                size = int(line.split('=')[1])
452            elif 'WAVELENGTH' in line:
453                wave = float(line.split('=')[1])
454            elif 'BIN' in line:
455                if line.split('=')[1] == '2x2':
456                    pixel=(102,102)
457                else:
458                    pixel = (51,51)
459            elif 'DISTANCE' in line:
460                distance = float(line.split('=')[1])
461            elif 'CENTER_X' in line:
462                center[0] = float(line.split('=')[1])
463            elif 'CENTER_Y' in line:
464                center[1] = float(line.split('=')[1])
465            head.append(line)
466    data = {'pixelSize':pixel,'wavelength':wave,'distance':distance,'center':center}
467    image = []
468    row = 0
469    pos = 512
470    image = np.zeros(shape=(size,size),dtype=np.int32)   
471    while row < size:
472        File.seek(pos)
473        line = ar.array('H',File.read(2*size))
474        image[row] = np.asarray(line)
475        if row == 3390: print len(line)
476        row += 1
477        pos += 2*size
478    File.close()
479    return lines[1:-2],data,size,image
480       
481def GetMAR345Data(filename):
482    import array as ar
483    import struct as st
484    try:
485        import pack_f as pf
486    except:
487        msg = wx.MessageDialog(None, message="Unable to load the GSAS MAR image decompression, pack_f",
488                               caption="Import Error",
489                               style=wx.ICON_ERROR | wx.OK | wx.STAY_ON_TOP)
490        msg.ShowModal()
491        return None,None,None,None
492
493    print 'Read Mar345 file: ',filename
494    File = open(filename,'rb')
495    head = File.read(4095)
496    numbers = st.unpack('<iiiiiiiiii',head[:40])
497    lines = head[128:].split('\n')
498    head = []
499    for line in lines:
500        line = line.strip()
501        if 'PIXEL' in line:
502            values = line.split()
503            pixel = (int(values[2]),int(values[4]))     #in microns
504        elif 'WAVELENGTH' in line:
505            wave = float(line.split()[1])
506        elif 'DISTANCE' in line:
507            distance = float(line.split()[1])           #in mm
508        elif 'CENTER' in line:
509            values = line.split()
510            center = [float(values[2])/10.,float(values[4])/10.]    #make in mm from pixels
511        if line: 
512            head.append(line)
513    data = {'pixelSize':pixel,'wavelength':wave,'distance':distance,'center':center}
514    for line in head:
515        if 'FORMAT' in line[0:6]:
516            items = line.split()
517            size = int(items[1])
518    pos = 4096
519    File.seek(pos)
520    line = File.read(8)
521    while 'CCP4' not in line:       #get past overflow list for now
522        line = File.read(8)
523        pos += 8
524    pos += 37
525    File.seek(pos)
526    raw = File.read()
527    File.close()
528    image = np.zeros(shape=(size,size),dtype=np.int32)
529    image = pf.pack_f(len(raw),raw,size,image)
530    return head,data,size,image.T
531   
532def GetTifData(filename):
533    # only works for APS Perkin-Elmer detector data files in "TIFF" format that are readable by Fit2D
534    import struct as st
535    import array as ar
536    print 'Read APS PE-detector tiff file: ',filename
537    File = open(filename,'Ur')
538    Meta = open(filename+'.metadata','Ur')
539    tag = File.read(3)
540    if tag != 'II*':
541        lines = ['not a APS PE-detector tiff file',]
542        return lines,0,0
543    head = Meta.readlines()
544    dataType = 0
545    for line in head:
546        line = line.strip()
547        if 'dataType' in line:
548            dataType = int(line.split('=')[1])
549    size = st.unpack('<i',File.read(4))[0]
550    image = np.zeros(shape=(size,size),dtype=np.int32)
551    row = 0
552    pos = 8
553    while row < size:
554        File.seek(pos)
555        if dataType == 5:
556            line = ar.array('f',File.read(4*size))
557        else:
558            line = ar.array('l',File.read(4*size))
559        image[row] = np.asarray(line)
560        row += 1
561        pos += 4*size
562    data = {'pixelSize':(200,200),'wavelength':0.10,'distance':100.0,'center':[204.8,204.8]} 
563    return head,data,size,image
564    File.close()   
565
566def ProjFileOpen(self):
567    file = open(self.GSASprojectfile,'rb')
568    print 'load from file: ',self.GSASprojectfile
569    wx.BeginBusyCursor()
570    try:
571        while True:
572            try:
573                data = cPickle.load(file)
574            except EOFError:
575                break
576            datum = data[0]
577            print 'load: ',datum[0]
578            if 'PKS' not in datum[0]:
579                if datum[0] not in ['Notebook','Controls','Phases'] and 'PWDR' not in datum[0]:            #temporary fix
580                    datum[0] = 'PWDR '+datum[0]
581            Id = self.PatternTree.AppendItem(parent=self.root,text=datum[0])
582            self.PatternTree.SetItemPyData(Id,datum[1])
583            for datus in data[1:]:
584                print '    load: ',datus[0]
585                sub = self.PatternTree.AppendItem(Id,datus[0])
586                self.PatternTree.SetItemPyData(sub,datus[1])
587            if 'PWDR' in datum[0] and not G2gd.GetPatternTreeItemId(self,Id, 'Comments'):
588                print 'no comments - add comments'
589                sub = self.PatternTree.AppendItem(Id,'Comments')
590                self.PatternTree.SetItemPyData(sub,['no comments'])
591               
592        file.close()
593        if not G2gd.GetPatternTreeItemId(self,self.root,'Notebook'):
594            sub = self.PatternTree.AppendItem(parent=self.root,text='Notebook')
595            self.PatternTree.SetItemPyData(sub,[''])
596            sub = self.PatternTree.AppendItem(parent=self.root,text='Controls')
597            self.PatternTree.SetItemPyData(sub,[0])
598           
599    finally:
600        wx.EndBusyCursor()
601    print 'project load successful'
602    self.NewPlot = True
603   
604def ProjFileSave(self):
605    if not self.PatternTree.IsEmpty():
606        file = open(self.GSASprojectfile,'wb')
607        print 'save to file: ',self.GSASprojectfile
608        wx.BeginBusyCursor()
609        try:
610            item, cookie = self.PatternTree.GetFirstChild(self.root)
611            while item:
612                data = []
613                name = self.PatternTree.GetItemText(item)
614                print 'save: ',name
615                data.append([name,self.PatternTree.GetItemPyData(item)])
616                item2, cookie2 = self.PatternTree.GetFirstChild(item)
617                while item2:
618                    name = self.PatternTree.GetItemText(item2)
619                    print '    save: ',name
620                    data.append([name,self.PatternTree.GetItemPyData(item2)])
621                    item2, cookie2 = self.PatternTree.GetNextChild(item, cookie2)                           
622                item, cookie = self.PatternTree.GetNextChild(self.root, cookie)                           
623                cPickle.dump(data,file,1)
624            file.close()
625        finally:
626            wx.EndBusyCursor()
627        print 'project save successful'
628       
629def PowderxyeSave(self):
630    file = open(self.powderfile,'wa')
631    print 'save powder pattern to file: ',self.powderfile
632    wx.BeginBusyCursor()
633    try:
634        x,y,w,yc,yb,yd = self.PatternTree.GetItemPyData(self.PickId)[1]
635        for i,X in enumerate(x):
636            file.write("%15.6g %15.6g %15.6g\n" % (X,y[i],1.0/math.sqrt(w[i])))
637        file.close()
638    finally:
639        wx.EndBusyCursor()
640    print 'powder pattern file written'
641   
642def PeakListSave(self,file,peaks):
643    print 'save peak list to file: ',self.peaklistfile
644    if not peaks:
645        dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
646        try:
647            result = dlg.ShowModal()
648        finally:
649            dlg.Destroy()
650        return
651    for peak in peaks:
652        file.write("%10.4f %12.2f %10.3f %10.3f \n" % \
653            (peak[0],peak[2],peak[4],peak[6]))
654    print 'peak list saved'
655             
656def IndexPeakListSave(self,peaks):
657    file = open(self.peaklistfile,'wa')
658    print 'save index peak list to file: ',self.peaklistfile
659    wx.BeginBusyCursor()
660    try:
661        if not peaks:
662            dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
663            try:
664                result = dlg.ShowModal()
665            finally:
666                dlg.Destroy()
667            return
668        for peak in peaks:
669            file.write("%12.6f\n" % (peak[7]))
670        file.close()
671    finally:
672        wx.EndBusyCursor()
673    print 'index peak list saved'
674   
675def ReadEXPPhase(filename):
676    import GSASIIspc as G2spc
677    import GSASIIcomp as G2cmp
678    import GSASIIElem as G2el
679    file = open(filename, 'Ur')
680    Phase = {}
681    S = 1
682    Expr = [{},{},{},{},{},{},{},{},{}]
683    while S:
684        S = file.readline()
685        if 'EXPR NPHAS' in S[:12]:
686            Num = S[12:-1].count('0')
687            NPhas = S[12:-1].split()
688        if 'CRS' in S[:3]:
689            N = int(S[3:4])-1
690            Expr[N][S[:12]] = S[12:-1]
691    file.close()
692    PNames = []
693    for n,N in enumerate(NPhas):
694        if N != '0':
695            result = n
696            key = 'CRS'+str(n+1)+'    PNAM'
697            PNames.append(Expr[n][key])
698    if Num < 8:
699        dlg = wx.SingleChoiceDialog(self, 'Which phase to read?', 'Read phase data', PNames, wx.CHOICEDLG_STYLE)
700        try:
701            if dlg.ShowModal() == wx.ID_OK:
702                result = dlg.GetSelection()
703        finally:
704            dlg.Destroy()       
705    EXPphase = Expr[result]
706    keyList = EXPphase.keys()
707    keyList.sort()
708    SGData = {}
709    if NPhas[result] == '1':
710        Ptype = 'nuclear'
711    elif NPhas[result] in ['2','3']:
712        Ptype = 'magnetic'
713    elif NPhas[result] == '4':
714        Ptype = 'macromolecular'
715    elif NPhas[result] == '10':
716        Ptype = 'Pawley'
717    for key in keyList:
718        if 'PNAM' in key:
719           PhaseName = EXPphase[key].strip()
720        elif 'ABC   ' in key:
721            abc = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                       
722        elif 'ANGLES' in key:
723            angles = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                                               
724        elif 'SG SYM' in key:
725            SpGrp = EXPphase[key][:15].strip()
726            E,SGData = G2spc.SpcGroup(SpGrp)
727    Atoms = []
728    if Ptype == 'nuclear':
729        for key in keyList:
730            if 'AT' in key:
731                if key[11:] == 'A':
732                    S = EXPphase[key]
733                elif key[11:] == 'B':
734                    S += EXPphase[key]
735                    Atom = [S[50:58].strip(),S[:10].strip(),'',
736                        float(S[10:20]),float(S[20:30]),float(S[30:40]),
737                        float(S[40:50]),'',int(S[60:62]),S[130:131]]
738                    if Atom[9] == 'I':
739                        Atom += [float(S[68:78]),0.,0.,0.,0.,0.,0.]
740                    elif Atom[9] == 'A':
741                        Atom += [0.0,float(S[68:78]),float(S[78:88]),
742                            float(S[88:98]),float(S[98:108]),
743                            float(S[108:118]),float(S[118:128])]
744                    XYZ = Atom[3:6]
745                    Atom[7],Atom[8] = G2spc.SytSym(XYZ,SGData)
746                    Atoms.append(Atom)
747    elif Ptype == 'macromolecular':
748        for key in keyList:
749            if 'AT' in key[6:8]:
750                S = EXPphase[key]
751                Atom = [int(S[56:60]),S[50:54].strip().upper(),S[54:56],
752                    S[46:51].strip(),S[:8].strip(),'',
753                    float(S[16:24]),float(S[24:32]),float(S[32:40]),
754                    float(S[8:16]),'1',1,'I',float(S[40:46]),0,0,0,0,0,0]
755                XYZ = Atom[6:9]
756                Atom[10],Atom[11] = G2spc.SytSym(XYZ,SGData)
757                Atoms.append(Atom)
758    Volume = G2cmp.calc_V(G2cmp.cell2A(abc+angles))
759    Phase['General'] = [PhaseName,Ptype,SGData,[False,]+abc+angles+[Volume,],[False,1.0],
760        0,0,0,0,0]
761    Phase['Atoms'] = Atoms
762    return Phase
763       
764def ReadPDBPhase(filename):
765    import GSASIIspc as G2spc
766    import GSASIIcomp as G2cmp
767    import GSASIIElem as G2el
768    EightPiSq = 8.*math.pi**2
769    file = open(filename, 'Ur')
770    Phase = {}
771    Title = ''
772    Compnd = ''
773    Atoms = []
774    A = np.zeros(shape=(3,3))
775    S = file.readline()
776    while S:
777        Atom = []
778        if 'TITLE' in S[:5]:
779            Title = S[10:72].strip()
780            S = file.readline()
781        elif 'COMPND    ' in S[:10]:
782            Compnd = S[10:72].strip()
783            S = file.readline()
784        elif 'CRYST' in S[:5]:
785            abc = S[7:34].split()
786            angles = S[34:55].split()
787            cell=[float(abc[0]),float(abc[1]),float(abc[2]),
788                float(angles[0]),float(angles[1]),float(angles[2])]
789            Volume = G2cmp.calc_V(G2cmp.cell2A(cell))
790            AA,AB = G2cmp.cell2AB(cell)
791            SpGrp = S[55:65]
792            E,SGData = G2spc.SpcGroup(SpGrp)
793            S = file.readline()
794        elif 'SCALE' in S[:5]:
795            V = (S[10:41].split())
796            A[int(S[5])-1] = [float(V[0]),float(V[1]),float(V[2])]
797            S = file.readline()
798        elif 'ATOM' in S[:4] or 'HETATM' in S[:6]:
799            XYZ = [float(S[31:39]),float(S[39:47]),float(S[47:55])]
800            XYZ = np.sum(AA*XYZ,axis=1)
801            SytSym,Mult = G2spc.SytSym(XYZ,SGData)
802            Uiso = float(S[61:67])/EightPiSq
803            Type = S[12:14].upper()
804            if Type[0] in '123456789':
805                Type = Type[1:]
806            Atom = [int(S[22:27]),S[17:20].upper(),S[20:22],
807                S[12:17].strip(),Type.strip(),'',XYZ[0],XYZ[1],XYZ[2],
808                float(S[55:61]),SytSym,Mult,'I',Uiso,0,0,0,0,0,0]
809            S = file.readline()
810            if 'ANISOU' in S[:6]:
811                Uij = S[30:72].split()
812                Uij = [float(Uij[0])/10000.,float(Uij[1])/10000.,float(Uij[2])/10000.,
813                    float(Uij[3])/10000.,float(Uij[4])/10000.,float(Uij[5])/10000.]
814                Atom = Atom[:14]+Uij
815                Atom[12] = 'A'
816                S = file.readline()
817            Atoms.append(Atom)
818        else:           
819            S = file.readline()
820    file.close()
821    if Title:
822        PhaseName = Title
823    elif Compnd:
824        PhaseName = Compnd
825    else:
826        PhaseName = 'None'
827    Phase['General'] = [PhaseName,'macromolecular',SGData,[False,]+cell+[Volume,],[False,1.0],
828        0,0,0,0,0]
829    Phase['Atoms'] = Atoms
830   
831    return Phase
832   
833def ReadCIFAtoms(self,data):
834    print data
Note: See TracBrowser for help on using the repository browser.