source: trunk/imports/G2pwd_fxye.py @ 3146

Last change on this file since 3146 was 3146, checked in by vondreele, 5 years ago

add time to G2img_MAR.py, fix open bug in G2pwd_fxye.py
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 16.5 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2017-11-01 14:05:28 +0000 (Wed, 01 Nov 2017) $
4# $Author: vondreele $
5# $Revision: 3146 $
6# $URL: trunk/imports/G2pwd_fxye.py $
7# $Id: G2pwd_fxye.py 3146 2017-11-01 14:05:28Z vondreele $
8########### SVN repository information ###################
9'''
10*Module G2pwd_fxye: GSAS data files*
11------------------------------------
12Routine to read in powder data in a variety of formats
13that are defined for GSAS.
14
15'''
16from __future__ import division, print_function
17import os.path as ospath
18import numpy as np
19import GSASIIobj as G2obj
20import GSASIIpath
21GSASIIpath.SetVersionNumber("$Revision: 3146 $")
22
23class GSAS_ReaderClass(G2obj.ImportPowderData):
24    'Routines to import powder data from a GSAS files'
25    def __init__(self):
26        super(self.__class__,self).__init__( # fancy way to self-reference
27            extensionlist=('.fxye','.raw','.gsas','.gda','.gsa','.gss','.RAW','.GSAS','.GDA','.GSA'),
28            strictExtension=False,
29            formatName = 'GSAS powder data',
30            longFormatName = 'GSAS powder data files (.fxye, .raw, .gsas...)'
31            )
32        self.clockWd = {}
33        self.TimeMap = {}
34        self.dnames = []
35        self.scriptable = True
36
37    # Validate the contents -- look for a bank line
38    def ContentsValidator(self, filename):
39        'Validrate by checking to see if the file has BANK lines & count them'
40        #print 'ContentsValidator: '+self.formatName
41        nBanks= 0
42        fp = open(filename,'r')
43        fname = ospath.basename(fp.name)
44        for i,line in enumerate(fp):
45            self.GSAS = True
46            if i==0: # first line is always a comment
47                continue
48            if i==1 and line[:4].lower() == 'inst':
49                # 2nd line is optional instrument parameter file
50                continue
51            if line[0] == '#': continue
52            if line[:4] == 'BANK':
53                self.dnames.append(fname+' '+' '.join(line.split()[:2]))
54                nBanks += 1
55                continue
56            elif line[:7] == 'Monitor': continue
57            elif line [:8] == 'TIME_MAP':          #LANSCE TOF data
58                continue
59#            else:
60#                if not all(c in ' 0123456789.eE+-' for c in line[:-1]):
61#                    self.errors = 'Unexpected information in line: '+str(i+1)
62#                    self.errors += '  '+str(line)
63#                    return False
64        if nBanks:
65            if not len(self.selections):
66                self.selections = range(nBanks)
67            fp.close()
68            return True
69        self.errors = 'No BANK records found'
70        fp.close()
71        return False # no bank records
72
73    def Reader(self,filename, ParentFrame=None, **kwarg):
74        '''Read a GSAS (old formats) file of type FXY, FXYE, ESD or STD types.
75        If multiple datasets are requested, use self.repeat and buffer caching.
76        '''
77        def GetFXYEdata(File,Pos,Bank):
78            File.seek(Pos)
79            x = []
80            y = []
81            w = []
82            S = File.readline()
83            while S and S[:4] != 'BANK' and S[0] != '#':
84                vals = S.split()
85                x.append(float(vals[0])/100.)               #CW: from centidegrees to degrees
86                f = float(vals[1])
87                s = float(vals[2])
88                if f <= 0.0 or s <= 0.0:
89                    y.append(0.0)
90                    w.append(0.0)
91                else:
92                    y.append(float(vals[1]))
93                    w.append(1.0/float(vals[2])**2)
94                S = File.readline()
95            N = len(x)
96            return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]   
97           
98        def GetFXYdata(File,Pos,Bank):
99            File.seek(Pos)
100            x = []
101            y = []
102            w = []
103            S = File.readline()
104            while S and S[:4] != 'BANK' and S[0] != '#':
105                vals = S.split()
106                x.append(float(vals[0])/100.)               #CW: from centidegrees to degrees
107                f = float(vals[1])
108                if f > 0.0:
109                    y.append(f)
110                    w.append(1.0/f)
111                else:             
112                    y.append(0.0)
113                    w.append(0.0)
114                S = File.readline()
115            N = len(x)
116            return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
117           
118        def GetESDdata(File,Pos,Bank):
119            File.seek(Pos)
120            cons = Bank.split()
121            if 'TIME_MAP' == cons[4]:
122                start = 0
123                step = 1
124            else:
125                start = float(cons[5])/100.0               #CW: from centidegrees to degrees
126                step = float(cons[6])/100.0
127            x = []
128            y = []
129            w = []
130            S = File.readline()
131            j = 0
132            while S and S[:4] != 'BANK' and S[0] != '#':
133                if 'TIME_MAP' in S:
134                    break
135                for i in range(0,80,16):
136                    if S[i:i+8] == 8*' ':
137                        break
138                    xi = start+step*j
139                    yi = sfloat(S[i:i+8])
140                    ei = sfloat(S[i+8:i+16])
141                    x.append(xi)
142                    if yi > 0.0:
143                        y.append(yi)
144                        w.append(1.0/ei**2)
145                    else:             
146                        y.append(0.0)
147                        w.append(0.0)
148                    j += 1
149                S = File.readline()
150            N = len(x)
151            if self.clockWd:
152                if len(self.clockWd) == 1:
153                    x = Tmap2TOF(self.TimeMap['1'],self.clockWd['1'])
154                else:
155                    x = Tmap2TOF(self.TimeMap[cons[5]],self.clockWd[cons[5]])
156            return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
157       
158        def GetSTDdata(File,Pos,Bank):
159            File.seek(Pos)
160            cons = Bank.split()
161            Nch = int(cons[2])
162            if self.clockWd:
163                start = 0
164                step = 1
165            else:
166                start = float(cons[5])/100.0               #CW: from centidegrees to degrees
167                step = float(cons[6])/100.0                 #NB TOF 0.1*ms!
168            x = []
169            y = []
170            w = []
171            S = File.readline()
172            j = 0
173            while S and S[:4] != 'BANK' and S[0] != '#':
174                if 'TIME_MAP' in S or '\x1a' in S:
175                    break
176                for i in range(0,80,8):
177                    if S[i:i+10] == 10*' ':
178                        break
179                    xi = start+step*j
180                    ni = max(sint(S[i:i+2]),1)
181                    yi = max(sfloat(S[i+2:i+8]),0.0)
182                    if yi:
183                        vi = yi/ni
184                    else:
185                        yi = 0.0
186                        vi = 0.0
187                    j += 1
188                    if j < Nch:
189                        x.append(xi)
190                        if vi <= 0.:
191                            y.append(0.)
192                            w.append(0.)
193                        else:
194                            y.append(yi)
195                            w.append(1.0/vi)
196                S = File.readline()
197            N = len(x)
198            if self.clockWd:
199                if len(self.clockWd) == 1:
200                    x = Tmap2TOF(self.TimeMap['1'],self.clockWd['1'])
201                else:
202                    x = Tmap2TOF(self.TimeMap[cons[5]],self.clockWd[cons[5]])
203            return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
204           
205        def GetALTdata(File,Pos,Bank):
206            File.seek(Pos)
207            cons = Bank.split()
208            x = []
209            y = []
210            w = []
211            S = File.readline()
212            j = 0
213            while S and S[:4] != 'BANK' and S[0] != '#':
214                for i in range(0,80,20):
215                    if S[i:i+8] == 9*' ':
216                        break
217                    xi = sfloat(S[i:i+9])/3200.
218                    yi = sfloat(S[i+9:i+16])/1000.
219                    ei = sfloat(S[i+16:i+21])/1000.
220                    x.append(xi)
221                    if yi > 0.0:
222                        y.append(yi)
223                        w.append(1.0/ei**2)
224                    else:             
225                        y.append(0.0)
226                        w.append(0.0)
227                    j += 1
228                S = File.readline()
229            N = len(x)
230            if self.clockWd:
231                x = Tmap2TOF(self.TimeMap[cons[5]],self.clockWd[cons[5]])
232            return [np.array(x),np.array(y),np.array(w),np.zeros(N),np.zeros(N),np.zeros(N)]
233           
234        def GetTimeMap(File,Pos,TimeMap):
235            File.seek(Pos)
236            cons = TimeMap[8:].split()
237            mapNo = cons[0]
238            if 'TIME_MAP10' in TimeMap: mapNo = '1'   #HIPD cluge!
239            Nch = int(cons[1])
240            Nrec = int(cons[2])
241            clockWd = float(cons[4])/1000.          #in mus
242            TMap = np.zeros(Nch+2,dtype=int)
243            ind = 0
244            for i in range(Nrec):
245                S = File.readline().rstrip('\n')
246                vals = S.split()
247                for val in vals:
248                    TMap[ind] = int(val)
249                    ind += 1
250            TMap = np.reshape(TMap,(-1,3))
251            TMax = TMap[-1][0]
252            Nch = TMap[-2][0]+(TMax-TMap[-2][1]+TMap[-2][2]-1)/TMap[-2][2]
253            TMap[-1] = [Nch+1,TMax,0]
254            TMap = TMap.T
255            TMap[0] -= 1
256            return TMap.T,clockWd,mapNo
257           
258        def Tmap2TOF(TMap,clockWd):
259            TOF = []
260            Tch,T,Step = TMap[0]
261            for tmap in TMap[1:]:
262                tch,t,step = tmap
263                TOF += [T+Step*(i-Tch) for i in range(Tch,tch)]
264                Tch,T,Step = tmap
265            TOF = np.array(TOF)*clockWd
266            return TOF
267
268        Banks = []
269        Pos = []
270        rdbuffer = kwarg.get('buffer')
271        title = ''
272        comments = None
273        fp = open(filename,'Ur')        #needs the 'U' to maintain py 2.7 compatibility
274        # reload previously saved values - used for multibank reads
275        if self.repeat and rdbuffer is not None:
276            Banks = rdbuffer.get('Banks')
277            Pos = rdbuffer.get('Pos')
278            self.selections = rdbuffer.get('selections')
279            comments = rdbuffer.get('comments')
280
281        # read through the file and find the beginning of each bank
282        # Save the offset (Pos), BANK line (Banks), comments for each bank
283        #
284        # This is going to need a fair amount of work to track line numbers
285        # in the input file.
286        if len(Banks) != len(Pos) or len(Banks) == 0:
287            i = -1
288            while True:
289                i += 1
290                S = fp.readline()
291                if len(S) == 0: break
292                   
293                if i==0: # first line is always a comment
294                    self.errors = 'Error reading title'
295                    title = S[:-1]
296                    comments = [[title,]]
297                    continue
298                if i==1 and S[:4].lower() == 'inst' and ':' in S:
299                    # 2nd line is instrument parameter file (optional)
300                    self.errors = 'Error reading instrument parameter filename'
301                    self.instparm = S.split(':')[1].strip('[]').strip()
302                    continue
303                if S[0] == '#': # allow comments anywhere in the file
304                    # comments in fact should only preceed BANK lines
305                    comments[-1].append(S[:-1])
306                    continue
307                if S[:4] == 'BANK':
308                    self.errors = 'Error reading bank:'
309                    self.errors += '  '+str(S)
310                    comments.append([title,])
311                    Banks.append(S)
312                    Pos.append(fp.tell())
313                if S[:8] == 'TIME_MAP':     #assumes one time map; HIPPO has multiple time maps
314                    if len(Banks) == 0:
315                        self.errors = 'Error reading time map before any bank lines'
316                    else:
317                        self.errors = 'Error reading time map after bank:\n  '+str(Banks[-1])
318                    timemap,clockwd,mapNo = GetTimeMap(fp,fp.tell(),S)
319                    self.TimeMap[mapNo] = timemap
320                    self.clockWd[mapNo] = clockwd
321                   
322
323        # Now select the bank to read
324        if not Banks: # use of ContentsValidator should prevent this error
325            print (self.formatName+' scan error: no BANK records')
326            selblk = None # no block to choose
327            self.errors = 'No BANK records found (strange!)'
328            return False
329        elif len(Banks) == 1: # only one Bank, don't ask
330            selblk = 0
331        elif self.repeat and self.selections is not None:
332            # we were called to repeat the read
333            #print 'debug: repeat #',self.repeatcount,'selection',self.selections[self.repeatcount]
334            selblk = self.selections[self.repeatcount]
335            self.repeatcount += 1
336            if self.repeatcount >= len(self.selections): self.repeat = False
337        else:                       # choose from options
338            selblk = self.selections[0] # select first in list
339            if len(self.selections) > 1: # prepare to loop through again
340                self.repeat = True
341                self.repeatcount = 1
342                if rdbuffer is not None:
343                    rdbuffer['Banks'] = Banks
344                    rdbuffer['Pos'] = Pos
345                    rdbuffer['selections'] = self.selections
346                    rdbuffer['comments'] = comments
347
348        # got a selection, now read it
349        Bank = Banks[selblk]
350        try:
351            bnkNo = int(Bank.split()[1])
352        except ValueError:
353            bnkNo = 1
354        if 'FXYE' in Bank:
355            self.errors = 'Error reading FXYE data in Bank\n  '+Banks[selblk]
356            self.powderdata = GetFXYEdata(fp,Pos[selblk],Banks[selblk])
357        elif 'FXY' in Bank:
358            self.errors = 'Error reading FXY data in Bank\n  '+Banks[selblk]
359            self.powderdata = GetFXYdata(fp,Pos[selblk],Banks[selblk])
360        elif 'ESD' in Bank:
361            self.errors = 'Error reading ESD data in Bank\n  '+Banks[selblk]
362            self.powderdata = GetESDdata(fp,Pos[selblk],Banks[selblk])
363        elif 'STD' in Bank:
364            self.errors = 'Error reading STD data in Bank\n  '+Banks[selblk]
365            self.powderdata = GetSTDdata(fp,Pos[selblk],Banks[selblk])
366        elif 'ALT' in Bank:
367            self.errors = 'Error reading ALT data in Bank\n  '+Banks[selblk]
368            self.powderdata = GetALTdata(fp,Pos[selblk],Banks[selblk])
369        else:
370            self.errors = 'Error reading STD data in Bank\n  '+Banks[selblk]
371            self.powderdata = GetSTDdata(fp,Pos[selblk],Banks[selblk])
372
373        self.errors = 'Error processing information after read complete'
374        if comments is not None:
375            self.comments = comments[selblk]
376        self.powderentry[0] = filename
377        self.powderentry[1] = Pos # position offset (never used, I hope)
378        self.powderentry[2] = bnkNo #selblk+1 # bank number
379        self.idstring = ospath.basename(filename) + ' Bank '+str(bnkNo) #selblk+1)
380        self.numbanks=len(Banks)
381        # scan comments for temperature & radius
382        Temperature = 300.
383        for S in self.comments:
384            if 'Temp' in S.split('=')[0]:
385                try:
386                    Temperature = float(S.split('=')[1])
387                except:
388                    pass
389            elif 'Gonio' in S.split('=')[0]:
390                try:
391                    self.Sample['Gonio. radius'] = float(S.split('=')[1])
392                except:
393                    pass
394            elif 'Omega' in S.split('=')[0] or 'Theta' in S.split('=')[0]:  #HIPD weirdness
395                try:
396                    self.Sample['Omega'] = 90.-float(S.split('=')[1])
397                except:
398                    pass
399            elif 'Chi' in S.split('=')[0]:
400                try:
401                    self.Sample['Chi'] = -float(S.split('=')[1])
402                except:
403                    pass                   
404            elif 'Phi' in S.split('=')[0]:
405                try:
406                    self.Sample['Phi'] = float(S.split('=')[1])
407                except:
408                    pass                   
409        self.Sample['Temperature'] = Temperature
410        fp.close()
411        return True       
412
413def sfloat(S):
414    'convert a string to a float, treating an all-blank string as zero'
415    if S.strip():
416        return float(S)
417    else:
418        return 0.0
419
420def sint(S):
421    'convert a string to an integer, treating an all-blank string as zero'
422    if S.strip():
423        return int(S)
424    else:
425        return 0
426
Note: See TracBrowser for help on using the repository browser.