API

βlvg(τ)

Expanding sphere, Large Velocity Gradient, or Sobolev case. Formula from De Jong, Boland, and Dalgarno (1980, A&A 91, 68)^[DeJong80]. Corrected by a factor of 2 in order to return 1 for τ=0. See Appendix B equation B-7.

βslab(τ)

Slab geometry (e.g. shocks) with power law approximations from De Jong, Dalgarno, and Chu (1975, ApJ 199, 69)^[DeJong75].

βsphere(τ)

Uniform sphere formula from Osterbrock (Astrophysics of Gaseous Nebulae and Active Galactic Nuclei) Appendix 2 with power law approximations for large and small optical depth.

parse_colliders(lines, levels, transitions)

Collider collision rate terms are of the form:

!NUMBER OF COLL PARTNERS
1
!COLLISIONS BETWEEN
1  HCO+ - H2 from Flower (1999)
!NUMBER OF COLL TRANS
210
!NUMBER OF COLL TEMPS
12
!COLL TEMPS
   10.0   20.0   30.0   50.0   70.0  100.0  150.0  200.0  250.0  300.0  350.0  400.0
!TRANS + UP + LOW + COLLRATES(cm^3 s^-1)
    1     2     1  2.6e-10 2.3e-10 2.1e-10 2.0e-10 1.9e-10 1.8e-10 2.0e-10 2.2e-10 2.3e-10 2.5e-10 2.7e-10 2.8e-10

parse_energies(lines)

Energy level terms begin on the 8th line and are of the form:

!LEVEL + ENERGIES(cm^-1) + WEIGHT + J
   1     0.000000000   1.0     0

parse_table(lines, types)

Parse a white-space delimited, fixed-width table into a set of columns with known data types.

parse_transitions(lines, levels::EnergyLevels)

Transition terms are of the form:

!TRANS + UP + LOW + EINSTEINA(s^-1) + FREQ(GHz) + E_u(K)
    1     2     1  4.251e-05          89.18839570     4.28

Specie(name::AbstractString; datadir=datadir)

Parse the LAMDA collision rate file for an atomic or molecular specie from the filename stem (i.e., without the ".dat" file extension). For example, the stem for the "hco+@xpol.dat" file would be "hco+@xpol".

blackbody_background(xnu; tbg::Real=TCMB)

Set the background radiation field from a blackbody with a given background temperature. The default temperature is set for the CMB.

calc_galactic_isrf(xnu; uv_scaling=1.0)

Galactic interstellar radiation field (ISRF) as parametrized by Hocuk et al. (2017; "H17") in Appendix B based on Zucconi et al. (2001). Six weighted black-body components are used corresponding to the values in H17 Table B1. The UV component of the ISRF is based on Draine (1978) and added according to the polynomial in H17 Equation B2.

The angularly-integrated intensity $J_ν$ is returned with units of erg / (s cm^2 sr Hz).

planck(ΔE, tex)

Planck function to compute black body intensity in [erg / (s cm^2 Hz sr)].

Arguments

• ΔE: Energy level difference in wavenumber [1/cm].
• tex: Excitation temperature in [K].

BackgroundField(trj, backi, totalb)

Arguments

• trj: Background brightness temperature [K] in the Rayleigh-Jeans limit.
• backi: Background intensity (integrated, $J_ν$) in [erg / (s cm^2 Hz sr)].
• totalb: Total background intensity. Note that because internal radiation fields are not implemented, totalb should be identical to backi.

get_collision_rates(mol::Specie, densities::Dict, tkin::Real)

Interpolate the collision rates for the kinetic temperature and sum over all colliders.

Returns

• crate::Matrix: Collision rate matrix (product of density and rate coeff.).
• ctot::Vector: Total collision rate for all transitions to a level.

thermal_h2_density(density, tkin)

Calculate densities for the ortho (symmetric) and para (anti-symmetric) spin-isomers of molecular hydrogen (H2) from a total molecule density. The ortho-to-para ratio assumes that the H2 was formed in thermal equillibrium at the give kinetic temperature. In the high-temperature limit (> 100 K) the ortho-to-para ratio is 3.

calc_radiation_temperature(xnu, tex, τl, intens_bg)

Calculate the radiation peak temperature or the background-subtracted line intensity in units of the equivalent radiation temperature in the Rayleigh-Jeans limit. This is equivalent to an observed "main beam" antenna temperature $T_\mathrm{mb}$ (antenna temperature corrected for the aperture efficiency, spill-over, and atmospheric opacity) if the emission completely fills the beam (i.e., a beam filling fraction of 1).

\[T_R = \frac{c^2}{2 k \nu^2} \left( I^\mathrm{em}_\nu - I^\mathrm{bg}_\nu \right)\]

On the first iteration, use the background intensity to estimate the level populations from optically thin statistical equillibrium. This assumes that the emission is optically thin by setting the escape probability β=1.

ng_accelerate!(sol::Solution)

Apply Ng-acceleration to level populations from previous three iterations. See S4.4.7 of "Radiative Transfer in Astrophysics: Theory, Numerical Methods, and Applications" lecture note series by C.P. Dullemond which is itself based on Olson, Auer, and Buchler (1985).

solve_rate_eq!(xpop, yrate, rhs)

Solve the system of linear equations for the excitation rates to and from each energy level. The result is computed using LU factorization and back-substition on the yrate matrix without pivoting.

The system of equations $Y x = r$ is inverted to $x = Y^{-1} r$ and solved, where $Y$ (yrate) is the rate matrix encoding the total rate of excitation/dexcitation from one level to another, $x$ (xpop) is the level populations, and $r$ is the right-hand side (rhs) expressing the time rate-of-change of the level populations. Statistical equillibrium assumes that the derivative rhs term is all zero except for the conservation equation.

Notes

These operations correspond to the LU factorization subroutines ludcmp (lower-upper decomposition) and lubksb (lower-upper back substition) Fortran functions in the F77 RADEX source code matrix.f and slatec.f. In practice these subroutines simply prepare (and copy) the arguments for the subroutine SGEIR, which does both the decomposition and solution through back-substitution. In lubksb, the equation for the last energy level is over-written with the conservation equation (i.e., ones) to make the rate matrix square and non-singular (otherwise the last "column" would be all zeros). This is also done here by over-writing the rates in-place with ones.

The LU factorization lu! replaces yrate in-place and returns a view. Some work products are still generated with LAPACK/BLAS, but benchmarking shows that these impose negligible overhead.

A more rigorous factorization would be to include the conservation equation as an additional row and solve the rectangular system of equations $(m+1,m)$. This is also what DESPOTIC does using the non-negative least squares (nnls) solver. Tests with QR factorization and Gaussian elimination however were less numerically stable than LU factorization as well as a factor of a few slower.

solve_rate_eq_reduced!(xpop, yrate, nr)

Solve the rate equation for a restricted/reduced number of levels. Levels higher than nr are assumed to be optically thin and only coupled by radiative processes. A cascade to levels less than or equal to nr is computed for higher-lying transitions.

This formalism may not be strictly accurate for complex molecules (e.g., $\mathrm{CH_3OH}$) where many transitions between levels may be radiatively forbidden.

step_collision!(sol::Solution, rdf::RunDef)

Compute contribution of collisional processes to the rate matrix. Modifies the Solution in-place.

step_radiative!(sol::Solution, rdf::RunDef)

Compute contribution of radiative processes to the rate matrix using the escape probability. Modifies the Solution in-place.

step_tau_tex!(sol::Solution, rdf::RunDef)

Compute the optical depth and excitation temperature for each line.

get_interp(A::Matrix, axes; B=<interpolation-type>)

Generate a 5D interpolation object for a given τ or Tₑₓ cube where the last dimension indexes the transition number. By default a second order (quadratic) interpolation is used with flat extrapolation for out-of-bounds access. The axis containing the transition number is not interpolated – an exactly matching value must be given.

Please ensure that the axes passed have a linear step size (nb. use log(n) or log(N) axes if using logarithmically spaced values). A gridded interpolation object can be created using the Gridded interpolation type.

Compare the interpolation to a uniform random sample and return an array of the deviations.

Note that this implementation currently only works for 5 parameter grids of (n, Tₖ, N, Δv, j) and where n is the volume density of the specified collision partner.

rungrid(mol, density, tkin, cdmol, deltav, transitions; ...)

Generate a grid of τ and Tₑₓ values over the set of physical conditions and selected transitions. Most parameters for a RunDef run definition can be set as keyword arguments (e.g., β, bg). Arguments can be passed as scalars or collections. Runs will be executed in parallel if Julia was started with multiple threads.

Resulting cubes are indexed with shape (n, Tₖ, N, δv, j) for volume density, kinetic temperature, column density, linewidth, and transition number. The transition index number corresponds to the TRANS value in the LAMDA data file or the index number in the output table.

runseq(sol, rdfs; ...)

Execute a series model definitions. All RunDef run definitions should be for the same Specie type (i.e., LAMDA file).

Wrap SLATEC

Wrap the slatec.f Fortran file distributed as part of RADEX. The slatec.f library is used to solve the system of equations for statistical equillibrium and return the populations for each level. The rate matrix is factored into lower and upper diagonal components (LU decomposition) and solved using partial-pivoting. The subroutines ludcmp (LU decompose) and lubksb (LU back-substition) format the inputs for the main SGEIR subroutine, which itself calls the LINPACK subroutines SGEFA and SGESL.

The Julia standard library includes LINPACK and uses LU-decomposition with partial-pivoting as the default algorithm for matrix inversion when calling A \ b, thus this module is not necessary but included for the purposes of verifying the output against the exact output of RADEX.

Usage

First compile the slatec.f file distributed in the src directory of RADEX into a shared library libslatec.so:

    gfortran -shared -O2 slatec.f -o libslatec.so -fPIC

The Fortran standard does not specify how symbol names are expressed in shared libraries and are mangled in different ways by different compilers and even versions of the same compiler. To check, use the nm program to list symbol names from the object file:

    nm libslatec.so

Look for entries in the output of nm that are similar to "ludcmp" and "lubksb". The gfortran compiler typically appends an underscore to the symbol name, thus ludcmp_ and lubksb_. If these are not the names, then this file will need to edited to reflect those names, as Julia's ccall must be given literals for the symbol, library name, output type, and input types (although this could be avoided with a macro that replaces these values at compile time).

Finally, ensure that the directory containing the library file is included in the shell environment variable LD_LIBRARY_PATH:

    # Run before starting Julia or include in shell configuration.
    export LD_LIBRARY_PATH=<path-to-lib>:$LD_LIBRARY_PATH
    # Or include when starting Julia.
    LD_LIBRARY_PATH=<path-to-lib> julia

To use, call the call_slatec! function to modify the rhs level population work array in-place.

call_slatec!(yrate, rhs)

Solve the system of equations in yrate using LU-decomposition and partial-pivoting to perform the matrix inversion operation yrate \ rhs. The results are stored by modifying the vector rhs in-place. The Fortran routines in RADEX are called natively and included for verification purposes.

For testing results with LU factorization as used by RADEX.