source: trunk/GSASIIIO.py @ 607

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

revise help to include tutorials on tree only; add kw args to all readers; cache cif for multiple passes; start on powder imports

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