source: trunk/GSASIIIO.py @ 573

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

add 'None' as one of the phase data plot options
trap failed charge flipping errors
add atom modify commands
add a new class - SingleFloatDialog? to GSASIIphsGUI.py
can now roll in place charge flip maps with u,d,l,r keys

  • Property svn:keywords set to Date Author Revision URL Id
File size: 53.2 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-25 16:54:23 +0000 (Wed, 25 Apr 2012) $
6# $Author: vondreele $
7# $Revision: 573 $
8# $URL: trunk/GSASIIIO.py $
9# $Id: GSASIIIO.py 573 2012-04-25 16:54:23Z 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    for i,export in enumerate(exports):
919        filename = ospath.join(head,name+'-%03d.'%(i)+ext)
920        prmname = filename.strip(ext)+'prm'
921        prm = open(prmname,'w')      #old style GSAS parm file
922        PickId = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, export)
923        Values,Names = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame, \
924            PickId, 'Instrument Parameters'))[1::2]     #get values & names
925        Inst = dict(zip(Names,Values))
926        print Inst['Type']
927        prm.write( '            123456789012345678901234567890123456789012345678901234567890        '+'\n')
928        prm.write( 'INS   BANK      1                                                               '+'\n')
929        prm.write(('INS   HTYPE   %sR                                                              '+'\n')%(Inst['Type']))
930        if 'Lam1' in Inst:              #Ka1 & Ka2
931            prm.write(('INS  1 ICONS%10.7f%10.7f    0.0000               0.990    0     0.500   '+'\n')%(Inst['Lam1'],Inst['Lam2']))
932        elif 'Lam' in Inst:             #single wavelength
933            prm.write(('INS  1 ICONS%10.7f%10.7f    0.0000               0.990    0     0.500   '+'\n')%(Inst['Lam'],0.0))
934        prm.write( 'INS  1 IRAD     0                                                               '+'\n')
935        prm.write( 'INS  1I HEAD                                                                    '+'\n')
936        prm.write( 'INS  1I ITYP    0    0.0000  180.0000         1                                 '+'\n')
937        prm.write(('INS  1DETAZM%10.3f                                                          '+'\n')%(Inst['Azimuth']))
938        prm.write( 'INS  1PRCF1     3    8   0.00100                                                '+'\n')
939        prm.write(('INS  1PRCF11     %15.6g%15.6g%15.6g%15.6g   '+'\n')%(Inst['U'],Inst['V'],Inst['W'],0.0))
940        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.))
941        prm.close()
942        file = open(filename,'w')
943        print 'save powder pattern to file: ',filename
944        x,y,w,yc,yb,yd = G2frame.PatternTree.GetItemPyData(PickId)[1]
945        file.write(powderfile+'\n')
946        file.write('Instrument parameter file:'+ospath.split(prmname)[1]+'\n')
947        file.write('BANK 1 %d %d CONS %.2f %.2f 0 0 FXYE\n'%(len(x),len(x),\
948            100.*x[0],100.*(x[1]-x[0])))
949        s = list(np.sqrt(1./np.array(w)))       
950        XYW = zip(x,y,s)
951        for X,Y,S in XYW:
952            file.write("%15.6g %15.6g %15.6g\n" % (100.*X,Y,max(S,1.0)))
953        file.close()
954        print 'powder pattern file '+filename+' written'
955       
956def powderXyeSave(G2frame,exports,powderfile):
957    head,tail = ospath.split(powderfile)
958    name,ext = tail.split('.')
959    for i,export in enumerate(exports):
960        filename = ospath.join(head,name+'-%03d.'%(i)+ext)
961        PickId = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, export)
962        file = open(filename,'w')
963        file.write('#%s\n'%(export))
964        print 'save powder pattern to file: ',filename
965        x,y,w,yc,yb,yd = G2frame.PatternTree.GetItemPyData(PickId)[1]
966        s = list(np.sqrt(1./np.array(w)))       
967        XYW = zip(x,y,s)
968        for X,Y,W in XYW:
969            file.write("%15.6g %15.6g %15.6g\n" % (X,Y,W))
970        file.close()
971        print 'powder pattern file '+filename+' written'
972       
973def PDFSave(G2frame,exports):   
974    for export in exports:
975        PickId = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, export)
976        SQname = 'S(Q)'+export[4:]
977        GRname = 'G(R)'+export[4:]
978        sqfilename = ospath.join(G2frame.dirname,export.replace(' ','_')[5:]+'.sq')
979        grfilename = ospath.join(G2frame.dirname,export.replace(' ','_')[5:]+'.gr')
980        sqId = G2gd.GetPatternTreeItemId(G2frame, PickId, SQname)
981        grId = G2gd.GetPatternTreeItemId(G2frame, PickId, GRname)
982        sqdata = np.array(G2frame.PatternTree.GetItemPyData(sqId)[1][:2]).T
983        grdata = np.array(G2frame.PatternTree.GetItemPyData(grId)[1][:2]).T
984        sqfile = open(sqfilename,'w')
985        grfile = open(grfilename,'w')
986        sqfile.write('#T S(Q) %s\n'%(export))
987        grfile.write('#T G(R) %s\n'%(export))
988        sqfile.write('#L Q     S(Q)\n')
989        grfile.write('#L R     G(R)\n')
990        for q,sq in sqdata:
991            sqfile.write("%15.6g %15.6g\n" % (q,sq))
992        sqfile.close()
993        for r,gr in grdata:
994            grfile.write("%15.6g %15.6g\n" % (r,gr))
995        grfile.close()
996   
997def PeakListSave(G2frame,file,peaks):
998    print 'save peak list to file: ',G2frame.peaklistfile
999    if not peaks:
1000        dlg = wx.MessageDialog(G2frame, 'No peaks!', 'Nothing to save!', wx.OK)
1001        try:
1002            result = dlg.ShowModal()
1003        finally:
1004            dlg.Destroy()
1005        return
1006    for peak in peaks:
1007        file.write("%10.4f %12.2f %10.3f %10.3f \n" % \
1008            (peak[0],peak[2],peak[4],peak[6]))
1009    print 'peak list saved'
1010             
1011def IndexPeakListSave(G2frame,peaks):
1012    file = open(G2frame.peaklistfile,'wa')
1013    print 'save index peak list to file: ',G2frame.peaklistfile
1014    wx.BeginBusyCursor()
1015    try:
1016        if not peaks:
1017            dlg = wx.MessageDialog(G2frame, 'No peaks!', 'Nothing to save!', wx.OK)
1018            try:
1019                result = dlg.ShowModal()
1020            finally:
1021                dlg.Destroy()
1022            return
1023        for peak in peaks:
1024            file.write("%12.6f\n" % (peak[7]))
1025        file.close()
1026    finally:
1027        wx.EndBusyCursor()
1028    print 'index peak list saved'
1029   
1030def SetNewPhase(Name='New Phase',SGData=G2spc.SpcGroup('P 1')[1],cell=[1.0,1.0,1.0,90.,90,90.,1.]):
1031    phaseData = {
1032        'General':{
1033            'Name':Name,
1034            'Type':'nuclear',
1035            'SGData':SGData,
1036            'Cell':[False,]+cell,
1037            'Pawley dmin':1.0,
1038            'Data plot type':'None',
1039            'SH Texture':{
1040                'Order':0,
1041                'Model':'cylindrical',
1042                'Sample omega':[False,0.0],
1043                'Sample chi':[False,0.0],
1044                'Sample phi':[False,0.0],
1045                'SH Coeff':[False,{}],
1046                'SHShow':False,
1047                'PFhkl':[0,0,1],
1048                'PFxyz':[0,0,1],
1049                'PlotType':'Pole figure'}},
1050        'Atoms':[],
1051        'Drawing':{},
1052        'Histograms':{},
1053        'Pawley ref':[],
1054        'Models':{},
1055        }
1056    return phaseData
1057   
1058def ReadEXPPhase(G2frame,filename):
1059    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1060    textureData = {'Order':0,'Model':'cylindrical','Sample omega':[False,0.0],
1061        'Sample chi':[False,0.0],'Sample phi':[False,0.0],'SH Coeff':[False,{}],
1062        'SHShow':False,'PFhkl':[0,0,1],'PFxyz':[0,0,1],'PlotType':'Pole figure'}
1063    shNcof = 0
1064    file = open(filename, 'Ur')
1065    S = 1
1066    Expr = [{},{},{},{},{},{},{},{},{}]
1067    while S:
1068        S = file.readline()
1069        if 'EXPR NPHAS' in S[:12]:
1070            Num = S[12:-1].count('0')
1071            NPhas = S[12:-1].split()
1072        if 'CRS' in S[:3]:
1073            N = int(S[3:4])-1
1074            Expr[N][S[:12]] = S[12:-1]
1075    file.close()
1076    PNames = []
1077    for n,N in enumerate(NPhas):
1078        if N != '0':
1079            result = n
1080            key = 'CRS'+str(n+1)+'    PNAM'
1081            PNames.append(Expr[n][key])
1082    if Num < 8:
1083        dlg = wx.SingleChoiceDialog(G2frame, 'Which phase to read?', 'Read phase data', PNames, wx.CHOICEDLG_STYLE)
1084        try:
1085            if dlg.ShowModal() == wx.ID_OK:
1086                result = dlg.GetSelection()
1087        finally:
1088            dlg.Destroy()       
1089    EXPphase = Expr[result]
1090    keyList = EXPphase.keys()
1091    keyList.sort()
1092    SGData = {}
1093    if NPhas[result] == '1':
1094        Ptype = 'nuclear'
1095    elif NPhas[result] in ['2','3']:
1096        Ptype = 'magnetic'
1097    elif NPhas[result] == '4':
1098        Ptype = 'macromolecular'
1099    elif NPhas[result] == '10':
1100        Ptype = 'Pawley'
1101    for key in keyList:
1102        if 'PNAM' in key:
1103           PhaseName = EXPphase[key].strip()
1104        elif 'ABC   ' in key:
1105            abc = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                       
1106        elif 'ANGLES' in key:
1107            angles = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                                               
1108        elif 'SG SYM' in key:
1109            SpGrp = EXPphase[key][:15].strip()
1110            E,SGData = G2spc.SpcGroup(SpGrp)
1111        elif 'OD    ' in key:
1112            SHdata = EXPphase[key].split() # may not have all 9 values
1113            SHvals = 9*[0]
1114            for i in range(9):
1115                try:
1116                    float(SHdata[i])
1117                    SHvals[i] = SHdata[i]
1118                except:
1119                    pass
1120            textureData['Order'] = int(SHvals[0])
1121            textureData['Model'] = shModels[int(SHvals[2])]
1122            textureData['Sample omega'] = [False,float(SHvals[6])]
1123            textureData['Sample chi'] = [False,float(SHvals[7])]
1124            textureData['Sample phi'] = [False,float(SHvals[8])]
1125            shNcof = int(SHvals[1])
1126    Atoms = []
1127    if Ptype == 'nuclear':
1128        for key in keyList:
1129            if 'AT' in key:
1130                if key[11:] == 'A':
1131                    S = EXPphase[key]
1132                elif key[11:] == 'B':
1133                    S += EXPphase[key]
1134                    Atom = [S[50:58].strip(),S[:10].strip(),'',
1135                        float(S[10:20]),float(S[20:30]),float(S[30:40]),
1136                        float(S[40:50]),'',int(S[60:62]),S[130:131]]
1137                    if Atom[9] == 'I':
1138                        Atom += [float(S[68:78]),0.,0.,0.,0.,0.,0.]
1139                    elif Atom[9] == 'A':
1140                        Atom += [0.0,float(S[68:78]),float(S[78:88]),
1141                            float(S[88:98]),float(S[98:108]),
1142                            float(S[108:118]),float(S[118:128])]
1143                    XYZ = Atom[3:6]
1144                    Atom[7],Atom[8] = G2spc.SytSym(XYZ,SGData)
1145                    Atom.append(ran.randint(0,sys.maxint))
1146                    Atoms.append(Atom)
1147    elif Ptype == 'macromolecular':
1148        for key in keyList:
1149            if 'AT' in key[6:8]:
1150                S = EXPphase[key]
1151                Atom = [S[56:60],S[50:54].strip().upper(),S[54:56],
1152                    S[46:51].strip(),S[:8].strip(),'',
1153                    float(S[16:24]),float(S[24:32]),float(S[32:40]),
1154                    float(S[8:16]),'1',1,'I',float(S[40:46]),0,0,0,0,0,0]
1155                XYZ = Atom[6:9]
1156                Atom[10],Atom[11] = G2spc.SytSym(XYZ,SGData)
1157                Atom.append(ran.randint(0,sys.maxint))
1158                Atoms.append(Atom)
1159    Volume = G2lat.calc_V(G2lat.cell2A(abc+angles))
1160    if shNcof:
1161        shCoef = {}
1162        nRec = [i+1 for i in range((shNcof-1)/6+1)]
1163        for irec in nRec:
1164            ODkey = keyList[0][:6]+'OD'+'%3dA'%(irec)
1165            indx = EXPphase[ODkey].split()
1166            ODkey = ODkey[:-1]+'B'
1167            vals = EXPphase[ODkey].split()
1168            for i,val in enumerate(vals):
1169                key = 'C(%s,%s,%s)'%(indx[3*i],indx[3*i+1],indx[3*i+2])
1170                shCoef[key] = float(val)
1171        textureData['SH Coeff'] = [False,shCoef]
1172       
1173    Phase = SetNewPhase(Name=PhaseName,SGData=SGData,cell=abc+angles+[Volume,])
1174    general = Phase['General']
1175    general['Type'] = Ptype
1176    general['SH Texture'] = textureData
1177    Phase['Atoms'] = Atoms
1178    return Phase
1179       
1180def ReadPDBPhase(filename):
1181    EightPiSq = 8.*math.pi**2
1182    file = open(filename, 'Ur')
1183    Phase = {}
1184    Title = ''
1185    Compnd = ''
1186    Atoms = []
1187    A = np.zeros(shape=(3,3))
1188    S = file.readline()
1189    while S:
1190        Atom = []
1191        if 'TITLE' in S[:5]:
1192            Title = S[10:72].strip()
1193            S = file.readline()
1194        elif 'COMPND    ' in S[:10]:
1195            Compnd = S[10:72].strip()
1196            S = file.readline()
1197        elif 'CRYST' in S[:5]:
1198            abc = S[7:34].split()
1199            angles = S[34:55].split()
1200            cell=[float(abc[0]),float(abc[1]),float(abc[2]),
1201                float(angles[0]),float(angles[1]),float(angles[2])]
1202            Volume = G2lat.calc_V(G2lat.cell2A(cell))
1203            AA,AB = G2lat.cell2AB(cell)
1204            SpGrp = S[55:65]
1205            E,SGData = G2spc.SpcGroup(SpGrp)
1206            while E:
1207                print G2spc.SGErrors(E)
1208                dlg = wx.TextEntryDialog(None,
1209                    SpGrp[:-1]+' is invalid \nN.B.: make sure spaces separate axial fields in symbol',
1210                    'ERROR in space group symbol','',style=wx.OK)
1211                if dlg.ShowModal() == wx.ID_OK:
1212                    SpGrp = dlg.GetValue()
1213                    E,SGData = G2spc.SpcGroup(SpGrp)
1214                else:
1215                    return None
1216                dlg.Destroy()               
1217            SGlines = G2spc.SGPrint(SGData)
1218            for line in SGlines: print line
1219            S = file.readline()
1220        elif 'SCALE' in S[:5]:
1221            V = (S[10:41].split())
1222            A[int(S[5])-1] = [float(V[0]),float(V[1]),float(V[2])]
1223            S = file.readline()
1224        elif 'ATOM' in S[:4] or 'HETATM' in S[:6]:
1225            XYZ = [float(S[31:39]),float(S[39:47]),float(S[47:55])]
1226            XYZ = np.inner(AB,XYZ)
1227            XYZ = np.where(abs(XYZ)<0.00001,0,XYZ)
1228            SytSym,Mult = G2spc.SytSym(XYZ,SGData)
1229            Uiso = float(S[61:67])/EightPiSq
1230            Type = S[12:14].upper()
1231            if Type[0] in '123456789':
1232                Type = Type[1:]
1233            Atom = [S[22:27].strip(),S[17:20].upper(),S[20:22],
1234                S[12:17].strip(),Type.strip(),'',XYZ[0],XYZ[1],XYZ[2],
1235                float(S[55:61]),SytSym,Mult,'I',Uiso,0,0,0,0,0,0]
1236            S = file.readline()
1237            if 'ANISOU' in S[:6]:
1238                Uij = S[30:72].split()
1239                Uij = [float(Uij[0])/10000.,float(Uij[1])/10000.,float(Uij[2])/10000.,
1240                    float(Uij[3])/10000.,float(Uij[4])/10000.,float(Uij[5])/10000.]
1241                Atom = Atom[:14]+Uij
1242                Atom[12] = 'A'
1243                S = file.readline()
1244            Atom.append(ran.randint(0,sys.maxint))
1245            Atoms.append(Atom)
1246        else:           
1247            S = file.readline()
1248    file.close()
1249    if Title:
1250        PhaseName = Title
1251    elif Compnd:
1252        PhaseName = Compnd
1253    else:
1254        PhaseName = 'None'
1255    Phase = SetNewPhase(Name=PhaseName,SGData=SGData,cell=cell+[Volume,])
1256    Phase['General']['Type'] = 'macromolecular'
1257    Phase['Atoms'] = Atoms
1258   
1259    return Phase
1260
1261######################################################################
1262# base classes for reading various types of data files
1263#   not used directly, only by subclassing
1264######################################################################
1265E,SGData = G2spc.SpcGroup('P 1') # data structure for default space group
1266class ImportPhase(object):
1267    '''Defines a base class for the reading of files with coordinates
1268    '''
1269    def __init__(self,
1270                 formatName,
1271                 longFormatName=None,
1272                 extensionlist=[],
1273                 strictExtension=False,
1274                 ):
1275        self.formatName = formatName # short string naming file type
1276        if longFormatName: # longer string naming file type
1277            self.longFormatName = longFormatName
1278        else:
1279            self.longFormatName = formatName
1280        # define extensions that are allowed for the file type
1281        # for windows, remove any extensions that are duplicate, as case is ignored
1282        if sys.platform == 'windows' and extensionlist:
1283            extensionlist = list(set([s.lower() for s in extensionlist]))
1284        self.extensionlist = extensionlist
1285        # If strictExtension is True, the file will not be read, unless
1286        # the extension matches one in the extensionlist
1287        self.strictExtension = strictExtension
1288        # define a default Phase structure
1289        self.Phase = SetNewPhase(Name='new phase',SGData=SGData)
1290        self.warnings = ''
1291        self.errors = ''
1292        #print 'created',self.__class__
1293
1294    def PhaseSelector(self, ChoiceList, ParentFrame=None,
1295                      title='Select a phase', size=None):
1296        ''' Provide a wx dialog to select a phase if the file contains more
1297        than one
1298        '''
1299        dlg = wx.SingleChoiceDialog(
1300            ParentFrame,
1301            title,
1302            'Phase Selection',
1303            ChoiceList,
1304            )
1305        if size: dlg.SetSize(size)
1306        if dlg.ShowModal() == wx.ID_OK:
1307            sel = dlg.GetSelection()
1308            dlg.Destroy()
1309            return sel
1310        else:
1311            dlg.Destroy()
1312            return None
1313
1314    def ShowBusy(self):
1315        wx.BeginBusyCursor()
1316
1317    def DoneBusy(self):
1318        wx.EndBusyCursor()
1319       
1320#    def Reader(self, filename, filepointer, ParentFrame=None):
1321#        '''This method must be supplied in the child class
1322#        it will read the file
1323#        '''
1324#        return True # if read OK
1325#        return False # if an error occurs
1326
1327    def ExtensionValidator(self, filename):
1328        '''This methods checks if the file has the correct extension
1329        Return False if this filename will not be supported by this reader
1330        Return True if the extension matches the list supplied by the reader
1331        Return None if the reader allows un-registered extensions
1332        '''
1333        if filename:
1334            ext = os.path.splitext(filename)[1]
1335            if sys.platform == 'windows': ext = ext.lower()
1336            if ext in self.extensionlist: return True
1337            if self.strictExtension: return False
1338        return None
1339
1340    def ContentsValidator(self, filepointer):
1341        '''This routine will attempt to determine if the file can be read
1342        with the current format.
1343        This will typically be overridden with a method that
1344        takes a quick scan of [some of]
1345        the file contents to do a "sanity" check if the file
1346        appears to match the selected format.
1347        Expected to be called via self.Validator()
1348        '''
1349        #filepointer.seek(0) # rewind the file pointer
1350        return True
1351
Note: See TracBrowser for help on using the repository browser.