source: trunk/GSASIIIO.py @ 473

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

Add FileDlgFixExt? to fix linux problem in FileDialog? where the extension is not taken from the wildcard choice.

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