source: branch/MPbranch/GSASIIIO.py

Last change on this file was 485, checked in by vondreele, 10 years ago
  • Property svn:keywords set to Date Author Revision URL Id
File size: 53.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: 2012-02-16 16:56:42 +0000 (Thu, 16 Feb 2012) $
6# $Author: vondreele $
7# $Revision: 485 $
8# $URL: branch/MPbranch/GSASIIIO.py $
9# $Id: GSASIIIO.py 485 2012-02-16 16:56:42Z 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 GSASIIpwdGUI as G2pdG
22import GSASIIElem as G2el
23import os
24import os.path as ospath
25
26def sfloat(S):
27    if S.strip():
28        return float(S)
29    else:
30        return 0.0
31
32def sint(S):
33    if S.strip():
34        return int(S)
35    else:
36        return 0
37
38def FileDlgFixExt(dlg,file):
39    #this is needed to fix a problem in linux wx.FileDialog
40    ext = dlg.GetWildcard().split('|')[2*dlg.GetFilterIndex()+1].strip('*')
41    if ext not in file:
42        file += ext
43    return file
44   
45def SelectPowderData(G2frame, filename):
46    """Selects banks of data from a filename of any GSAS powder data format
47    Input - filename: any GSAS powder data formatted file (currently STD, FXYE, FXY & ESD)
48    Returns - a list of banks to be read; each entry in list is a tuple containing:
49    filename: same as input filename
50    Pos: position for start of data; record just after BANK record
51    Bank: the BANK record
52    """
53    File = open(filename,'Ur')
54    Title = '''
55First line of this file:
56'''+File.readline()
57    dlg = wx.MessageDialog(G2frame, Title, 'Is this the file you want?', 
58        wx.YES_NO | wx.ICON_QUESTION)
59    try:
60        result = dlg.ShowModal()
61    finally:
62        dlg.Destroy()
63    if result == wx.ID_NO: return (0,0)
64    Temperature = 300
65   
66    if '.xye' in filename:      #Topas style xye file (e.g. 2-th, I, sig) - no iparm file/no BANK record
67        dlg = wx.MessageDialog(G2frame,'''Is this laboratory Cu Ka1/Ka2 data?
68(No = 0.6A wavelength synchrotron data)
69Change wavelength in Instrument Parameters if needed''','Data type?',
70            wx.YES_NO | wx.ICON_QUESTION)
71        try:
72            result = dlg.ShowModal()
73        finally:
74            dlg.Destroy()
75        print result
76        if result == wx.ID_YES:
77            Iparm = {}                                               #Assume CuKa lab data
78            Iparm['INS   HTYPE '] = 'PXC '
79            Iparm['INS  1 ICONS'] = '  1.540500  1.544300       0.0         0       0.7    0       0.5   '
80            Iparm['INS  1PRCF1 '] = '    3    8      0.01                                                '
81            Iparm['INS  1PRCF11'] = '   2.000000E+00  -2.000000E+00   5.000000E+00   0.000000E+00        '
82            Iparm['INS  1PRCF12'] = '   0.000000E+00   0.000000E+00   0.150000E-01   0.150000E-01        '
83        else:
84            Iparm = {}                                               #Assume 0.6A synchrotron data
85            Iparm['INS   HTYPE '] = 'PXC '
86            Iparm['INS  1 ICONS'] = '  0.600000  0.000000       0.0         0      0.99    0       0.5   '
87            Iparm['INS  1PRCF1 '] = '    3    8      0.01                                                '
88            Iparm['INS  1PRCF11'] = '   1.000000E+00  -1.000000E+00   0.300000E+00   0.000000E+00        '
89            Iparm['INS  1PRCF12'] = '   0.000000E+00   0.000000E+00   0.100000E-01   0.100000E-01        '
90                       
91       
92    else:                       #GSAS style fxye or fxy file (e.g. 100*2-th, I, sig)
93        G2frame.IparmName = GetInstrumentFile(G2frame,filename)
94        if G2frame.IparmName:
95            Iparm = GetInstrumentData(G2frame.IparmName)
96        else:
97            Iparm = {}                                               #Assume CuKa lab data if no iparm file
98            Iparm['INS   HTYPE '] = 'PXC '
99            Iparm['INS  1 ICONS'] = '  1.540500  1.544300       0.0         0       0.7    0       0.5   '
100            Iparm['INS  1PRCF1 '] = '    3    8      0.01                                                '
101            Iparm['INS  1PRCF11'] = '   2.000000E+00  -2.000000E+00   5.000000E+00   0.000000E+00        '
102            Iparm['INS  1PRCF12'] = '   0.000000E+00   0.000000E+00   0.150000E-01   0.150000E-01        '
103    S = 1
104    Banks = []
105    Pos = []
106    FoundData = []
107    Comments = []
108    wx.BeginBusyCursor()
109    try:
110        while S:
111            S = File.readline()
112            if S[:1] != '#':
113                if S[:4] == 'BANK':
114                    Banks.append(S)
115                    Pos.append(File.tell())
116                elif '.xye' in filename:    #No BANK in a xye file
117                    Banks.append('BANK 1 XYE')
118                    Pos.append(File.tell())
119                    break
120            else:
121                Comments.append(S[:-1])
122                if 'Temp' in S.split('=')[0]:
123                    Temperature = float(S.split('=')[1])
124        File.close()
125    finally:
126        wx.EndBusyCursor()
127    if Comments:
128       print 'Comments on file:'
129       for Comment in Comments: print Comment
130    if Banks:
131        result = [0]
132        if len(Banks) >= 2:
133            dlg = wx.MultiChoiceDialog(G2frame, 'Which scans do you want?', 'Select scans', Banks, wx.CHOICEDLG_STYLE)
134            try:
135                if dlg.ShowModal() == wx.ID_OK:
136                    result = dlg.GetSelections()
137                else:
138                    result = []
139            finally:
140                dlg.Destroy()
141        for i in result:
142            FoundData.append((filename,Pos[i],Banks[i]))
143    else:
144        dlg = wx.MessageDialog(G2frame, 'ERROR - this is not a GSAS powder data file', 'No BANK records', wx.OK | wx.ICON_ERROR)
145        try:
146            result = dlg.ShowModal()
147        finally:
148            dlg.Destroy()
149    return FoundData,Iparm,Comments,Temperature
150
151def GetInstrumentFile(G2frame,filename):
152    import os.path as op
153    dlg = wx.FileDialog(G2frame,'Choose an instrument file','.', '', 'GSAS iparm file (*.prm)|*.prm|All files(*.*)|*.*', 
154        wx.OPEN|wx.CHANGE_DIR)
155    Tname = filename[:filename.index('.')]+'.prm'
156    if op.exists(Tname):
157        G2frame.IparmName = Tname       
158    if G2frame.IparmName: dlg.SetFilename(G2frame.IparmName)
159    filename = ''
160    try:
161        if dlg.ShowModal() == wx.ID_OK:
162            filename = dlg.GetPath()
163    finally:
164        dlg.Destroy()
165    return filename
166
167def GetInstrumentData(IparmName):
168    file = open(IparmName, 'Ur')
169    S = 1
170    Iparm = {}
171    while S:
172        S = file.readline()
173        Iparm[S[:12]] = S[12:-1]
174    return Iparm
175   
176def GetPowderPeaks(fileName):
177    sind = lambda x: math.sin(x*math.pi/180.)
178    asind = lambda x: 180.*math.asin(x)/math.pi
179    Cuka = 1.54052
180    File = open(fileName,'Ur')
181    Comments = []
182    peaks = []
183    S = File.readline()
184    while S:
185        if S[:1] == '#':
186            Comments.append(S[:-1])
187        else:
188            item = S.split()
189            if len(item) == 1:
190                peaks.append([float(item[0]),1.0])
191            elif len(item) > 1:
192                peaks.append([float(item[0]),float(item[0])])
193        S = File.readline()
194    File.close()
195    if Comments:
196       print 'Comments on file:'
197       for Comment in Comments: print Comment
198    Peaks = []
199    if peaks[0][0] > peaks[-1][0]:          # d-spacings - assume CuKa
200        for peak in peaks:
201            dsp = peak[0]
202            sth = Cuka/(2.0*dsp)
203            if sth < 1.0:
204                tth = 2.0*asind(sth)
205            else:
206                break
207            Peaks.append([tth,peak[1],True,False,0,0,0,dsp,0.0])
208    else:                                   #2-thetas - assume Cuka (for now)
209        for peak in peaks:
210            tth = peak[0]
211            dsp = Cuka/(2.0*sind(tth/2.0))
212            Peaks.append([tth,peak[1],True,False,0,0,0,dsp,0.0])
213    return Comments,Peaks
214
215def GetPawleyPeaks(filename):
216    rt2ln2x2 = 2.35482
217    File = open(filename,'Ur')
218    PawleyPeaks = []
219    S = File.readline()         #skip header
220    S = File.readline()
221    item = S.split()
222    while S:
223        h,k,l = int(item[0]),int(item[1]),int(item[2])
224        mult = int(item[3])
225        tth = float(item[5])
226        sig = float(item[6])/rt2ln2x2
227        Iobs = float(item[7])*mult
228        PawleyPeaks.append([h,k,l,mult,tth,False,Iobs,0.0])
229        S = File.readline()
230        item = S.split()
231        if item[3] == '-100.0000':       #find trailer
232            break
233    File.close()
234    return PawleyPeaks
235   
236def GetHKLData(filename):
237    print 'Reading: '+filename
238    File = open(filename,'Ur')
239    HKLref = []
240    HKLmin = [1000,1000,1000]
241    HKLmax = [0,0,0]
242    FoMax = 0
243    ifFc = False
244    S = File.readline()
245    while '#' in S[0]:        #get past comments if any
246        S = File.readline()       
247    if '_' in S:         #cif style .hkl file
248        while 'loop_' not in S:         #skip preliminaries if any - can't have 'loop_' in them!
249            S = File.readline()       
250        S = File.readline()             #get past 'loop_' line
251        pos = 0
252        hpos = kpos = lpos = Fosqpos = Fcsqpos = sigpos = -1
253        while S:
254            if '_' in S:
255                if 'index_h' in S:
256                    hpos = pos
257                elif 'index_k' in S:
258                    kpos = pos
259                elif 'index_l' in S:
260                    lpos = pos
261                elif 'F_squared_meas' in S:
262                    Fosqpos = pos
263                elif 'F_squared_calc' in S:
264                    Fcsqpos = pos
265                elif 'F_squared_sigma' in S:
266                    sigpos = pos
267                pos += 1
268            else:
269                data = S.split()
270                if data:                    #avoid blank lines
271                    HKL = np.array([int(data[hpos]),int(data[kpos]),int(data[lpos])])
272                    h,k,l = HKL
273                    Fosq = float(data[Fosqpos])
274                    if sigpos != -1:
275                        sigFosq = float(data[sigpos])
276                    else:
277                        sigFosq = 1.
278                    if Fcsqpos != -1:
279                        Fcsq = float(data[Fcsqpos])
280                        if Fcsq:
281                            ifFc = True
282                    else:
283                        Fcsq = 0.
284                       
285                    HKLmin = [min(h,HKLmin[0]),min(k,HKLmin[1]),min(l,HKLmin[2])]
286                    HKLmax = [max(h,HKLmax[0]),max(k,HKLmax[1]),max(l,HKLmax[2])]
287                    FoMax = max(FoMax,Fosq)
288                    HKLref.append([HKL,Fosq,sigFosq,Fcsq,0,0,0])                 #room for Fcp, Fcpp & phase
289            S = File.readline()
290    else:                   #dumb h,k,l,Fo,sigFo .hkl file
291        while S:
292            h,k,l,Fo,sigFo = S.split()
293            HKL = np.array([int(h),int(k),int(l)])
294            h,k,l = HKL
295            Fo = float(Fo)
296            sigFo = float(sigFo)
297            HKLmin = [min(h,HKLmin[0]),min(k,HKLmin[1]),min(l,HKLmin[2])]
298            HKLmax = [max(h,HKLmax[0]),max(k,HKLmax[1]),max(l,HKLmax[2])]
299            FoMax = max(FoMax,Fo)
300            HKLref.append([HKL,Fo**2,2.*Fo*sigFo,0,0,0,0])                 #room for Fc, Fcp, Fcpp & phase
301            S = File.readline()
302    File.close()
303    return HKLref,HKLmin,HKLmax,FoMax,ifFc
304
305def GetPowderData(filename,Pos,Bank,DataType):
306    '''Reads one BANK of data from GSAS raw powder data file
307    input:
308    filename: GSAS raw powder file dataname
309    Pos: start of data in file just after BANK record
310    Bank: the BANK record
311    DataType: powder data type, e.g. "PXC" for Powder X-ray CW data
312    returns: list [x,y,e,yc,yb]
313    x: np.array of x-axis values
314    y: np.array of powder pattern intensities
315    w: np.array of w=sig(intensity)^2 values
316    yc: np.array of calc. intensities (zero)
317    yb: np.array of calc. background (zero)
318    yd: np.array of obs-calc profiles
319    '''
320    print 'Reading: '+filename
321    print 'Bank:    '+Bank[:-1]
322    if 'FXYE' in Bank:
323        return GetFXYEdata(filename,Pos,Bank,DataType)
324    elif ' XYE' in Bank:
325        return GetXYEdata(filename,Pos,Bank,DataType)
326    elif 'FXY' in Bank:
327        return GetFXYdata(filename,Pos,Bank,DataType)
328    elif 'ESD' in Bank:
329        return GetESDdata(filename,Pos,Bank,DataType)
330    elif 'STD' in Bank:
331        return GetSTDdata(filename,Pos,Bank,DataType)
332    else:
333        return GetSTDdata(filename,Pos,Bank,DataType)
334    return []
335
336def GetFXYEdata(filename,Pos,Bank,DataType):
337    File = open(filename,'Ur')
338    File.seek(Pos)
339    x = []
340    y = []
341    w = []
342    S = File.readline()
343    while S and S[:4] != 'BANK':
344        vals = S.split()
345        if DataType[2] == 'C':
346            x.append(float(vals[0])/100.)               #CW: from centidegrees to degrees
347        elif DataType[2] == 'T':
348            x.append(float(vals[0])/1000.0)             #TOF: from musec to millisec
349        f = float(vals[1])
350        if f <= 0.0:
351            y.append(0.0)
352            w.append(1.0)
353        else:
354            y.append(float(vals[1]))
355            w.append(1.0/float(vals[2])**2)
356        S = File.readline()
357    File.close()
358    N = len(x)
359    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
360   
361def GetXYEdata(filename,Pos,Bank,DataType):
362    File = open(filename,'Ur')
363    File.seek(Pos)
364    x = []
365    y = []
366    w = []
367    S = File.readline()
368    while S:
369        vals = S.split()
370        try:
371            x.append(float(vals[0]))
372            f = float(vals[1])
373            if f <= 0.0:
374                y.append(0.0)
375                w.append(1.0)
376            else:
377                y.append(float(vals[1]))
378                w.append(1.0/float(vals[2])**2)
379            S = File.readline()
380        except ValueError:
381            break
382    File.close()
383    N = len(x)
384    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
385   
386   
387def GetFXYdata(filename,Pos,Bank,DataType):
388    File = open(filename,'Ur')
389    File.seek(Pos)
390    x = []
391    y = []
392    w = []
393    S = File.readline()
394    while S and S[:4] != 'BANK':
395        vals = S.split()
396        if DataType[2] == 'C':
397            x.append(float(vals[0])/100.)               #CW: from centidegrees to degrees
398        elif DataType[2] == 'T':
399            x.append(float(vals[0])/1000.0)             #TOF: from musec to millisec
400        f = float(vals[1])
401        if f > 0.0:
402            y.append(f)
403            w.append(1.0/f)
404        else:             
405            y.append(0.0)
406            w.append(1.0)
407        S = File.readline()
408    File.close()
409    N = len(x)
410    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
411   
412def GetESDdata(filename,Pos,Bank,DataType):
413    File = open(filename,'Ur')
414    cons = Bank.split()
415    if DataType[2] == 'C':
416        start = float(cons[5])/100.0               #CW: from centidegrees to degrees
417        step = float(cons[6])/100.0
418    elif DataType[2] == 'T':
419        start = float(cons[5])/1000.0              #TOF: from musec to millisec
420        step = float(cons[6])/1000.0
421    File.seek(Pos)
422    x = []
423    y = []
424    w = []
425    S = File.readline()
426    j = 0
427    while S and S[:4] != 'BANK':
428        for i in range(0,80,16):
429            xi = start+step*j
430            yi = sfloat(S[i:i+8])
431            ei = sfloat(S[i+8:i+16])
432            x.append(xi)
433            if yi > 0.0:
434                y.append(yi)
435                w.append(1.0/ei**2)
436            else:             
437                y.append(0.0)
438                w.append(1.0)
439            j += 1
440        S = File.readline()
441    File.close()
442    N = len(x)
443    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
444
445def GetSTDdata(filename,Pos,Bank,DataType):
446    File = open(filename,'Ur')
447    cons = Bank.split()
448    Nch = int(cons[2])
449    if DataType[2] == 'C':
450        start = float(cons[5])/100.0               #CW: from centidegrees to degrees
451        step = float(cons[6])/100.0
452    elif DataType[2] == 'T':
453        start = float(cons[5])/1000.0              #TOF: from musec to millisec - not likely!
454        step = float(cons[6])/1000.0
455    File.seek(Pos)
456    x = []
457    y = []
458    w = []
459    S = File.readline()
460    j = 0
461    while S and S[:4] != 'BANK':
462        for i in range(0,80,8):
463            xi = start+step*j
464            ni = max(sint(S[i:i+2]),1)
465            yi = max(sfloat(S[i+2:i+8]),0.0)
466            if yi:
467                vi = yi/ni
468            else:
469                yi = 0.0
470                vi = 1.0
471            j += 1
472            if j < Nch:
473                x.append(xi)
474                y.append(yi)
475                w.append(1.0/vi)
476        S = File.readline()
477    File.close()
478    N = len(x)
479    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
480   
481def CheckImageFile(G2frame,imagefile):
482    if not ospath.exists(imagefile):
483        dlg = wx.FileDialog(G2frame, 'Bad image file name; choose name', '.', '',\
484        'Any image file (*.tif;*.tiff;*.mar*;*.avg;*.sum;*.img)\
485        |*.tif;*.tiff;*.mar*;*.avg;*.sum;*.img|\
486        Any detector tif (*.tif;*.tiff)|*.tif;*.tiff|\
487        MAR file (*.mar*)|*.mar*|\
488        GE Image (*.avg;*.sum)|*.avg;*.sum|\
489        ADSC Image (*.img)|*.img|\
490        All files (*.*)|*.*',wx.OPEN|wx.CHANGE_DIR)
491        try:
492            dlg.SetFilename(ospath.split(imagefile)[1])
493            if dlg.ShowModal() == wx.ID_OK:
494                imagefile = dlg.GetPath()
495            else:
496                imagefile = False
497        finally:
498            dlg.Destroy()
499    return imagefile
500       
501def GetImageData(G2frame,imagefile,imageOnly=False):       
502    ext = ospath.splitext(imagefile)[1]
503    Comments = []
504    if ext == '.tif' or ext == '.tiff':
505        Comments,Data,Npix,Image = GetTifData(imagefile)
506    elif ext == '.img':
507        Comments,Data,Npix,Image = GetImgData(imagefile)
508        Image[0][0] = 0
509    elif ext == '.mar3450' or ext == '.mar2300':
510        Comments,Data,Npix,Image = GetMAR345Data(imagefile)
511    elif ext in ['.sum','.avg','']:
512        Comments,Data,Npix,Image = GetGEsumData(imagefile)
513    elif ext == '.G2img':
514        Comments,Data,Npix,Image = GetG2Image(imagefile)
515    if imageOnly:
516        return Image
517    else:
518        return Comments,Data,Npix,Image
519       
520def PutG2Image(filename,Comments,Data,Npix,image):
521    File = open(filename,'wb')
522    cPickle.dump([Comments,Data,Npix,image],File,1)
523    File.close()
524    return
525   
526def GetG2Image(filename):
527    File = open(filename,'rb')
528    Comments,Data,Npix,image = cPickle.load(File)
529    File.close()
530    return Comments,Data,Npix,image
531   
532def GetGEsumData(filename,imageOnly=False):
533    import struct as st
534    import array as ar
535    if not imageOnly:
536        print 'Read GE sum file: ',filename   
537    File = open(filename,'rb')
538    if '.sum' in filename:
539        head = ['GE detector sum data from APS 1-ID',]
540        sizexy = [2048,2048]
541    elif '.avg' in filename:
542        head = ['GE detector avg data from APS 1-ID',]
543        sizexy = [2048,2048]
544    else:
545        head = ['GE detector raw data from APS 1-ID',]
546        File.seek(18)
547        size,nframes = st.unpack('<ih',File.read(6))
548        sizexy = [2048,2048]
549        pos = 8192
550        File.seek(pos)
551    Npix = sizexy[0]*sizexy[1]
552    if '.sum' in filename:
553        image = np.array(ar.array('f',File.read(4*Npix)),dtype=np.int32)
554    elif '.avg' in filename:
555        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
556    else:
557        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
558        while nframes > 1:
559            image += np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
560            nframes -= 1
561    image = np.reshape(image,(sizexy[1],sizexy[0]))
562    data = {'pixelSize':(200,200),'wavelength':0.15,'distance':250.0,'center':[204.8,204.8],'size':sizexy} 
563    File.close()   
564    if imageOnly:
565        return image
566    else:
567        return head,data,Npix,image
568       
569def GetImgData(filename,imageOnly=False):
570    import struct as st
571    import array as ar
572    if not imageOnly:
573        print 'Read ADSC img file: ',filename
574    File = open(filename,'rb')
575    head = File.read(511)
576    lines = head.split('\n')
577    head = []
578    center = [0,0]
579    for line in lines[1:-2]:
580        line = line.strip()[:-1]
581        if line:
582            if 'SIZE1' in line:
583                size = int(line.split('=')[1])
584                Npix = size*size
585            elif 'WAVELENGTH' in line:
586                wave = float(line.split('=')[1])
587            elif 'BIN' in line:
588                if line.split('=')[1] == '2x2':
589                    pixel=(102,102)
590                else:
591                    pixel = (51,51)
592            elif 'DISTANCE' in line:
593                distance = float(line.split('=')[1])
594            elif 'CENTER_X' in line:
595                center[0] = float(line.split('=')[1])
596            elif 'CENTER_Y' in line:
597                center[1] = float(line.split('=')[1])
598            head.append(line)
599    data = {'pixelSize':pixel,'wavelength':wave,'distance':distance,'center':center,'size':[size,size]}
600    image = []
601    row = 0
602    pos = 512
603    File.seek(pos)
604    image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
605    image = np.reshape(image,(sizexy[1],sizexy[0]))
606#    image = np.zeros(shape=(size,size),dtype=np.int32)   
607#    while row < size:
608#        File.seek(pos)
609#        line = ar.array('H',File.read(2*size))
610#        image[row] = np.asarray(line)
611#        row += 1
612#        pos += 2*size
613    File.close()
614    if imageOnly:
615        return image
616    else:
617        return lines[1:-2],data,Npix,image
618       
619def GetMAR345Data(filename,imageOnly=False):
620    import array as ar
621    import struct as st
622    try:
623        import pack_f as pf
624    except:
625        msg = wx.MessageDialog(None, message="Unable to load the GSAS MAR image decompression, pack_f",
626                               caption="Import Error",
627                               style=wx.ICON_ERROR | wx.OK | wx.STAY_ON_TOP)
628        msg.ShowModal()
629        return None,None,None,None
630
631    if not imageOnly:
632        print 'Read Mar345 file: ',filename
633    File = open(filename,'rb')
634    head = File.read(4095)
635    numbers = st.unpack('<iiiiiiiiii',head[:40])
636    lines = head[128:].split('\n')
637    head = []
638    for line in lines:
639        line = line.strip()
640        if 'PIXEL' in line:
641            values = line.split()
642            pixel = (int(values[2]),int(values[4]))     #in microns
643        elif 'WAVELENGTH' in line:
644            wave = float(line.split()[1])
645        elif 'DISTANCE' in line:
646            distance = float(line.split()[1])           #in mm
647        elif 'CENTER' in line:
648            values = line.split()
649            center = [float(values[2])/10.,float(values[4])/10.]    #make in mm from pixels
650        if line: 
651            head.append(line)
652    data = {'pixelSize':pixel,'wavelength':wave,'distance':distance,'center':center}
653    for line in head:
654        if 'FORMAT' in line[0:6]:
655            items = line.split()
656            size = int(items[1])
657            Npix = size*size
658    pos = 4096
659    data['size'] = [size,size]
660    File.seek(pos)
661    line = File.read(8)
662    while 'CCP4' not in line:       #get past overflow list for now
663        line = File.read(8)
664        pos += 8
665    pos += 37
666    File.seek(pos)
667    raw = File.read()
668    File.close()
669    image = np.zeros(shape=(size,size),dtype=np.int32)
670    image = pf.pack_f(len(raw),raw,size,image)
671    if imageOnly:
672        return image.T              #transpose to get it right way around
673    else:
674        return head,data,Npix,image.T
675       
676def GetTifData(filename,imageOnly=False):
677    import struct as st
678    import array as ar
679    File = open(filename,'rb')
680    dataType = 5
681    try:
682        Meta = open(filename+'.metadata','Ur')
683        head = Meta.readlines()
684        for line in head:
685            line = line.strip()
686            if 'dataType=' in line:
687                dataType = int(line.split('=')[1])
688        Meta.close()
689    except IOError:
690        print 'no metadata file found - will try to read file anyway'
691        head = ['no metadata file found',]
692       
693    tag = File.read(2)
694    byteOrd = '<'
695    if tag == 'II' and int(st.unpack('<h',File.read(2))[0]) == 42:     #little endian
696        IFD = int(st.unpack(byteOrd+'i',File.read(4))[0])
697    elif tag == 'MM' and int(st.unpack('>h',File.read(2))[0]) == 42:   #big endian
698        byteOrd = '>'
699        IFD = int(st.unpack(byteOrd+'i',File.read(4))[0])       
700    else:
701        lines = ['not a detector tiff file',]
702        return lines,0,0,0
703    File.seek(IFD)                                                  #get number of directory entries
704    NED = int(st.unpack(byteOrd+'h',File.read(2))[0])
705    IFD = {}
706    for ied in range(NED):
707        Tag,Type = st.unpack(byteOrd+'Hh',File.read(4))
708        nVal = st.unpack(byteOrd+'i',File.read(4))[0]
709        if Type == 1:
710            Value = st.unpack(byteOrd+nVal*'b',File.read(nVal))
711        elif Type == 2:
712            Value = st.unpack(byteOrd+'i',File.read(4))
713        elif Type == 3:
714            Value = st.unpack(byteOrd+nVal*'h',File.read(nVal*2))
715            x = st.unpack(byteOrd+nVal*'h',File.read(nVal*2))
716        elif Type == 4:
717            Value = st.unpack(byteOrd+nVal*'i',File.read(nVal*4))
718        elif Type == 5:
719            Value = st.unpack(byteOrd+nVal*'i',File.read(nVal*4))
720        elif Type == 11:
721            Value = st.unpack(byteOrd+nVal*'f',File.read(nVal*4))
722        IFD[Tag] = [Type,nVal,Value]
723#    for key in IFD:
724#        print key,IFD[key]
725    sizexy = [IFD[256][2][0],IFD[257][2][0]]
726    [nx,ny] = sizexy
727    Npix = nx*ny
728    if 272 in IFD:
729        ifd = IFD[272]
730        File.seek(ifd[2][0])
731        S = File.read(ifd[1])
732        if 'PILATUS' in S:
733            tifType = 'Pilatus'
734            dataType = 0
735            pixy = (172,172)
736            File.seek(4096)
737            if not imageOnly:
738                print 'Read Pilatus tiff file: ',filename
739            image = ar.array('L',File.read(4*Npix))
740            image = np.array(np.asarray(image),dtype=np.int32)
741    elif 262 in IFD and IFD[262][2][0] > 4:
742        tifType = 'DND'
743        pixy = (158,158)
744        File.seek(512)
745        if not imageOnly:
746            print 'Read DND SAX/WAX-detector tiff file: ',filename
747        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
748    elif sizexy == [1536,1536]:
749        tifType = 'APS Gold'
750        pixy = (150,150)
751        File.seek(64)
752        if not imageOnly:
753            print 'Read Gold tiff file:',filename
754        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
755    elif sizexy == [2048,2048] or sizexy == [1024,1024]:
756        if IFD[273][2][0] == 8:
757            if IFD[258][2][0] == 32:
758                tifType = 'PE'
759                pixy = (200,200)
760                File.seek(8)
761                if not imageOnly:
762                    print 'Read APS PE-detector tiff file: ',filename
763                if dataType == 5:
764                    image = np.array(ar.array('f',File.read(4*Npix)),dtype=np.float32)
765                else:
766                    image = np.array(ar.array('I',File.read(4*Npix)),dtype=np.int32)
767        elif IFD[273][2][0] == 4096:
768            tifType = 'MAR'
769            pixy = (158,158)
770            File.seek(4096)
771            if not imageOnly:
772                print 'Read MAR CCD tiff file: ',filename
773            image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
774    elif sizexy == [4096,4096]:
775        if IFD[273][2][0] == 8:
776            if IFD[258][2][0] == 16:
777                tifType = 'scanCCD'
778                pixy = (9,9)
779                File.seek(8)
780                if not imageOnly:
781                    print 'Read APS scanCCD tiff file: ',filename
782                image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
783#    elif sizexy == [960,960]:
784#        tiftype = 'PE-BE'
785#        pixy = (200,200)
786#        File.seek(8)
787#        if not imageOnly:
788#            print 'Read Gold tiff file:',filename
789#        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
790           
791    else:
792        lines = ['not a known detector tiff file',]
793        return lines,0,0,0
794       
795    image = np.reshape(image,(sizexy[1],sizexy[0]))
796    center = [pixy[0]*sizexy[0]/2000,pixy[1]*sizexy[1]/2000]
797    data = {'pixelSize':pixy,'wavelength':0.10,'distance':100.0,'center':center,'size':sizexy}
798    File.close()   
799    if imageOnly:
800        return image
801    else:
802        return head,data,Npix,image
803   
804def ProjFileOpen(G2frame):
805    file = open(G2frame.GSASprojectfile,'rb')
806    print 'load from file: ',G2frame.GSASprojectfile
807    wx.BeginBusyCursor()
808    try:
809        while True:
810            try:
811                data = cPickle.load(file)
812            except EOFError:
813                break
814            datum = data[0]
815           
816            Id = G2frame.PatternTree.AppendItem(parent=G2frame.root,text=datum[0])
817            if 'PWDR' in datum[0]:               
818                G2frame.PatternTree.SetItemPyData(Id,datum[1][:3])     #temp. trim off junk
819            else:
820                G2frame.PatternTree.SetItemPyData(Id,datum[1])
821            for datus in data[1:]:
822                sub = G2frame.PatternTree.AppendItem(Id,datus[0])
823                G2frame.PatternTree.SetItemPyData(sub,datus[1])
824            if 'IMG' in datum[0]:                   #retreive image default flag & data if set
825                Data = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Image Controls'))
826                if Data['setDefault']:
827                    G2frame.imageDefault = Data               
828        file.close()
829       
830    finally:
831        wx.EndBusyCursor()
832    print 'project load successful'
833    G2frame.NewPlot = True
834   
835def ProjFileSave(G2frame):
836    if not G2frame.PatternTree.IsEmpty():
837        file = open(G2frame.GSASprojectfile,'wb')
838        print 'save to file: ',G2frame.GSASprojectfile
839        wx.BeginBusyCursor()
840        try:
841            item, cookie = G2frame.PatternTree.GetFirstChild(G2frame.root)
842            while item:
843                data = []
844                name = G2frame.PatternTree.GetItemText(item)
845                data.append([name,G2frame.PatternTree.GetItemPyData(item)])
846                item2, cookie2 = G2frame.PatternTree.GetFirstChild(item)
847                while item2:
848                    name = G2frame.PatternTree.GetItemText(item2)
849                    data.append([name,G2frame.PatternTree.GetItemPyData(item2)])
850                    item2, cookie2 = G2frame.PatternTree.GetNextChild(item, cookie2)                           
851                item, cookie = G2frame.PatternTree.GetNextChild(G2frame.root, cookie)                           
852                cPickle.dump(data,file,1)
853            file.close()
854        finally:
855            wx.EndBusyCursor()
856        print 'project save successful'
857
858def SaveIntegration(G2frame,PickId,data):
859    azms = G2frame.Integrate[1]
860    X = G2frame.Integrate[2][:-1]
861    Xminmax = [X[0],X[-1]]
862    N = len(X)
863    Id = G2frame.PatternTree.GetItemParent(PickId)
864    name = G2frame.PatternTree.GetItemText(Id)
865    name = name.replace('IMG ','PWDR ')
866    Comments = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id, 'Comments'))
867    names = ['Type','Lam','Zero','Polariz.','U','V','W','X','Y','SH/L','Azimuth'] 
868    codes = [0 for i in range(11)]
869    LRazm = data['LRazimuth']
870    Azms = []
871    if data['fullIntegrate'] and data['outAzimuths'] == 1:
872        Azms = [45.0,]                              #a poor man's average?
873    else:
874        for i,azm in enumerate(azms[:-1]):
875            Azms.append((azms[i+1]+azm)/2.)
876    for i,azm in enumerate(azms[:-1]):
877        item, cookie = G2frame.PatternTree.GetFirstChild(G2frame.root)
878        Id = 0
879        while item:
880            Name = G2frame.PatternTree.GetItemText(item)
881            if name == Name:
882                Id = item
883            item, cookie = G2frame.PatternTree.GetNextChild(G2frame.root, cookie)
884        parms = ['PXC',data['wavelength'],0.0,0.99,1.0,-0.10,0.4,0.30,1.0,0.0001,Azms[i]]    #set polarization for synchrotron radiation!
885        Y = G2frame.Integrate[0][i]
886        W = 1./Y                    #probably not true
887        Sample = G2pdG.SetDefaultSample()
888        if Id:
889            G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id, 'Comments'),Comments)                   
890            G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Limits'),[tuple(Xminmax),Xminmax])
891            G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Background'),[['chebyschev',1,3,1.0,0.0,0.0],
892                            {'nDebye':0,'debyeTerms':[],'nPeaks':0,'peaksList':[]}])
893            G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Instrument Parameters'),[tuple(parms),parms[:],codes,names])
894            G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Peak List'),[])
895            G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Index Peak List'),[])
896            G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Unit Cells List'),[])             
897            G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Reflection Lists'),{})             
898        else:
899            Id = G2frame.PatternTree.AppendItem(parent=G2frame.root,text=name+" Azm= %.2f"%(Azms[i]))
900            G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Comments'),Comments)                   
901            G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Limits'),[tuple(Xminmax),Xminmax])
902            G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Background'),[['chebyschev',1,3,1.0,0.0,0.0],
903                            {'nDebye':0,'debyeTerms':[],'nPeaks':0,'peaksList':[]}])
904            G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Instrument Parameters'),[tuple(parms),parms[:],codes,names])
905            G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Sample Parameters'),Sample)
906            G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Peak List'),[])
907            G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Index Peak List'),[])
908            G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Unit Cells List'),[])
909            G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Reflection Lists'),{})             
910        G2frame.PatternTree.SetItemPyData(Id,[[''],[np.array(X),np.array(Y),np.array(W),np.zeros(N),np.zeros(N),np.zeros(N)]])
911    G2frame.PatternTree.SelectItem(Id)
912    G2frame.PatternTree.Expand(Id)
913    G2frame.PatternId = Id
914           
915def powderFxyeSave(G2frame,exports,powderfile):
916    head,tail = ospath.split(powderfile)
917    name,ext = tail.split('.')
918    wx.BeginBusyCursor()
919    for i,export in enumerate(exports):
920        filename = ospath.join(head,name+'-%03d.'%(i)+ext)
921        prmname = filename.strip(ext)+'prm'
922        prm = open(prmname,'w')      #old style GSAS parm file
923        PickId = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, export)
924        Values,Names = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame, \
925            PickId, 'Instrument Parameters'))[1::2]     #get values & names
926        Inst = dict(zip(Names,Values))
927        print Inst['Type']
928        prm.write( '            123456789012345678901234567890123456789012345678901234567890        '+'\n')
929        prm.write( 'INS   BANK      1                                                               '+'\n')
930        prm.write(('INS   HTYPE   %sR                                                              '+'\n')%(Inst['Type']))
931        if 'Lam1' in Inst:              #Ka1 & Ka2
932            prm.write(('INS  1 ICONS%10.7f%10.7f    0.0000               0.990    0     0.500   '+'\n')%(Inst['Lam1'],Inst['Lam2']))
933        elif 'Lam' in Inst:             #single wavelength
934            prm.write(('INS  1 ICONS%10.7f%10.7f    0.0000               0.990    0     0.500   '+'\n')%(Inst['Lam'],0.0))
935        prm.write( 'INS  1 IRAD     0                                                               '+'\n')
936        prm.write( 'INS  1I HEAD                                                                    '+'\n')
937        prm.write( 'INS  1I ITYP    0    0.0000  180.0000         1                                 '+'\n')
938        prm.write(('INS  1DETAZM%10.3f                                                          '+'\n')%(Inst['Azimuth']))
939        prm.write( 'INS  1PRCF1     3    8   0.00100                                                '+'\n')
940        prm.write(('INS  1PRCF11     %15.6g%15.6g%15.6g%15.6g   '+'\n')%(Inst['U'],Inst['V'],Inst['W'],0.0))
941        prm.write(('INS  1PRCF12     %15.6g%15.6g%15.6g%15.6g   '+'\n')%(Inst['X'],Inst['Y'],Inst['SH/L']/2.,Inst['SH/L']/2.))
942        prm.close()
943        file = open(filename,'w')
944        print 'save powder pattern to file: ',filename
945        try:
946            x,y,w,yc,yb,yd = G2frame.PatternTree.GetItemPyData(PickId)[1]
947            file.write(powderfile+'\n')
948            file.write('BANK 1 %d %d CONS %.2f %.2f 0 0 FXYE\n'%(len(x),len(x),\
949                100.*x[0],100.*(x[1]-x[0])))
950            s = list(np.sqrt(1./np.array(w)))       
951            XYW = zip(x,y,s)
952            for X,Y,S in XYW:
953                file.write("%15.6g %15.6g %15.6g\n" % (100.*X,Y,max(S,1.0)))
954            file.close()
955        finally:
956            wx.EndBusyCursor()
957        print 'powder pattern file written'
958       
959def powderXyeSave(G2frame,exports,powderfile):
960    head,tail = ospath.split(powderfile)
961    name,ext = tail.split('.')
962    for i,export in enumerate(exports):
963        filename = ospath.join(head,name+'-%03d.'%(i)+ext)
964        PickId = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, export)
965        file = open(filename,'w')
966        file.write('#%s\n'%(export))
967        print 'save powder pattern to file: ',filename
968        wx.BeginBusyCursor()
969        try:
970            x,y,w,yc,yb,yd = G2frame.PatternTree.GetItemPyData(PickId)[1]
971            s = list(np.sqrt(1./np.array(w)))       
972            XYW = zip(x,y,s)
973            for X,Y,W in XYW:
974                file.write("%15.6g %15.6g %15.6g\n" % (X,Y,W))
975            file.close()
976        finally:
977            wx.EndBusyCursor()
978        print 'powder pattern file written'
979       
980def PDFSave(G2frame,exports):   
981    for export in exports:
982        PickId = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, export)
983        SQname = 'S(Q)'+export[4:]
984        GRname = 'G(R)'+export[4:]
985        sqfilename = ospath.join(G2frame.dirname,export.replace(' ','_')[5:]+'.sq')
986        grfilename = ospath.join(G2frame.dirname,export.replace(' ','_')[5:]+'.gr')
987        sqId = G2gd.GetPatternTreeItemId(G2frame, PickId, SQname)
988        grId = G2gd.GetPatternTreeItemId(G2frame, PickId, GRname)
989        sqdata = np.array(G2frame.PatternTree.GetItemPyData(sqId)[1][:2]).T
990        grdata = np.array(G2frame.PatternTree.GetItemPyData(grId)[1][:2]).T
991        sqfile = open(sqfilename,'w')
992        grfile = open(grfilename,'w')
993        sqfile.write('#T S(Q) %s\n'%(export))
994        grfile.write('#T G(R) %s\n'%(export))
995        sqfile.write('#L Q     S(Q)\n')
996        grfile.write('#L R     G(R)\n')
997        for q,sq in sqdata:
998            sqfile.write("%15.6g %15.6g\n" % (q,sq))
999        sqfile.close()
1000        for r,gr in grdata:
1001            grfile.write("%15.6g %15.6g\n" % (r,gr))
1002        grfile.close()
1003   
1004def PeakListSave(G2frame,file,peaks):
1005    print 'save peak list to file: ',G2frame.peaklistfile
1006    if not peaks:
1007        dlg = wx.MessageDialog(G2frame, 'No peaks!', 'Nothing to save!', wx.OK)
1008        try:
1009            result = dlg.ShowModal()
1010        finally:
1011            dlg.Destroy()
1012        return
1013    for peak in peaks:
1014        file.write("%10.4f %12.2f %10.3f %10.3f \n" % \
1015            (peak[0],peak[2],peak[4],peak[6]))
1016    print 'peak list saved'
1017             
1018def IndexPeakListSave(G2frame,peaks):
1019    file = open(G2frame.peaklistfile,'wa')
1020    print 'save index peak list to file: ',G2frame.peaklistfile
1021    wx.BeginBusyCursor()
1022    try:
1023        if not peaks:
1024            dlg = wx.MessageDialog(G2frame, 'No peaks!', 'Nothing to save!', wx.OK)
1025            try:
1026                result = dlg.ShowModal()
1027            finally:
1028                dlg.Destroy()
1029            return
1030        for peak in peaks:
1031            file.write("%12.6f\n" % (peak[7]))
1032        file.close()
1033    finally:
1034        wx.EndBusyCursor()
1035    print 'index peak list saved'
1036   
1037def ReadEXPPhase(G2frame,filename):
1038    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1039    textureData = {'Order':0,'Model':'cylindrical','Sample omega':[False,0.0],
1040        'Sample chi':[False,0.0],'Sample phi':[False,0.0],'SH Coeff':[False,{}],
1041        'SHShow':False,'PFhkl':[0,0,1],'PFxyz':[0,0,1],'PlotType':'Pole figure'}
1042    shNcof = 0
1043    file = open(filename, 'Ur')
1044    S = 1
1045    Expr = [{},{},{},{},{},{},{},{},{}]
1046    while S:
1047        S = file.readline()
1048        if 'EXPR NPHAS' in S[:12]:
1049            Num = S[12:-1].count('0')
1050            NPhas = S[12:-1].split()
1051        if 'CRS' in S[:3]:
1052            N = int(S[3:4])-1
1053            Expr[N][S[:12]] = S[12:-1]
1054    file.close()
1055    PNames = []
1056    for n,N in enumerate(NPhas):
1057        if N != '0':
1058            result = n
1059            key = 'CRS'+str(n+1)+'    PNAM'
1060            PNames.append(Expr[n][key])
1061    if Num < 8:
1062        dlg = wx.SingleChoiceDialog(G2frame, 'Which phase to read?', 'Read phase data', PNames, wx.CHOICEDLG_STYLE)
1063        try:
1064            if dlg.ShowModal() == wx.ID_OK:
1065                result = dlg.GetSelection()
1066        finally:
1067            dlg.Destroy()       
1068    EXPphase = Expr[result]
1069    keyList = EXPphase.keys()
1070    keyList.sort()
1071    SGData = {}
1072    if NPhas[result] == '1':
1073        Ptype = 'nuclear'
1074    elif NPhas[result] in ['2','3']:
1075        Ptype = 'magnetic'
1076    elif NPhas[result] == '4':
1077        Ptype = 'macromolecular'
1078    elif NPhas[result] == '10':
1079        Ptype = 'Pawley'
1080    for key in keyList:
1081        if 'PNAM' in key:
1082           PhaseName = EXPphase[key].strip()
1083        elif 'ABC   ' in key:
1084            abc = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                       
1085        elif 'ANGLES' in key:
1086            angles = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                                               
1087        elif 'SG SYM' in key:
1088            SpGrp = EXPphase[key][:15].strip()
1089            E,SGData = G2spc.SpcGroup(SpGrp)
1090        elif 'OD    ' in key:
1091            SHdata = EXPphase[key].split() # may not have all 9 values
1092            SHvals = 9*[0]
1093            for i in range(9):
1094                try:
1095                    float(SHdata[i])
1096                    SHvals[i] = SHdata[i]
1097                except:
1098                    pass
1099            textureData['Order'] = int(SHvals[0])
1100            textureData['Model'] = shModels[int(SHvals[2])]
1101            textureData['Sample omega'] = [False,float(SHvals[6])]
1102            textureData['Sample chi'] = [False,float(SHvals[7])]
1103            textureData['Sample phi'] = [False,float(SHvals[8])]
1104            shNcof = int(SHvals[1])
1105    Atoms = []
1106    if Ptype == 'nuclear':
1107        for key in keyList:
1108            if 'AT' in key:
1109                if key[11:] == 'A':
1110                    S = EXPphase[key]
1111                elif key[11:] == 'B':
1112                    S += EXPphase[key]
1113                    Atom = [S[50:58].strip(),S[:10].strip(),'',
1114                        float(S[10:20]),float(S[20:30]),float(S[30:40]),
1115                        float(S[40:50]),'',int(S[60:62]),S[130:131]]
1116                    if Atom[9] == 'I':
1117                        Atom += [float(S[68:78]),0.,0.,0.,0.,0.,0.]
1118                    elif Atom[9] == 'A':
1119                        Atom += [0.0,float(S[68:78]),float(S[78:88]),
1120                            float(S[88:98]),float(S[98:108]),
1121                            float(S[108:118]),float(S[118:128])]
1122                    XYZ = Atom[3:6]
1123                    Atom[7],Atom[8] = G2spc.SytSym(XYZ,SGData)
1124                    Atom.append(ran.randint(0,sys.maxint))
1125                    Atoms.append(Atom)
1126    elif Ptype == 'macromolecular':
1127        for key in keyList:
1128            if 'AT' in key[6:8]:
1129                S = EXPphase[key]
1130                Atom = [S[56:60],S[50:54].strip().upper(),S[54:56],
1131                    S[46:51].strip(),S[:8].strip(),'',
1132                    float(S[16:24]),float(S[24:32]),float(S[32:40]),
1133                    float(S[8:16]),'1',1,'I',float(S[40:46]),0,0,0,0,0,0]
1134                XYZ = Atom[6:9]
1135                Atom[10],Atom[11] = G2spc.SytSym(XYZ,SGData)
1136                Atom.append(ran.randint(0,sys.maxint))
1137                Atoms.append(Atom)
1138    Volume = G2lat.calc_V(G2lat.cell2A(abc+angles))
1139    if shNcof:
1140        shCoef = {}
1141        nRec = [i+1 for i in range((shNcof-1)/6+1)]
1142        for irec in nRec:
1143            ODkey = keyList[0][:6]+'OD'+'%3dA'%(irec)
1144            indx = EXPphase[ODkey].split()
1145            ODkey = ODkey[:-1]+'B'
1146            vals = EXPphase[ODkey].split()
1147            for i,val in enumerate(vals):
1148                key = 'C(%s,%s,%s)'%(indx[3*i],indx[3*i+1],indx[3*i+2])
1149                shCoef[key] = float(val)
1150        textureData['SH Coeff'] = [False,shCoef]
1151       
1152    Phase = {
1153            'General':{
1154                'Name':PhaseName,
1155                'Type':Ptype,
1156                'SGData':SGData,
1157                'Cell':[False,]+abc+angles+[Volume,],
1158                'Pawley dmin':1.0},
1159            'Atoms':Atoms,
1160            'Drawing':{},
1161            'Histograms':{},
1162            'Pawley ref':[],
1163            'Models':{},
1164            'SH Texture':textureData
1165            }
1166           
1167    return Phase
1168       
1169def ReadPDBPhase(filename):
1170    EightPiSq = 8.*math.pi**2
1171    file = open(filename, 'Ur')
1172    Phase = {}
1173    Title = ''
1174    Compnd = ''
1175    Atoms = []
1176    A = np.zeros(shape=(3,3))
1177    S = file.readline()
1178    while S:
1179        Atom = []
1180        if 'TITLE' in S[:5]:
1181            Title = S[10:72].strip()
1182            S = file.readline()
1183        elif 'COMPND    ' in S[:10]:
1184            Compnd = S[10:72].strip()
1185            S = file.readline()
1186        elif 'CRYST' in S[:5]:
1187            abc = S[7:34].split()
1188            angles = S[34:55].split()
1189            cell=[float(abc[0]),float(abc[1]),float(abc[2]),
1190                float(angles[0]),float(angles[1]),float(angles[2])]
1191            Volume = G2lat.calc_V(G2lat.cell2A(cell))
1192            AA,AB = G2lat.cell2AB(cell)
1193            SpGrp = S[55:65]
1194            E,SGData = G2spc.SpcGroup(SpGrp)
1195            while E:
1196                print G2spc.SGErrors(E)
1197                dlg = wx.TextEntryDialog(None,
1198                    SpGrp[:-1]+' is invalid \nN.B.: make sure spaces separate axial fields in symbol',
1199                    'ERROR in space group symbol','',style=wx.OK)
1200                if dlg.ShowModal() == wx.ID_OK:
1201                    SpGrp = dlg.GetValue()
1202                    E,SGData = G2spc.SpcGroup(SpGrp)
1203                else:
1204                    return None
1205                dlg.Destroy()               
1206            SGlines = G2spc.SGPrint(SGData)
1207            for line in SGlines: print line
1208            S = file.readline()
1209        elif 'SCALE' in S[:5]:
1210            V = (S[10:41].split())
1211            A[int(S[5])-1] = [float(V[0]),float(V[1]),float(V[2])]
1212            S = file.readline()
1213        elif 'ATOM' in S[:4] or 'HETATM' in S[:6]:
1214            XYZ = [float(S[31:39]),float(S[39:47]),float(S[47:55])]
1215            XYZ = np.inner(AB,XYZ)
1216            XYZ = np.where(abs(XYZ)<0.00001,0,XYZ)
1217            SytSym,Mult = G2spc.SytSym(XYZ,SGData)
1218            Uiso = float(S[61:67])/EightPiSq
1219            Type = S[12:14].upper()
1220            if Type[0] in '123456789':
1221                Type = Type[1:]
1222            Atom = [S[22:27].strip(),S[17:20].upper(),S[20:22],
1223                S[12:17].strip(),Type.strip(),'',XYZ[0],XYZ[1],XYZ[2],
1224                float(S[55:61]),SytSym,Mult,'I',Uiso,0,0,0,0,0,0]
1225            S = file.readline()
1226            if 'ANISOU' in S[:6]:
1227                Uij = S[30:72].split()
1228                Uij = [float(Uij[0])/10000.,float(Uij[1])/10000.,float(Uij[2])/10000.,
1229                    float(Uij[3])/10000.,float(Uij[4])/10000.,float(Uij[5])/10000.]
1230                Atom = Atom[:14]+Uij
1231                Atom[12] = 'A'
1232                S = file.readline()
1233            Atom.append(ran.randint(0,sys.maxint))
1234            Atoms.append(Atom)
1235        else:           
1236            S = file.readline()
1237    file.close()
1238    if Title:
1239        PhaseName = Title
1240    elif Compnd:
1241        PhaseName = Compnd
1242    else:
1243        PhaseName = 'None'
1244    Phase['General'] = {'Name':PhaseName,'Type':'macromolecular','SGData':SGData,
1245        'Cell':[False,]+cell+[Volume,]}
1246    Phase['Atoms'] = Atoms
1247    Phase['Drawing'] = {}
1248    Phase['Histograms'] = {}
1249   
1250    return Phase
1251
1252######################################################################
1253# base classes for reading various types of data files
1254#   not used directly, only by subclassing
1255######################################################################
1256E,SGData = G2spc.SpcGroup('P 1') # data structure for default space group
1257class ImportPhase(object):
1258    '''Defines a base class for the reading of files with coordinates
1259    '''
1260    def __init__(self,
1261                 formatName,
1262                 longFormatName=None,
1263                 extensionlist=[],
1264                 strictExtension=False,
1265                 ):
1266        self.formatName = formatName # short string naming file type
1267        if longFormatName: # longer string naming file type
1268            self.longFormatName = longFormatName
1269        else:
1270            self.longFormatName = formatName
1271        # define extensions that are allowed for the file type
1272        # for windows, remove any extensions that are duplicate, if case is ignored
1273        if sys.platform == 'windows' and extensionlist:
1274            extensionlist = list(set([s.lower() for s in extensionlist]))
1275        self.extensionlist = extensionlist
1276        # If strictExtension is True, the file will not be read, unless
1277        # the extension matches one in the extensionlist
1278        self.strictExtension = strictExtension
1279        # define a default Phase structure
1280        self.Phase = {}
1281        for i in 'General', 'Atoms', 'Drawing', 'Histograms':
1282            self.Phase[i] = {}
1283        self.Phase['General']['Name'] = 'default'
1284        self.Phase['General']['Type'] = 'nuclear'
1285        self.Phase['General']['SGData'] = SGData
1286        self.Phase['General']['Cell'] = [
1287            False, # refinement flag
1288            1.,1.,1.,    # a,b,c
1289            90.,90.,90., # alpha, beta, gamma
1290            1.           # volume
1291            ]
1292        self.warnings = ''
1293        self.errors = ''
1294        #print 'created',self.__class__
1295
1296    def PhaseSelector(self, ChoiceList, ParentFrame=None,
1297                      title='Select a phase', size=None):
1298        ''' Provide a wx dialog to select a phase if the file contains more
1299        than one
1300        '''
1301        dlg = wx.SingleChoiceDialog(
1302            ParentFrame,
1303            title,
1304            'Phase Selection',
1305            ChoiceList,
1306            )
1307        if size: dlg.SetSize(size)
1308        if dlg.ShowModal() == wx.ID_OK:
1309            sel = dlg.GetSelection()
1310            dlg.Destroy()
1311            return sel
1312        else:
1313            dlg.Destroy()
1314            return None
1315
1316    def ShowBusy(self):
1317        wx.BeginBusyCursor()
1318
1319    def DoneBusy(self):
1320        wx.EndBusyCursor()
1321       
1322#    def Reader(self, filename, filepointer, ParentFrame=None):
1323#        '''This method must be supplied in the child class
1324#        it will read the file
1325#        '''
1326#        return True # if read OK
1327#        return False # if an error occurs
1328
1329    def ExtensionValidator(self, filename):
1330        '''This methods checks if the file has the correct extension
1331        Return False if this filename will not be supported by this reader
1332        Return True if the extension matches the list supplied by the reader
1333        Return None if the reader allows un-registered extensions
1334        '''
1335        if filename:
1336            ext = os.path.splitext(filename)[1]
1337            if sys.platform == 'windows': ext = ext.lower()
1338            if ext in self.extensionlist: return True
1339            if self.strictExtension: return False
1340        return None
1341
1342    def ContentsValidator(self, filepointer):
1343        '''This routine will attempt to determine if the file can be read
1344        with the current format.
1345        This will typically be overridden with a method that
1346        takes a quick scan of [some of]
1347        the file contents to do a "sanity" check if the file
1348        appears to match the selected format.
1349        Expected to be called via self.Validator()
1350        '''
1351        #filepointer.seek(0) # rewind the file pointer
1352        return True
1353
Note: See TracBrowser for help on using the repository browser.