Source code for pyspeckit.spectrum.models.hcn

HCN Hyperfine Fitter
This is an HCN fitter...
ref for line params:

Module API

from __future__ import print_function
import numpy as np
from six import iteritems
import matplotlib.cbook as mpcb
import copy
    from import fits as pyfits
except ImportError:
    import pyfits
    import scipy.interpolate
    import scipy.ndimage
    scipyOK = True
except ImportError:

from .. import units
from . import fitter,model,modelgrid
from . import hyperfine

freq_dict={# splatalogue frequencies
           # frequencies from Ahrens 2002 & Loughnane 2013
           '11-01':88.63041376e9, # F=1-1
           '12-01':88.63184666e9, # F=2-1
           '10-01':88.63393544e9, # F=0-1
aval_dict = { # not used
Line strengths of the 15 hyperfine components in J = 1 -  0  transition. The
thickness of the lines indicates their relative weight compared to the others.
Line strengths are normalized in such a way that summing over all initial J = 1
levels gives the degeneracy of the  J = 0 levels, i.e.,  for JF1F = 012,
three for JF1F = 011, and one for JF1F = 010. Thus, the sum over all 15
transitions gives the total spin degeneracy
line_strength_dict = { # effectively the degeneracy per rotation state...

relative_strength_total_degeneracy = {

line_names = freq_dict.keys()

ckms = units.speedoflight_ms / 1e3 #2.99792458e5
voff_lines_dict = dict([(k,(v-88.6318470e9)/88.6318470e9*ckms) for k,v in iteritems(freq_dict)])

hcn_vtau = hyperfine.hyperfinemodel(line_names, voff_lines_dict, freq_dict, line_strength_dict, relative_strength_total_degeneracy)
hcn_amp = hcn_vtau.ampfitter
hcn_vtau_fitter = hcn_vtau.fitter
hcn_vtau_vheight_fitter = hcn_vtau.vheight_fitter
hcn_varyhf_fitter = hcn_vtau.varyhf_fitter 
hcn_varyhf_amp_fitter = hcn_vtau.varyhf_amp_fitter 
hcn_varyhf_amp_width_fitter = hcn_vtau.varyhf_amp_width_fitter 

[docs]def hcn_radex(xarr, density=4, column=13, xoff_v=0.0, width=1.0, grid_vwidth=1.0, grid_vwidth_scale=False, texgrid=None, taugrid=None, hdr=None, path_to_texgrid='', path_to_taugrid='', temperature_gridnumber=3, debug=False, verbose=False, **kwargs): """ Use a grid of RADEX-computed models to make a model line spectrum The RADEX models have to be available somewhere. OR they can be passed as arrays. If as arrays, the form should be: texgrid = ((minfreq1,maxfreq1,texgrid1),(minfreq2,maxfreq2,texgrid2)) xarr must be a SpectroscopicAxis instance xoff_v, width are both in km/s grid_vwidth is the velocity assumed when computing the grid in km/s this is important because tau = modeltau / width (see, e.g., Draine 2011 textbook pgs 219-230) grid_vwidth_scale is True or False: False for LVG, True for Sphere """ if texgrid is None and taugrid is None: if path_to_texgrid == '' or path_to_taugrid=='': raise IOError("Must specify model grids to use.") else: taugrid = [pyfits.getdata(path_to_taugrid)] texgrid = [pyfits.getdata(path_to_texgrid)] hdr = pyfits.getheader(path_to_taugrid) yinds,xinds = np.indices(taugrid[0].shape[1:]) densityarr = (xinds+hdr['CRPIX1']-1)*hdr['CD1_1']+hdr['CRVAL1'] # log density columnarr = (yinds+hdr['CRPIX2']-1)*hdr['CD2_2']+hdr['CRVAL2'] # log column minfreq = (4.8,) maxfreq = (5.0,) elif len(taugrid)==len(texgrid) and hdr is not None: minfreq,maxfreq,texgrid = zip(*texgrid) minfreq,maxfreq,taugrid = zip(*taugrid) yinds,xinds = np.indices(taugrid[0].shape[1:]) densityarr = (xinds+hdr['CRPIX1']-1)*hdr['CD1_1']+hdr['CRVAL1'] # log density columnarr = (yinds+hdr['CRPIX2']-1)*hdr['CD2_2']+hdr['CRVAL2'] # log column else: raise Exception # Convert X-units to frequency in GHz xarr = copy.copy(xarr) xarr.convert_to_unit('Hz', quiet=True) tau_nu_cumul = np.zeros(len(xarr)) gridval1 = np.interp(density, densityarr[0,:], xinds[0,:]) gridval2 = np.interp(column, columnarr[:,0], yinds[:,0]) if np.isnan(gridval1) or np.isnan(gridval2): raise ValueError("Invalid column/density") if scipyOK: tau = [scipy.ndimage.map_coordinates(tg[temperature_gridnumber,:,:],np.array([[gridval2],[gridval1]]),order=1) for tg in taugrid] tex = [scipy.ndimage.map_coordinates(tg[temperature_gridnumber,:,:],np.array([[gridval2],[gridval1]]),order=1) for tg in texgrid] else: raise ImportError("Couldn't import scipy, therefore cannot interpolate") #tau = modelgrid.line_params_2D(gridval1,gridval2,densityarr,columnarr,taugrid[temperature_gridnumber,:,:]) #tex = modelgrid.line_params_2D(gridval1,gridval2,densityarr,columnarr,texgrid[temperature_gridnumber,:,:]) if verbose: print("density %20.12g column %20.12g: tau %20.12g tex %20.12g" % (density, column, tau, tex)) if debug: import pdb; pdb.set_trace() return hcn_vtau(xarr,Tex=tex,tau=tau,xoff_v=xoff_v,width=width,**kwargs)