LTE Molecule Model¶
There is a tool for modeling the full LTE spectrum of a molecule. It uses the JPL or CDMS databases through their astroquery interfaces.
A very simple example looks like this:
import astropy.units as u
from pyspeckit.spectrum.models.lte_molecule import get_molecular_parameters, generate_fitter, generate_model
freqs, aij, deg, EU, partfunc = get_molecular_parameters('CH3OH',
fmin=90*u.GHz,
fmax=100*u.GHz)
def modfunc(xarr, vcen, width, tex, column):
return generate_model(xarr, vcen, width, tex, column, freqs=freqs, aij=aij,
deg=deg, EU=EU, partfunc=partfunc)
fitter = generate_fitter(modfunc, name="CH3OH")
Details can be found in the API documentation:
LTE Molecule Modeling Tool¶
Uses astroquery to obtain molecular parameters.
Equations are based on Mangum & Shirley 2015 (2015PASP..127..266M)
Module API¶
- pyspeckit.spectrum.models.lte_molecule.Jnu(nu, T)[source] [github] [bitbucket]¶
RJ equivalent temperature (MS15 eqn 24)
- pyspeckit.spectrum.models.lte_molecule.Jnu_cgs(nu, T)[source] [github] [bitbucket]¶
RJ equivalent temperature (MS15 eqn 24) (use cgs constants for speed)
- pyspeckit.spectrum.models.lte_molecule.generate_fitter(model_func, name)[source] [github] [bitbucket]¶
Generator for hnco fitter class
- pyspeckit.spectrum.models.lte_molecule.generate_model(xarr, vcen, width, tex, column, freqs, aij, deg, EU, partfunc, background=None, tbg=2.73, get_tau=False, get_tau_sticks=False)[source] [github] [bitbucket]¶
Model Generator
- Parameters:
xarr : Quantity array [Hz]
The X-axis (frequency)
vcen : Quantity [km/s]
The central velocity
width : Quantity [km/s]
The 1-sigma Gaussian line width [not FWHM]
tex : Quantity [K]
Excitation temperature
column : Quantity [cm^-2]
The total column density of the molecule. If specified as a value < 30, it will be assumed to be the log10 of the column density.
freqs : Quantity array [Hz]
The central frequency of the lines to model. Subsequent parameters, aij, deg, EU, must have the same shape as
freqs.aij : array [log s^-1]
The Einstein A coefficients, assumed to be in units of log10 of the rate in 1/s
deg : array
The transition degeneracy
EU : Quantity array [erg]
The upper state energy in ergs
partfunc : function
The partition function. Can also be specified as an array of values in Kelvin.
background : None
An additional background field to include in the radiative transfer model. Must be in Kelvin. [Presently not correctly implemented]
tbg : Quantity [K]
The background brightness temperature, defaulting to 2.73 for the CMB temperature. This background will be treated as uniform with frequency and subtracted.
get_tau : bool, default False
If specified, the optical depth in each frequency bin will be returned
get_tau_sticks : bool, default False
If specified, the optical depth in each transition will be returned
Examples
Example case to produce a model:
from pyspeckit.spectrum.models.lte_molecule import get_molecular_parameters, generate_model, generate_fitter freqs, aij, deg, eu, partfunc = get_molecular_parameters('CH3OH') def modfunc(xarr, vcen, width, tex, column): return generate_model(xarr, vcen, width, tex, column, freqs=freqs, aij=aij, deg=deg, EU=EU, partfunc=partfunc) fitter = generate_fitter(modfunc, name="CH3OH")
- pyspeckit.spectrum.models.lte_molecule.get_molecular_parameters(molecule_name, tex=50, fmin=<Quantity 1. GHz>, fmax=<Quantity 1. THz>, catalog='JPL', molecule_tag=None, parse_name_locally=True, flags=0, return_table=False, **kwargs)[source] [github] [bitbucket]¶
Get the molecular parameters for a molecule from the JPL or CDMS catalog
(this version should, in principle, be entirely self-consistent)
- Parameters:
molecule_name : string
The string name of the molecule (normal name, like CH3OH or CH3CH2OH, but it has to match the JPL catalog spec). Will use partial string matching if the full string name fails.
tex : float
Optional excitation temperature (basically checks if the partition function calculator works)
catalog : ‘JPL’ or ‘CDMS’
Which catalog to pull from
fmin : quantity with frequency units
fmax : quantity with frequency units
The minimum and maximum frequency to search over
molecule_tag : int, optional
If specified, this will override the molecule name. You can specify molecules based on the ‘TAG’ column in the JPL table
parse_name_locally : bool
Option passed to the query tool to specify whether to use regex to search for the molecule name
flags : int
Regular expression flags to pass to the regex search
return_table : bool
Also return the parameter table?
Examples
>>> from pyspeckit.spectrum.models.lte_molecule import get_molecular_parameters >>> freqs, aij, deg, EU, partfunc = get_molecular_parameters('CH2CHCN', ... fmin=220*u.GHz, ... fmax=222*u.GHz, ) >>> freqs, aij, deg, EU, partfunc = get_molecular_parameters('CH3OH', ... fmin=90*u.GHz, ... fmax=100*u.GHz)
- pyspeckit.spectrum.models.lte_molecule.line_tau(tex, total_column, partition_function, degeneracy, frequency, energy_upper, einstein_A=None)[source] [github] [bitbucket]¶
Given the excitation temperature of the state, total column density of the molecule, the partition function, the degeneracy of the state, the frequency of the state, and the upper-state energy level, return the optical depth of that transition.
This is a helper function for the LTE molecule calculations. It implements the equations
\[\int \tau_\nu d\nu = \frac{c^2}{8 \pi \nu^2} A_{ij} N_u \left[ \exp\left( \frac{h \nu}{k_B T_{ex}} \right) - 1 \right]\]or
\[\tau_\nu = \frac{c^2}{8 \pi \nu^2} A_{ij} N_u \phi_\nu \left[ \exp\left( \frac{h \nu}{k_B T_{ex}} \right) - 1 \right]\]where
\[N_{u} = N_{tot} \frac{g_u}{Q} \exp\left(\frac{-E_u}{k_B T_{ex}} \right)\]based on Equations 11 and 29 of Mangum & Shirley 2015 (2015PASP..127..266M)
The return value is therefore
\[\tau_\nu / \phi_\nu\]The line profile function is, for a Gaussian, given by eqn A1:
\[\phi_\nu = \frac{1}{\sqrt{2 \pi} \sigma} \exp \left[ -\frac{(\nu-\nu_0)^2}{2 \sigma^2} \right]\]
- pyspeckit.spectrum.models.lte_molecule.line_tau_cgs(tex, total_column, partition_function, degeneracy, frequency, energy_upper, einstein_A)[source] [github] [bitbucket]¶
Given the excitation temperature of the state, total column density of the molecule, the partition function, the degeneracy of the state, the frequency of the state, and the upper-state energy level, return the optical depth of that transition.
Unlike
line_tau(), this function requires inputs in CGS units.This is a helper function for the LTE molecule calculations. It implements the equations
\[\int \tau_\nu d\nu = \frac{c^2}{8 \pi \nu^2} A_{ij} N_u \left[ \exp\left( \frac{h \nu}{k_B T_{ex}} \right) - 1 \right]\]or
\[\tau_\nu = \frac{c^2}{8 \pi \nu^2} A_{ij} N_u \phi_\nu \left[ \exp\left( \frac{h \nu}{k_B T_{ex}} \right) - 1 \right]\]where
\[N_{u} = N_{tot} \frac{g_u}{Q} \exp\left(\frac{-E_u}{k_B T_{ex}} \right)\]based on Equations 11 and 29 of Mangum & Shirley 2015 (2015PASP..127..266M)
The return value is therefore
\[\tau_\nu / \phi_\nu\]The line profile function is, for a Gaussian, given by eqn A1:
\[\phi_\nu = \frac{1}{\sqrt{2 \pi} \sigma} \exp \left[ -\frac{(\nu-\nu_0)^2}{2 \sigma^2} \right]\]
- pyspeckit.spectrum.models.lte_molecule.ntot_of_nupper(nupper, eupper, tex, Q_rot, degeneracy=1)[source] [github] [bitbucket]¶
Given an N_upper, E_upper, tex, Q_rot, and degeneracy for a single state, give N_tot
Mangum & Shirley 2015 eqn 31
\[N_{tot}/N_u = Q_{rot} / g_u \exp(E_u/k T_{ex})\]- Example:
>>> import astropy.units as u >>> tex = 50*u.K >>> kkms = 100*u.K*u.km/u.s >>> from pyspeckit.spectrum.models import lte_molecule >>> freqs, aij, deg, EU, partfunc = lte_molecule.get_molecular_parameters(molecule_name='HNCO v=0', molecule_name_jpl='HNCO', fmin=87*u.GHz, fmax=88*u.GHz) >>> nupper = lte_molecule.nupper_of_kkms(kkms, freqs, 10**aij) >>> ntot = lte_molecule.ntot_of_nupper(nupper, EU*u.erg, tex, Q_rot=partfunc(tex), degeneracy=deg)
- pyspeckit.spectrum.models.lte_molecule.nupper_of_kkms(kkms, freq, Aul, replace_bad=None)[source] [github] [bitbucket]¶
Mangum & Shirley 2015 eqn 82 gives, for the optically thin, Rayleigh-Jeans, negligible background approximation:
\[N_{tot} = (3 k) / (8 \pi^3 \nu S \mu^2 R_i) (Q/g) \exp(E_u/k T_{ex}) \int(T_R/f dv)\]Eqn 31:
\[ \begin{align}\begin{aligned}N_{tot}/N_u = Q_{rot} / g_u \exp(E_u/k T_{ex})\\-> N_{tot} = N_u Q_{rot} / g_u \exp(E_u/k T_{ex}) -> Nu = N_{tot} g / Q_{rot} \exp(-E_u / k T_{ex})\end{aligned}\end{align} \]To get Nu of an observed line, then:
\[Nu Q_{rot} / g_u \exp(E_u/k T_{ex}) = (3 k) / (8 \pi^3 \nu S \mu^2 R_i) (Q/g) \exp(E_u/k T_{ex}) \int(T_R/f dv)\]This term cancels:
\[Q_{rot} / g_u \exp(E_u/k T_{ex})\]Leaving:
\[N_u = (3 k) / (8 \pi^3 \nu S \mu^2 R_i) \int(T_R/f dv)\]$int(T_R/f dv)$ is the optically thin integrated intensity in K km/s dnu/nu = dv/c [doppler eqn], so to get $int(T_R dnu)$, sub in $dv = c/nu dnu$
\[N_u = (3 k c) / (8 \pi^3 \nu^2 S \mu^2 R_i) \int(T_R/f d\nu)\]We then need to deal with the S mu^2 R_i term. We assume R_i = 1, since we are not measuring any hyperfine transitions (R_i is the hyperfine degeneracy; eqn 75) Equation 11:
\[A_{ul} = 64 \pi^4 \nu^3 / (3 h c^3) |\mu_{ul}|^2\]Equation 62:
\[ \begin{align}\begin{aligned}|\mu_{ul}|^2 = S \mu^2\\-> S \mu^2 = (3 h c^3 A_{ul}) / (64 \pi^4 \nu^3)\end{aligned}\end{align} \]Plugging that in gives
\[Nu = (3 k c) / (8 \pi^3 \nu^2 ((3 h c^3 A_{ul}) / (64 \pi^4 \nu^3))) \int(T_R/f d\nu) = (3 k c 64 \pi^4 \nu^3) / (8 \pi^3 \nu^2 3 h c^3 A_{ul}) \int(T_R/f d\nu) = (8 \pi \nu k / (A_{ul} c^2 h)) \int(T_R/f d\nu)\]which is the equation implemented below. We could also have left this in dv units by substituting du = nu/c dv:
\[= (8 \pi \nu^2 k / (A_{ul} c^3 h)) \int (T_R/f dv)\]