source: trunk/GSASIIIO.py @ 457

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

begin distance angle calcs
move gpxfile routines from GSASIIstruct.py to GSASIIIO.py
move getVCov & ValEsd? to GSASIImath.py
add some text to help/gsasII.html

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