source: trunk/GSASIIIO.py @ 526

Last change on this file since 526 was 526, checked in by vondreele, 10 years ago

make SetNewPhase? routine in GSASIIIO.py to initialize phase info
begin Fourier map search routine
fix atom plotting if no map

  • Property svn:keywords set to Date Author Revision URL Id
File size: 53.3 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-04-02 19:30:28 +0000 (Mon, 02 Apr 2012) $
6# $Author: vondreele $
7# $Revision: 526 $
8# $URL: trunk/GSASIIIO.py $
9# $Id: GSASIIIO.py 526 2012-04-02 19:30:28Z 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 SetNewPhase(Name='New Phase',SGData=G2spc.SpcGroup('P 1')[1],cell=[1.0,1.0,1.0,90.,90,90.,1.]):
1038    phaseData = {
1039        'General':{
1040            'Name':Name,
1041            'Type':'nuclear',
1042            'SGData':SGData,
1043            'Cell':[False,]+cell,
1044            'Pawley dmin':1.0,
1045            'Data plot type':'Mustrain',
1046            'SH Texture':{
1047                'Order':0,
1048                'Model':'cylindrical',
1049                'Sample omega':[False,0.0],
1050                'Sample chi':[False,0.0],
1051                'Sample phi':[False,0.0],
1052                'SH Coeff':[False,{}],
1053                'SHShow':False,
1054                'PFhkl':[0,0,1],
1055                'PFxyz':[0,0,1],
1056                'PlotType':'Pole figure'}},
1057        'Atoms':[],
1058        'Drawing':{},
1059        'Histograms':{},
1060        'Pawley ref':[],
1061        'Models':{},
1062        }
1063    return phaseData
1064   
1065def ReadEXPPhase(G2frame,filename):
1066    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1067    textureData = {'Order':0,'Model':'cylindrical','Sample omega':[False,0.0],
1068        'Sample chi':[False,0.0],'Sample phi':[False,0.0],'SH Coeff':[False,{}],
1069        'SHShow':False,'PFhkl':[0,0,1],'PFxyz':[0,0,1],'PlotType':'Pole figure'}
1070    shNcof = 0
1071    file = open(filename, 'Ur')
1072    S = 1
1073    Expr = [{},{},{},{},{},{},{},{},{}]
1074    while S:
1075        S = file.readline()
1076        if 'EXPR NPHAS' in S[:12]:
1077            Num = S[12:-1].count('0')
1078            NPhas = S[12:-1].split()
1079        if 'CRS' in S[:3]:
1080            N = int(S[3:4])-1
1081            Expr[N][S[:12]] = S[12:-1]
1082    file.close()
1083    PNames = []
1084    for n,N in enumerate(NPhas):
1085        if N != '0':
1086            result = n
1087            key = 'CRS'+str(n+1)+'    PNAM'
1088            PNames.append(Expr[n][key])
1089    if Num < 8:
1090        dlg = wx.SingleChoiceDialog(G2frame, 'Which phase to read?', 'Read phase data', PNames, wx.CHOICEDLG_STYLE)
1091        try:
1092            if dlg.ShowModal() == wx.ID_OK:
1093                result = dlg.GetSelection()
1094        finally:
1095            dlg.Destroy()       
1096    EXPphase = Expr[result]
1097    keyList = EXPphase.keys()
1098    keyList.sort()
1099    SGData = {}
1100    if NPhas[result] == '1':
1101        Ptype = 'nuclear'
1102    elif NPhas[result] in ['2','3']:
1103        Ptype = 'magnetic'
1104    elif NPhas[result] == '4':
1105        Ptype = 'macromolecular'
1106    elif NPhas[result] == '10':
1107        Ptype = 'Pawley'
1108    for key in keyList:
1109        if 'PNAM' in key:
1110           PhaseName = EXPphase[key].strip()
1111        elif 'ABC   ' in key:
1112            abc = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                       
1113        elif 'ANGLES' in key:
1114            angles = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                                               
1115        elif 'SG SYM' in key:
1116            SpGrp = EXPphase[key][:15].strip()
1117            E,SGData = G2spc.SpcGroup(SpGrp)
1118        elif 'OD    ' in key:
1119            SHdata = EXPphase[key].split() # may not have all 9 values
1120            SHvals = 9*[0]
1121            for i in range(9):
1122                try:
1123                    float(SHdata[i])
1124                    SHvals[i] = SHdata[i]
1125                except:
1126                    pass
1127            textureData['Order'] = int(SHvals[0])
1128            textureData['Model'] = shModels[int(SHvals[2])]
1129            textureData['Sample omega'] = [False,float(SHvals[6])]
1130            textureData['Sample chi'] = [False,float(SHvals[7])]
1131            textureData['Sample phi'] = [False,float(SHvals[8])]
1132            shNcof = int(SHvals[1])
1133    Atoms = []
1134    if Ptype == 'nuclear':
1135        for key in keyList:
1136            if 'AT' in key:
1137                if key[11:] == 'A':
1138                    S = EXPphase[key]
1139                elif key[11:] == 'B':
1140                    S += EXPphase[key]
1141                    Atom = [S[50:58].strip(),S[:10].strip(),'',
1142                        float(S[10:20]),float(S[20:30]),float(S[30:40]),
1143                        float(S[40:50]),'',int(S[60:62]),S[130:131]]
1144                    if Atom[9] == 'I':
1145                        Atom += [float(S[68:78]),0.,0.,0.,0.,0.,0.]
1146                    elif Atom[9] == 'A':
1147                        Atom += [0.0,float(S[68:78]),float(S[78:88]),
1148                            float(S[88:98]),float(S[98:108]),
1149                            float(S[108:118]),float(S[118:128])]
1150                    XYZ = Atom[3:6]
1151                    Atom[7],Atom[8] = G2spc.SytSym(XYZ,SGData)
1152                    Atom.append(ran.randint(0,sys.maxint))
1153                    Atoms.append(Atom)
1154    elif Ptype == 'macromolecular':
1155        for key in keyList:
1156            if 'AT' in key[6:8]:
1157                S = EXPphase[key]
1158                Atom = [S[56:60],S[50:54].strip().upper(),S[54:56],
1159                    S[46:51].strip(),S[:8].strip(),'',
1160                    float(S[16:24]),float(S[24:32]),float(S[32:40]),
1161                    float(S[8:16]),'1',1,'I',float(S[40:46]),0,0,0,0,0,0]
1162                XYZ = Atom[6:9]
1163                Atom[10],Atom[11] = G2spc.SytSym(XYZ,SGData)
1164                Atom.append(ran.randint(0,sys.maxint))
1165                Atoms.append(Atom)
1166    Volume = G2lat.calc_V(G2lat.cell2A(abc+angles))
1167    if shNcof:
1168        shCoef = {}
1169        nRec = [i+1 for i in range((shNcof-1)/6+1)]
1170        for irec in nRec:
1171            ODkey = keyList[0][:6]+'OD'+'%3dA'%(irec)
1172            indx = EXPphase[ODkey].split()
1173            ODkey = ODkey[:-1]+'B'
1174            vals = EXPphase[ODkey].split()
1175            for i,val in enumerate(vals):
1176                key = 'C(%s,%s,%s)'%(indx[3*i],indx[3*i+1],indx[3*i+2])
1177                shCoef[key] = float(val)
1178        textureData['SH Coeff'] = [False,shCoef]
1179       
1180    Phase = SetNewPhase(Name=PhaseName,SGData=SGData,cell=abc+angles+[Volume,])
1181    general = Phase['General']
1182    general['Type'] = Ptype
1183    general['SH Texture'] = textureData
1184    Phase['Atoms'] = Atoms
1185    return Phase
1186       
1187def ReadPDBPhase(filename):
1188    EightPiSq = 8.*math.pi**2
1189    file = open(filename, 'Ur')
1190    Phase = {}
1191    Title = ''
1192    Compnd = ''
1193    Atoms = []
1194    A = np.zeros(shape=(3,3))
1195    S = file.readline()
1196    while S:
1197        Atom = []
1198        if 'TITLE' in S[:5]:
1199            Title = S[10:72].strip()
1200            S = file.readline()
1201        elif 'COMPND    ' in S[:10]:
1202            Compnd = S[10:72].strip()
1203            S = file.readline()
1204        elif 'CRYST' in S[:5]:
1205            abc = S[7:34].split()
1206            angles = S[34:55].split()
1207            cell=[float(abc[0]),float(abc[1]),float(abc[2]),
1208                float(angles[0]),float(angles[1]),float(angles[2])]
1209            Volume = G2lat.calc_V(G2lat.cell2A(cell))
1210            AA,AB = G2lat.cell2AB(cell)
1211            SpGrp = S[55:65]
1212            E,SGData = G2spc.SpcGroup(SpGrp)
1213            while E:
1214                print G2spc.SGErrors(E)
1215                dlg = wx.TextEntryDialog(None,
1216                    SpGrp[:-1]+' is invalid \nN.B.: make sure spaces separate axial fields in symbol',
1217                    'ERROR in space group symbol','',style=wx.OK)
1218                if dlg.ShowModal() == wx.ID_OK:
1219                    SpGrp = dlg.GetValue()
1220                    E,SGData = G2spc.SpcGroup(SpGrp)
1221                else:
1222                    return None
1223                dlg.Destroy()               
1224            SGlines = G2spc.SGPrint(SGData)
1225            for line in SGlines: print line
1226            S = file.readline()
1227        elif 'SCALE' in S[:5]:
1228            V = (S[10:41].split())
1229            A[int(S[5])-1] = [float(V[0]),float(V[1]),float(V[2])]
1230            S = file.readline()
1231        elif 'ATOM' in S[:4] or 'HETATM' in S[:6]:
1232            XYZ = [float(S[31:39]),float(S[39:47]),float(S[47:55])]
1233            XYZ = np.inner(AB,XYZ)
1234            XYZ = np.where(abs(XYZ)<0.00001,0,XYZ)
1235            SytSym,Mult = G2spc.SytSym(XYZ,SGData)
1236            Uiso = float(S[61:67])/EightPiSq
1237            Type = S[12:14].upper()
1238            if Type[0] in '123456789':
1239                Type = Type[1:]
1240            Atom = [S[22:27].strip(),S[17:20].upper(),S[20:22],
1241                S[12:17].strip(),Type.strip(),'',XYZ[0],XYZ[1],XYZ[2],
1242                float(S[55:61]),SytSym,Mult,'I',Uiso,0,0,0,0,0,0]
1243            S = file.readline()
1244            if 'ANISOU' in S[:6]:
1245                Uij = S[30:72].split()
1246                Uij = [float(Uij[0])/10000.,float(Uij[1])/10000.,float(Uij[2])/10000.,
1247                    float(Uij[3])/10000.,float(Uij[4])/10000.,float(Uij[5])/10000.]
1248                Atom = Atom[:14]+Uij
1249                Atom[12] = 'A'
1250                S = file.readline()
1251            Atom.append(ran.randint(0,sys.maxint))
1252            Atoms.append(Atom)
1253        else:           
1254            S = file.readline()
1255    file.close()
1256    if Title:
1257        PhaseName = Title
1258    elif Compnd:
1259        PhaseName = Compnd
1260    else:
1261        PhaseName = 'None'
1262    Phase = SetNewPhase(Name=PhaseName,SGData=SGData,cell=cell+[Volume,])
1263    Phase['General']['Type'] = 'macromolecular'
1264    Phase['Atoms'] = Atoms
1265   
1266    return Phase
1267
1268######################################################################
1269# base classes for reading various types of data files
1270#   not used directly, only by subclassing
1271######################################################################
1272E,SGData = G2spc.SpcGroup('P 1') # data structure for default space group
1273class ImportPhase(object):
1274    '''Defines a base class for the reading of files with coordinates
1275    '''
1276    def __init__(self,
1277                 formatName,
1278                 longFormatName=None,
1279                 extensionlist=[],
1280                 strictExtension=False,
1281                 ):
1282        self.formatName = formatName # short string naming file type
1283        if longFormatName: # longer string naming file type
1284            self.longFormatName = longFormatName
1285        else:
1286            self.longFormatName = formatName
1287        # define extensions that are allowed for the file type
1288        # for windows, remove any extensions that are duplicate, as case is ignored
1289        if sys.platform == 'windows' and extensionlist:
1290            extensionlist = list(set([s.lower() for s in extensionlist]))
1291        self.extensionlist = extensionlist
1292        # If strictExtension is True, the file will not be read, unless
1293        # the extension matches one in the extensionlist
1294        self.strictExtension = strictExtension
1295        # define a default Phase structure
1296        self.Phase = SetNewPhase(Name='new phase',SGData=SGData)
1297        self.warnings = ''
1298        self.errors = ''
1299        #print 'created',self.__class__
1300
1301    def PhaseSelector(self, ChoiceList, ParentFrame=None,
1302                      title='Select a phase', size=None):
1303        ''' Provide a wx dialog to select a phase if the file contains more
1304        than one
1305        '''
1306        dlg = wx.SingleChoiceDialog(
1307            ParentFrame,
1308            title,
1309            'Phase Selection',
1310            ChoiceList,
1311            )
1312        if size: dlg.SetSize(size)
1313        if dlg.ShowModal() == wx.ID_OK:
1314            sel = dlg.GetSelection()
1315            dlg.Destroy()
1316            return sel
1317        else:
1318            dlg.Destroy()
1319            return None
1320
1321    def ShowBusy(self):
1322        wx.BeginBusyCursor()
1323
1324    def DoneBusy(self):
1325        wx.EndBusyCursor()
1326       
1327#    def Reader(self, filename, filepointer, ParentFrame=None):
1328#        '''This method must be supplied in the child class
1329#        it will read the file
1330#        '''
1331#        return True # if read OK
1332#        return False # if an error occurs
1333
1334    def ExtensionValidator(self, filename):
1335        '''This methods checks if the file has the correct extension
1336        Return False if this filename will not be supported by this reader
1337        Return True if the extension matches the list supplied by the reader
1338        Return None if the reader allows un-registered extensions
1339        '''
1340        if filename:
1341            ext = os.path.splitext(filename)[1]
1342            if sys.platform == 'windows': ext = ext.lower()
1343            if ext in self.extensionlist: return True
1344            if self.strictExtension: return False
1345        return None
1346
1347    def ContentsValidator(self, filepointer):
1348        '''This routine will attempt to determine if the file can be read
1349        with the current format.
1350        This will typically be overridden with a method that
1351        takes a quick scan of [some of]
1352        the file contents to do a "sanity" check if the file
1353        appears to match the selected format.
1354        Expected to be called via self.Validator()
1355        '''
1356        #filepointer.seek(0) # rewind the file pointer
1357        return True
1358
Note: See TracBrowser for help on using the repository browser.