source: trunk/GSASIIIO.py @ 467

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

more on file directories
fix for image controls & mask save
modify image controls data
  • Property svn:keywords set to Date Author Revision URL Id
File size: 64.5 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-01 22:03:51 +0000 (Wed, 01 Feb 2012) $
6# $Author: vondreele $
7# $Revision: 467 $
8# $URL: trunk/GSASIIIO.py $
9# $Id: GSASIIIO.py 467 2012-02-01 22:03:51Z 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(*.*)|*.*', 
147        wx.OPEN|wx.CHANGE_DIR)
148    Tname = filename[:filename.index('.')]+'.prm'
149    if op.exists(Tname):
150        self.IparmName = Tname       
151    if self.IparmName: dlg.SetFilename(self.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(self,imagefile):
475    if not ospath.exists(imagefile):
476        dlg = wx.FileDialog(self, '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(self,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(self):
798    file = open(self.GSASprojectfile,'rb')
799    print 'load from file: ',self.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 = self.PatternTree.AppendItem(parent=self.root,text=datum[0])
810            if 'PWDR' in datum[0]:               
811                self.PatternTree.SetItemPyData(Id,datum[1][:3])     #temp. trim off junk
812            else:
813                self.PatternTree.SetItemPyData(Id,datum[1])
814            for datus in data[1:]:
815                sub = self.PatternTree.AppendItem(Id,datus[0])
816                self.PatternTree.SetItemPyData(sub,datus[1])
817            if 'IMG' in datum[0]:                   #retreive image default flag & data if set
818                Data = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Image Controls'))
819                if Data['setDefault']:
820                    self.imageDefault = Data               
821        file.close()
822       
823    finally:
824        wx.EndBusyCursor()
825    print 'project load successful'
826    self.NewPlot = True
827   
828def ProjFileSave(self):
829    if not self.PatternTree.IsEmpty():
830        file = open(self.GSASprojectfile,'wb')
831        print 'save to file: ',self.GSASprojectfile
832        wx.BeginBusyCursor()
833        try:
834            item, cookie = self.PatternTree.GetFirstChild(self.root)
835            while item:
836                data = []
837                name = self.PatternTree.GetItemText(item)
838                data.append([name,self.PatternTree.GetItemPyData(item)])
839                item2, cookie2 = self.PatternTree.GetFirstChild(item)
840                while item2:
841                    name = self.PatternTree.GetItemText(item2)
842                    data.append([name,self.PatternTree.GetItemPyData(item2)])
843                    item2, cookie2 = self.PatternTree.GetNextChild(item, cookie2)                           
844                item, cookie = self.PatternTree.GetNextChild(self.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[1]:
1050            if makeBack:
1051                last = max(last,int(name[1].strip('bak'))+1)
1052            else:
1053                last = max(last,int(name[1].strip('bak')))
1054    GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.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   
1134                   
1135def GetPWDRdata(GPXfile,PWDRname):
1136    ''' Returns powder data from GSASII gpx file
1137    input:
1138        GPXfile = .gpx full file name
1139        PWDRname = powder histogram name as obtained from GetHistogramNames
1140    return:
1141        PWDRdata = powder data dictionary with:
1142            Data - powder data arrays, Limits, Instrument Parameters, Sample Parameters
1143       
1144    '''
1145    file = open(GPXfile,'rb')
1146    PWDRdata = {}
1147    while True:
1148        try:
1149            data = cPickle.load(file)
1150        except EOFError:
1151            break
1152        datum = data[0]
1153        if datum[0] == PWDRname:
1154            PWDRdata['Data'] = datum[1][1]          #powder data arrays
1155            PWDRdata[data[2][0]] = data[2][1]       #Limits
1156            PWDRdata[data[3][0]] = data[3][1]       #Background
1157            PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
1158            PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
1159            try:
1160                PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
1161            except IndexError:
1162                PWDRdata['Reflection lists'] = {}
1163    file.close()
1164    return PWDRdata
1165   
1166def GetHKLFdata(GPXfile,HKLFname):
1167    ''' Returns single crystal data from GSASII gpx file
1168    input:
1169        GPXfile = .gpx full file name
1170        HKLFname = single crystal histogram name as obtained from GetHistogramNames
1171    return:
1172        HKLFdata = single crystal data list of reflections: for each reflection:
1173            HKLF = [np.array([h,k,l]),FoSq,sigFoSq,FcSq,Fcp,Fcpp,phase]
1174    '''
1175    file = open(GPXfile,'rb')
1176    HKLFdata = []
1177    while True:
1178        try:
1179            data = cPickle.load(file)
1180        except EOFError:
1181            break
1182        datum = data[0]
1183        if datum[0] == HKLFname:
1184            HKLFdata = datum[1:][0]
1185    file.close()
1186    return HKLFdata
1187   
1188def SaveIntegration(self,PickId,data):
1189    azms = self.Integrate[1]
1190    X = self.Integrate[2][:-1]
1191    Xminmax = [X[0],X[-1]]
1192    N = len(X)
1193    Id = self.PatternTree.GetItemParent(PickId)
1194    name = self.PatternTree.GetItemText(Id)
1195    name = name.replace('IMG ','PWDR ')
1196    Comments = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Comments'))
1197    names = ['Type','Lam','Zero','Polariz.','U','V','W','X','Y','SH/L','Azimuth'] 
1198    codes = [0 for i in range(11)]
1199    LRazm = data['LRazimuth']
1200    Azms = []
1201    if data['fullIntegrate'] and data['outAzimuths'] == 1:
1202        Azms = [45.0,]                              #a poor man's average?
1203    else:
1204        for i,azm in enumerate(azms[:-1]):
1205            Azms.append((azms[i+1]+azm)/2.)
1206    for i,azm in enumerate(azms[:-1]):
1207        item, cookie = self.PatternTree.GetFirstChild(self.root)
1208        Id = 0
1209        while item:
1210            Name = self.PatternTree.GetItemText(item)
1211            if name == Name:
1212                Id = item
1213            item, cookie = self.PatternTree.GetNextChild(self.root, cookie)
1214        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!
1215        Y = self.Integrate[0][i]
1216        W = 1./Y                    #probably not true
1217        Sample = G2pdG.SetDefaultSample()
1218        if Id:
1219            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id, 'Comments'),Comments)                   
1220            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Limits'),[tuple(Xminmax),Xminmax])
1221            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Background'),[['chebyschev',1,3,1.0,0.0,0.0],
1222                            {'nDebye':0,'debyeTerms':[],'nPeaks':0,'peaksList':[]}])
1223            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Instrument Parameters'),[tuple(parms),parms[:],codes,names])
1224            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Peak List'),[])
1225            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Index Peak List'),[])
1226            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Unit Cells List'),[])             
1227            self.PatternTree.SetItemPyData(G2gd.GetPatternTreeItemId(self,Id,'Reflection Lists'),{})             
1228        else:
1229            Id = self.PatternTree.AppendItem(parent=self.root,text=name+" Azm= %.2f"%(Azms[i]))
1230            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Comments'),Comments)                   
1231            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Limits'),[tuple(Xminmax),Xminmax])
1232            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Background'),[['chebyschev',1,3,1.0,0.0,0.0],
1233                            {'nDebye':0,'debyeTerms':[],'nPeaks':0,'peaksList':[]}])
1234            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Instrument Parameters'),[tuple(parms),parms[:],codes,names])
1235            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Sample Parameters'),Sample)
1236            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Peak List'),[])
1237            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Index Peak List'),[])
1238            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Unit Cells List'),[])
1239            self.PatternTree.SetItemPyData(self.PatternTree.AppendItem(Id,text='Reflection Lists'),{})             
1240        self.PatternTree.SetItemPyData(Id,[[''],[np.array(X),np.array(Y),np.array(W),np.zeros(N),np.zeros(N),np.zeros(N)]])
1241    self.PatternTree.SelectItem(Id)
1242    self.PatternTree.Expand(Id)
1243    self.PatternId = Id
1244           
1245def powderFxyeSave(self,exports,powderfile):
1246    head,tail = ospath.split(powderfile)
1247    name,ext = tail.split('.')
1248    wx.BeginBusyCursor()
1249    for i,export in enumerate(exports):
1250        filename = ospath.join(head,name+'-%03d.'%(i)+ext)
1251        prmname = filename.strip(ext)+'prm'
1252        prm = open(prmname,'w')      #old style GSAS parm file
1253        PickId = G2gd.GetPatternTreeItemId(self, self.root, export)
1254        Values,Names = self.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(self, \
1255            PickId, 'Instrument Parameters'))[1::2]     #get values & names
1256        Inst = dict(zip(Names,Values))
1257        print Inst['Type']
1258        prm.write( '            123456789012345678901234567890123456789012345678901234567890        '+'\n')
1259        prm.write( 'INS   BANK      1                                                               '+'\n')
1260        prm.write(('INS   HTYPE   %sR                                                              '+'\n')%(Inst['Type']))
1261        if 'Lam1' in Inst:              #Ka1 & Ka2
1262            prm.write(('INS  1 ICONS%10.7f%10.7f    0.0000               0.990    0     0.500   '+'\n')%(Inst['Lam1'],Inst['Lam2']))
1263        elif 'Lam' in Inst:             #single wavelength
1264            prm.write(('INS  1 ICONS%10.7f%10.7f    0.0000               0.990    0     0.500   '+'\n')%(Inst['Lam'],0.0))
1265        prm.write( 'INS  1 IRAD     0                                                               '+'\n')
1266        prm.write( 'INS  1I HEAD                                                                    '+'\n')
1267        prm.write( 'INS  1I ITYP    0    0.0000  180.0000         1                                 '+'\n')
1268        prm.write(('INS  1DETAZM%10.3f                                                          '+'\n')%(Inst['Azimuth']))
1269        prm.write( 'INS  1PRCF1     3    8   0.00100                                                '+'\n')
1270        prm.write(('INS  1PRCF11     %15.6g%15.6g%15.6g%15.6g   '+'\n')%(Inst['U'],Inst['V'],Inst['W'],0.0))
1271        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.))
1272        prm.close()
1273        file = open(filename,'w')
1274        print 'save powder pattern to file: ',filename
1275        try:
1276            x,y,w,yc,yb,yd = self.PatternTree.GetItemPyData(PickId)[1]
1277            file.write(powderfile+'\n')
1278            file.write('BANK 1 %d %d CONS %.2f %.2f 0 0 FXYE\n'%(len(x),len(x),\
1279                100.*x[0],100.*(x[1]-x[0])))
1280            s = list(np.sqrt(1./np.array(w)))       
1281            XYW = zip(x,y,s)
1282            for X,Y,S in XYW:
1283                file.write("%15.6g %15.6g %15.6g\n" % (100.*X,Y,max(S,1.0)))
1284            file.close()
1285        finally:
1286            wx.EndBusyCursor()
1287        print 'powder pattern file written'
1288       
1289def powderXyeSave(self,exports,powderfile):
1290    head,tail = ospath.split(powderfile)
1291    name,ext = tail.split('.')
1292    for i,export in enumerate(exports):
1293        filename = ospath.join(head,name+'-%03d.'%(i)+ext)
1294        PickId = G2gd.GetPatternTreeItemId(self, self.root, export)
1295        file = open(filename,'w')
1296        file.write('#%s\n'%(export))
1297        print 'save powder pattern to file: ',filename
1298        wx.BeginBusyCursor()
1299        try:
1300            x,y,w,yc,yb,yd = self.PatternTree.GetItemPyData(PickId)[1]
1301            s = list(np.sqrt(1./np.array(w)))       
1302            XYW = zip(x,y,s)
1303            for X,Y,W in XYW:
1304                file.write("%15.6g %15.6g %15.6g\n" % (X,Y,W))
1305            file.close()
1306        finally:
1307            wx.EndBusyCursor()
1308        print 'powder pattern file written'
1309       
1310def PDFSave(self,exports):   
1311    for export in exports:
1312        PickId = G2gd.GetPatternTreeItemId(self, self.root, export)
1313        SQname = 'S(Q)'+export[4:]
1314        GRname = 'G(R)'+export[4:]
1315        sqfilename = ospath.join(self.dirname,export.replace(' ','_')[5:]+'.sq')
1316        grfilename = ospath.join(self.dirname,export.replace(' ','_')[5:]+'.gr')
1317        sqId = G2gd.GetPatternTreeItemId(self, PickId, SQname)
1318        grId = G2gd.GetPatternTreeItemId(self, PickId, GRname)
1319        sqdata = np.array(self.PatternTree.GetItemPyData(sqId)[1][:2]).T
1320        grdata = np.array(self.PatternTree.GetItemPyData(grId)[1][:2]).T
1321        sqfile = open(sqfilename,'w')
1322        grfile = open(grfilename,'w')
1323        sqfile.write('#T S(Q) %s\n'%(export))
1324        grfile.write('#T G(R) %s\n'%(export))
1325        sqfile.write('#L Q     S(Q)\n')
1326        grfile.write('#L R     G(R)\n')
1327        for q,sq in sqdata:
1328            sqfile.write("%15.6g %15.6g\n" % (q,sq))
1329        sqfile.close()
1330        for r,gr in grdata:
1331            grfile.write("%15.6g %15.6g\n" % (r,gr))
1332        grfile.close()
1333   
1334def PeakListSave(self,file,peaks):
1335    print 'save peak list to file: ',self.peaklistfile
1336    if not peaks:
1337        dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
1338        try:
1339            result = dlg.ShowModal()
1340        finally:
1341            dlg.Destroy()
1342        return
1343    for peak in peaks:
1344        file.write("%10.4f %12.2f %10.3f %10.3f \n" % \
1345            (peak[0],peak[2],peak[4],peak[6]))
1346    print 'peak list saved'
1347             
1348def IndexPeakListSave(self,peaks):
1349    file = open(self.peaklistfile,'wa')
1350    print 'save index peak list to file: ',self.peaklistfile
1351    wx.BeginBusyCursor()
1352    try:
1353        if not peaks:
1354            dlg = wx.MessageDialog(self, 'No peaks!', 'Nothing to save!', wx.OK)
1355            try:
1356                result = dlg.ShowModal()
1357            finally:
1358                dlg.Destroy()
1359            return
1360        for peak in peaks:
1361            file.write("%12.6f\n" % (peak[7]))
1362        file.close()
1363    finally:
1364        wx.EndBusyCursor()
1365    print 'index peak list saved'
1366   
1367def ReadEXPPhase(self,filename):
1368    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1369    textureData = {'Order':0,'Model':'cylindrical','Sample omega':[False,0.0],
1370        'Sample chi':[False,0.0],'Sample phi':[False,0.0],'SH Coeff':[False,{}],
1371        'SHShow':False,'PFhkl':[0,0,1],'PFxyz':[0,0,1],'PlotType':'Pole figure'}
1372    shNcof = 0
1373    file = open(filename, 'Ur')
1374    S = 1
1375    Expr = [{},{},{},{},{},{},{},{},{}]
1376    while S:
1377        S = file.readline()
1378        if 'EXPR NPHAS' in S[:12]:
1379            Num = S[12:-1].count('0')
1380            NPhas = S[12:-1].split()
1381        if 'CRS' in S[:3]:
1382            N = int(S[3:4])-1
1383            Expr[N][S[:12]] = S[12:-1]
1384    file.close()
1385    PNames = []
1386    for n,N in enumerate(NPhas):
1387        if N != '0':
1388            result = n
1389            key = 'CRS'+str(n+1)+'    PNAM'
1390            PNames.append(Expr[n][key])
1391    if Num < 8:
1392        dlg = wx.SingleChoiceDialog(self, 'Which phase to read?', 'Read phase data', PNames, wx.CHOICEDLG_STYLE)
1393        try:
1394            if dlg.ShowModal() == wx.ID_OK:
1395                result = dlg.GetSelection()
1396        finally:
1397            dlg.Destroy()       
1398    EXPphase = Expr[result]
1399    keyList = EXPphase.keys()
1400    keyList.sort()
1401    SGData = {}
1402    if NPhas[result] == '1':
1403        Ptype = 'nuclear'
1404    elif NPhas[result] in ['2','3']:
1405        Ptype = 'magnetic'
1406    elif NPhas[result] == '4':
1407        Ptype = 'macromolecular'
1408    elif NPhas[result] == '10':
1409        Ptype = 'Pawley'
1410    for key in keyList:
1411        if 'PNAM' in key:
1412           PhaseName = EXPphase[key].strip()
1413        elif 'ABC   ' in key:
1414            abc = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                       
1415        elif 'ANGLES' in key:
1416            angles = [float(EXPphase[key][:10]),float(EXPphase[key][10:20]),float(EXPphase[key][20:30])]                                               
1417        elif 'SG SYM' in key:
1418            SpGrp = EXPphase[key][:15].strip()
1419            E,SGData = G2spc.SpcGroup(SpGrp)
1420        elif 'OD    ' in key:
1421            SHdata = EXPphase[key].split()
1422            textureData['Order'] = int(SHdata[0])
1423            textureData['Model'] = shModels[int(SHdata[2])]
1424            textureData['Sample omega'] = [False,float(SHdata[6])]
1425            textureData['Sample chi'] = [False,float(SHdata[7])]
1426            textureData['Sample phi'] = [False,float(SHdata[8])]
1427            shNcof = int(SHdata[1])
1428    Atoms = []
1429    if Ptype == 'nuclear':
1430        for key in keyList:
1431            if 'AT' in key:
1432                if key[11:] == 'A':
1433                    S = EXPphase[key]
1434                elif key[11:] == 'B':
1435                    S += EXPphase[key]
1436                    Atom = [S[50:58].strip(),S[:10].strip(),'',
1437                        float(S[10:20]),float(S[20:30]),float(S[30:40]),
1438                        float(S[40:50]),'',int(S[60:62]),S[130:131]]
1439                    if Atom[9] == 'I':
1440                        Atom += [float(S[68:78]),0.,0.,0.,0.,0.,0.]
1441                    elif Atom[9] == 'A':
1442                        Atom += [0.0,float(S[68:78]),float(S[78:88]),
1443                            float(S[88:98]),float(S[98:108]),
1444                            float(S[108:118]),float(S[118:128])]
1445                    XYZ = Atom[3:6]
1446                    Atom[7],Atom[8] = G2spc.SytSym(XYZ,SGData)
1447                    Atom.append(ran.randint(0,sys.maxint))
1448                    Atoms.append(Atom)
1449    elif Ptype == 'macromolecular':
1450        for key in keyList:
1451            if 'AT' in key[6:8]:
1452                S = EXPphase[key]
1453                Atom = [int(S[56:60]),S[50:54].strip().upper(),S[54:56],
1454                    S[46:51].strip(),S[:8].strip(),'',
1455                    float(S[16:24]),float(S[24:32]),float(S[32:40]),
1456                    float(S[8:16]),'1',1,'I',float(S[40:46]),0,0,0,0,0,0]
1457                XYZ = Atom[6:9]
1458                Atom[10],Atom[11] = G2spc.SytSym(XYZ,SGData)
1459                Atom.append(ran.randint(0,sys.maxint))
1460                Atoms.append(Atom)
1461    Volume = G2lat.calc_V(G2lat.cell2A(abc+angles))
1462    if shNcof:
1463        shCoef = {}
1464        nRec = [i+1 for i in range((shNcof-1)/6+1)]
1465        for irec in nRec:
1466            ODkey = keyList[0][:6]+'OD'+'%3dA'%(irec)
1467            indx = EXPphase[ODkey].split()
1468            ODkey = ODkey[:-1]+'B'
1469            vals = EXPphase[ODkey].split()
1470            for i,val in enumerate(vals):
1471                key = 'C(%s,%s,%s)'%(indx[3*i],indx[3*i+1],indx[3*i+2])
1472                shCoef[key] = float(val)
1473        textureData['SH Coeff'] = [False,shCoef]
1474       
1475    Phase = {
1476            'General':{
1477                'Name':PhaseName,
1478                'Type':Ptype,
1479                'SGData':SGData,
1480                'Cell':[False,]+abc+angles+[Volume,],
1481                'Pawley dmin':1.0},
1482            'Atoms':Atoms,
1483            'Drawing':{},
1484            'Histograms':{},
1485            'Pawley ref':[],
1486            'Models':{},
1487            'SH Texture':textureData
1488            }
1489           
1490    return Phase
1491       
1492def ReadPDBPhase(filename):
1493    EightPiSq = 8.*math.pi**2
1494    file = open(filename, 'Ur')
1495    Phase = {}
1496    Title = ''
1497    Compnd = ''
1498    Atoms = []
1499    A = np.zeros(shape=(3,3))
1500    S = file.readline()
1501    while S:
1502        Atom = []
1503        if 'TITLE' in S[:5]:
1504            Title = S[10:72].strip()
1505            S = file.readline()
1506        elif 'COMPND    ' in S[:10]:
1507            Compnd = S[10:72].strip()
1508            S = file.readline()
1509        elif 'CRYST' in S[:5]:
1510            abc = S[7:34].split()
1511            angles = S[34:55].split()
1512            cell=[float(abc[0]),float(abc[1]),float(abc[2]),
1513                float(angles[0]),float(angles[1]),float(angles[2])]
1514            Volume = G2lat.calc_V(G2lat.cell2A(cell))
1515            AA,AB = G2lat.cell2AB(cell)
1516            SpGrp = S[55:65]
1517            E,SGData = G2spc.SpcGroup(SpGrp)
1518            if E: 
1519                print ' ERROR in space group symbol ',SpGrp,' in file ',filename
1520                print ' N.B.: make sure spaces separate axial fields in symbol' 
1521                print G2spc.SGErrors(E)
1522                return None
1523            SGlines = G2spc.SGPrint(SGData)
1524            for line in SGlines: print line
1525            S = file.readline()
1526        elif 'SCALE' in S[:5]:
1527            V = (S[10:41].split())
1528            A[int(S[5])-1] = [float(V[0]),float(V[1]),float(V[2])]
1529            S = file.readline()
1530        elif 'ATOM' in S[:4] or 'HETATM' in S[:6]:
1531            XYZ = [float(S[31:39]),float(S[39:47]),float(S[47:55])]
1532            XYZ = np.inner(AB,XYZ)
1533            SytSym,Mult = G2spc.SytSym(XYZ,SGData)
1534            Uiso = float(S[61:67])/EightPiSq
1535            Type = S[12:14].upper()
1536            if Type[0] in '123456789':
1537                Type = Type[1:]
1538            Atom = [S[22:27].strip(),S[17:20].upper(),S[20:22],
1539                S[12:17].strip(),Type.strip(),'',XYZ[0],XYZ[1],XYZ[2],
1540                float(S[55:61]),SytSym,Mult,'I',Uiso,0,0,0,0,0,0]
1541            S = file.readline()
1542            if 'ANISOU' in S[:6]:
1543                Uij = S[30:72].split()
1544                Uij = [float(Uij[0])/10000.,float(Uij[1])/10000.,float(Uij[2])/10000.,
1545                    float(Uij[3])/10000.,float(Uij[4])/10000.,float(Uij[5])/10000.]
1546                Atom = Atom[:14]+Uij
1547                Atom[12] = 'A'
1548                S = file.readline()
1549            Atom.append(ran.randint(0,sys.maxint))
1550            Atoms.append(Atom)
1551        else:           
1552            S = file.readline()
1553    file.close()
1554    if Title:
1555        PhaseName = Title
1556    elif Compnd:
1557        PhaseName = Compnd
1558    else:
1559        PhaseName = 'None'
1560    Phase['General'] = {'Name':PhaseName,'Type':'macromolecular','SGData':SGData,
1561        'Cell':[False,]+cell+[Volume,]}
1562    Phase['Atoms'] = Atoms
1563    Phase['Drawing'] = {}
1564    Phase['Histograms'] = {}
1565   
1566    return Phase
1567   
1568def ReadCIFPhase(filename):
1569    anisoNames = ['aniso_u_11','aniso_u_22','aniso_u_33','aniso_u_12','aniso_u_13','aniso_u_23']
1570    file = open(filename, 'Ur')
1571    Phase = {}
1572    Title = ospath.split(filename)[-1]
1573    print '\n Reading cif file: ',Title
1574    Compnd = ''
1575    Atoms = []
1576    A = np.zeros(shape=(3,3))
1577    S = file.readline()
1578    while S:
1579        if '_symmetry_space_group_name_H-M' in S:
1580            SpGrp = S.split("_symmetry_space_group_name_H-M")[1].strip().strip('"').strip("'")
1581            E,SGData = G2spc.SpcGroup(SpGrp)
1582            if E:
1583                print ' ERROR in space group symbol ',SpGrp,' in file ',filename
1584                print ' N.B.: make sure spaces separate axial fields in symbol' 
1585                print G2spc.SGErrors(E)
1586                return None
1587            S = file.readline()
1588        elif '_cell' in S:
1589            if '_cell_length_a' in S:
1590                a = S.split('_cell_length_a')[1].strip().strip('"').strip("'").split('(')[0]
1591            elif '_cell_length_b' in S:
1592                b = S.split('_cell_length_b')[1].strip().strip('"').strip("'").split('(')[0]
1593            elif '_cell_length_c' in S:
1594                c = S.split('_cell_length_c')[1].strip().strip('"').strip("'").split('(')[0]
1595            elif '_cell_angle_alpha' in S:
1596                alp = S.split('_cell_angle_alpha')[1].strip().strip('"').strip("'").split('(')[0]
1597            elif '_cell_angle_beta' in S:
1598                bet = S.split('_cell_angle_beta')[1].strip().strip('"').strip("'").split('(')[0]
1599            elif '_cell_angle_gamma' in S:
1600                gam = S.split('_cell_angle_gamma')[1].strip().strip('"').strip("'").split('(')[0]
1601            S = file.readline()
1602        elif 'loop_' in S:
1603            labels = {}
1604            i = 0
1605            while S:
1606                S = file.readline()
1607                if '_atom_site' in S.strip()[:10]:
1608                    labels[S.strip().split('_atom_site_')[1].lower()] = i
1609                    i += 1
1610                else:
1611                    break
1612            if labels:
1613                if 'aniso_label' not in labels:
1614                    while S:
1615                        atom = ['','','',0,0,0,1.0,'','','I',0.01,0,0,0,0,0,0]
1616                        S.strip()
1617                        if len(S.split()) != len(labels):
1618                            if 'loop_' in S:
1619                                break
1620                            S += file.readline().strip()
1621                        data = S.split()
1622                        if len(data) != len(labels):
1623                            break
1624                        for key in labels:
1625                            if key == 'type_symbol':
1626                                atom[1] = data[labels[key]]
1627                            elif key == 'label':
1628                                atom[0] = data[labels[key]]
1629                            elif key == 'fract_x':
1630                                atom[3] = float(data[labels[key]].split('(')[0])
1631                            elif key == 'fract_y':
1632                                atom[4] = float(data[labels[key]].split('(')[0])
1633                            elif key == 'fract_z':
1634                                atom[5] = float(data[labels[key]].split('(')[0])
1635                            elif key == 'occupancy':
1636                                atom[6] = float(data[labels[key]].split('(')[0])
1637                            elif key == 'thermal_displace_type':
1638                                if data[labels[key]].lower() == 'uiso':
1639                                    atom[9] = 'I'
1640                                    atom[10] = float(data[labels['u_iso_or_equiv']].split('(')[0])
1641                                else:
1642                                    atom[9] = 'A'
1643                                    atom[10] = 0.0
1644                                   
1645                        atom[7],atom[8] = G2spc.SytSym(atom[3:6],SGData)
1646                        atom.append(ran.randint(0,sys.maxint))
1647                        Atoms.append(atom)
1648                        S = file.readline()
1649                else:
1650                    while S:
1651                        S.strip()
1652                        data = S.split()
1653                        if len(data) != len(labels):
1654                            break
1655                        name = data[labels['aniso_label']]
1656                        for atom in Atoms:
1657                            if name == atom[0]:
1658                                for i,uname in enumerate(anisoNames):
1659                                    atom[i+11] = float(data[labels[uname]].split('(')[0])
1660                        S = file.readline()
1661                                                                       
1662        else:           
1663            S = file.readline()
1664    file.close()
1665    if Title:
1666        PhaseName = Title
1667    else:
1668        PhaseName = 'None'
1669    SGlines = G2spc.SGPrint(SGData)
1670    for line in SGlines: print line
1671    cell = [float(a),float(b),float(c),float(alp),float(bet),float(gam)]
1672    Volume = G2lat.calc_V(G2lat.cell2A(cell))
1673    Phase['General'] = {'Name':PhaseName,'Type':'nuclear','SGData':SGData,
1674        'Cell':[False,]+cell+[Volume,]}
1675    Phase['Atoms'] = Atoms
1676    Phase['Drawing'] = {}
1677    Phase['Histograms'] = {}
1678   
1679    return Phase
1680
Note: See TracBrowser for help on using the repository browser.