source: trunk/GSASIIIO.py @ 630

Last change on this file since 630 was 630, checked in by vondreele, 11 years ago

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