1
2
3 """
4 Toolbox for Transitions, used in various applications concerning GASTRoNOoM.
5
6 Author: R. Lombaert
7
8 """
9
10 import os
11 import re
12 import subprocess
13 import copy
14 from glob import glob
15 from scipy import pi, exp, linspace, argmin, array, diff, mean, isfinite
16 from scipy.interpolate import interp1d
17 from scipy.integrate import trapz
18 import numpy as np
19 from astropy import units as u
20 import types
21
22 from cc.ivs.sigproc import funclib
23
24 import cc.path
25 from cc.ivs.sigproc import funclib
26 from cc.modeling.objects import Molecule
27 from cc.tools.io import Database, DataIO
28 from cc.tools.readers import SphinxReader
29 from cc.tools.readers import FitsReader, TxtReader
30 from cc.tools.numerical import Interpol
31 from cc.statistics import BasicStats as bs
32 from cc.data import LPTools
33 from cc.tools.units import Equivalency as eq
34
35
37
38
39 '''
40 Get the line strengths from Transition() objects, either from data or from
41 models defined by the mode.
42
43 Additional keywords for the functions that return LSs from the Transition()
44 object can be passed along as **kwargs.
45
46 Modes:
47 - dint: Integrated line strength of an observed line in spectral mode
48 - mint: Intrinsic integrated line strength of a modeled line
49 - dtmb: Integrated line strength in main-beam temperature (Kelvin*km/s)
50 of a heterodyne instrument
51 - mtmb: Modeled line strength after convolution with beam profile and
52 conversion to main-beam temperature.
53 - cint: A combination of dint and mint. Data objects are assumed to be
54 listed first.
55 - ctmb: A combination of dtmb and mtmb. Data objects are assumed to be
56 listed first.
57 For data: Blended lines are returned as negative. If a line is in a blend,
58 but not attributed a line strength, the line strength of the 'main
59 component' is returned, also as a negative value.
60
61 For now no errors for any mode, except mode==dint or cint.
62
63 @param trl: The Transition() objects.
64 @type trl: list[Transition()]
65
66 @keyword mode: The mode in which the method is called. Either 'dint',
67 'mint', 'mtmb', 'dtmb', 'cint' or 'ctmb' values.
68
69 (default: 'dint')
70 @type mode: str
71 @keyword nans: Set undefined line strengths as nans. Errors are set as a
72 nan if it concerns mode==dint. Otherwise, they are not set.
73
74 (default: 1)
75 @type nans: bool
76 @keyword n_data: The number of data Star() objects, assuming they are the
77 first in the star_grid. Only required if mode == 'combo'.
78
79 (default: 0)
80 @type n_data: int
81 @keyword scale: Scale data to antenna of 1 m**2, necessary if comparing data
82 from different telescopes
83
84 (default: 0)
85 @type scale: bool
86
87 @return: The requested line strengths in W/m2, or K*km/s, as well as errors
88 if applicable. If a combo mode is requested, errors are given when
89 available, and listed as None if not available (Plotting2 module
90 knows how to deal with this).
91 @rtype: (list[float],list[float])
92
93 '''
94
95 modes = {'dint': 'getIntIntUnresolved',\
96 'mint': 'getIntIntIntSphinx',\
97 'dtmb': 'getIntTmbData',\
98 'mtmb': 'getIntTmbSphinx'}
99
100 n_data = int(n_data)
101 if mode[0] == 'c':
102 mode = 'd%s'%mode[1:]
103 elif mode[0] == 'd':
104 n_data = len(trl)
105 else:
106 n_data = 0
107 mode = mode.lower()
108 if mode not in modes.keys():
109 print 'Mode unrecognized. Can be dint, mint, dtmb, mtmb, cint or ctmb.'
110 return
111
112 allints = []
113 allerrs = []
114 for it,t in enumerate(trl):
115 scaling = t.telescope_size**2
116 if n_data and it == n_data:
117 mode = 'm%s'%mode[1:]
118 if t is None:
119 allints.append(nans and float('nan') or None)
120 continue
121
122 nls = getattr(t,modes[mode])(**kwargs)
123
124 if mode == 'dint':
125 if nls[0] == 'inblend':
126 for tb in nls[2]:
127 bls,ebls,blends = getattr(tb,modes[mode])(**kwargs)
128 if bls != 'inblend':
129 nls = (abs(bls)*(-1),ebls,blends)
130 break
131 nint = (nans and nls[0] is None) and float('nan') or nls[0]
132 nerr = (nans and nls[0] is None) and float('nan') or nls[1]
133 allints.append(nint)
134 allerrs.append(nerr)
135
136
137 elif mode == 'dtmb':
138 if scale == 1:
139 allints.append(((nans and nls is None) and float('nan') or nls[0])/scaling)
140 allerrs.append((nans and nls is None) and float('nan') or nls[1])
141 else:
142 allints.append((nans and nls is None) and float('nan') or nls[0])
143 allerrs.append((nans and nls is None) and float('nan') or nls[1])
144
145 else:
146 if scale == 1:
147 allints.append(((nans and nls is None) and float('nan') or nls)/scaling)
148 allerrs.append(None)
149 else:
150 allints.append((nans and nls is None) and float('nan') or nls)
151 allerrs.append(None)
152
153
154 allints, allerrs = array(allints), array(allerrs)
155 return (allints,allerrs)
156
157
159
160 '''
161 Select a transition from a list of Star() objects based on a given
162 criterion.
163
164 Different modes are possible: 'index' and 'sample'. The former returns a
165 transition with given list-index in the Star()['GAS_LINES'] which should
166 always return the same transition for every object. The latter returns all
167 transitions that are equal to the a given sample transition (ie following
168 the equality rules of a Transition() object).
169
170 @param sg: The grid of Star() objects.
171 @type sg: list[Star()]
172 @param criterion: The selection criterion, based on the mode.
173 @type criterion: int/Transition()
174
175 @keyword mode: The selection mode. For now either 'index' or 'sample'.
176
177 (default: 'index')
178 @type mode: string
179
180 @return: The selected Transition() objects are returned.
181 @rtype: list[Transition()]
182
183 '''
184
185 if mode.lower() == 'index':
186 criterion = int(criterion)
187 elif mode.lower() == 'sample':
188 pass
189 else:
190 print 'Keyword mode was not recognized. Set to "index" or "sample".'
191 return
192
193 if mode == 'index':
194 trl = [s['GAS_LINES'][criterion] for s in sg]
195 elif mode == 'sample':
196 trl = [s.getTransition(criterion) for s in sg]
197 return trl
198
199
200
203
204 '''
205 Extract a list a of unique transitions included in a list of Star() objects
206 and sort them.
207
208 A selection can be made on data type, which is either a telescope or
209 instrument, or all, or all unresolved, or all resolved.
210
211 The list of transitions is copied to make sure no funky references mess
212 with the original transitions.
213
214 The data can be reset, since data do not uniquely identify a transition
215 object. However, in some cases it is useful to keep the data in the object
216 if the extracted list will always work for the same Star() object.
217
218 @param star_grid: The list of Star() objects from which all 'GAS_LINES'
219 keys are taken and collected in a set.
220 @type star_grid: list[Star()]
221
222 @keyword sort_freq: Sort the transitions according to their frequencies.
223 Otherwise, they are sorted by their wavelengths.
224
225 (default: 1)
226 @type sort_freq: bool
227 @keyword sort_molec: Sort the transitions by molecule first.
228
229 (default: 1)
230 @type sort_molec: bool
231 @keyword dtype: 'all': return all lines,
232 'resolved': return all resolved lines,
233 'unresolved': return all unresolved lines,
234 'pacs'/'spire'/'apex'/'jcmt'/...: return all telescope or
235 instrument specific lines selection is based on presence of
236 string in telescope name.
237 Invalid definitions return an empty list.
238
239 (default: 'all')
240 @type dtype: str
241 @keyword reset_data: Reset the data properties to default, so they have to
242 be read again.
243
244 (default: 1)
245 @type reset_data: bool
246
247 @return: a list of unique transitions included in all Star() objects in
248 star_grid
249 @rtype: list[Transition()]
250
251 '''
252
253 dtype = dtype.upper()
254 selection = sorted(set([trans
255 for star in star_grid
256 for trans in star['GAS_LINES']
257 if trans]),\
258 key=lambda x:sort_freq \
259 and (sort_molec and x.molecule.molecule or '',\
260 x.frequency) \
261 or (sort_molec and x.molecule.molecule or '',\
262 x.wavelength))
263
264 if dtype == 'UNRESOLVED':
265 selection = [trans for trans in selection if trans.unresolved]
266 elif dtype == 'RESOLVED':
267 selection = [trans for trans in selection if not trans.unresolved]
268 elif not dtype == 'ALL':
269 selection = [trans for trans in selection if dtype in trans.telescope]
270
271 selection = [copy.deepcopy(trans) for trans in selection]
272 if reset_data:
273 for t in selection: t.resetData()
274
275 return selection
276
277
278
280
281 '''
282 Update telescope.spec file.
283
284 This method checks for requested transitions that are not yet present in
285 .spec file of the telescope, and adds them.
286
287 @param trans_list: The transitions in the CC input
288 @type trans_list: list[Transition()]
289
290 '''
291
292 telescopes = list(set([trans.telescope for trans in trans_list]))
293 for telescope in telescopes:
294 if not ('HIFI' in telescope):
295 print 'Warning! Telescope beam efficiencies for %s'%telescope + \
296 ' are added arbitrarily and thus are not the correct values.'
297 old_spec = DataIO.readFile(os.path.join(cc.path.gdata,telescope+'.spec'))
298 line_spec_list = [line
299 for line in old_spec
300 if line.find('LINE_SPEC') >-1 and line[0] != '#']
301 line_spec_quant = [line.split()[0:9]
302 for line in line_spec_list]
303
304
305
306
307
308
309 new_lsl = []
310 new_lsq = []
311 for lsq,lsl in zip(line_spec_quant,line_spec_list):
312 if line_spec_quant.count(lsq) == 1:
313 new_lsq.append(lsq)
314 new_lsl.append(lsl)
315 else:
316 if lsq not in new_lsq:
317 if float(lsl.split()[10]) != 1.:
318 new_lsq.append(lsq)
319 new_lsl.append(lsl)
320
321 these_trans = [tr
322 for tr in trans_list
323 if tr.telescope == telescope]
324 these_trans = [tr.getLineSpec()
325 for tr in these_trans
326 if tr.getLineSpec().split()[0:9] not in new_lsq]
327
328 new_lsl = sorted(new_lsl + these_trans,\
329 key=lambda x: [x.split()[0],float(x.split()[9])])
330 new_lsl = [line.replace(' ','\t') for line in new_lsl]
331 try:
332 line_spec_index = [line[0:9]
333 for line in old_spec].index('LINE_SPEC')
334 new_spec = old_spec[0:line_spec_index] + new_lsl
335 except ValueError:
336 new_spec = old_spec + new_lsl
337
338 DataIO.writeFile(os.path.join(cc.path.gdata,telescope+'.spec'),\
339 new_spec+['\n######################################'])
340
341
342
344
345 '''
346 Create a Transition instance based on a Star object and a standard CC input
347 line, with 11 or 12 entries in a list. This method assumes the Transition
348 definition has been screened by DataIO.checkEntryInfo. If the 12th entry is
349 not present, n_quad is taken from the Star() object. It is set to 100 if no
350 Star() object is given.
351
352 If a star object is not given, the method creates Molecule() objects itself
353 For now only 12C16O, 13C16O, 1H1H16O and p1H1H16O are possible.
354
355 @param trans: the input line with 11 or 12 entries, the first being the
356 molecule, followed by all 10 or 11 CC input parameters for a
357 transition. It is possible also to give the GASTRoNOoM syntax
358 for a transition, without n_quad in the end. In that case,
359 n_quad is taken from Star() or equal to 100. Any parameters
360 given after the 11th, respectively the 10th, parameter are
361 ignored.
362 @type trans: list[string]
363
364 @keyword star: The star object providing basic info
365
366 (default: None)
367 @type star: Star()
368 @keyword def_molecs: Default molecules needed for the requested transitions
369 None is returned if the requested molecule is not
370 present. If both this and star are not given, a few
371 default molecules are loaded.
372
373 (default: None)
374 @type def_molecs: dict(string: Molecule())
375 @keyword kwargs: Any additional keywords given here are added to the
376 Transition object creation. If not valid there, an error
377 will be thrown. The keys given here will overwrite what is
378 in star! Usually, only use this option when working with
379 defaults (and star is not given)
380
381 (default: {})
382 @type kwargs: dict
383
384 @return: The transition object is returned with all info included
385 @rtype: Transition()
386
387 '''
388
389
390 if trans[0].find('.') != -1:
391 imolec = DataIO.getInputData(keyword='MOLEC_TYPE',\
392 filename='Molecule.dat',\
393 make_float=0).index(trans[0])
394 molec_short = DataIO.getInputData(keyword='TYPE_SHORT',\
395 filename='Molecule.dat')[imolec]
396 else:
397 molec_short = trans[0]
398
399
400 if star is None and def_molecs is None:
401 def_molecs = {'12C16O':Molecule.Molecule('12C16O',61,61,240),\
402 '13C16O':Molecule.Molecule('13C16O',61,61,240),\
403 '1H1H16O':Molecule.Molecule('1H1H16O',39,90,1157),\
404 'p1H1H16O':Molecule.Molecule('p1H1H16O',32,90,1029)}
405
406
407 if star is None:
408 star = Star.Star()
409 molec = def_molecs.get(molec_short,None)
410 path_gastronoom = None
411 else:
412 molec = star.getMolecule(molec_short)
413 path_gastronoom = star.path_gastronoom
414
415
416 n_quad = star['N_QUAD']
417
418
419
420
421
422 if len(trans) > 11:
423 try:
424 n_quad = int(trans[11])
425 except ValueError:
426 pass
427
428
429 extra_pars = {'fraction_tau_step':star['FRACTION_TAU_STEP'],\
430 'min_tau_step':star['MIN_TAU_STEP'],\
431 'write_intensities':star['WRITE_INTENSITIES'],\
432 'tau_max':star['TAU_MAX'],'tau_min':star['TAU_MIN'],\
433 'check_tau_step':star['CHECK_TAU_STEP']}
434 extra_pars.update(kwargs)
435
436 if not molec is None:
437 return Transition(molecule=molec,n_quad=n_quad,\
438 vup=int(trans[1]),jup=int(trans[2]),\
439 kaup=int(trans[3]),kcup=int(trans[4]),\
440 vlow=int(trans[5]),jlow=int(trans[6]),\
441 kalow=int(trans[7]),kclow=int(trans[8]),\
442 telescope=trans[9],offset=float(trans[10]),\
443 path_gastronoom=path_gastronoom,\
444 **extra_pars)
445 else:
446 return None
447
448
449
451
452 '''
453 Return the modelids for a given model folder, as well as path_gastronoom.
454
455 @param filepath: The path to the model folder.
456 @type filepath: str
457
458 @return: the 3 ids and the path_gastronoom are returned
459 @rtype: (model_id,molec_id,trans_id,path_gastronoom)
460
461 '''
462
463
464 filepath,trans_id = os.path.split(filepath)
465
466 filepath = os.path.split(filepath)[0]
467
468 path_gastronoom = os.path.split(filepath)[1]
469
470
471 cc.path.gout = os.path.join(cc.path.gastronoom,path_gastronoom)
472
473 cooling_log = os.path.join(cc.path.gout,'models',trans_id,'cooling_id.log')
474 if os.path.isfile(cooling_log):
475 model_id = DataIO.readFile(cooling_log)[0]
476
477 mline_log = os.path.join(cc.path.gout,'models',trans_id,'mline_id.log')
478 if os.path.isfile(mline_log):
479 molec_id = DataIO.readFile(mline_log)[0]
480
481 else:
482 molec_id = trans_id
483
484 else:
485 model_id = trans_id
486 molec_id = trans_id
487
488 return (model_id,molec_id,trans_id,path_gastronoom)
489
490
491
493
494 '''
495 Make a Transition() based on the filename of a Sphinx file.
496
497 For this, information is taken from the Molecule() model database and the
498 sphinx file has to be associated with a molecule available there.
499
500 Information is also be pulled from the sph database, but can be turned off.
501
502 @param filename: The sphinx file name, including path
503 @type filename: string
504
505 @keyword mline_db: The mline database, which can be passed in case one
506 wants to reduce overhead. Not required though.
507
508 (default: None)
509 @type mline_db: Database()
510
511 @return: The transition
512 @rtype: Transition()
513
514 '''
515
516 filepath,fn = os.path.split(filename)
517 if not filepath:
518 raise IOError('Please include the full filepath of the Sphinx ' + \
519 'filename, needed to determine the name of the ' + \
520 'database.')
521
522 model_id,molec_id,trans_id,path_gastronoom = getModelIds(filepath)
523
524
525 file_components = fn.rstrip('.dat').split('_')[2:]
526 molec_name = file_components.pop(0)
527
528
529
530 molec = Molecule.makeMoleculeFromDb(molecule=molec_name,molec_id=molec_id,\
531 path_gastronoom=path_gastronoom,\
532 mline_db=mline_db)
533
534
535 if molec is None:
536 return None
537
538 telescope = file_components.pop(-2)
539
540
541 pattern = re.compile(r'^(\D+)(\d*.?\d*)$')
542 numbers = [pattern.search(string).groups() for string in file_components]
543 numbers = dict([(s.lower(),n) for s,n in numbers])
544
545
546 fn = os.path.join(filepath,'sphinx_parameters.log')
547 if not os.path.isfile(fn):
548 print 'File sphinx_parameters.log does not exist. Using default ' + \
549 'values for various sphinx parameters. THIS MAY BE INCORRECT. ' +\
550 'Currently only extracting n_quad. Talk to Robin Lombaert.'
551 sph2file = DataIO.readFile(filename=filename,delimiter=None)
552 i1 = sph2file[3].find('frequency points and ')
553 i2 = sph2file[3].find(' impact parameters')
554 ddict = {}
555 ddict['n_quad'] = int(sph2file[3][i1+len('frequency points and '):i2])
556 else:
557 ddict = DataIO.readDict(filename=fn,convert_floats=1,convert_ints=1)
558 numbers.update(ddict)
559
560
561 trans = Transition(molecule=molec,telescope=telescope,\
562 path_gastronoom=path_gastronoom,**numbers)
563
564
565 trans.setModelId(trans_id)
566
567 return trans
568
569
570
572
573 '''
574 Reconstruct a sphinx database based on existing model_ids, presence of sph2
575 files and inclusion in the mline database.
576
577 Not based at first on the cooling database because it is possible different
578 mline or sphinx models with different ids exist for the same cooling id.
579
580 @param path_gastronoom: The path_gastronoom to the output folder
581 @type path_gastronoom: string
582
583 '''
584
585
586 cc.path.gout = os.path.join(cc.path.gastronoom,path_gastronoom)
587
588
589 sph_db_path = os.path.join(cc.path.gout,'GASTRoNOoM_sphinx_models.db')
590 if os.path.isfile(sph_db_path):
591 i = 0
592 backup_file = '%s_backupSphDbRetrieval_%i'%(sph_db_path,i)
593 while os.path.isfile(backup_file):
594 i += 1
595 backup_file = '%s_backupSphDbRetrieval_%i'%(sph_db_path,i)
596 subprocess.call(['mv %s %s'%(sph_db_path,backup_file)],shell=True)
597
598
599 all_sph2 = glob(os.path.join(cc.path.gout,'models','*','sph2*'))
600
601
602 mline_db = Database.Database(os.path.join(cc.path.gout,\
603 'GASTRoNOoM_mline_models.db'))
604 trans_db = Database.Database(os.path.join(cc.path.gout,\
605 'GASTRoNOoM_sphinx_models.db'))
606
607 for i,sph2 in enumerate(all_sph2):
608 if i%1000 == 0:
609 print('Saving sphinx result %i, with filename %s.'%(i,sph2))
610 fp,filename = os.path.split(sph2)
611 model_id,molec_id,trans_id,pg = getModelIds(fp)
612 trans = makeTransitionFromSphinx(filename=sph2,\
613 mline_db=mline_db)
614
615
616 if trans is None:
617 continue
618 if not trans_db.has_key(model_id):
619 trans_db[model_id] = dict()
620 if not trans_db[model_id].has_key(molec_id):
621 trans_db[model_id][molec_id] = dict()
622 if not trans_db[model_id][molec_id].has_key(trans_id):
623 trans_db[model_id][molec_id][trans_id] = dict()
624 if not trans_db[model_id][molec_id][trans_id].has_key(str(trans)):
625 trans_db[model_id][molec_id][trans_id][str(trans)] \
626 = trans.makeDict()
627 else:
628 print "There's already an entry in the sphinx db for"
629 print "model_id: %s"%model_id
630 print "mline_id: %s"%molec_id
631 print "sphinx_id: %s"%trans_id
632 print "transition: %s"%str(trans)
633 print "Check what is going on!"
634
635 trans_db.sync()
636
637
638
639 -def makeTransitionsFromRadiat(molec,telescope,ls_min,ls_max,ls_unit='GHz',\
640 n_quad=100,offset=0.0,path_gastronoom=None,\
641 no_vib=0,fraction_tau_step=1e-2,\
642 write_intensities=0,tau_max=12.,tau_min=-6.,\
643 check_tau_step=1e-2,min_tau_step=1e-4):
644
645 '''
646 Make Transition() objects from a Radiat file of a molecule, within a given
647 wavelength/frequency range.
648
649 Requires a Molecule() object to work!
650
651 @param molec: The molecule for which the line list is made.
652 @type molec: Molecule()
653 @param telescope: The telescope for which the Transition() list is made.
654 @type telescope: string
655 @param ls_min: The minimum allowed wavelength/frequency for the transitions
656 @type ls_min: float
657 @param ls_max: The maximum allowed wavelength/frequency for the transitions
658 @type ls_max: float
659
660 @keyword ls_unit: The unit of the wavelength/frequency range. Can be: GHz,
661 MHz, Hz, MICRON, MM, CM, M
662
663 (default: 'GHz')
664 @type ls_unit: string
665 @keyword n_quad: The N_QUAD value for GASTRoNOoM sphinx calculations.
666
667 (default: 100)
668 @type n_quad: int
669 @keyword fraction_tau_step: tau_total*fraction_tau_step gives min.
670 delta_tau in strahl.f. If too low,
671 min_tau_step will be used.
672
673 (default: 1e-2)
674 @type fraction_tau_step: float
675 @keyword min_tau_step: minimum of delta_tau in strahl.f
676
677 (default: 1e-4)
678 @type min_tau_step: float
679 @keyword write_intensities: set to 1 to write the intensities of first
680 50 impact-parameters at the end of sphinx
681
682 (default: 0)
683 @type write_intensities: bool
684 @keyword tau_max: maximum optical depth used for the calculation of the
685 formal integral
686
687 (default: 12.)
688 @type tau_max: float
689 @keyword tau_min: maximum optical depth used for the calculation of the
690 formal integral
691
692 (default: -6.)
693 @type tau_min: float
694 @keyword check_tau_step: check.par.in sphinx if step in tau not too
695 large
696
697 (default: 0.01)
698 @type check_tau_step: float
699 @keyword offset: The offset from center position for calculations in sphinx
700
701 (default: 0.0)
702 @type offset: float
703 @keyword path_gastronoom: model output folder in the GASTRoNOoM home
704
705 (default: None)
706 @type path_gastronoom: string
707 @keyword no_vib: Do not include vibrational states in the output list.
708
709 (default: 0)
710 @type no_vib: bool
711
712 @return: The newly made Transition() objects for all transitions in range
713 @rtype: list[Transition()]
714
715 '''
716
717 trkeys = {'fraction_tau_step':fraction_tau_step,\
718 'min_tau_step':min_tau_step,'tau_max':tau_max,\
719 'write_intensities':write_intensities,'tau_min':tau_min,\
720 'check_tau_step':check_tau_step,'n_quad':n_quad,\
721 "offset":offset,'path_gastronoom':path_gastronoom,\
722 'telescope':telescope,'molecule':molec}
723 radiat = molec.radiat
724 wave = radiat.getTFrequency(unit=ls_unit)
725 low = radiat.getTLower()
726 up = radiat.getTUpper()
727 if not molec.spec_indices:
728
729
730
731 ny_low = molec.ny_low
732 n_vib = int(molec.ny_up/ny_low) + 1
733 indices = [[i+1,int(i/ny_low),i-ny_low*int(i/ny_low)]
734 for i in xrange(molec.ny_low*n_vib)]
735 nl = [Transition(vup=int(indices[u-1][1]),\
736 jup=int(indices[u-1][2]),\
737 vlow=int(indices[l-1][1]),\
738 jlow=int(indices[l-1][2]),\
739 **trkeys)
740 for l,u,w in zip(low,up,wave)
741 if w > ls_min and w < ls_max]
742 else:
743 indices = molec.radiat_indices
744 quantum = ['v','j','ka','kc']
745 nl = []
746 for l,u,w in zip(low,up,wave):
747 if w > ls_min and w < ls_max:
748 quantum_dict = dict()
749
750 for i in xrange(1,len(indices[0])):
751 quantum_dict[quantum[i-1]+'up'] = \
752 int(indices[u-1][i])
753 quantum_dict[quantum[i-1]+'low'] = \
754 int(indices[l-1][i])
755 trkeys.update(quantum_dict)
756 nl.append(Transition(**trkeys))
757 if no_vib:
758 nl = [line for line in nl if line.vup == 0]
759 return nl
760
761
762
764
765 '''
766 Check uniqueness of a list of transitions.
767
768 Same transitions are replaced by a single transition with all datafiles
769 from the originals.
770
771 Based on the parameters of the transitions. If they are the same, only one
772 transition is added to the output list, but the datafiles are all included.
773
774 Datafiles do not identify a transition!
775
776 @param trans_list: The transitions to be checked for uniqueness
777 @type trans_list: list[Transition()]
778
779 @return: The merged transitions with all datafiles included.
780 @rtype: tuple[Transition()]
781
782 '''
783
784 merged = []
785 for trans in trans_list:
786 if trans not in merged:
787 merged.append(trans)
788 else:
789
790 if not trans.datafiles is None:
791 ddict = dict(zip(trans.datafiles,trans.fittedlprof))
792 merged[merged.index(trans)].addDatafile(ddict)
793 return merged
794
795
796
798
799 '''
800 A class to deal with transitions in GASTRoNOoM.
801
802 '''
803
804 - def __init__(self,molecule,telescope=None,vup=0,jup=0,kaup=0,kcup=0,\
805 nup=None,vlow=0,jlow=0,kalow=0,kclow=0,nlow=None,offset=0.0,\
806 frequency=None,exc_energy=None,int_intensity_log=None,\
807 n_quad=100,fraction_tau_step=1e-2,min_tau_step=1e-4,\
808 write_intensities=0,tau_max=12.,tau_min=-6.,\
809 check_tau_step=1e-2,vibrational='',path_gastronoom=None):
810
811 '''
812
813 Initiate a Transition instance, setting all values for the
814 quantummechanical parameters to zero (by default).
815
816 @param molecule: The molecule to which this transition belongs
817 @type molecule: Molecule()
818
819 @keyword telescope: The telescope with which the transition is observed
820
821 (default: None)
822 @type telescope: string
823 @keyword vup: The upper vibrational level
824
825 (default: 0)
826 @type vup: int
827 @keyword jup: The upper rotational level
828
829 (default: 0)
830 @type jup: int
831 @keyword kaup: The upper level of the first projection of the ang mom
832
833 (default: 0)
834 @type kaup: int
835 @keyword kcup: The upper level of the second projection of the ang mom
836
837 (default: 0)
838 @type kcup: int
839 @keyword nup: if not None it is equal to Kaup and only relevant for
840 SO/hcn-type molecules
841
842 (default: 0)
843 @type nup: int
844 @keyword vlow: The lower vibrational level
845
846 (default: 0)
847 @type vlow: int
848 @keyword jlow: The lower rotational level
849
850 (default: 0)
851 @type jlow: int
852 @keyword kalow: The lower level of the first projection of the ang mom
853
854 (default: 0)
855 @type kalow: int
856 @keyword kclow: The lower level of the second projection of the ang mom
857
858 (default: 0)
859 @type kclow: int
860 @keyword nlow: if not None it is equal to Kalow, and only relevant for
861 SO-type molecules
862
863 (default: None)
864 @type nlow: int
865 @keyword offset: The offset of the radiation peak with respect to the
866 central pixel of the observing instrument. Only
867 relevant for non-point corrected PACS/SPIRE or
868 resolved data (so not relevant when the intrinsic line
869 strengths are used from the models)
870
871 (default: 0.0)
872 @type offset: float
873 @keyword n_quad: Number of impact par. in quadrature int(Inu pdp).
874 Relevant for the ray-tracing of the line profiles
875
876 (default: 100)
877 @type n_quad: int
878 @keyword fraction_tau_step: tau_total*fraction_tau_step gives min.
879 delta_tau in strahl.f. If too low,
880 min_tau_step will be used.
881
882 (default: 1e-2)
883 @type fraction_tau_step: float
884 @keyword min_tau_step: minimum of delta_tau in strahl.f
885
886 (default: 1e-4)
887 @type min_tau_step: float
888 @keyword write_intensities: set to 1 to write the intensities of first
889 50 impact-parameters at the end of sphinx
890
891 (default: 0)
892 @type write_intensities: bool
893 @keyword tau_max: maximum optical depth used for the calculation of the
894 formal integral
895
896 (default: 12.)
897 @type tau_max: float
898 @keyword tau_min: maximum optical depth used for the calculation of the
899 formal integral
900
901 (default: -6.)
902 @type tau_min: float
903 @keyword check_tau_step: check.par.in sphinx if step in tau not too
904 large
905
906 (default: 0.01)
907 @type check_tau_step: float
908 @keyword frequency: if not None the frequency of the transition is
909 taken to be this parameter in Hz, if None the
910 frequency of this transition is derived from the
911 GASTRoNOoM radiative input file
912
913 No indices are searched for if frequency is given
914 here. Usually used only for line listing.
915
916 (default: None)
917 @type frequency: float
918 @keyword exc_energy: the excitation energy (cm-1, from CDMS/JPL) if you
919 want it included in the Transition() object,
920 not mandatory!
921
922 (default: None)
923 @type exc_energy: float
924 @keyword int_intensity_log: the integrated intensity of the line in
925 logscale from CDMS/JPL, if you want it
926 included in the Transition() object,
927 not mandatory!
928
929 (default: None)
930 @type int_intensity_log: float
931 @keyword vibrational: In the case of line lists, this keyword indicates
932 the type of vibrational excitation, relevant fi
933 for H2O. Empty string if vibrational groundstate
934
935 (default: '')
936 @type vibrational: string
937 @keyword path_gastronoom: model output folder in the GASTRoNOoM home
938
939 (default: None)
940 @type path_gastronoom: string
941
942 '''
943
944 self.molecule = molecule
945
946
947 if telescope is None:
948 self.telescope = 'N.A.'
949 else:
950 telescope = telescope.upper()
951 if telescope.find('H2O') != -1 and telescope.find('PACS') != -1 \
952 and not self.molecule.isWater():
953 telescope = telescope.replace('-H2O','')
954 elif telescope.find('H2O') == -1 and telescope.find('PACS') != -1 \
955 and self.molecule.isWater():
956
957
958
959
960 telescope = '%s-H2O'%telescope
961 self.telescope = telescope
962 props = LPTools.readTelescopeProperties(self.telescope)
963 self.telescope_size = props[0]
964 self.tel_abs_err = props[1]
965
966
967
968 self.vup = int(vup)
969 self.jup = int(jup)
970 self.kaup = int(kaup)
971 self.kcup = int(kcup)
972 self.vlow = int(vlow)
973 self.jlow = int(jlow)
974 self.kalow = int(kalow)
975 self.kclow = int(kclow)
976 self.offset = float(offset)
977
978
979 self.n_quad = int(n_quad)
980
981
982 self.fraction_tau_step = fraction_tau_step
983 self.min_tau_step = min_tau_step
984 self.write_intensities = write_intensities
985 self.tau_max = tau_max
986 self.tau_min = tau_min
987 self.check_tau_step = check_tau_step
988
989 self.__model_id = None
990 if nup is None or nlow is None:
991 self.nup = self.kaup
992 self.nlow = self.kalow
993
994
995 self.exc_energy = exc_energy
996 self.int_intensity_log = int_intensity_log
997 self.vibrational = vibrational
998 self.sphinx = None
999 self.path_gastronoom = path_gastronoom
1000 self.unresolved = 'PACS' in self.telescope or 'SPIRE' in self.telescope
1001
1002
1003 self.datafiles = None
1004 self.lpdata = None
1005 self.fittedlprof = None
1006
1007 if frequency is None:
1008
1009 self.__setIndices()
1010 else:
1011 self.frequency = frequency
1012 self.c = 2.99792458e10
1013 self.h = 6.62606957e-27
1014 self.k = 1.3806488e-16
1015 self.wavelength = self.c/self.frequency
1016
1017
1018
1019
1020
1021 self.vlsr = None
1022
1023
1024 self.best_vlsr = None
1025 self.best_mtmb = None
1026 self.chi2_best_vlsr = None
1027
1028
1029 if self.unresolved:
1030 self.unreso = dict()
1031 self.unreso_err = dict()
1032 self.unreso_blends = dict()
1033 else:
1034 self.unreso = None
1035 self.unreso_err = None
1036 self.unreso_blends = None
1037
1038
1039
1041
1042 '''
1043 Printing a transition as it should appear in the GASTRoNOoM input file.
1044
1045 @return: The Transition string
1046 @rtype: string
1047
1048 '''
1049
1050 ll = [self.molecule.molecule_full,self.vup,self.jup,self.kaup,\
1051 self.kcup,self.vlow,self.jlow,self.kalow,self.kclow,\
1052 self.telescope,self.offset]
1053 tstr = 'TRANSITION={} {:d} {:d} {:d} {:d} {:d} {:d} {:d} {:d} {} {:.2f}'
1054 return tstr.format(*ll)
1055
1056
1057
1059
1060 '''
1061 Compare two transitions and return true if equal.
1062
1063 The condition is the string representation of this Transition(). Note
1064 that the other properties included in the self.makeDict() dictionary are
1065 not compared. Those properties do not determine whether a transition is
1066 equal or not.
1067
1068 In this sense equal refers to the quantum numbers, telescope and offset
1069 properties, ie the tags that go into the GASTRoNOoM model.
1070
1071 @return: The comparison
1072 @rtype: bool
1073
1074 '''
1075
1076 try:
1077 if str(self) == str(other):
1078 return True
1079 else:
1080 return False
1081 except AttributeError:
1082 return False
1083
1084
1085
1087
1088 '''
1089 Compare two transitions and return true if not equal.
1090
1091 The condition is the string representation of this Transition(). Note
1092 that the other properties included in the self.makeDict() dictionary are
1093 not compared. Those properties do not determine whether a transition is
1094 equal or not.
1095
1096 In this sense equal refers to the quantum numbers, telescope and offset
1097 properties, ie the tags that go into the GASTRoNOoM model.
1098
1099 @return: The negative comparison
1100 @rtype: bool
1101
1102 '''
1103
1104 try:
1105 if str(self) != str(other):
1106 return True
1107 else:
1108 return False
1109 except AttributeError:
1110 return True
1111
1112
1113
1115
1116 '''
1117 Return a hash number based on the string of the transition.
1118
1119 The condition is the string representation of this Transition(). Note
1120 that the other properties included in the self.makeDict() dictionary are
1121 not compared. Those properties do not determine whether a transition is
1122 equal or not.
1123
1124 In this sense equal refers to the quantum numbers, telescope and offset
1125 properties, ie the tags that go into the GASTRoNOoM model.
1126
1127 @return: The hash number:
1128 @rtype: int
1129
1130 '''
1131
1132 return hash(str(self))
1133
1134
1135
1165
1166
1168
1169 '''
1170 Make a string that suits telescope.spec inputfile.
1171
1172 @return: The Line spec string
1173 @rtype: string
1174
1175 '''
1176
1177 return 'LINE_SPEC=%s %i %i %i %i %i %i %i %i %.2f %.2f ! %s'\
1178 %(self.molecule.molecule_full,self.vup,self.jup,self.kaup,\
1179 self.kcup,self.vlow,self.jlow,self.kalow,self.kclow,\
1180 self.getBeamwidth(),self.getEfficiency(),\
1181 self.molecule.molecule)
1182
1183
1184
1186
1187 '''
1188 Calculate the beamwidth specific for the telescope.
1189
1190 @return: The beamwidth
1191 @rtype: float
1192
1193 '''
1194
1195
1196 if 'PACS' in self.telescope.upper():
1197 data = DataIO.readCols(os.path.join(cc.path.aux,\
1198 'Pacs_beamsize_v4.dat'))
1199 wav = data[0]/10000.
1200 beam = data[1]
1201
1202
1203 interp = interp1d(wav,beam)
1204 try:
1205 return interp(self.wavelength)
1206 except ValueError:
1207 print 'WARNING! The wavelength of %s falls '%str(self) +\
1208 'outside of the interpolation range when determining '+\
1209 'the beamwidth. Extrapolating linearly...'
1210 if self.wavelength < wav[0]:
1211 return Interpol.linInterpol(wav[:2],beam[:2],\
1212 self.wavelength)
1213 else:
1214 return Interpol.linInterpol(wav[-2:],beam[-2:],\
1215 self.wavelength)
1216
1217
1218
1219 telescope_diameter = self.telescope_size * 100
1220 return 1.22*self.wavelength/telescope_diameter*60.*60.*360./(2.*pi)
1221
1222
1223
1225
1226 '''
1227 Calculate telescope beam efficiency.
1228
1229 Telescope specific!
1230
1231 The beam efficiency is included in the .spec files for GASTRoNOoM.
1232
1233 This number however is not used in any calculations of the line profile
1234 and is included here for reference only.
1235
1236 The numbers currently are correct only for HIFI.
1237
1238 @return: The beam efficiency
1239 @rtype: float
1240
1241 '''
1242
1243 if self.unresolved:
1244
1245 return 0.60
1246 if self.telescope.find('HIFI') > -1:
1247 if self.frequency > 1116.2e9 and self.frequency < 1400.e9:
1248
1249 eff_mb0 = 0.66
1250 else:
1251 eff_mb0 = 0.76
1252 eff_f = 0.96
1253 sigma = 3.8e-4
1254 eff_mb = eff_mb0 * exp(-1*(4*pi*sigma/self.wavelength)**2)
1255 return eff_mb/eff_f
1256 else:
1257 return 1.00
1258
1259
1260
1262
1263 '''
1264 Return a dict with transition string, and other relevant parameters.
1265
1266 @keyword in_progress: add an extra dict entry "IN_PROGRESS" if the
1267 transition is still being calculated.
1268
1269 (default: 0)
1270 @type in_progress: bool
1271
1272 @return: The transition dictionary including all relevant, defining
1273 information
1274 @rtype: dict()
1275
1276 '''
1277
1278 dd = dict([('TRANSITION',str(self).replace('TRANSITION=','')),\
1279 ('N_QUAD',self.n_quad),\
1280 ('FRACTION_TAU_STEP',self.fraction_tau_step),\
1281 ('MIN_TAU_STEP',self.min_tau_step),\
1282 ('WRITE_INTENSITIES',self.write_intensities),\
1283 ('TAU_MAX',self.tau_max),('TAU_MIN',self.tau_min),\
1284 ('CHECK_TAU_STEP',self.check_tau_step)])
1285
1286 if int(in_progress):
1287 dd['IN_PROGRESS'] = 1
1288
1289 return dd
1290
1291
1292
1294
1295 '''
1296 Return a string in the sphinx filename format.
1297
1298 @keyword number: the number in the filename (sph*, sph1 or sph2, hence
1299 can be *, 1, 2)
1300
1301 (default: '*')
1302 @type number: string
1303 @keyword include_path: Include the full filepath.
1304
1305 (default: 0)
1306 @type include_path: bool
1307
1308 @return: The sphinx filename for this transition
1309 @rtype: string
1310
1311 '''
1312
1313 try:
1314 number = str(int(number))
1315 except ValueError:
1316 number = str(number)
1317 quantum_dict = dict()
1318 if not self.molecule.spec_indices:
1319 for i,attr in enumerate(['vup','jup','vlow','jlow']):
1320 quantum_dict[i] = (attr,getattr(self,attr))
1321 elif self.molecule.spec_indices == 2:
1322 for i,attr in enumerate(['vup','Jup','Nup','vlow','jlow','Nlow']):
1323 quantum_dict[i] = (attr,getattr(self,attr.lower()))
1324 elif self.molecule.spec_indices == 3:
1325 for i,attr in enumerate(['vup','jup','vlow','jlow']):
1326 quantum_dict[i] = (attr,getattr(self,attr))
1327 else:
1328 for i,attr in enumerate(['vup','jup','Kaup','Kcup',\
1329 'vlow','jlow','Kalow','Kclow']):
1330 quantum_dict[i] = (attr,getattr(self,attr.lower()))
1331
1332 fn = '_'.join(['sph%s%s'%(number,self.getModelId()),\
1333 self.molecule.molecule] + \
1334 ['%s%i'%(quantum_dict[k][0],quantum_dict[k][1])
1335 for k in sorted(quantum_dict.keys())] + \
1336 [self.telescope,'OFFSET%.2f.dat'%(self.offset)])
1337
1338 if include_path:
1339 fn = os.path.join(cc.path.gastronoom,self.path_gastronoom,'models',\
1340 self.getModelId(),fn)
1341
1342 return fn
1343
1344
1345
1347
1348 '''
1349 Get frequency of the transition from the radiat.dat files.
1350
1351 @return: frequency in Hz
1352 @rtype: float
1353
1354 '''
1355
1356 return self.frequency
1357
1358
1359
1361
1362 '''
1363 Return the energy level of the upper state of this transition.
1364
1365 @keyword unit: The unit of the returned values. Can be any valid units
1366 str from the astropy units module (energy), or the unit
1367 itself. 'cm-1' and 'cm^-1' are accepted as well.
1368
1369 (default: 1./u.cm)
1370 @type unit: string/unit
1371
1372 @return: energy level in chosen unit
1373 @rtype: float
1374
1375 '''
1376
1377 if self.molecule.radiat is None:
1378 print '%s_radiat.dat not found. Cannot find energy levels.'\
1379 %self.molecule.molecule
1380 return
1381
1382 return self.molecule.radiat.getLEnergy(index=self.lup,unit=unit)[0]
1383
1384
1385
1387
1388 '''
1389 Return the energy level of the lower state of this transition.
1390
1391 @keyword unit: The unit of the returned values. Can be any valid units
1392 str from the astropy units module (energy), or the unit
1393 itself. 'cm-1' and 'cm^-1' are accepted as well.
1394
1395 (default: 1./u.cm)
1396 @type unit: string/unit
1397
1398 @return: energy level in chosen unit
1399 @rtype: float
1400
1401 '''
1402
1403 if self.molecule.radiat is None:
1404 print '%s_radiat.dat not found. Cannot find energy levels.'\
1405 %self.molecule.molecule
1406 return
1407
1408 return self.molecule.radiat.getLEnergy(index=self.llow,unit=unit)[0]
1409
1410
1411
1413
1414 '''
1415 Set the transition index of this transition from the spectroscopy file.
1416
1417 The level indices for lower and upper level are also set.
1418
1419 '''
1420
1421
1422
1423
1424
1425 if not self.molecule.spec_indices:
1426 self.lup = self.jup + self.vup*self.molecule.ny_low + 1
1427 self.llow = self.jlow + self.vlow*self.molecule.ny_low + 1
1428 else:
1429 indices = self.molecule.radiat_indices
1430
1431 quantum_up = [q
1432 for i,q in enumerate([self.vup,self.jup,\
1433 self.kaup,self.kcup])
1434 if i<len(indices[0])-1]
1435 quantum_low = [q
1436 for i,q in enumerate([self.vlow,self.jlow,\
1437 self.kalow,self.kclow])
1438 if i<len(indices[0])-1]
1439
1440
1441
1442
1443 self.lup = indices[[i[1:]
1444 for i in indices].index(quantum_up)][0]
1445 self.llow = indices[[i[1:]
1446 for i in indices].index(quantum_low)][0]
1447
1448
1449
1450
1451 self.tindex = self.molecule.radiat.getTI(lup=self.lup,llow=self.llow)
1452 if not isinstance(self.tindex,int):
1453 msg = 'Something fishy is going on in Transition.py... '+\
1454 'non-unique or invalid transition indices for %s!'\
1455 %self.getInputString(include_nquad=0)
1456 raise(IndexError(msg))
1457
1458
1459 freq = self.molecule.radiat.getTFrequency(index=self.tindex,unit='Hz')
1460 self.frequency = freq
1461
1462
1463
1465
1466 '''
1467 Make a label suitable to be used on an axis in a plot.
1468
1469 At the moment always returns a label typical for integrated line
1470 strength.
1471
1472 @keyword include_mol: Include the molecule name in the label.
1473
1474 (default: 1)
1475 @type include_mol: bool
1476
1477 @return: The label
1478 @rtype: str
1479
1480 '''
1481
1482 t = self.makeLabel(inc_vib=0).replace('$','',2)
1483 if include_mol:
1484 m = self.molecule.molecule_plot.replace('$','',4)
1485 axislabel = 'I_{\mathrm{%s} (%s)}'%(m,t)
1486 else:
1487 axislabel = 'I_{%s}'%(t)
1488 return axislabel
1489
1490
1491
1492 - def makeLabel(self,inc_vib=1,return_vib=0):
1493
1494 '''
1495 Return a short-hand label for this particular transition.
1496
1497 These labels can be used for plot line identifications for instance.
1498
1499 If self.vibrational is not None, it always concerns a line list and is
1500 included as well.
1501
1502 @keyword inc_vib: Include the vibrational state in the label. Is always
1503 True is self.vibrational is not None.
1504
1505 (default: 1)
1506 @type inc_vib: bool
1507 @keyword return_vib: Only return a label for the vibrational state.
1508 Does not work if vibrational is not None.
1509
1510 (default: 0)
1511 @type return_vib: bool
1512
1513 @return: The transition label
1514 @rtype: string
1515
1516 '''
1517
1518 if self.vibrational:
1519 if not self.molecule.spec_indices:
1520 return '%s: %i,%i - %i,%i' \
1521 %(self.vibrational,self.vup,self.jup,self.vlow,\
1522 self.jlow)
1523 elif self.molecule.spec_indices == 2:
1524 return '%s: %i,%i$_{%i}$ - %i,%i$_{%i}$'\
1525 %(self.vibrational,self.vup,self.jup,self.nup,\
1526 self.vlow,self.jlow,self.nlow)
1527 elif self.molecule.spec_indices == 3:
1528 return '%s: %i$_{%i}$ - %i$_{%i}$' \
1529 %(self.vibrational,self.jup,self.nup,self.jlow,\
1530 self.nlow)
1531 else:
1532 return '%s: %i,%i$_{%i,%i}$ - %i,%i$_{%i,%i}$'\
1533 %(self.vibrational,self.vup,self.jup,self.kaup,\
1534 self.kcup,self.vlow,self.jlow,self.kalow,self.kclow)
1535 else:
1536 if not self.molecule.spec_indices:
1537 if ((self.vup == 0 and self.vlow ==0) or not inc_vib)\
1538 and not return_vib:
1539 return r'$J=%i - %i$' %(self.jup,self.jlow)
1540 elif return_vib:
1541 return r'$\nu=%i$' %(self.vup)
1542 else:
1543 return r'$\nu=%i$, $J=%i-%i$' \
1544 %(self.vup,self.jup,self.jlow)
1545 elif self.molecule.spec_indices == 1:
1546 ugly = r'J_{\mathrm{K}_\mathrm{a}, \mathrm{K}_\mathrm{c}}'
1547 if ((self.vup == 0 and self.vlow ==0) or not inc_vib) \
1548 and not return_vib:
1549 return r'$%s=%i_{%i,%i} - %i_{%i,%i}$'\
1550 %(ugly,self.jup,self.kaup,self.kcup,\
1551 self.jlow,self.kalow,self.kclow)
1552 elif return_vib:
1553 if self.vup == 0:
1554 return r'$\nu=0$'
1555 else:
1556 dvup = {1:2, 2:3}
1557 return r'$\nu_%i=1$'%(dvup[self.vup])
1558 else:
1559 dvup = {1:2, 2:3}
1560 return r'$\nu_%i=1$, $%s=%i_{%i,%i} - %i_{%i,%i}$'\
1561 %(dvup[self.vup],ugly,self.jup,self.kaup,\
1562 self.kcup,self.jlow,self.kalow,self.kclow)
1563 elif (self.molecule.spec_indices == 2 \
1564 or self.molecule.spec_indices == 3):
1565 ugly = r'J_{\mathrm{N}}'
1566 if ((self.vup == 0 and self.vlow ==0) or not inc_vib)\
1567 and not return_vib:
1568 return r'$%s=%i_{%i} - %i_{%i}$'\
1569 %(ugly,self.jup,self.nup,self.jlow,self.nlow)
1570 elif return_vib:
1571 return r'$\nu=%i$' %(self.vup)
1572 else:
1573 return r'$\nu=%i$, $%s=%i_{%i} - %i_{%i}$'\
1574 %(self.vup,ugly,self.jup,self.nup,\
1575 self.jlow,self.nlow)
1576
1577 else:
1578 return '%i,%i$_{%i,%i}$ - %i,%i$_{%i,%i}$'\
1579 %(self.vup,self.jup,self.kaup,self.kcup,\
1580 self.vlow,self.jlow,self.kalow,self.kclow)
1581
1582
1583
1585
1586 '''
1587 Set a model_id for the transition, identifying the model for SPHINX!
1588
1589 May or may not be equal to the molecule's id, depending on the input
1590 parameters.
1591
1592 @param model_id: The model id
1593 @type model_id: string
1594
1595 '''
1596
1597 self.__model_id = model_id
1598
1599
1600
1602
1603 '''
1604 Return the model_id associated with this transition.
1605
1606 None if not yet set.
1607
1608 Empty string if the model calculation failed.
1609
1610 @return: The model id of this transition
1611 @rtype: string
1612
1613 '''
1614
1615 return self.__model_id
1616
1617
1618
1620
1621 '''
1622 Is this a Molecule() object?
1623
1624 @return: Molecule()?
1625 @rtype: bool
1626
1627 '''
1628
1629 return False
1630
1631
1632
1634
1635 '''
1636 Is this a Transition() object?
1637
1638 @return: Transition()?
1639 @rtype: bool
1640
1641 '''
1642
1643 return True
1644
1645
1646
1648
1649 '''
1650 Return True if successfully calculated by sphinx, False otherwise.
1651
1652 @return: Finished calculation?
1653 @rtype: bool
1654
1655 '''
1656
1657 return self.__model_id
1658
1659
1660
1671
1672
1673
1675
1676 '''
1677 Reset the data read for this object. Datafiles will have to be added
1678 anew through addDataFile or setData.
1679
1680 '''
1681
1682 self.datafiles = None
1683 self.lpdata = None
1684 self.fittedlprof = None
1685 self.best_vlsr = None
1686 self.best_mtmb = None
1687 self.chi2_best_vlsr = None
1688
1689
1690
1692
1693 '''
1694 Add a datadict for this transition. Includes filenames as keys, and
1695 fit results as values (can be None, in which case the filename is
1696 excluded)
1697
1698 The filenames are saved in self.datafiles, the fitresults in
1699 self.fittedlprof.
1700
1701 If datadict has been given a valid value, the self.lpdata list is set
1702 back to None, so the datafiles can all be read again.
1703
1704 When called for a SPIRE or PACS transition, no datafiles are added.
1705
1706 The fit results include:
1707 The gas terminal velocity, its error, the soft parabola function and
1708 possibly the extra gaussian will be set. It is possible the SP is a
1709 gaussian instead, if a soft parabola did not work well (the no gamma
1710 parameter is included in the dictionary).
1711
1712 The fit results are taken from the radio database which has the option
1713 to (re-)run the LPTools.fitLP method. If no fit results are available
1714 there, no data can be used for this instance of Transition().
1715
1716 Typically, the entry db[star_name][transition] is what is given here.
1717
1718 @param datadict: the filenames and associated fit results
1719 @type datadict: dict()
1720
1721 @keyword replace: Replace data if they had already been added before.
1722
1723 (default: 0)
1724 @type replace: bool
1725
1726 '''
1727
1728
1729
1730 if not type(datadict) is types.DictType or not datadict:
1731 return
1732
1733
1734 if self.unresolved:
1735 return
1736
1737
1738 if replace:
1739 self.resetData()
1740
1741
1742
1743
1744 if self.datafiles is None:
1745 self.datafiles = []
1746 self.fittedlprof = []
1747
1748 for k in sorted(datadict.keys()):
1749 if not datadict[k] is None:
1750 if not os.path.split(k)[0]:
1751 self.datafiles.append(os.path.join(cc.path.dradio,k))
1752 else:
1753 self.datafiles.append(k)
1754 self.fittedlprof.append(datadict[k])
1755
1756
1757 if not self.datafiles:
1758 self.datafiles, self.fittedlprof = None,None
1759 print 'No data found for %s. If there should be data: '%str(self)+\
1760 'Did you run fitLP() in the Radio db?'
1761
1762
1763
1764 self.lpdata = None
1765
1766
1767
1769
1770 '''
1771 Read the datafiles associated with this transition if available.
1772
1773 '''
1774
1775 if self.unresolved:
1776 return
1777 if self.lpdata is None:
1778 if not self.datafiles is None:
1779 self.lpdata = []
1780 for idf,df in enumerate(self.datafiles):
1781 if df[-5:] == '.fits':
1782 lprof = FitsReader.FitsReader(fn=df)
1783 else:
1784 lprof = TxtReader.TxtReader(fn=df)
1785 lprof.setNoise(self.getVexp(idf))
1786 self.lpdata.append(lprof)
1787
1788
1789
1790 - def setData(self,trans,replace=0):
1791
1792 """
1793 Set data equal to the data in a different Transition object, but for
1794 the same transition. Does not erase data if they have already been set
1795 in this object, unless explicitly requested.
1796
1797 A check is ran to see if the transition in both objects really is the
1798 same. Sphinx model id may differ!
1799
1800 Avoids too much overhead when reading data from files.
1801
1802 The same is done for the fitresults from LPTools.fitLP() taken out of
1803 the Radio database.
1804
1805 @param trans: The other Transition() object, assumes the transition
1806 is the same, but possibly different sphinx models.
1807 @type trans: Transition()
1808
1809 @keyword replace: Replace data if they had already been added before.
1810
1811 (default: 0)
1812 @type replace: bool
1813
1814 """
1815
1816
1817 if replace:
1818 self.resetData()
1819
1820
1821 if self == trans:
1822 if self.datafiles is None:
1823 self.datafiles = trans.datafiles
1824 self.fittedlprof = trans.fittedlprof
1825 self.lpdata = trans.lpdata
1826
1827
1828
1830
1831 """
1832 Return the vlsr read from the fits file of a resolved-data object, or
1833 taken from the Star.dat file in case of unresolved data. Note that
1834 resolved data may also return vlsr from Star.dat if the vlsr in the
1835 data file is significantly different from the value in Star.dat.
1836
1837 This is taken from the first lpdata object available by default, but
1838 can be chosen through the index keyword..
1839
1840 Returns 0.0 if not an unresolved line, and there are no data available.
1841
1842 This is different from the getBestVlsr() method, which determines the
1843 best matching vlsr between data and sphinx, if both are available.
1844
1845 @keyword index: The data list index of the requested noise value
1846
1847 (default: 0)
1848 @type index: int
1849
1850 @return: the source velocity taken from the fits file OR Star.dat. 0
1851 if data are not available. In km/s!
1852 @rtype: float
1853
1854 """
1855
1856 self.readData()
1857 if self.lpdata:
1858 return self.lpdata[index].getVlsr()
1859 elif not self.unresolved and not self.lpdata:
1860 return 0.0
1861 else:
1862 return self.vlsr
1863
1864
1865
1867
1868 '''
1869 Return the noise value of the FIRST data object in this transition, if
1870 available.
1871
1872 Note that None is returned if no data are available.
1873
1874 A different index than the default allows access to the other data
1875 objects.
1876
1877 @keyword index: The data list index of the requested noise value
1878
1879 (default: 0)
1880 @type index: int
1881
1882 @return: The noise value
1883 @rtype: float
1884
1885 '''
1886
1887 self.readData()
1888 if self.lpdata:
1889 return self.lpdata[index].getNoise(vexp=self.getVexp(index=index))
1890 else:
1891 return None
1892
1893
1894
1896
1897 '''
1898 Get the gas terminal velocity as estimated from a line profile fitting
1899 routine for the FIRST data object.
1900 A different index than the default allows access to the other data
1901 objects.
1902
1903 @keyword index: The data list index of the requested noise value
1904
1905 (default: 0)
1906 @type index: int
1907
1908 @return: vexp
1909 @rtype: float
1910
1911 '''
1912
1913 if self.fittedlprof:
1914 return self.fittedlprof[index]['vexp']
1915 else:
1916 return None
1917
1918
1920
1921 """
1922 If self.best_vlsr is None, the best source velocity will be guessed by
1923 comparing chi^2 values between different values for the source velocity
1924
1925 May be different from input value vlsr, the original expected
1926 source velocity (from Star.dat)!
1927
1928 By default, based on the first dataset in self.lpdata. Note that
1929 multiple datasets might be included. If so, a different index can be
1930 given to do the scaling based on a different data object. Scaling is
1931 done for only one datase. Since different datasets are for the same
1932 telescope, the same line and (usually) the same conditions, the source
1933 velocity is not expected to be different for the different datasets
1934 anyway.
1935
1936 Method will attempt to call self.readData and readSphinx if either are
1937 not available. If data are still not available, the initial guess is
1938 returned. If data are available, but sphinx isn't, vlsr from the fits
1939 files is returned, and the initial guess is returned in case of txt
1940 files.
1941
1942 @keyword index: The data list index of the requested noise value
1943
1944 (default: 0)
1945 @type index: int
1946
1947 @return: the best guess vlsr, or the initial guess if no sphinx or data
1948 are available [will return vlsr included in fitsfiles if
1949 available].
1950 @rtype: float
1951
1952 """
1953
1954
1955
1956
1957 if not self.best_vlsr is None:
1958 return self.best_vlsr
1959
1960
1961
1962 self.readData()
1963
1964
1965
1966
1967 self.readSphinx()
1968 if self.unresolved or not self.lpdata or not self.sphinx:
1969 return self.getVlsr(index=index)
1970
1971
1972 noise = self.getNoise(index=index)
1973 dvel = self.lpdata[index].getVelocity()
1974 dtmb = self.lpdata[index].getFlux()
1975 mvel = self.sphinx.getVelocity()
1976 mtmb = self.sphinx.getLPTmb()
1977
1978
1979
1980
1981
1982 interpolator = interp1d(x=mvel+self.getVlsr(index=index),\
1983 y=mtmb,fill_value=0.0,bounds_error=False)
1984
1985
1986
1987
1988
1989
1990 res = dvel[1]-dvel[0]
1991 nstep = int(0.5*self.getVexp(index=index)/res+1)
1992
1993
1994 mtmb_grid = [interpolator(dvel+i*res) for i in range(-nstep,nstep+1)]
1995
1996
1997
1998
1999 vlsr_grid = [self.getVlsr(index)-i*res for i in range(-nstep,nstep+1)]
2000
2001
2002 chisquared = [bs.calcChiSquared(data=dtmb[dtmb>=-3*noise],\
2003 model=mtmbi[dtmb>=-3*noise],\
2004 noise=noise)
2005 for mtmbi in mtmb_grid]
2006
2007
2008
2009 imin = argmin(chisquared)
2010 self.chi2_best_vlsr = chisquared[imin]
2011 self.best_vlsr = vlsr_grid[imin]
2012
2013
2014 self.best_mtmb = mtmb_grid[imin]
2015
2016 return self.best_vlsr
2017
2018
2020
2021 """
2022 Calculate the integrated intrinsic intensity of the sphinx line profile
2023 in SI or cgs units. Velocity is converted to frequency before
2024 integration.
2025
2026 Returns None if no sphinx profile is available yet!
2027
2028 @keyword cont_subtract: Subtract the continuum value outside the line
2029 from the whole line profile.
2030
2031 (default: 1)
2032 @type cont_subtract: bool
2033
2034 @keyword units: The unit system in which the integrated intensity is
2035 returned. Can be 'si' or 'cgs'.
2036
2037 (default: 'si')
2038 @type units: string
2039
2040 @return: The integrated intrinsic intensity of the line profile in
2041 SI or cgs units. (W/m2 or erg/s/cm2)
2042 @rtype: float
2043
2044 """
2045
2046 units = units.lower()
2047 if self.sphinx is None:
2048 self.readSphinx()
2049 if self.sphinx is None:
2050 return
2051
2052 mvel = self.sphinx.getVelocityIntrinsic()*10**5
2053
2054 mint = self.sphinx.getLPIntrinsic(cont_subtract=cont_subtract)
2055
2056
2057
2058
2059
2060 vel_to_freq = u.doppler_radio(self.frequency*u.Hz)
2061 freqgrid = (mvel*u.cm/u.s).to(u.Hz,equivalencies=vel_to_freq).value
2062
2063
2064
2065
2066 intint_cgs = -1*trapz(x=freqgrid[isfinite(mint)],\
2067 y=mint[isfinite(mint)])
2068 intint_si = intint_cgs*10**-3
2069 if intint_cgs < 0:
2070 print 'WARNING! Negative integrated flux found for %s with id %s!'\
2071 %(str(self),self.getModelId())
2072
2073 if units == 'si':
2074 return intint_si
2075 else:
2076 return intint_cgs
2077
2078
2079
2081
2082 """
2083 Calculate the integrated convolved intensity of the sphinx line profile
2084 over velocity.
2085
2086 Returns None if no sphinx profile is available yet!
2087
2088 @keyword cont_subtract: Subtract the continuum value outside the line
2089 from the whole line profile.
2090
2091 (default: 1)
2092 @type cont_subtract: bool
2093
2094 @return: The integrated convolved intensity of the line profile in
2095 erg km/s/s/cm2
2096 @rtype: float
2097
2098 """
2099
2100 if self.sphinx is None:
2101 self.readSphinx()
2102 if self.sphinx is None:
2103 return
2104 mvel = self.sphinx.getVelocity()
2105 mcon = self.sphinx.getLPConvolved(cont_subtract=cont_subtract)
2106
2107 return trapz(x=mvel,y=mcon)
2108
2109
2110
2112
2113 """
2114 Calculate the integrated Tmb of the sphinx line profile over velocity.
2115
2116 Returns None if no sphinx profile is available yet!
2117
2118 @keyword cont_subtract: Subtract the continuum value outside the line
2119 from the whole line profile.
2120
2121 (default: 1)
2122 @type cont_subtract: bool
2123
2124 @return: The integrated model Tmb profile in K km/s
2125 @rtype: float
2126
2127 """
2128
2129 if self.sphinx is None:
2130 self.readSphinx()
2131 if self.sphinx is None:
2132 return
2133 mvel = self.sphinx.getVelocity()
2134 mtmb = self.sphinx.getLPTmb(cont_subtract=cont_subtract)
2135
2136 return trapz(x=mvel,y=mtmb)
2137
2138
2139
2141
2142 """
2143 Get the peak Tmb of the sphinx line profile.
2144
2145 Returns None if no sphinx profile is available yet.
2146
2147 Is equal to the mean of the 5 points around the center of the profile.
2148
2149 @return: The peak Tmb of the sphinx line profile
2150 @rtype: float
2151
2152 """
2153
2154 if self.sphinx is None:
2155 self.readSphinx()
2156 if self.sphinx is None:
2157 return
2158 mtmb = self.sphinx.getLPTmb(cont_subtract=1)
2159 imid = len(mtmb)/2
2160 return mean(mtmb[imid-2:imid+3])
2161
2162
2163
2165
2166 """
2167 Calculate the integrated Tmb of the data line profile over velocity.
2168
2169 Note that by default only the first of data profiles is used for this,
2170 if there are multiple profiles available for this transition. (i.e.
2171 multiple observations of the same transition with the same telescope)
2172
2173 A different index than the default allows access to the other data
2174 objects.
2175
2176 Makes use of the results from the LPTools.fitLP method taken from the
2177 Radio database upon adding datafiles. If no extra
2178 gaussian is used, the integrated data profile is returned. Otherwise,
2179 the soft parabola fit is integrated instead to avoid taking into
2180 account an absorption feature in the profile.
2181
2182 The fitted line profile can be forced to be used for the integrated line
2183 strength.
2184
2185 The uncertainty is also returned. Three options:
2186 - The default absolute flux calibration uncertainty from
2187 Telescope.dat
2188 - The above + the fitting uncertainty [TBI]
2189 - The abs flux cal uncertainty set in Radio.py + the fitting
2190 uncertainty [TBI]
2191 The fitting uncertainty is currently not yet implemented, nor the option
2192 to add the the flux calibration uncertainty to Radio.py. [TBI]
2193
2194 Returns None if no data are available.
2195
2196 CGS or SI units can be requested as well, where the profile is converted
2197 to Fnu before integration via Fnu = Tmb/(pi*tel_diameter**2/8/k_B).
2198
2199 This does not work for PACS or SPIRE data.
2200
2201 @keyword index: The data list index of the requested noise value
2202
2203 (default: 0)
2204 @type index: int
2205 @keyword use_fit: Force the use of the fitted line profile.
2206
2207 (default: 0)
2208 @type use_fit: bool
2209 @keyword units: The unit of the integrated line strength. Can be returned
2210 in K Km/s (tmb), erg/s/cm2 (cgs) or w/m2 (si).
2211
2212 (default: 'tmb')
2213 @type units: str
2214
2215 @return: The integrated line strength and the relative uncertainty in
2216 the requested units.
2217 @rtype: (float,float)
2218
2219 """
2220
2221 units = units.lower()
2222
2223
2224 self.readData()
2225
2226
2227 if self.fittedlprof is None:
2228 return (None, None)
2229
2230
2231 if self.fittedlprof[index].has_key('abs_err'):
2232 abs_err = self.fittedlprof[index]['abs_err']
2233 else:
2234 abs_err = self.tel_abs_err
2235
2236
2237 if not units in ['si','cgs']:
2238 if use_fit or not self.fittedlprof[index]['fitabs'] is None:
2239
2240
2241 return (self.fittedlprof[index]['fgintint'],abs_err)
2242 else:
2243
2244
2245 return (self.fittedlprof[index]['dintint'],abs_err)
2246
2247
2248
2249
2250
2251 dvel = self.lpdata[index].getVelocity()
2252 vlsr = self.lpdata[index].getVlsr()
2253 vexp = self.fittedlprof[index]['vexp']
2254
2255
2256 window = self.fittedlprof[index]['intwindow']
2257 keep = np.abs(dvel-vlsr)<=(window*vexp)
2258 dvel = dvel[keep]
2259
2260
2261 if use_fit or not self.fittedlprof[index]['fitabs'] is None:
2262 functype = self.fittedlprof[index]['fitprof'][0]
2263 pars = self.fittedlprof[index]['fitprof'][1]
2264 dflux = funclib.evaluate(functype,dvel,pars)
2265 else:
2266 dflux = self.lpdata[index].getFlux()
2267 dflux = dflux[keep]
2268 dvel = dvel[-np.isnan(dflux)]
2269 dflux = dflux[-np.isnan(dflux)]
2270
2271
2272 tmb = eq.Tmb(self.telescope_size)
2273 fcgs = (dflux*u.K).to(u.erg/u.s/u.cm/u.cm/u.Hz,equivalencies=tmb).value
2274 dop = u.doppler_radio(self.frequency*u.Hz)
2275 freq = (dvel*u.km/u.s).to(u.Hz,equivalencies=dop).value
2276
2277
2278
2279
2280
2281 fint = -1*trapz(y=fcgs,x=freq)
2282
2283
2284 if units == 'si':
2285 fint = fint*10**-3
2286
2287 return (fint,abs_err)
2288
2289
2290
2292
2293 """
2294 Get the peak Tmb of the data line profile.
2295
2296 Returns None if no sphinx profile or data are available yet.
2297
2298 Is equal to the mean of the 5 points in the data profile around the
2299 center of the sphinx profile (ie in the same velocity bin as
2300 getPeakTmbSphinx). Makes use of the best_vlsr, so that is ran first.
2301
2302 Note that by default only the first of data profiles is used for this,
2303 if there are multiple profiles available for this transition. (i.e.
2304 multiple observations of the same transition with the same telescope)
2305
2306 A different index than the default allows access to the other data
2307 objects.
2308
2309 This does not work for PACS or SPIRE data.
2310
2311 @keyword index: The data list index of the requested noise value
2312
2313 (default: 0)
2314 @type index: int
2315 @keyword method: The method applied: either 'mean' or 'fit'. Mean
2316 derives the peak value from the mean of the npoints
2317 flux points around the vlsr. Fit takes the central peak
2318 flux at from the fit.
2319
2320 (default: 'mean')
2321 @type method: str
2322 @keyword kwargs: Extra keywords passed on to the LPTools method for peak
2323 determination.
2324
2325 (default: {})
2326 @type kwargs: dict
2327
2328 @return: The peak Tmb of the sphinx line profile
2329 @rtype: float
2330
2331 """
2332
2333 method = method.lower()
2334 if method not in ['mean','fit']:
2335 method = 'mean'
2336
2337 self.readData()
2338 if self.fittedlprof is None:
2339 return None
2340
2341
2342
2343 if method == 'fit': return self.fittedlprof[index]['peak']
2344
2345
2346
2347
2348 return LPTools.getPeakLPData(lprof=self.lpdata[index],**kwargs)
2349
2350
2351
2353
2354 """
2355 Set the integrated intensity for this transition measured in given
2356 filename. (in SI units of W/m2)
2357
2358 A negative value is given for those lines suspected of being in a blend
2359 both by having at least 2 model transitions in a fitted line's breadth
2360 or by having a fitted_fwhm/intrinsic_fwhm of ~ 1.2 or more.
2361
2362 The vlsr is also passed through this function as that is only available
2363 from the data objects (in this case the Spire or Pacs class). For
2364 unresolved lines, vlsr is read from Star.dat.
2365
2366 @param fn: The data filename that contains the measured
2367 integrated intensity.
2368 @type fn: string
2369 @param dint: The value for the integrated intensity in W/m2. If the
2370 line is part of a blend that has already been added, this
2371 may also say 'inblend'. All transitions involved in the
2372 blend are then given by st_blends.
2373 @type dint: float or string
2374 @param dint_err: The fitting uncertainty on the intensity + absolute
2375 flux calibration uncertainty of 20%. Relative value!
2376 @type dint_err: float
2377 @param vlsr: The source velocity with respect to the local standard of
2378 rest in cm/s.
2379 @type vlsr: float
2380
2381 @keyword st_blends: List of sample transitions involved in a line blend
2382 detected by finding multiple central wavs in a
2383 wavelength resolution bin
2384
2385 (default: None)
2386 @type st_blends: list[Transition()]
2387
2388 """
2389
2390 if dint == 'inblend':
2391 self.unreso[fn] = dint
2392 else:
2393 self.unreso[fn] = float(dint)
2394 self.unreso_err[fn] = not dint_err is None and float(dint_err) or None
2395 self.unreso_blends[fn] = st_blends
2396
2397 self.vlsr = vlsr * 10**(-5)
2398
2399
2400
2402
2403 '''
2404 If already set, the integrated intensity can be accessed here based on
2405 filename (multiple measurements can be available for a single
2406 transition).
2407
2408 Always set and returned in W/m2!
2409
2410 Must have been set through setIntIntUnresolved()! Otherwise returns
2411 None.
2412
2413 A negative value is given for those lines suspected of being in a blend
2414 both by having at least 2 model transitions in a fitted line's breadth
2415 or by having a fitted_fwhm/intrinsic_fwhm of ~ 1.2 or more.
2416
2417 @keyword fn: The data filename that contains the measured integrated
2418 intensity. Can be set to '' or None if simply the
2419 first entry in the keys() list is to be used. Mostly only
2420 one line strength is associated with the object anyway.
2421
2422 (default: '')
2423 @type fn: string
2424
2425 @return: The integrated intensity measured in unresolved data for this
2426 filename, in SI units of W/m2, and the fitting uncertainty +
2427 absolute flux calibration uncertainty (relative value!), and
2428 the blends if applicable.
2429 @rtype: (float,float,list)
2430
2431 '''
2432
2433 if self.unreso.has_key(fn):
2434 return (self.unreso[fn],\
2435 self.unreso_err[fn],\
2436 self.unreso_blends[fn])
2437 elif not fn and self.unreso.keys():
2438 k = sorted(self.unreso.keys())[0]
2439 return (self.unreso[k],\
2440 self.unreso_err[k],\
2441 self.unreso_blends[k])
2442 else:
2443 return (None,None,None)
2444
2445
2446
2447 - def getLoglikelihood(self,use_bestvlsr=1,index=0,normalise=1,\
2448 vmin=0.0,vmax=0.0,use_fit=0):
2449
2450 """
2451 Calculate the loglikelihood of comparison between sphinx and dataset.
2452
2453 Gives a measure for the goodness of the fit of the SHAPE of the
2454 profiles.
2455
2456 Note that by default only the first of data profiles is used for this,
2457 if there are multiple profiles available for this transition. (i.e.
2458 multiple observations of the same transition with the same telescope)
2459
2460 A different index than the default allows access to the other data
2461 objects.
2462
2463 Done for the dataset with given index! Makes use of the interpolated
2464 sphinx profile for the best vlsr, see self.getBestVlsr() if use_bestvlsr
2465 is True. If this keyword is False, interpolates the sphinx model for the
2466 vlsr from Star.dat or the fits file.
2467
2468 Returns None if sphinx or data profile are not available.
2469
2470 Rescales the sphinx profile according to the difference in integrated
2471 Tmb between dataset and sphinx profile.
2472
2473 @keyword use_bestvlsr: Use the fitted best-guess for the v_lsr when
2474 determining the velocity grid for the model. If
2475 not, the vlsr from the Star.dat file or the fits
2476 file is used.
2477
2478 (default: 1)
2479 @type use_bestvlsr: bool
2480 @keyword index: The data list index of the requested noise value
2481
2482 (default: 0)
2483 @type index: int
2484 @keyword normalise: Normalise the data and model lines to the integrated
2485 line strength of the observed line.
2486
2487 (default: 1)
2488 @type normalise: bool
2489 @keyword vmin: The minimum value in km/s of the spectral line window.
2490 Ideally this is the same for all lines under
2491 consideration. If an invalid spectral window is given
2492 (vmax==vmin, or vmin>vmax), the spectral window is taken
2493 from the line fit results. This leads to a different
2494 window for each line under consideration, and is not
2495 recommended for calculating the loglikelihood.
2496
2497 (default: 0.0)
2498 @type vmin: float
2499 @keyword vmax: The maximum value in km/s of the spectral line window.
2500 Ideally this is the same for all lines under
2501 consideration. If an invalid spectral window is given
2502 (vmax==vmin, or vmin>vmax), the spectral window is taken
2503 from the line fit results. This leads to a different
2504 window for each line under consideration, and is not
2505 recommended for calculating the loglikelihood.
2506
2507 (default: 0.0)
2508 @type vmax: float
2509 @keyword use_fit: Force the use of the fitted line profile.
2510
2511 (default: 0)
2512 @type use_fit: bool
2513
2514 @return: The loglikelihood
2515 @rtype: float
2516
2517 """
2518
2519 if use_bestvlsr and self.best_vlsr is None:
2520 self.getBestVlsr(index=index)
2521 if use_bestvlsr and self.best_vlsr is None:
2522 print 'Using standard v_lsr from Star.dat or fits file for LLL.'
2523 use_bestvlsr = 0
2524
2525 vel = self.lpdata[index].getVelocity()
2526 noise = self.getNoise(index=index)
2527 window = self.fittedlprof[index]['intwindow']
2528 vexp = self.getVexp(index=index)
2529 vlsr = self.getVlsr(index=index)
2530
2531
2532
2533 if vmin != vmax and vmin < vmax:
2534 selection = np.logical_and(vel>=vmin, vel<=vmax)
2535
2536
2537
2538 else:
2539 print "WARNING: Invalid window for loglikelihood. Check vmin/vmax."
2540 selection = abs(vel-vlsr)<=window*vexp
2541
2542
2543
2544
2545
2546
2547
2548
2549 dsel = self.lpdata[index].getFlux()[selection]
2550
2551
2552
2553
2554
2555
2556 line_strength = self.getIntTmbData(index=index,use_fit=use_fit)[0]
2557
2558
2559
2560 if normalise:
2561 norm_factor = 1./line_strength
2562 else:
2563 norm_factor = 1.
2564 dsel = dsel*norm_factor
2565
2566
2567
2568
2569
2570
2571
2572
2573 shift_factor = line_strength/self.getIntTmbSphinx()
2574
2575
2576 if use_bestvlsr:
2577 msel = self.best_mtmb[selection]
2578 msel = msel*shift_factor
2579 else:
2580 mvel = self.sphinx.getVelocity()
2581 mtmb = self.sphinx.getLPTmb()
2582 interpolator = interp1d(x=mvel+self.getVlsr(index=index),y=mtmb,\
2583 fill_value=0.0,bounds_error=False)
2584 mtmb_interp = interpolator(vel[selection])
2585 msel = mtmb_interp*shift_factor
2586
2587
2588 msel = msel*norm_factor
2589
2590 return bs.calcLoglikelihood(data=dsel,model=msel,noise=noise)
2591