1
2
3 """
4 A class for reading and managing collisional rate data.
5
6 Author: R. Lombaert
7
8 """
9
10 import os, collections
11 import numpy as np
12
13 from cc.tools.readers.SpectroscopyReader import SpectroscopyReader
14 from cc.tools.io import DataIO
15
16 import matplotlib.pyplot as p
17
18 from scipy.interpolate import interp1d
19 from scipy.interpolate import InterpolatedUnivariateSpline as spline1d
20
21
22
24
25 '''
26 A Reader for collisional rate data.
27
28 This class inherits from SpectroscopyReader.
29
30 Methods for reading and managing collisional rate data are provided. In case
31 the basic methods are not sufficient, subclasses provide code-specific
32 methods. The basic methods apply to the collis files for GASTRoNOoM.
33
34 Inheriting from CollisReader are:
35 - LamdaReader for reading Lamda molecular data files, used in MCP/ALI.
36
37 Other molecular spectroscopic data are typically handled by MolReader. Level
38 populations are handled by PopReader. For GASTRoNOoM this is combined in
39 MlineReader for rad trans output, and available through RadiatReader for
40 spectroscopic input.
41
42 '''
43
45
46 '''
47 Initialize an CollisReader object. The filename and additional
48 args/kwargs are passed to the Reader parent class, where the dict is
49 made and the filename is stored.
50
51 Additional args/kwargs are used for the dict creation.
52
53 @param fn: The filename of collisional rate data.
54 @type fn: str
55
56 '''
57
58 super(CollisReader,self).__init__(fn,*args,**kwargs)
59
60
61 if self.__class__ == CollisReader:
62 self.read()
63
64
65
67
68 '''
69 Read the collision rates file. Assumes GASTRoNOoM format.
70
71 To read ALI/MCP collision rates (which are in the Lamda format), make
72 use of the LamdaReader, which redefines this method. The other retrieval
73 methods remain valid.
74
75 Each transition is stored as an index, and gives upper and lower level
76 index as well as the collision rate.
77
78 '''
79
80
81
82 collis = np.loadtxt(self.fn)
83 self['pars'] = dict()
84
85
86
87 ntrans = DataIO.findZero(0,collis)
88 n0 = DataIO.findNumber(ntrans,collis)-ntrans
89
90
91 ntemp = (len(collis)-2*(ntrans+n0))/(ntrans+n0+1)
92
93
94 dtype = [('index',int),('lup',int),('llow',int),('rates',np.ndarray)]
95 self['coll_trans'] = np.empty(shape=(ntrans,),dtype=dtype)
96 self['coll_trans']['index'] = range(1,ntrans+1)
97
98
99 self['coll_trans']['lup'] = collis[:ntrans]
100 self['coll_trans']['llow'] = collis[ntrans+n0:2*ntrans+n0]
101
102
103 start_i = 2*(ntrans+n0)
104 Tgrid = [collis[start_i+i*(ntrans + n0 + 1)] for i in range(ntemp)]
105
106
107
108 Tgrid_real = [Ti for Ti in Tgrid if Ti != 0.]
109 ntemp_real = len(Tgrid_real)
110
111
112
113 rates = np.empty(shape=(ntrans,ntemp_real))
114 start_i = start_i + 1
115 for i,Ti in enumerate(Tgrid):
116 if Ti == 0.0: continue
117 this_i = start_i+i*(ntrans + n0 + 1)
118 rates[:,i] = collis[this_i:this_i+ntrans]
119
120
121 for i in range(ntrans):
122 self['coll_trans']['rates'][i] = rates[i,:]
123
124
125 self['pars']['ncoll_trans'] = ntrans
126 self['pars']['ncoll_temp'] = ntemp_real
127 self['coll_temp'] = np.array(Tgrid_real)
128
129
130
131 - def getTI(self,itype='coll_trans',lup=None,llow=None):
132
133 '''
134 Return the indices of the transitions read from the collisional rate
135 data.
136
137 A specific index (or array of indices) can be requested by specifying
138 the lower and upper level index.
139
140 LamdaReader first inherits from MolReader, so the default will be itype
141 'trans' for LamdaReader.
142
143 @keyword itype: The type of index. Default is for CollisReader. Needs to
144 be here in case LamdaReader calls this method, where the
145 type has to be specified. It will call both this method
146 after the version in MolReader and pass on the call.
147
148 (default: 'coll_trans')
149 @type itype: str
150 @keyword lup: The index of the upper energy level in the transition. If
151 both this and llow are None, all transition indices are
152 returned.
153
154 (default: None)
155 @type lup: int
156 @keyword llow: The index of the lower energy level in the transition. If
157 both this and lup are None, all transition indices are
158 returned.
159
160 (default: None)
161 @type llow: int
162
163 @return: The transition indices
164 @rtype: array
165
166 '''
167
168 return super(CollisReader,self).getTI(itype=itype,lup=lup,llow=llow)
169
170
171
172 - def getTUpper(self,index=None,itype='coll_trans'):
173
174 '''
175 Return the indices of the upper states of the <nline> included
176 transitions.
177
178 These are NOT the quantum numbers! Especially for CO-type molecules,
179 the J-level is equal to level_index-1.
180
181 In case a single upper level is requested via index, the level index is
182 extracted from the array.
183
184 @keyword index: The index. In case of default, all are returned. Can be
185 any array-like object that includes indices
186
187 (default: None)
188 @type index: int/array
189 @keyword itype: The type of index. Default is for CollisReader.
190 MolReader needs this to be 'trans'.
191
192 (default: 'coll_trans')
193 @type itype: str
194
195 @return: The upper level indices/x ordered by transition index.
196 @rtype: array
197
198 '''
199
200 return super(CollisReader,self).getTUpper(itype=itype,index=index)
201
202
203
204 - def getTLower(self,index=None,itype='coll_trans'):
205
206 '''
207 Return the indices of the lower states of the <nline> included
208 transitions.
209
210 These are NOT the quantum numbers! Especially for CO-type molecules,
211 the J-level is equal to level_index-1.
212
213 In case a single lower level is requested via index, the level index is
214 extracted from the array.
215
216 @keyword index: The index. In case of default, all are returned. Can be
217 any array-like object that includes indices
218
219 (default: None)
220 @type index: int/array
221 @keyword itype: The type of index. Default is for CollisReader.
222 MolReader needs this to be 'trans'.
223
224 (default: 'coll_trans')
225 @type itype: str
226
227 @return: The lower level indices/x ordered by transition index.
228 @rtype: array
229
230 '''
231
232 return super(CollisReader,self).getTLower(itype=itype,index=index)
233
234
235
236 - def getRates(self,index=None,lup=None,llow=None):
237
238 '''
239 Retrieve the collision rates, sorted by collisional transition index.
240
241 An index can be specified to retrieve a specific set of collision rates.
242
243 Alternatively, the lup and llow can be specified. index takes
244 precedence
245
246 The rates are given as a function of the temperature, accessible by
247 getTemp.
248
249 In case a single set of rates is requested via index, the rates array is
250 extracted from the encompassing array.
251
252 @keyword index: The index. In case of default, all are returned. Can be
253 any array-like object that includes indices
254
255 (default: None)
256 @type index: int/array
257 @keyword lup: The index of the upper energy level in the transition. If
258 both this and llow are None, all transition indices are
259 returned.
260
261 (default: None)
262 @type lup: int
263 @keyword llow: The index of the lower energy level in the transition. If
264 both this and lup are None, all transition indices are
265 returned.
266
267 (default: None)
268 @type llow: int
269
270 @return: The collision rates (in cm^3 s^-1) are returned as arrays
271 as a function of temperature for given coll_trans indices.
272 @rtype: array
273
274 '''
275
276
277
278 if not index is None or llow is None or lup is None:
279 return self.get('coll_trans','rates',index)
280
281
282 index = self.getTI(llow=llow,lup=lup,itype='coll_trans')
283 return self.get('coll_trans','rates',index)
284
285
286
288
289 '''
290 Retrieve the temperature grid for the collision rates.
291
292 @return: The temperature grid in K
293 @rtype: array
294
295 '''
296
297 return self['coll_temp']
298
299
300
301 - def setInterp(self,index=None,itype='spline',*args,**kwargs):
302
303 '''
304 Set the interpolator for the collision rates versus temperature.
305
306 Additional arguments can be passed to the interpolator object.
307
308 @keyword index: The transition index. If default, all collision rates
309 are interpolated. If an iterable object (such as a list)
310 the method iterates over the indices. If a single value,
311 only one set of rates is interpolated.
312
313 (default: None)
314 @type index: int/list
315 @keyword itype: The type of interpolator. Either spline or linear.
316
317 (default: 'spline')
318 @type itype: str
319
320 '''
321
322
323 itype = itype.lower()
324 if itype == 'linear':
325 interp = interp1d
326 else:
327 interp = spline1d
328
329
330 if index is None:
331 index = range(1,self['pars']['ncoll_trans']+1)
332 elif not isinstance(indices,collections.Iterable) \
333 or isinstance(indices,str):
334 index = [index]
335
336 self['icoll'] = dict()
337
338
339 for i in index:
340 self['icoll'][i] = interp(x=self['coll_temp'],\
341 y=self.get('coll_trans','rates',i),\
342 *args,**kwargs)
343
344
345
347
348 '''
349 Get the interpolator for a given transition index.
350
351 @param index: The level index
352 @type index: int
353
354 @return: The colision rates (in cm^3 s^-1) interpolator as a function of
355 temperature in K.
356 @rtype: interpolator
357
358 '''
359
360 return self['icoll'][index]
361
362
363
365
366 '''
367 Plot the collision rates for given transitions as a function of
368 temperature.
369
370 @keyword fn: Filename of the plot, including path and extension. If None
371 the plot is shown in the python session.
372
373 (default: None)
374 @type fn: str
375 @keyword indices: The coll_trans indices to be plotted. Default in case
376 all transitions are to be plotted.
377
378 (default: None)
379 @type indices: array
380
381 '''
382
383 p.clf()
384 fig1 = p.figure(1)
385 if indices is None:
386 indices = range(1,self['pars']['ncoll_trans']+1)
387 for ii in indices:
388 p.semilogy(self.getTemp(),self.getRates(index=ii))
389 p.xlabel('T (K)', fontsize = 14)
390 p.ylabel('Collision Rates (cm$^3$ s$^{-1}$)', fontsize = 14)
391 p.xlim(xmin=min(self.getTemp()))
392 p.xlim(xmax=max(self.getTemp()))
393 if not fn is None:
394 p.savefig(fn,bbox_inches='tight')
395 else:
396 p.show()
397