source: trunk/GSASIIIO.py @ 580

Last change on this file since 580 was 580, checked in by toby, 10 years ago

finish import structure factor; refactor import classes

  • Property svn:keywords set to Date Author Revision URL Id
File size: 57.1 KB
Line 
1"""GSASIIIO: functions for IO of data
2   Copyright: 2008, Robert B. Von Dreele (Argonne National Laboratory)
3"""
4########### SVN repository information ###################
5# $Date: 2012-04-30 15:57:52 +0000 (Mon, 30 Apr 2012) $
6# $Author: toby $
7# $Revision: 580 $
8# $URL: trunk/GSASIIIO.py $
9# $Id: GSASIIIO.py 580 2012-04-30 15:57:52Z toby $
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   
236# this will be removed eventually
237def GetHKLData(filename):
238    print 'Reading: '+filename
239    File = open(filename,'Ur')
240    HKLref = []
241    HKLmin = [1000,1000,1000]
242    HKLmax = [0,0,0]
243    FoMax = 0
244    ifFc = False
245    S = File.readline()
246    while '#' in S[0]:        #get past comments if any
247        S = File.readline()       
248    if '_' in S:         #cif style .hkl file
249        while 'loop_' not in S:         #skip preliminaries if any - can't have 'loop_' in them!
250            S = File.readline()       
251        S = File.readline()             #get past 'loop_' line
252        pos = 0
253        hpos = kpos = lpos = Fosqpos = Fcsqpos = sigpos = -1
254        while S:
255            if '_' in S:
256                if 'index_h' in S:
257                    hpos = pos
258                elif 'index_k' in S:
259                    kpos = pos
260                elif 'index_l' in S:
261                    lpos = pos
262                elif 'F_squared_meas' in S:
263                    Fosqpos = pos
264                elif 'F_squared_calc' in S:
265                    Fcsqpos = pos
266                elif 'F_squared_sigma' in S:
267                    sigpos = pos
268                pos += 1
269            else:
270                data = S.split()
271                if data:                    #avoid blank lines
272                    HKL = np.array([int(data[hpos]),int(data[kpos]),int(data[lpos])])
273                    h,k,l = HKL
274                    Fosq = float(data[Fosqpos])
275                    if sigpos != -1:
276                        sigFosq = float(data[sigpos])
277                    else:
278                        sigFosq = 1.
279                    if Fcsqpos != -1:
280                        Fcsq = float(data[Fcsqpos])
281                        if Fcsq:
282                            ifFc = True
283                    else:
284                        Fcsq = 0.
285                       
286                    HKLmin = [min(h,HKLmin[0]),min(k,HKLmin[1]),min(l,HKLmin[2])]
287                    HKLmax = [max(h,HKLmax[0]),max(k,HKLmax[1]),max(l,HKLmax[2])]
288                    FoMax = max(FoMax,Fosq)
289                    HKLref.append([HKL,Fosq,sigFosq,Fcsq,0,0,0])                 #room for Fcp, Fcpp & phase
290            S = File.readline()
291    else:                   #dumb h,k,l,Fo,sigFo .hkl file
292        while S:
293            h,k,l,Fo,sigFo = S.split()
294            HKL = np.array([int(h),int(k),int(l)])
295            h,k,l = HKL
296            Fo = float(Fo)
297            sigFo = float(sigFo)
298            HKLmin = [min(h,HKLmin[0]),min(k,HKLmin[1]),min(l,HKLmin[2])]
299            HKLmax = [max(h,HKLmax[0]),max(k,HKLmax[1]),max(l,HKLmax[2])]
300            FoMax = max(FoMax,Fo)
301            HKLref.append([HKL,Fo**2,2.*Fo*sigFo,0,0,0,0])                 #room for Fc, Fcp, Fcpp & phase
302            S = File.readline()
303    File.close()
304    return HKLref,HKLmin,HKLmax,FoMax,ifFc
305
306def GetPowderData(filename,Pos,Bank,DataType):
307    '''Reads one BANK of data from GSAS raw powder data file
308    input:
309    filename: GSAS raw powder file dataname
310    Pos: start of data in file just after BANK record
311    Bank: the BANK record
312    DataType: powder data type, e.g. "PXC" for Powder X-ray CW data
313    returns: list [x,y,e,yc,yb]
314    x: np.array of x-axis values
315    y: np.array of powder pattern intensities
316    w: np.array of w=sig(intensity)^2 values
317    yc: np.array of calc. intensities (zero)
318    yb: np.array of calc. background (zero)
319    yd: np.array of obs-calc profiles
320    '''
321    print 'Reading: '+filename
322    print 'Bank:    '+Bank[:-1]
323    if 'FXYE' in Bank:
324        return GetFXYEdata(filename,Pos,Bank,DataType)
325    elif ' XYE' in Bank:
326        return GetXYEdata(filename,Pos,Bank,DataType)
327    elif 'FXY' in Bank:
328        return GetFXYdata(filename,Pos,Bank,DataType)
329    elif 'ESD' in Bank:
330        return GetESDdata(filename,Pos,Bank,DataType)
331    elif 'STD' in Bank:
332        return GetSTDdata(filename,Pos,Bank,DataType)
333    else:
334        return GetSTDdata(filename,Pos,Bank,DataType)
335    return []
336
337def GetFXYEdata(filename,Pos,Bank,DataType):
338    File = open(filename,'Ur')
339    File.seek(Pos)
340    x = []
341    y = []
342    w = []
343    S = File.readline()
344    while S and S[:4] != 'BANK':
345        vals = S.split()
346        if DataType[2] == 'C':
347            x.append(float(vals[0])/100.)               #CW: from centidegrees to degrees
348        elif DataType[2] == 'T':
349            x.append(float(vals[0])/1000.0)             #TOF: from musec to millisec
350        f = float(vals[1])
351        if f <= 0.0:
352            y.append(0.0)
353            w.append(1.0)
354        else:
355            y.append(float(vals[1]))
356            w.append(1.0/float(vals[2])**2)
357        S = File.readline()
358    File.close()
359    N = len(x)
360    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
361   
362def GetXYEdata(filename,Pos,Bank,DataType):
363    File = open(filename,'Ur')
364    File.seek(Pos)
365    x = []
366    y = []
367    w = []
368    S = File.readline()
369    while S:
370        vals = S.split()
371        try:
372            x.append(float(vals[0]))
373            f = float(vals[1])
374            if f <= 0.0:
375                y.append(0.0)
376                w.append(1.0)
377            else:
378                y.append(float(vals[1]))
379                w.append(1.0/float(vals[2])**2)
380            S = File.readline()
381        except ValueError:
382            break
383    File.close()
384    N = len(x)
385    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
386   
387   
388def GetFXYdata(filename,Pos,Bank,DataType):
389    File = open(filename,'Ur')
390    File.seek(Pos)
391    x = []
392    y = []
393    w = []
394    S = File.readline()
395    while S and S[:4] != 'BANK':
396        vals = S.split()
397        if DataType[2] == 'C':
398            x.append(float(vals[0])/100.)               #CW: from centidegrees to degrees
399        elif DataType[2] == 'T':
400            x.append(float(vals[0])/1000.0)             #TOF: from musec to millisec
401        f = float(vals[1])
402        if f > 0.0:
403            y.append(f)
404            w.append(1.0/f)
405        else:             
406            y.append(0.0)
407            w.append(1.0)
408        S = File.readline()
409    File.close()
410    N = len(x)
411    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
412   
413def GetESDdata(filename,Pos,Bank,DataType):
414    File = open(filename,'Ur')
415    cons = Bank.split()
416    if DataType[2] == 'C':
417        start = float(cons[5])/100.0               #CW: from centidegrees to degrees
418        step = float(cons[6])/100.0
419    elif DataType[2] == 'T':
420        start = float(cons[5])/1000.0              #TOF: from musec to millisec
421        step = float(cons[6])/1000.0
422    File.seek(Pos)
423    x = []
424    y = []
425    w = []
426    S = File.readline()
427    j = 0
428    while S and S[:4] != 'BANK':
429        for i in range(0,80,16):
430            xi = start+step*j
431            yi = sfloat(S[i:i+8])
432            ei = sfloat(S[i+8:i+16])
433            x.append(xi)
434            if yi > 0.0:
435                y.append(yi)
436                w.append(1.0/ei**2)
437            else:             
438                y.append(0.0)
439                w.append(1.0)
440            j += 1
441        S = File.readline()
442    File.close()
443    N = len(x)
444    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
445
446def GetSTDdata(filename,Pos,Bank,DataType):
447    File = open(filename,'Ur')
448    cons = Bank.split()
449    Nch = int(cons[2])
450    if DataType[2] == 'C':
451        start = float(cons[5])/100.0               #CW: from centidegrees to degrees
452        step = float(cons[6])/100.0
453    elif DataType[2] == 'T':
454        start = float(cons[5])/1000.0              #TOF: from musec to millisec - not likely!
455        step = float(cons[6])/1000.0
456    File.seek(Pos)
457    x = []
458    y = []
459    w = []
460    S = File.readline()
461    j = 0
462    while S and S[:4] != 'BANK':
463        for i in range(0,80,8):
464            xi = start+step*j
465            ni = max(sint(S[i:i+2]),1)
466            yi = max(sfloat(S[i+2:i+8]),0.0)
467            if yi:
468                vi = yi/ni
469            else:
470                yi = 0.0
471                vi = 1.0
472            j += 1
473            if j < Nch:
474                x.append(xi)
475                y.append(yi)
476                w.append(1.0/vi)
477        S = File.readline()
478    File.close()
479    N = len(x)
480    return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
481   
482def CheckImageFile(G2frame,imagefile):
483    if not ospath.exists(imagefile):
484        dlg = wx.FileDialog(G2frame, 'Bad image file name; choose name', '.', '',\
485        'Any image file (*.tif;*.tiff;*.mar*;*.avg;*.sum;*.img)\
486        |*.tif;*.tiff;*.mar*;*.avg;*.sum;*.img|\
487        Any detector tif (*.tif;*.tiff)|*.tif;*.tiff|\
488        MAR file (*.mar*)|*.mar*|\
489        GE Image (*.avg;*.sum)|*.avg;*.sum|\
490        ADSC Image (*.img)|*.img|\
491        All files (*.*)|*.*',wx.OPEN|wx.CHANGE_DIR)
492        try:
493            dlg.SetFilename(''+ospath.split(imagefile)[1])
494            if dlg.ShowModal() == wx.ID_OK:
495                imagefile = dlg.GetPath()
496            else:
497                imagefile = False
498        finally:
499            dlg.Destroy()
500    return imagefile
501       
502def GetImageData(G2frame,imagefile,imageOnly=False):       
503    ext = ospath.splitext(imagefile)[1]
504    Comments = []
505    if ext == '.tif' or ext == '.tiff':
506        Comments,Data,Npix,Image = GetTifData(imagefile)
507    elif ext == '.img':
508        Comments,Data,Npix,Image = GetImgData(imagefile)
509        Image[0][0] = 0
510    elif ext == '.mar3450' or ext == '.mar2300':
511        Comments,Data,Npix,Image = GetMAR345Data(imagefile)
512    elif ext in ['.sum','.avg','']:
513        Comments,Data,Npix,Image = GetGEsumData(imagefile)
514    elif ext == '.G2img':
515        Comments,Data,Npix,Image = GetG2Image(imagefile)
516    if imageOnly:
517        return Image
518    else:
519        return Comments,Data,Npix,Image
520       
521def PutG2Image(filename,Comments,Data,Npix,image):
522    File = open(filename,'wb')
523    cPickle.dump([Comments,Data,Npix,image],File,1)
524    File.close()
525    return
526   
527def GetG2Image(filename):
528    File = open(filename,'rb')
529    Comments,Data,Npix,image = cPickle.load(File)
530    File.close()
531    return Comments,Data,Npix,image
532   
533def GetGEsumData(filename,imageOnly=False):
534    import struct as st
535    import array as ar
536    if not imageOnly:
537        print 'Read GE sum file: ',filename   
538    File = open(filename,'rb')
539    if '.sum' in filename:
540        head = ['GE detector sum data from APS 1-ID',]
541        sizexy = [2048,2048]
542    elif '.avg' in filename:
543        head = ['GE detector avg data from APS 1-ID',]
544        sizexy = [2048,2048]
545    else:
546        head = ['GE detector raw data from APS 1-ID',]
547        File.seek(18)
548        size,nframes = st.unpack('<ih',File.read(6))
549        sizexy = [2048,2048]
550        pos = 8192
551        File.seek(pos)
552    Npix = sizexy[0]*sizexy[1]
553    if '.sum' in filename:
554        image = np.array(ar.array('f',File.read(4*Npix)),dtype=np.int32)
555    elif '.avg' in filename:
556        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
557    else:
558        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
559        while nframes > 1:
560            image += np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
561            nframes -= 1
562    image = np.reshape(image,(sizexy[1],sizexy[0]))
563    data = {'pixelSize':(200,200),'wavelength':0.15,'distance':250.0,'center':[204.8,204.8],'size':sizexy} 
564    File.close()   
565    if imageOnly:
566        return image
567    else:
568        return head,data,Npix,image
569       
570def GetImgData(filename,imageOnly=False):
571    import struct as st
572    import array as ar
573    if not imageOnly:
574        print 'Read ADSC img file: ',filename
575    File = open(filename,'rb')
576    head = File.read(511)
577    lines = head.split('\n')
578    head = []
579    center = [0,0]
580    for line in lines[1:-2]:
581        line = line.strip()[:-1]
582        if line:
583            if 'SIZE1' in line:
584                size = int(line.split('=')[1])
585                Npix = size*size
586            elif 'WAVELENGTH' in line:
587                wave = float(line.split('=')[1])
588            elif 'BIN' in line:
589                if line.split('=')[1] == '2x2':
590                    pixel=(102,102)
591                else:
592                    pixel = (51,51)
593            elif 'DISTANCE' in line:
594                distance = float(line.split('=')[1])
595            elif 'CENTER_X' in line:
596                center[0] = float(line.split('=')[1])
597            elif 'CENTER_Y' in line:
598                center[1] = float(line.split('=')[1])
599            head.append(line)
600    data = {'pixelSize':pixel,'wavelength':wave,'distance':distance,'center':center,'size':[size,size]}
601    image = []
602    row = 0
603    pos = 512
604    File.seek(pos)
605    image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
606    image = np.reshape(image,(sizexy[1],sizexy[0]))
607#    image = np.zeros(shape=(size,size),dtype=np.int32)   
608#    while row < size:
609#        File.seek(pos)
610#        line = ar.array('H',File.read(2*size))
611#        image[row] = np.asarray(line)
612#        row += 1
613#        pos += 2*size
614    File.close()
615    if imageOnly:
616        return image
617    else:
618        return lines[1:-2],data,Npix,image
619       
620def GetMAR345Data(filename,imageOnly=False):
621    import array as ar
622    import struct as st
623    try:
624        import pack_f as pf
625    except:
626        msg = wx.MessageDialog(None, message="Unable to load the GSAS MAR image decompression, pack_f",
627                               caption="Import Error",
628                               style=wx.ICON_ERROR | wx.OK | wx.STAY_ON_TOP)
629        msg.ShowModal()
630        return None,None,None,None
631
632    if not imageOnly:
633        print 'Read Mar345 file: ',filename
634    File = open(filename,'rb')
635    head = File.read(4095)
636    numbers = st.unpack('<iiiiiiiiii',head[:40])
637    lines = head[128:].split('\n')
638    head = []
639    for line in lines:
640        line = line.strip()
641        if 'PIXEL' in line:
642            values = line.split()
643            pixel = (int(values[2]),int(values[4]))     #in microns
644        elif 'WAVELENGTH' in line:
645            wave = float(line.split()[1])
646        elif 'DISTANCE' in line:
647            distance = float(line.split()[1])           #in mm
648        elif 'CENTER' in line:
649            values = line.split()
650            center = [float(values[2])/10.,float(values[4])/10.]    #make in mm from pixels
651        if line: 
652            head.append(line)
653    data = {'pixelSize':pixel,'wavelength':wave,'distance':distance,'center':center}
654    for line in head:
655        if 'FORMAT' in line[0:6]:
656            items = line.split()
657            size = int(items[1])
658            Npix = size*size
659    pos = 4096
660    data['size'] = [size,size]
661    File.seek(pos)
662    line = File.read(8)
663    while 'CCP4' not in line:       #get past overflow list for now
664        line = File.read(8)
665        pos += 8
666    pos += 37
667    File.seek(pos)
668    raw = File.read()
669    File.close()
670    image = np.zeros(shape=(size,size),dtype=np.int32)
671    image = pf.pack_f(len(raw),raw,size,image)
672    if imageOnly:
673        return image.T              #transpose to get it right way around
674    else:
675        return head,data,Npix,image.T
676       
677def GetTifData(filename,imageOnly=False):
678    import struct as st
679    import array as ar
680    File = open(filename,'rb')
681    dataType = 5
682    try:
683        Meta = open(filename+'.metadata','Ur')
684        head = Meta.readlines()
685        for line in head:
686            line = line.strip()
687            if 'dataType=' in line:
688                dataType = int(line.split('=')[1])
689        Meta.close()
690    except IOError:
691        print 'no metadata file found - will try to read file anyway'
692        head = ['no metadata file found',]
693       
694    tag = File.read(2)
695    byteOrd = '<'
696    if tag == 'II' and int(st.unpack('<h',File.read(2))[0]) == 42:     #little endian
697        IFD = int(st.unpack(byteOrd+'i',File.read(4))[0])
698    elif tag == 'MM' and int(st.unpack('>h',File.read(2))[0]) == 42:   #big endian
699        byteOrd = '>'
700        IFD = int(st.unpack(byteOrd+'i',File.read(4))[0])       
701    else:
702        lines = ['not a detector tiff file',]
703        return lines,0,0,0
704    File.seek(IFD)                                                  #get number of directory entries
705    NED = int(st.unpack(byteOrd+'h',File.read(2))[0])
706    IFD = {}
707    for ied in range(NED):
708        Tag,Type = st.unpack(byteOrd+'Hh',File.read(4))
709        nVal = st.unpack(byteOrd+'i',File.read(4))[0]
710        if Type == 1:
711            Value = st.unpack(byteOrd+nVal*'b',File.read(nVal))
712        elif Type == 2:
713            Value = st.unpack(byteOrd+'i',File.read(4))
714        elif Type == 3:
715            Value = st.unpack(byteOrd+nVal*'h',File.read(nVal*2))
716            x = st.unpack(byteOrd+nVal*'h',File.read(nVal*2))
717        elif Type == 4:
718            Value = st.unpack(byteOrd+nVal*'i',File.read(nVal*4))
719        elif Type == 5:
720            Value = st.unpack(byteOrd+nVal*'i',File.read(nVal*4))
721        elif Type == 11:
722            Value = st.unpack(byteOrd+nVal*'f',File.read(nVal*4))
723        IFD[Tag] = [Type,nVal,Value]
724#    for key in IFD:
725#        print key,IFD[key]
726    sizexy = [IFD[256][2][0],IFD[257][2][0]]
727    [nx,ny] = sizexy
728    Npix = nx*ny
729    if 272 in IFD:
730        ifd = IFD[272]
731        File.seek(ifd[2][0])
732        S = File.read(ifd[1])
733        if 'PILATUS' in S:
734            tifType = 'Pilatus'
735            dataType = 0
736            pixy = (172,172)
737            File.seek(4096)
738            if not imageOnly:
739                print 'Read Pilatus tiff file: ',filename
740            image = ar.array('L',File.read(4*Npix))
741            image = np.array(np.asarray(image),dtype=np.int32)
742    elif 262 in IFD and IFD[262][2][0] > 4:
743        tifType = 'DND'
744        pixy = (158,158)
745        File.seek(512)
746        if not imageOnly:
747            print 'Read DND SAX/WAX-detector tiff file: ',filename
748        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
749    elif sizexy == [1536,1536]:
750        tifType = 'APS Gold'
751        pixy = (150,150)
752        File.seek(64)
753        if not imageOnly:
754            print 'Read Gold tiff file:',filename
755        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
756    elif sizexy == [2048,2048] or sizexy == [1024,1024]:
757        if IFD[273][2][0] == 8:
758            if IFD[258][2][0] == 32:
759                tifType = 'PE'
760                pixy = (200,200)
761                File.seek(8)
762                if not imageOnly:
763                    print 'Read APS PE-detector tiff file: ',filename
764                if dataType == 5:
765                    image = np.array(ar.array('f',File.read(4*Npix)),dtype=np.float32)
766                else:
767                    image = np.array(ar.array('I',File.read(4*Npix)),dtype=np.int32)
768        elif IFD[273][2][0] == 4096:
769            tifType = 'MAR'
770            pixy = (158,158)
771            File.seek(4096)
772            if not imageOnly:
773                print 'Read MAR CCD tiff file: ',filename
774            image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
775    elif sizexy == [4096,4096]:
776        if IFD[273][2][0] == 8:
777            if IFD[258][2][0] == 16:
778                tifType = 'scanCCD'
779                pixy = (9,9)
780                File.seek(8)
781                if not imageOnly:
782                    print 'Read APS scanCCD tiff file: ',filename
783                image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
784#    elif sizexy == [960,960]:
785#        tiftype = 'PE-BE'
786#        pixy = (200,200)
787#        File.seek(8)
788#        if not imageOnly:
789#            print 'Read Gold tiff file:',filename
790#        image = np.array(ar.array('H',File.read(2*Npix)),dtype=np.int32)
791           
792    else:
793        lines = ['not a known detector tiff file',]
794        return lines,0,0,0
795       
796    image = np.reshape(image,(sizexy[1],sizexy[0]))
797    center = [pixy[0]*sizexy[0]/2000,pixy[1]*sizexy[1]/2000]
798    data = {'pixelSize':pixy,'wavelength':0.10,'distance':100.0,'center':center,'size':sizexy}
799    File.close()   
800    if imageOnly:
801        return image
802    else:
803        return head,data,Npix,image
804   
805def ProjFileOpen(G2frame):
806    file = open(G2frame.GSASprojectfile,'rb')
807    print 'load from file: ',G2frame.GSASprojectfile
808    wx.BeginBusyCursor()
809    try:
810        while True:
811            try:
812                data = cPickle.load(file)
813            except EOFError:
814                break
815            datum = data[0]
816           
817            Id = G2frame.PatternTree.AppendItem(parent=G2frame.root,text=datum[0])
818            if 'PWDR' in datum[0]:               
819                G2frame.PatternTree.SetItemPyData(Id,datum[1][:3])     #temp. trim off junk
820            else:
821                G2frame.PatternTree.SetItemPyData(Id,datum[1])
822            for datus in data[1:]:
823                sub = G2frame.PatternTree.AppendItem(Id,datus[0])
824                G2frame.PatternTree.SetItemPyData(sub,datus[1])
825            if 'IMG' in datum[0]:                   #retreive image default flag & data if set
826                Data = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Image Controls'))
827                if Data['setDefault']:
828                    G2frame.imageDefault = Data               
829        file.close()
830       
831    finally:
832        wx.EndBusyCursor()
833    print 'project load successful'
834    G2frame.NewPlot = True
835   
836def ProjFileSave(G2frame):
837    if not G2frame.PatternTree.IsEmpty():
838        file = open(G2frame.GSASprojectfile,'wb')
839        print 'save to file: ',G2frame.GSASprojectfile
840        wx.BeginBusyCursor()
841        try:
842            item, cookie = G2frame.PatternTree.GetFirstChild(G2frame.root)
843            while item:
844                data = []
845                name = G2frame.PatternTree.GetItemText(item)
846                data.append([name,G2frame.PatternTree.GetItemPyData(item)])
847                item2, cookie2 = G2frame.PatternTree.GetFirstChild(item)
848                while item2:
849                    name = G2frame.PatternTree.GetItemText(item2)
850                    data.append([name,G2frame.PatternTree.GetItemPyData(item2)])
851                    item2, cookie2 = G2frame.PatternTree.GetNextChild(item, cookie2)                           
852                item, cookie = G2frame.PatternTree.GetNextChild(G2frame.root, cookie)                           
853                cPickle.dump(data,file,1)
854            file.close()
855        finally:
856            wx.EndBusyCursor()
857        print 'project save successful'
858
859def SaveIntegration(G2frame,PickId,data):
860    azms = G2frame.Integrate[1]
861    X = G2frame.Integrate[2][:-1]
862    Xminmax = [X[0],X[-1]]
863    N = len(X)
864    Id = G2frame.PatternTree.GetItemParent(PickId)
865    name = G2frame.PatternTree.GetItemText(Id)
866    name = name.replace('IMG ','PWDR ')
867    Comments = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id, 'Comments'))
868    names = ['Type','Lam','Zero','Polariz.','U','V','W','X','Y','SH/L','Azimuth'] 
869    codes = [0 for i in range(11)]
870    LRazm = data['LRazimuth']
871    Azms = []
872    if data['fullIntegrate'] and data['outAzimuths'] == 1:
873        Azms = [45.0,]                              #a poor man's average?
874    else:
875        for i,azm in enumerate(azms[:-1]):
876            Azms.append((azms[i+1]+azm)/2.)
877    for i,azm in enumerate(azms[:-1]):
878        item, cookie = G2frame.PatternTree.GetFirstChild(G2frame.root)
879        Id = 0
880        while item:
881            Name = G2frame.PatternTree.GetItemText(item)
882            if name == Name:
883                Id = item
884            item, cookie = G2frame.PatternTree.GetNextChild(G2frame.root, cookie)
885        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!
886        Y = G2frame.Integrate[0][i]
887        W = 1./Y                    #probably not true
888        Sample = G2pdG.SetDefaultSample()
889        if Id:
890            G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id, 'Comments'),Comments)                   
891            G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Limits'),[tuple(Xminmax),Xminmax])
892            G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Background'),[['chebyschev',1,3,1.0,0.0,0.0],
893                            {'nDebye':0,'debyeTerms':[],'nPeaks':0,'peaksList':[]}])
894            G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Instrument Parameters'),[tuple(parms),parms[:],codes,names])
895            G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Peak List'),[])
896            G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Index Peak List'),[])
897            G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Unit Cells List'),[])             
898            G2frame.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(G2frame,Id,'Reflection Lists'),{})             
899        else:
900            Id = G2frame.PatternTree.AppendItem(parent=G2frame.root,text=name+" Azm= %.2f"%(Azms[i]))
901            G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Comments'),Comments)                   
902            G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Limits'),[tuple(Xminmax),Xminmax])
903            G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Background'),[['chebyschev',1,3,1.0,0.0,0.0],
904                            {'nDebye':0,'debyeTerms':[],'nPeaks':0,'peaksList':[]}])
905            G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Instrument Parameters'),[tuple(parms),parms[:],codes,names])
906            G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Sample Parameters'),Sample)
907            G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Peak List'),[])
908            G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Index Peak List'),[])
909            G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Unit Cells List'),[])
910            G2frame.PatternTree.SetItemPyData(G2frame.PatternTree.AppendItem(Id,text='Reflection Lists'),{})             
911        G2frame.PatternTree.SetItemPyData(Id,[[''],[np.array(X),np.array(Y),np.array(W),np.zeros(N),np.zeros(N),np.zeros(N)]])
912    G2frame.PatternTree.SelectItem(Id)
913    G2frame.PatternTree.Expand(Id)
914    G2frame.PatternId = Id
915           
916def powderFxyeSave(G2frame,exports,powderfile):
917    head,tail = ospath.split(powderfile)
918    name,ext = tail.split('.')
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        x,y,w,yc,yb,yd = G2frame.PatternTree.GetItemPyData(PickId)[1]
946        file.write(powderfile+'\n')
947        file.write('Instrument parameter file:'+ospath.split(prmname)[1]+'\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        print 'powder pattern file '+filename+' written'
956       
957def powderXyeSave(G2frame,exports,powderfile):
958    head,tail = ospath.split(powderfile)
959    name,ext = tail.split('.')
960    for i,export in enumerate(exports):
961        filename = ospath.join(head,name+'-%03d.'%(i)+ext)
962        PickId = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, export)
963        file = open(filename,'w')
964        file.write('#%s\n'%(export))
965        print 'save powder pattern to file: ',filename
966        x,y,w,yc,yb,yd = G2frame.PatternTree.GetItemPyData(PickId)[1]
967        s = list(np.sqrt(1./np.array(w)))       
968        XYW = zip(x,y,s)
969        for X,Y,W in XYW:
970            file.write("%15.6g %15.6g %15.6g\n" % (X,Y,W))
971        file.close()
972        print 'powder pattern file '+filename+' written'
973       
974def PDFSave(G2frame,exports):   
975    for export in exports:
976        PickId = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, export)
977        SQname = 'S(Q)'+export[4:]
978        GRname = 'G(R)'+export[4:]
979        sqfilename = ospath.join(G2frame.dirname,export.replace(' ','_')[5:]+'.sq')
980        grfilename = ospath.join(G2frame.dirname,export.replace(' ','_')[5:]+'.gr')
981        sqId = G2gd.GetPatternTreeItemId(G2frame, PickId, SQname)
982        grId = G2gd.GetPatternTreeItemId(G2frame, PickId, GRname)
983        sqdata = np.array(G2frame.PatternTree.GetItemPyData(sqId)[1][:2]).T
984        grdata = np.array(G2frame.PatternTree.GetItemPyData(grId)[1][:2]).T
985        sqfile = open(sqfilename,'w')
986        grfile = open(grfilename,'w')
987        sqfile.write('#T S(Q) %s\n'%(export))
988        grfile.write('#T G(R) %s\n'%(export))
989        sqfile.write('#L Q     S(Q)\n')
990        grfile.write('#L R     G(R)\n')
991        for q,sq in sqdata:
992            sqfile.write("%15.6g %15.6g\n" % (q,sq))
993        sqfile.close()
994        for r,gr in grdata:
995            grfile.write("%15.6g %15.6g\n" % (r,gr))
996        grfile.close()
997   
998def PeakListSave(G2frame,file,peaks):
999    print 'save peak list to file: ',G2frame.peaklistfile
1000    if not peaks:
1001        dlg = wx.MessageDialog(G2frame, 'No peaks!', 'Nothing to save!', wx.OK)
1002        try:
1003            result = dlg.ShowModal()
1004        finally:
1005            dlg.Destroy()
1006        return
1007    for peak in peaks:
1008        file.write("%10.4f %12.2f %10.3f %10.3f \n" % \
1009            (peak[0],peak[2],peak[4],peak[6]))
1010    print 'peak list saved'
1011             
1012def IndexPeakListSave(G2frame,peaks):
1013    file = open(G2frame.peaklistfile,'wa')
1014    print 'save index peak list to file: ',G2frame.peaklistfile
1015    wx.BeginBusyCursor()
1016    try:
1017        if not peaks:
1018            dlg = wx.MessageDialog(G2frame, 'No peaks!', 'Nothing to save!', wx.OK)
1019            try:
1020                result = dlg.ShowModal()
1021            finally:
1022                dlg.Destroy()
1023            return
1024        for peak in peaks:
1025            file.write("%12.6f\n" % (peak[7]))
1026        file.close()
1027    finally:
1028        wx.EndBusyCursor()
1029    print 'index peak list saved'
1030   
1031def SetNewPhase(Name='New Phase',SGData=G2spc.SpcGroup('P 1')[1],cell=[1.0,1.0,1.0,90.,90,90.,1.]):
1032    phaseData = {
1033        'General':{
1034            'Name':Name,
1035            'Type':'nuclear',
1036            'SGData':SGData,
1037            'Cell':[False,]+cell,
1038            'Pawley dmin':1.0,
1039            'Data plot type':'None',
1040            'SH Texture':{
1041                'Order':0,
1042                'Model':'cylindrical',
1043                'Sample omega':[False,0.0],
1044                'Sample chi':[False,0.0],
1045                'Sample phi':[False,0.0],
1046                'SH Coeff':[False,{}],
1047                'SHShow':False,
1048                'PFhkl':[0,0,1],
1049                'PFxyz':[0,0,1],
1050                'PlotType':'Pole figure'}},
1051        'Atoms':[],
1052        'Drawing':{},
1053        'Histograms':{},
1054        'Pawley ref':[],
1055        'Models':{},
1056        }
1057    return phaseData
1058   
1059def ReadEXPPhase(G2frame,filename):
1060    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1061    textureData = {'Order':0,'Model':'cylindrical','Sample omega':[False,0.0],
1062        'Sample chi':[False,0.0],'Sample phi':[False,0.0],'SH Coeff':[False,{}],
1063        'SHShow':False,'PFhkl':[0,0,1],'PFxyz':[0,0,1],'PlotType':'Pole figure'}
1064    shNcof = 0
1065    file = open(filename, 'Ur')
1066    S = 1
1067    Expr = [{},{},{},{},{},{},{},{},{}]
1068    while S:
1069        S = file.readline()
1070        if 'EXPR NPHAS' in S[:12]:
1071            Num = S[12:-1].count('0')
1072            NPhas = S[12:-1].split()
1073        if 'CRS' in S[:3]:
1074            N = int(S[3:4])-1
1075            Expr[N][S[:12]] = S[12:-1]
1076    file.close()
1077    PNames = []
1078    for n,N in enumerate(NPhas):
1079        if N != '0':
1080            result = n
1081            key = 'CRS'+str(n+1)+'    PNAM'
1082            PNames.append(Expr[n][key])
1083    if Num < 8:
1084        dlg = wx.SingleChoiceDialog(G2frame, 'Which phase to read?', 'Read phase data', PNames, wx.CHOICEDLG_STYLE)
1085        try:
1086            if dlg.ShowModal() == wx.ID_OK:
1087                result = dlg.GetSelection()
1088        finally:
1089            dlg.Destroy()       
1090    EXPphase = Expr[result]
1091    keyList = EXPphase.keys()
1092    keyList.sort()
1093    SGData = {}
1094    if NPhas[result] == '1':
1095        Ptype = 'nuclear'
1096    elif NPhas[result] in ['2','3']:
1097        Ptype = 'magnetic'
1098    elif NPhas[result] == '4':
1099        Ptype = 'macromolecular'
1100    elif NPhas[result] == '10':
1101        Ptype = 'Pawley'
1102    for key in keyList:
1103        if 'PNAM' in key:
1104           PhaseName = EXPphase[key].strip()
1105        elif 'ABC   ' in key:
1106            abc = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                       
1107        elif 'ANGLES' in key:
1108            angles = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                                               
1109        elif 'SG SYM' in key:
1110            SpGrp = EXPphase[key][:15].strip()
1111            E,SGData = G2spc.SpcGroup(SpGrp)
1112        elif 'OD    ' in key:
1113            SHdata = EXPphase[key].split() # may not have all 9 values
1114            SHvals = 9*[0]
1115            for i in range(9):
1116                try:
1117                    float(SHdata[i])
1118                    SHvals[i] = SHdata[i]
1119                except:
1120                    pass
1121            textureData['Order'] = int(SHvals[0])
1122            textureData['Model'] = shModels[int(SHvals[2])]
1123            textureData['Sample omega'] = [False,float(SHvals[6])]
1124            textureData['Sample chi'] = [False,float(SHvals[7])]
1125            textureData['Sample phi'] = [False,float(SHvals[8])]
1126            shNcof = int(SHvals[1])
1127    Atoms = []
1128    if Ptype == 'nuclear':
1129        for key in keyList:
1130            if 'AT' in key:
1131                if key[11:] == 'A':
1132                    S = EXPphase[key]
1133                elif key[11:] == 'B':
1134                    S += EXPphase[key]
1135                    Atom = [S[50:58].strip(),S[:10].strip(),'',
1136                        float(S[10:20]),float(S[20:30]),float(S[30:40]),
1137                        float(S[40:50]),'',int(S[60:62]),S[130:131]]
1138                    if Atom[9] == 'I':
1139                        Atom += [float(S[68:78]),0.,0.,0.,0.,0.,0.]
1140                    elif Atom[9] == 'A':
1141                        Atom += [0.0,float(S[68:78]),float(S[78:88]),
1142                            float(S[88:98]),float(S[98:108]),
1143                            float(S[108:118]),float(S[118:128])]
1144                    XYZ = Atom[3:6]
1145                    Atom[7],Atom[8] = G2spc.SytSym(XYZ,SGData)
1146                    Atom.append(ran.randint(0,sys.maxint))
1147                    Atoms.append(Atom)
1148    elif Ptype == 'macromolecular':
1149        for key in keyList:
1150            if 'AT' in key[6:8]:
1151                S = EXPphase[key]
1152                Atom = [S[56:60],S[50:54].strip().upper(),S[54:56],
1153                    S[46:51].strip(),S[:8].strip(),'',
1154                    float(S[16:24]),float(S[24:32]),float(S[32:40]),
1155                    float(S[8:16]),'1',1,'I',float(S[40:46]),0,0,0,0,0,0]
1156                XYZ = Atom[6:9]
1157                Atom[10],Atom[11] = G2spc.SytSym(XYZ,SGData)
1158                Atom.append(ran.randint(0,sys.maxint))
1159                Atoms.append(Atom)
1160    Volume = G2lat.calc_V(G2lat.cell2A(abc+angles))
1161    if shNcof:
1162        shCoef = {}
1163        nRec = [i+1 for i in range((shNcof-1)/6+1)]
1164        for irec in nRec:
1165            ODkey = keyList[0][:6]+'OD'+'%3dA'%(irec)
1166            indx = EXPphase[ODkey].split()
1167            ODkey = ODkey[:-1]+'B'
1168            vals = EXPphase[ODkey].split()
1169            for i,val in enumerate(vals):
1170                key = 'C(%s,%s,%s)'%(indx[3*i],indx[3*i+1],indx[3*i+2])
1171                shCoef[key] = float(val)
1172        textureData['SH Coeff'] = [False,shCoef]
1173       
1174    Phase = SetNewPhase(Name=PhaseName,SGData=SGData,cell=abc+angles+[Volume,])
1175    general = Phase['General']
1176    general['Type'] = Ptype
1177    general['SH Texture'] = textureData
1178    Phase['Atoms'] = Atoms
1179    return Phase
1180       
1181def ReadPDBPhase(filename):
1182    EightPiSq = 8.*math.pi**2
1183    file = open(filename, 'Ur')
1184    Phase = {}
1185    Title = ''
1186    Compnd = ''
1187    Atoms = []
1188    A = np.zeros(shape=(3,3))
1189    S = file.readline()
1190    while S:
1191        Atom = []
1192        if 'TITLE' in S[:5]:
1193            Title = S[10:72].strip()
1194            S = file.readline()
1195        elif 'COMPND    ' in S[:10]:
1196            Compnd = S[10:72].strip()
1197            S = file.readline()
1198        elif 'CRYST' in S[:5]:
1199            abc = S[7:34].split()
1200            angles = S[34:55].split()
1201            cell=[float(abc[0]),float(abc[1]),float(abc[2]),
1202                float(angles[0]),float(angles[1]),float(angles[2])]
1203            Volume = G2lat.calc_V(G2lat.cell2A(cell))
1204            AA,AB = G2lat.cell2AB(cell)
1205            SpGrp = S[55:65]
1206            E,SGData = G2spc.SpcGroup(SpGrp)
1207            while E:
1208                print G2spc.SGErrors(E)
1209                dlg = wx.TextEntryDialog(None,
1210                    SpGrp[:-1]+' is invalid \nN.B.: make sure spaces separate axial fields in symbol',
1211                    'ERROR in space group symbol','',style=wx.OK)
1212                if dlg.ShowModal() == wx.ID_OK:
1213                    SpGrp = dlg.GetValue()
1214                    E,SGData = G2spc.SpcGroup(SpGrp)
1215                else:
1216                    return None
1217                dlg.Destroy()               
1218            SGlines = G2spc.SGPrint(SGData)
1219            for line in SGlines: print line
1220            S = file.readline()
1221        elif 'SCALE' in S[:5]:
1222            V = (S[10:41].split())
1223            A[int(S[5])-1] = [float(V[0]),float(V[1]),float(V[2])]
1224            S = file.readline()
1225        elif 'ATOM' in S[:4] or 'HETATM' in S[:6]:
1226            XYZ = [float(S[31:39]),float(S[39:47]),float(S[47:55])]
1227            XYZ = np.inner(AB,XYZ)
1228            XYZ = np.where(abs(XYZ)<0.00001,0,XYZ)
1229            SytSym,Mult = G2spc.SytSym(XYZ,SGData)
1230            Uiso = float(S[61:67])/EightPiSq
1231            Type = S[12:14].upper()
1232            if Type[0] in '123456789':
1233                Type = Type[1:]
1234            Atom = [S[22:27].strip(),S[17:20].upper(),S[20:22],
1235                S[12:17].strip(),Type.strip(),'',XYZ[0],XYZ[1],XYZ[2],
1236                float(S[55:61]),SytSym,Mult,'I',Uiso,0,0,0,0,0,0]
1237            S = file.readline()
1238            if 'ANISOU' in S[:6]:
1239                Uij = S[30:72].split()
1240                Uij = [float(Uij[0])/10000.,float(Uij[1])/10000.,float(Uij[2])/10000.,
1241                    float(Uij[3])/10000.,float(Uij[4])/10000.,float(Uij[5])/10000.]
1242                Atom = Atom[:14]+Uij
1243                Atom[12] = 'A'
1244                S = file.readline()
1245            Atom.append(ran.randint(0,sys.maxint))
1246            Atoms.append(Atom)
1247        else:           
1248            S = file.readline()
1249    file.close()
1250    if Title:
1251        PhaseName = Title
1252    elif Compnd:
1253        PhaseName = Compnd
1254    else:
1255        PhaseName = 'None'
1256    Phase = SetNewPhase(Name=PhaseName,SGData=SGData,cell=cell+[Volume,])
1257    Phase['General']['Type'] = 'macromolecular'
1258    Phase['Atoms'] = Atoms
1259   
1260    return Phase
1261
1262######################################################################
1263# base classes for reading various types of data files
1264#   not used directly, only by subclassing
1265######################################################################
1266E,SGData = G2spc.SpcGroup('P 1') # data structure for default space group
1267class ImportBaseclass(object):
1268    '''Defines a base class for the importing of data files (diffraction
1269    data, coordinates,...
1270    '''
1271    def __init__(self,
1272                 formatName,
1273                 longFormatName=None,
1274                 extensionlist=[],
1275                 strictExtension=False,
1276                 ):
1277        self.formatName = formatName # short string naming file type
1278        if longFormatName: # longer string naming file type
1279            self.longFormatName = longFormatName
1280        else:
1281            self.longFormatName = formatName
1282        # define extensions that are allowed for the file type
1283        # for windows, remove any extensions that are duplicate, as case is ignored
1284        if sys.platform == 'windows' and extensionlist:
1285            extensionlist = list(set([s.lower() for s in extensionlist]))
1286        self.extensionlist = extensionlist
1287        # If strictExtension is True, the file will not be read, unless
1288        # the extension matches one in the extensionlist
1289        self.strictExtension = strictExtension
1290        self.warnings = ''
1291        self.errors = ''
1292        #print 'created',self.__class__
1293
1294    def BlockSelector(self, ChoiceList, ParentFrame=None,
1295                      title='Select a block',
1296                      size=None, header='Block Selector'):
1297        ''' Provide a wx dialog to select a block if the file contains more
1298        than one set of data and one must be selected
1299        '''
1300        dlg = wx.SingleChoiceDialog(
1301            ParentFrame,
1302            title, header,
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
1352class ImportPhase(ImportBaseclass):
1353    '''Defines a base class for the reading of files with coordinates
1354    '''
1355    def __init__(self,
1356                 formatName,
1357                 longFormatName=None,
1358                 extensionlist=[],
1359                 strictExtension=False,
1360                 ):
1361        # call parent __init__
1362        ImportBaseclass.__init__(self,formatName,
1363                                            longFormatName,
1364                                            extensionlist,
1365                                            strictExtension)
1366        # define a default Phase structure
1367        self.Phase = SetNewPhase(Name='new phase',SGData=SGData)
1368
1369    def PhaseSelector(self, ChoiceList, ParentFrame=None,
1370                      title='Select a phase', size=None,
1371                      header='Phase Selector'):
1372        ''' Provide a wx dialog to select a phase if the file contains more
1373        than one phase
1374        '''
1375        self.BlockSelector(ChoiceList, ParentFrame, title, size, header)
1376
1377######################################################################
1378class ImportStructFactor(ImportBaseclass):
1379    '''Defines a base class for the reading of files with tables
1380    of structure factors
1381    '''
1382    def __init__(self,
1383                 formatName,
1384                 longFormatName=None,
1385                 extensionlist=[],
1386                 strictExtension=False,
1387                 ):
1388        ImportBaseclass.__init__(self,formatName,
1389                                            longFormatName,
1390                                            extensionlist,
1391                                            strictExtension)
1392
1393        # define contents of Structure Factor entry
1394        self.Controls = { # dictionary with plotting controls
1395            'Type' : 'Fosq',
1396            'ifFc' : False,    #
1397            'HKLmax' : [None,None,None],
1398            'HKLmin' : [None,None,None],
1399            'FoMax' : None,   # maximum observed structure factor as Fo
1400            'Zone' : '001',
1401            'Layer' : 0,
1402            'Scale' : 1.0,
1403            'log-lin' : 'lin',
1404            }
1405        self.Parameters = [ # list with data collection parameters
1406            ('SXC',1.5428),
1407            ['SXC',1.5428],
1408            ['Type','Lam']
1409            ]
1410        self.RefList = []
1411
1412    def UpdateParameters(self,Type=None,Wave=None):
1413        HistType = self.Parameters[0][0]
1414        HistWave = self.Parameters[0][1]
1415        if Type is not None:
1416            HistType = Type
1417        if Wave is not None:
1418            HistWave = Wave
1419        self.Parameters = [ # overwrite entire list
1420            (HistType,HistWave),
1421            [HistType,HistWave],
1422            ['Type','Lam']
1423            ]
1424           
1425    def UpdateControls(self,Type='Fosq',FcalcPresent=False):
1426        '''Scan through the reflections to update the Controls dictionary
1427        '''
1428        self.Controls['Type'] = Type
1429        self.Controls['iffc'] = FcalcPresent
1430        HKLmax = [None,None,None]
1431        HKLmin = [None,None,None]
1432        Fo2max = None
1433        for HKL,Fo2,SFo2,Fc,Fcp,Fcpp,phase in self.RefList:
1434            if Fo2max is None:
1435                Fo2max = Fo2
1436            else:
1437                Fo2max = max(Fo2max,Fo2)
1438            for i,hkl in enumerate(HKL):
1439                if HKLmax[i] is None:
1440                    HKLmax[i] = hkl
1441                    HKLmin[i] = hkl
1442                else:
1443                    HKLmax[i] = max(HKLmax[i],hkl)
1444                    HKLmin[i] = min(HKLmin[i],hkl)
1445        self.Controls['HKLmax'] = HKLmax
1446        self.Controls['HKLmin'] = HKLmin
1447        if Type ==  'Fosq':
1448            self.Controls['FoMax'] = np.sqrt(Fo2max)
1449        elif Type ==  'Fo':
1450            self.Controls['FoMax'] = Fo2max
1451        else:
1452            print "Unsupported Stract Fact type in ImportStructFactor.UpdateControls"
1453            raise Exception,"Unsupported Stract Fact type in ImportStructFactor.UpdateControls"
Note: See TracBrowser for help on using the repository browser.