source: trunk/imports/G2pwd_fxye.py @ 3266

Last change on this file since 3266 was 3266, checked in by vondreele, 4 years ago

change modulated wave type assignments so one for each kind of wave (pos, frac, etc).
implement optional display of PNT data set positions on pole figure plots. Picker displays histo. name.

  • 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: 2018-02-08 19:31:55 +0000 (Thu, 08 Feb 2018) $
4# $Author: vondreele $
5# $Revision: 3266 $
6# $URL: trunk/imports/G2pwd_fxye.py $
7# $Id: G2pwd_fxye.py 3266 2018-02-08 19:31:55Z 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: 3266 $")
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'] = 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.