1
2
3 """
4 Toolbox for reading, parsing, producing line lists for numerous molecules, and
5 from different databases.
6
7 Author: R. Lombaert
8
9 """
10
11 import os
12 import re
13 import string
14 from astropy import units as u
15
16 import cc.path
17 from cc.tools.io import DataIO
18 from cc.modeling.objects import Transition
19
20
21
23
24 '''
25 Reading, managing and parsing of line lists taken from online databases.
26
27 '''
28
29 - def __init__(self,fn,x_min=1e-10,x_max=1e10,unit='micron',\
30 min_strength=None,max_exc=None):
31 '''
32 Creating a LineList object.
33
34 @param fn: The full path and filename to the spectroscopy file. Must
35 contain either the CDMS or the JPL string.
36 @type fn: str
37
38 @keyword x_min: Minimum value for reading data in frequency/wavelength
39 If default, the lowest frequency is the one in the file.
40
41 (default: 1e-10)
42 @type x_min: float
43 @keyword x_max: Maximum value for reading data in frequency/wavelength
44 If default, the highest frequency is the one in the file
45
46 (default: 1e10)
47 @type x_max: float
48 @keyword unit: Unit of x_min/x_max (micron,GHz,MHz,Hz). Any
49 astropy.constants unit works. This can also be a Unit()
50 class object.
51
52 (default: 'micron')
53 @type unit: string
54 @keyword min_strength: if None all are included, else only lines with
55 strengths above this value are included
56 (log scale, nm2*Mhz)
57
58 (default: None)
59 @type min_strength: float
60 @keyword max_exc: if None all are included, else only lines with
61 excitation energies below this value are included
62 (cm-1)
63
64 (default: None)
65 @type max_exc: float
66
67 '''
68
69 self.fn = fn
70 self.min_strength = min_strength
71 self.max_exc = max_exc
72
73
74 if self.fn.upper().find('JPL') != -1:
75 self.catstring = 'JPL'
76 else:
77 self.catstring = 'CDMS'
78
79
80 if isinstance(unit,str) and unit.lower() in ['1 / cm','cm-1','cm^-1']:
81 unit = u.Unit("1 / cm")
82 elif isinstance(unit,str):
83 unit = getattr(u,unit)
84
85
86
87
88 if unit.is_equivalent(u.micron):
89 self.x_min = (x_max*unit).to(u.MHz,equivalencies=u.spectral())
90 self.x_max = (x_min*unit).to(u.MHz,equivalencies=u.spectral())
91 else:
92 self.x_max = (x_max*unit).to(u.MHz,equivalencies=u.spectral())
93 self.x_min = (x_min*unit).to(u.MHz,equivalencies=u.spectral())
94
95
96 self.line_list = []
97 self.read()
98
99
100
102
103 '''
104 Return an integer matching the given string from the catalog
105 (to deal with 'A#' numerals).
106
107 @param numeral: The numeral to be converted to integer
108 @type numeral: string
109
110 '''
111
112 try:
113 return int(numeral)
114 except ValueError:
115 alpha = string.letters[:26]
116 alphanumerals = [100+10*i for i,letter in enumerate(alpha)]
117 return alphanumerals[alpha.index(numeral[0])] + int(numeral[1])
118
119
120
122
123 '''
124 Parse Line Lists of standard format, such as for JPL or CDMS.
125
126 @param data: The content of the inputfile, no replaced spaces or
127 delimiter used when reading with DataIO.readFile!
128 @type data: list[string]
129
130 '''
131
132 data = [line
133 for line in data
134 if float(line[0:13]) <= self.x_max.value \
135 and float(line[0:13]) >= self.x_min.value]
136 if self.min_strength:
137 data = [line
138 for line in data
139 if float(line[21:29]) >= self.min_strength]
140 if self.max_exc:
141 data = [line
142 for line in data
143 if float(line[31:41]) <= self.max_exc]
144 vib_pattern = re.compile(r'(v\d?=\d)')
145
146 data = [[float(line[0:13]),\
147 line[61:63].strip() \
148 and self.makeCatInt(line[61:63]) \
149 or 0,self.makeCatInt(line[55:57]),\
150 line[57:59].strip() \
151 and self.makeCatInt(line[57:59]) \
152 or 0,\
153 line[59:61].strip() \
154 and self.makeCatInt(line[59:61]) \
155 or 0,\
156 line[73:75].strip() \
157 and self.makeCatInt(line[73:75]) \
158 or 0,\
159 self.makeCatInt(line[67:69]),\
160 line[69:71].strip() \
161 and self.makeCatInt(line[69:71]) \
162 or 0,\
163 line[71:73].strip() \
164 and self.makeCatInt(line[71:73]) \
165 or 0,\
166 vib_pattern.search(line[81:len(line)]) \
167 and vib_pattern.search(line[81:len(line)]).groups()[0]\
168 or '',\
169 self.catstring,\
170 float(line[21:29]),\
171 float(line[31:41])]
172 for line in data]
173 return data
174
175
176
178
179 '''
180 Read data from CDMS line list catalogs for a specific molecule.
181
182 '''
183
184 data = DataIO.readFile(self.fn,\
185 replace_spaces=0)
186 print 'Reading data from CDMS database for'
187 print self.fn
188
189
190
191 uncertainties = [float(line[13:21]) for line in data]
192 if min(uncertainties) < 0 and max(uncertainties) == 0:
193 self.x_min = self.x_min.to(1./u.cm,equivalencies=u.spectral())
194 self.x_max = self.x_max.to(1./u.cm,equivalencies=u.spectral())
195 elif min(uncertainties) < 0 and max(uncertainties) > 0:
196 raise ValueError('Uncertainties in CDMS input file for ' + \
197 'file %s are ambiguous.'\
198 %self.fn)
199
200 data = self.__parseCatalog(data)
201
202
203 rcm = u.Unit("1 / cm")
204 if self.x_min.unit == rcm:
205 data = sorted([[(entry*rcm).to(u.MHz,equivalencies=u.spectral())
206 if not i else entry
207 for i,entry in enumerate(line)]
208 for line in data])
209 self.line_list = data
210
211
212
214
215 '''
216 Read data from JPL line list catalogs for a specific molecule.
217
218 '''
219
220 data = DataIO.readFile(self.fn,\
221 replace_spaces=0)
222 print 'Reading data from JPL database for'
223 print self.fn
224 data = self.__parseCatalog(data)
225 self.line_list = data
226
227
228
230
231 '''
232 Start reading all requested line list data.
233
234 '''
235
236 self.line_list = []
237 if self.catstring == 'CDMS':
238 self.__readCDMS()
239 elif self.catstring == 'JPL':
240 self.__readJPL()
241
242
243
245
246 '''
247 Return the loaded line list, with frequencies in MHz.
248
249 '''
250
251 return self.line_list
252
253
254
255 - def makeTransitions(self,molecule,telescope=None,offset=0.0,n_quad=100,\
256 fraction_tau_step=1e-2,min_tau_step=1e-4,\
257 write_intensities=0,tau_max=12.,tau_min=-6.,\
258 check_tau_step=1e-2,):
259
260 '''
261 Make a list of transitions from the line list that was read.
262
263 Default transitions parameters are the same as those used in
264 Transition.Transition.__init__().
265
266 Requires a Molecule() object to be passed to the method call.
267
268 @param molecule: The molecule for these transitions
269 @type molecule: Molecule()
270
271 @keyword telescope: The telescope that observed given line
272
273 (default: None)
274 @type telescope: string
275 @keyword offset: The offset from center pixel of observed transition
276
277 (default: None)
278 @type offset: float
279 @keyword n_quad: The number of grid points in the formal integration
280 when calculating the line profile in sphinx.
281
282 (default: 100)
283 @type n_quad: int
284 @keyword fraction_tau_step: tau_total*fraction_tau_step gives min.
285 delta_tau in strahl.f. If too low,
286 min_tau_step will be used.
287
288 (default: 1e-2)
289 @type fraction_tau_step: float
290 @keyword min_tau_step: minimum of delta_tau in strahl.f
291
292 (default: 1e-4)
293 @type min_tau_step: float
294 @keyword write_intensities: set to 1 to write the intensities of first
295 50 impact-parameters at the end of sphinx
296
297 (default: 0)
298 @type write_intensities: bool
299 @keyword tau_max: maximum optical depth used for the calculation of the
300 formal integral
301
302 (default: 12.)
303 @type tau_max: float
304 @keyword tau_min: maximum optical depth used for the calculation of the
305 formal integral
306
307 (default: -6.)
308 @type tau_min: float
309 @keyword check_tau_step: check.par.in sphinx if step in tau not too
310 large
311
312 (default: 0.01)
313 @type check_tau_step: float
314
315 @return: The transitions
316 @rtype: list[Transition]
317
318 '''
319
320 pars = {'fraction_tau_step':fraction_tau_step,\
321 'min_tau_step':min_tau_step,'tau_max':tau_max,\
322 'write_intensities':write_intensities,'tau_min':tau_min,\
323 'check_tau_step':check_tau_step,'n_quad':n_quad,\
324 "offset":offset,'telescope':telescope}
325 if molecule.spec_indices == 0:
326 trans_list = [Transition.Transition(molecule=molecule,\
327 frequency=float(trans[0])*10**6,\
328 exc_energy=float(trans[12]),\
329 int_intensity_log=float(trans[11]),\
330 vup=int(trans[3]),\
331 jup=int(trans[2]),\
332 vlow=int(trans[7]),\
333 jlow=int(trans[6]),\
334 vibrational=trans[9],\
335 **pars)
336 for trans in self.getLineList()]
337 else:
338 trans_list = [Transition.Transition(molecule=molecule,\
339 frequency=float(trans[0])*10**6,\
340 exc_energy=float(trans[12]),\
341 int_intensity_log=float(trans[11]),\
342 vup=int(trans[1]),\
343 jup=int(trans[2]),\
344 kaup=int(trans[3]),\
345 kcup=int(trans[4]),\
346 vlow=int(trans[5]),\
347 jlow=int(trans[6]),\
348 kalow=int(trans[7]),\
349 kclow=int(trans[8]),\
350 vibrational=trans[9],\
351 **pars)
352 for trans in self.getLineList()]
353
354 return trans_list
355