source: trunk/GSASIIIO.py @ 464

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

revamp directory behavior

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