Main modules

Cosmological calculator

toolslyman.cosmology_calc.cdist_to_z(cdist, cosmo=None, zmin=0.0001, zmax=1100, n_bins=1000)[source]

Convert comoving distance(s) to redshift(s) using an interpolated inverse of the comoving distance-redshift relation.

Parameters:

  • cdist (Quantity or array-like) – Comoving distance(s) to convert (can be with or without units).

  • cosmo (astropy.cosmology instance, optional) – Cosmological model to use. Defaults to toolslyman.cosmology.cosmo if None.

  • zmin (float, optional) – Minimum redshift to consider for interpolation (default: 1e-4).

  • zmax (float, optional) – Maximum redshift to consider for interpolation (default: 1100).

  • n_bins (int, optional) – Number of redshift samples to build the interpolation grid (default: 1000).

Returns:

z_arr (ndarray) – Redshift(s) corresponding to the input comoving distance(s).

Notes

  • The function uses a spline interpolation in log-log space for better accuracy across a wide range of redshifts and distances.

  • If any values fall outside the initial interpolation range, the grid is automatically extended and the interpolation recalculated.

Lyman-alpha forest transmission

toolslyman.lyman_alpha_transmission.approximation_voigt_faddeeva(x, a)[source]

Compute the Voigt profile using the Faddeeva function.

The Voigt profile combines Doppler and pressure broadening effects and can be computed using the Faddeeva function (complex error function). This approximation is particularly useful for calculating the Voigt profile efficiently.

Parameters:
xfloat or ndarray

The frequency or wavelength offset from the line center, scaled by the Doppler width. For example, ( x =

rac{
u -
u_0}{Delta
u_d} ), where (
u ) is the frequency, (
u_0 ) is the center frequency, and ( Delta
u_d ) is the Doppler width.
afloat

The line a-factor, which characterizes the ratio of the Lorentzian width to the Doppler width. It is defined as ( a =

rac{gamma}{Delta
u_d} ), where ( gamma ) is the Lorentzian width and ( Delta
u_d ) is the Doppler width.
Returns:
Hfloat or ndarray

The value of the Voigt profile at the given x values. The output shape matches the input shape.

toolslyman.lyman_alpha_transmission.direct_voigt_profile(x, a)[source]

Compute the Voigt profile using direct numerical integration of the line profile function.

The Voigt profile is a combination of a Lorentzian and a Gaussian profile, used to model spectral lines broadened by both Doppler and pressure broadening effects. This function calculates the Voigt profile by directly integrating the line profile function.

Parameters:
xarray_like

The array of positions (e.g., frequency or wavelength offsets) at which to evaluate the Voigt profile.

afloat or array_like

The line a-factor, which characterizes the ratio of the Lorentzian to Gaussian broadening. It is defined as ( a =

rac{gamma}{sigma sqrt{2}} ), where ( gamma ) is the Lorentzian width and ( sigma ) is the Gaussian width.
Returns:
Vndarray

The Voigt profile evaluated at the input positions x. The shape of the output matches the shape of the input x.

toolslyman.lyman_alpha_transmission.lyA_scattering_cross_section(nu, Tk)[source]

Calculate the Lyman-alpha scattering cross-section for a given frequency and temperature.

Parameters:

  • nu (float or astropy.units.quantity.Quantity) – The frequency at which to calculate the cross-section. Can be given as a float in Hz or as an astropy Quantity.

  • Tk (float or astropy.units.quantity.Quantity) – The temperature of the gas. Can be given as a float in K or as an astropy Quantity.

Returns:

sigmaA (astropy.units.quantity.Quantity) – The Lyman-alpha scattering cross-section at the given frequency and temperature, in units of cm^2.

Notes

This function uses parameters from e.g. Dijkstra (2014, arXiv:1406.7292). The Voigt profile is approximated using the Faddeeva function.

toolslyman.lyman_alpha_transmission.tau_nu_lyA(z, xHI, Tk, dr, nu=<Quantity 2.46e+15 Hz>, X_H=0.76, cosmo=None)[source]

Estimate the Lyman-alpha optical depth on a simulation grid.

Parameters:

  • z (float) – The redshift at which the optical depth is being calculated.

  • xHI (ndarray) – The neutral hydrogen fraction at each grid point.

  • Tk (float or astropy.units.quantity.Quantity) – The temperature of the gas at each grid point. Can be given as a float in K or as an astropy Quantity.

  • dr (astropy.units.quantity.Quantity) – The comoving path length through the grid cells. Should be an astropy Quantity.

  • nu (float or astropy.units.quantity.Quantity, optional) – The frequency at which to calculate the optical depth, default is the Lyman-alpha frequency (2.46e15 Hz).

  • X_H (float, optional) – The hydrogen mass fraction, default is 0.76.

  • cosmo (astropy.cosmology instance, optional) – The cosmology to use for the calculation. If None, Planck18 cosmology is assumed.

Returns:

tau_nu (ndarray) – The Lyman-alpha optical depth at each grid point.

Notes

This function uses Eq. 7 from Smith, A. et al. (2015, arXiv:1409.4480). The cross-section and Voigt profile parameters are taken from e.g. Dijkstra (2014, arXiv:1406.7292).

Example

>>> from astropy import units as u
>>> from astropy.cosmology import Planck18
>>> z = 6.0
>>> xHI = np.random.rand(100, 100, 100)
>>> Tk = 1e4 * u.K
>>> dr = 1.0 * u.Mpc
>>> tau = tau_nu_lyA(z, xHI, Tk, dr)
toolslyman.lyman_alpha_transmission.tau_nu_lyA_FGPA(z, xHI, dn, f_alpha=0.416, lambda_alpha=<Quantity 1216. Angstrom>, k_res=1.0, X_H=0.76, cosmo=None)[source]

Estimate the Lyman-alpha optical depth on a simulation grid.

Parameters:

  • z (float) – The redshift at which the optical depth is being calculated.

  • xHI (ndarray) – The neutral hydrogen fraction at each grid point.

  • dn (ndarray) – The overdensity at each grid point.

  • f_alpha (float) – The Lyman alpha oscillator strength, default is 0.416.

  • lambda_alpha (float or astropy.units.quantity.Quantity, optional) – The Lyman-alpha wavelength (1216 A).

  • f_rescale (float) – The normalization factor to account for the unresolved small-scale fluctuations in the density and velocity fields, default is 1.0.

  • X_H (float, optional) – The hydrogen mass fraction, default is 0.76.

  • cosmo (astropy.cosmology instance, optional) – The cosmology to use for the calculation. If None, Planck18 cosmology is assumed.

Returns:

tau_nu (ndarray) – The Lyman-alpha optical depth at each grid point.

Notes

This function uses equation 21 in Choudhury, Paranjape & Bosman (2021).

Example

>>> from astropy import units as u
>>> from astropy.cosmology import Planck18
>>> z = 6.0
>>> xHI = np.radom.rand(100, 100, 100)
>>> Tk = 1e4 * u.K
>>> dr = 1.0 * u.Mpc
>>> tau = tau_nu_lyA(z, xHI, Tk, dr)
toolslyman.lyman_alpha_transmission.tau_nu_lyA_los(z, xHI, dn, los_len=<Quantity 30. Mpc>, box_len=None, los_axis=2, method='FGPA', X_H=0.76, cosmo=None, **kwargs)[source]

Estimate the effective Lyman-alpha optical depth along the line of sight, smoothed over different scales.

Parameters:

  • z (float) – The redshift at which the optical depth is being calculated.

  • xHI (ndarray) – The neutral hydrogen fraction at each grid point.

  • dn (ndarray) – The overdensity at each grid point.

  • los_len (astropy.units.quantity.Quantity or int) – The length along the line of sight over which the effective optical depth is estimated. If provided as an integer, it is assumed to be the number of grid cells.

  • method (str) – The formalism used to estimate the Lyman-alpha optical depth, default is ‘FGPA’.

Returns:

tau_nu (ndarray) – The Lyman-alpha optical depth at each grid point.

Example

>>> from astropy import units as u
>>> from astropy.cosmology import Planck18
>>> z = 6.0
>>> xHI = np.random.rand(100, 100, 100)
>>> tau = tau_nu_lyA_los(z, xHI)

Damped Lyman-alpha transmission

toolslyman.lyman_alpha_damping_wings.column_density_along_skewer(z_source, xHI, dn, dr, X_H=0.76, cosmo=None)[source]

Compute the cumulative neutral hydrogen column density along a skewer.

Parameters:

  • z_source (float) – Redshift of the background source.

  • xHI (ndarray) – Neutral hydrogen fraction along the skewer.

  • dn (ndarray) – Overdensity field (δ = ρ/ρ̄ - 1).

  • dr (Quantity or float) – Comoving cell length (can be with or without units).

  • X_H (float, optional) – Hydrogen mass fraction. Default is 0.76.

  • cosmo (astropy.cosmology, optional) – Cosmology instance. If None, uses default from toolslyman.

Returns:

N_HI (ndarray) – Comoving neutral hydrogen column density along the skewer (in cm^-2).

toolslyman.lyman_alpha_damping_wings.compute_spectrum(xvel_in, xvel_out, cdens, temp, lambda0, fvalue, mass, damped, periodic)[source]

Returns optical depth array

xvel_in : velocity (km/s) of array element (x-coordinate of spectrum) xvel_out : velocity (km/s) of array element (x-coordinate of spectrum) cdens : column density (particles/cm^2) temp : temperature (K) lambda0 : rest wavelength (Å) fvalue : oscillator strength mass : mass of atom (g)

toolslyman.lyman_alpha_damping_wings.optical_depth_lyA_along_skewer(z_source, xHI, dn, dr=None, z_arr=None, temp=<Quantity 10000. K>, X_H=0.76, cosmo=None, f_alpha=0.4164, damped=True, verbose=False, use_compute_spectrum=True)[source]

Compute the Lyman-alpha optical depth (τ) along one or more cosmological skewers.

This function supports both single and multiple skewers (2D arrays). The source is assumed to be at the origin or index 0 of the skewer array if dr is provided. Use np.roll to adjust the line-of-sight skewer if necessary.

Parameters:

  • z_source (float) – Redshift of the background source (e.g., a quasar).

  • xHI (ndarray) – Neutral hydrogen fraction along the skewer. Shape can be (N,) or (N_skewers, N).

  • dn (ndarray) – Baryon overdensity field (δ = ρ/ρ̄ - 1). Shape must match xHI.

  • dr (Quantity or float) – Comoving length of each cell along the skewer (e.g., in Mpc or cm).

  • z_arr (ndarray) – Redshift array corresponding to the skewer.

  • temp (Quantity or ndarray, optional) – Gas temperature in Kelvin. Can be scalar, 1D, or 2D array matching shape of xHI. Default is 1e4 K.

  • X_H (float, optional) – Hydrogen mass fraction. Default is 0.76.

  • cosmo (astropy.cosmology.Cosmology, optional) – Cosmology instance to use. If None, uses the default from toolslyman.

  • f_alpha (float, optional) – Oscillator strength for the Lyman-alpha transition. Default is 0.4164.

  • damped (bool, optional) – If True, include the damping wing using a Voigt profile. Otherwise, use only the Doppler core (Gaussian).

  • verbose (bool, optional) – If True, print progress every 10 steps. Also prints skewer number if multiple skewers are processed.

Returns:

  • tau_lambda (ndarray) – Optical depth as a function of observed wavelength. Shape is (N_lambda,) for single skewer or (N_skewers, N_lambda) for multiple skewers.

  • lambda_obs (Quantity) – Observed wavelength grid corresponding to tau_lambda, in Ångström.