NonLinear Module Source file:halofit

Classes

Dependencies

Subroutines
  • allocate_LUT(lut, n)
    Allocates memory for the HMcode look-up HM_tables
  • assign_HM_cosmology(this, State, cosm)
  • calculate_nowiggle(k, z, Pk, Pk_nw, cosm)
    Calculate the normalised no wiggle power spectrum at a range of k and a Comes from the Eisenstein & Hu approximation
    • real(dl) intent(in) :: k(:)
    • real(dl) intent(in) :: z
    • real(dl) intent(in) :: Pk(:)
    • real(dl) allocatable, intent(out) :: Pk_nw(:)
    • HM_cosmology intent(in)  :: cosm
  • calculate_psmooth(k, z, Pk, Pk_smt, cosm)
    Calculate the normalised smoothed power spectrum at a range of k
    • real(dl) intent(in) :: k(:)
    • real(dl) intent(in) :: z
    • real(dl) intent(in) :: Pk(:)
    • real(dl) allocatable, intent(out) :: Pk_smt(:)
    • HM_cosmology intent(in)  :: cosm
  • fill_conc(this, z, cosm, lut)
  • fill_growtab(cosm)
    Fills a table of values of the scale-independent growth function
  • fill_plintab(CAMB_PK, iz, cosm)
    Fills internal HMcode HM_tables for the linear power spectrum at z = 0
  • fill_sigtab(this, cosm)
  • fill_table(min, max, arr, n)
    Fills array 'arr' in equally spaced intervals
    • real(dl) intent(in) :: min
    • real(dl) intent(in) :: max
    • real(dl) allocatable :: arr(:)
    • integer intent(in) :: n
  • fit_cubic(a, b, c, d, x1, y1, x2, y2, x3, y3, x4, y4)
    Given xi, yi i = 1,2,3,4 fits a cubic between these points
    • real(dl) intent(out) :: a
    • real(dl) intent(out) :: b
    • real(dl) intent(out) :: c
    • real(dl) intent(out) :: d
    • real(dl) intent(in) :: x1
    • real(dl) intent(in) :: y1
    • real(dl) intent(in) :: x2
    • real(dl) intent(in) :: y2
    • real(dl) intent(in) :: x3
    • real(dl) intent(in) :: y3
    • real(dl) intent(in) :: x4
    • real(dl) intent(in) :: y4
  • fit_line(a0, a1, x1, y1, x2, y2)
    Given xi, yi i = 1,2 fits a line between these points
    • real(dl) intent(out) :: a0
    • real(dl) intent(out) :: a1
    • real(dl) intent(in) :: x1
    • real(dl) intent(in) :: y1
    • real(dl) intent(in) :: x2
    • real(dl) intent(in) :: y2
  • fit_quadratic(a0, a1, a2, x1, y1, x2, y2, x3, y3)
    Given xi, yi i = 1,2,3 fits a quadratic between these points
    • real(dl) intent(out) :: a0
    • real(dl) intent(out) :: a1
    • real(dl) intent(out) :: a2
    • real(dl) intent(in) :: x1
    • real(dl) intent(in) :: y1
    • real(dl) intent(in) :: x2
    • real(dl) intent(in) :: y2
    • real(dl) intent(in) :: x3
    • real(dl) intent(in) :: y3
  • halomod_tables(this, z, cosm, lut)
  • HMcode_redshift(this, CAMB_Pk, j, npass, imead_pass, cosm, lut, p_den, p_num)
    Everything HMcode does for one redshift; cosm and lut are the caller's working space, so this can be called either serially or from a thread with its own private copies of them
    • THalofit   :: this
    • MatterPowerData   :: CAMB_Pk
    • integer intent(in) :: j
    • integer intent(in) :: npass
    • integer intent(in) :: imead_pass(:)
    • HM_cosmology   :: cosm
    • HM_tables   :: lut
    • real(dl) allocatable :: p_den(:, :)
    • real(dl) allocatable :: p_num(:, :)
  • init_wiggle(cosm)
    Isolate the power spectrum wiggle
  • initialise_HM_cosmology(this, CAMB_PK, cosm, iz)
  • Mead_growth_terms(z, cosm, x, y, Om_m)
    Growth terms shared by the Mead (2017; 1606.05345) delta_c and Delta_v fitting functions
    • real(dl) intent(in) :: z
    • HM_cosmology intent(in)  :: cosm
    • real(dl) intent(out) :: x
    • real(dl) intent(out) :: y
    • real(dl) intent(out) :: Om_m
  • ode_growth(xi, ti, tf, acc, vi, x, v, t, imeth, cosm)
    Solves 2nd order ODE x''(t) from ti to tf and writes out array of x, v, t values
    • real(dl)  :: xi
    • real(dl)  :: ti
    • real(dl)  :: tf
    • real(dl)  :: acc
    • real(dl)  :: vi
    • real(dl) allocatable :: x(:)
    • real(dl) allocatable :: v(:)
    • real(dl) allocatable :: t(:)
    • integer  :: imeth
    • HM_cosmology intent(in)  :: cosm
  • omegas_hf(aa, om_m0, om_v0, wval, waval, om_m, om_v)
    Evolution of omega_m and omega_lambda with expansion factor, generalized to w0-wa
    • real(dl) intent(in) :: aa
    • real(dl) intent(in) :: om_m0
    • real(dl) intent(in) :: om_v0
    • real(dl) intent(in) :: wval
    • real(dl) intent(in) :: waval
    • real(dl) intent(out) :: om_m
    • real(dl) intent(out) :: om_v
  • PKequal(State, redshift, w_eff, wa_eff, w_hf, wa_hf)
    used by halofit_casarini: arXiv:0810.0190, arXiv:1601.07230 Solve for the constant-w model (w_hf, wa_hf=0) that has the same comoving distance between redshift and last scattering as the actual dark energy model, keeping all the other densities (and the present dark energy density) fixed. The Takahashi halofit fit is then evaluated for that equivalent constant-w model. w_eff, wa_eff are the effective w0-wa of the actual model, used only to bracket the root.
    • CAMBdata target  :: State
    • real(dl) intent(in) :: redshift
    • real(dl) intent(in) :: w_eff
    • real(dl) intent(in) :: wa_eff
    • real(dl) intent(out) :: w_hf
    • real(dl) intent(out) :: wa_hf
  • reverse(n, arry, output)
    This reverses the contents of arry!
    • integer intent(in) :: n
    • real(dl) intent(in) :: arry(n)
    • real(dl) allocatable, intent(out) :: output(:)
  • SiCi(x, Six, Cix)
    Calculates the 'sine integral' Si(x) and 'cosine integral' Ci(x) together The large-x expansions share the same auxiliary functions f and g, so are done once
    • real(dl) intent(in) :: x
    • real(dl) intent(out) :: Six
    • real(dl) intent(out) :: Cix
  • smooth_array_Gaussian(x, f, sigma)
    Smooth an array f(x) using a Gaussian kernel If x is equally spaced (as it is here) the kernel depends only on |i-j| and can be tabulated once, which avoids the n^2 exponentials of the general case
    • real(dl) intent(in) :: x(:) x coordinates
    • real(dl) intent(inout) :: f(:) Array to smooth
    • real(dl) intent(in) :: sigma Width of smoothing Gaussian
  • Tk_nw_init(cosm, s, alpha)
    The k-independent parameters of the no-wiggle transfer function; only needed once per cosmology
    • HM_cosmology intent(in)  :: cosm
    • real(dl) intent(out) :: s
    • real(dl) intent(out) :: alpha
  • wint(fac, r, sig, d1, d2)
    • real(dl) intent(in) :: fac(nint_wint)
    • real(dl) intent(in) :: r
    • real(dl) intent(out) :: sig
    • real(dl) intent(out) :: d1
    • real(dl) intent(out) :: d2
  • wint_pk_table(CAMB_Pk, itf, fac)
    Tabulates the part of the wint integrand that does not depend on the filter scale r, so that the bisection for r_nl does not repeat the (expensive) power spectrum look-ups
    • MatterPowerData   :: CAMB_Pk
    • integer intent(in) :: itf
    • real(dl) intent(out) :: fac(nint_wint)
Functions
  • real(dl)
    AH(z, cosm)
    The Hubble acceleration function \ddot{a}/a
  • real(dl)
    cached_grow(z, cosm)
    The growth factor, using the value cached by fill_plintab where possible
  • real(dl)
    cosmic_density(cosm)
    The z = 0 cosmological matter density
  • real(dl)
    dc_Mead(z, cosm)
    delta_c fitting function from Mead (2017; 1606.05345)
  • real(dl)
    distance_error(obj, w)
    Fractional difference between the comoving distance from a_z to a_star in a model with constant dark energy equation of state w and the same distance in the actual model
    • class(*)  :: obj
    • real(dl) intent(in) :: w
  • real(dl)
    dtauda_wconst(obj, a)
    d tau/d a for a model with constant dark energy equation of state, and all other densities as in the actual model. Agrees with dtauda to the last bit when w is the actual (constant) equation of state.
    • class(*)  :: obj
    • real(dl) intent(in) :: a
  • real(dl)
    Dv_Mead(z, cosm)
    Delta_v fitting function from Mead (2017; 1606.05345)
  • real(dl)
    f_Mead(x, y, p0, p1, p2, p3)
    Equation A3 in Mead (2017)
    • real(dl) intent(in) :: x
    • real(dl) intent(in) :: y
    • real(dl) intent(in) :: p0
    • real(dl) intent(in) :: p1
    • real(dl) intent(in) :: p2
    • real(dl) intent(in) :: p3
  • real(dl)
    f_star(lut, cosm)
  • real(dl)
    fd(v)
    d'=f(d) in ODE solver
    • real(dl) intent(in) :: v
  • real(dl)
    find(n, x, xtab, ytab, iorder, ifind, imeth)
    Given two arrays x and y this routine interpolates to find the y_i value at position x_i
    • integer intent(in) :: n
    • real(dl) intent(in) :: x
    • real(dl) intent(in) :: xtab(n)
    • real(dl) intent(in) :: ytab(n)
    • integer intent(in) :: iorder
    • integer intent(in) :: ifind
    • integer intent(in) :: imeth
  • real(dl)
    find_pk(k, itype, cosm)
    Look-up and interpolation for P(k,z=0)
    • real(dl) intent(in) :: k
    • integer intent(in) :: itype
    • HM_cosmology intent(in)  :: cosm
  • real(dl)
    fv(d, v, a, cosm)
    v'=f(v) in ODE solver
    • real(dl) intent(in) :: d
    • real(dl) intent(in) :: v
    • real(dl) intent(in) :: a
    • HM_cosmology intent(in)  :: cosm
  • real(dl)
    gnu(nu)
    Select the mass function
    • real(dl) intent(in) :: nu
  • real(dl)
    grow(z, cosm)
    Finds the scale-independent growth function at redshift z
  • real(dl)
    growint(z, cosm)
    Approximate growth function from integrating Omega_m(a)^gamma dln(a) up to a = 1
  • real(dl)
    growint_integrand(a, y, z, itype, cosm)
    Integrand for the approximate growth integral y, z and itype are unused, but are needed to match the integrate() interface
    • real(dl) intent(in) :: a
    • real(dl) intent(in) :: y
    • real(dl) intent(in) :: z
    • integer intent(in) :: itype
    • HM_cosmology intent(in)  :: cosm
  • real(dl)
    gst(nu)
    Sheth & Tormen (1999) mass function!
    • real(dl) intent(in) :: nu
  • logical
    HM_par_inner()
    The loops over redshift are the outermost parallel regions, and give much better scaling than the small regions used to set up each redshift. Those inner regions are therefore only used when the redshift loop is not itself running in parallel (e.g. for a single redshift). $ use omp_lib, only: omp_in_parallel
  • real(dl)
    Hubble2(z, cosm)
    This calculates the dimensionless squared hubble parameter at redshift z (H/H_0)^2! Ignores contributions from radiation (not accurate at high z, but consistent with simulations)!
  • integer
    int_split(n, x, xtab)
    Finds the position of the value in the table by continually splitting it in half
    • integer intent(in) :: n
    • real(dl) intent(in) :: x
    • real(dl) intent(in) :: xtab(n)
  • real(dl)
    integrate(a, b, y, z, itype, cosm, acc, iorder, f)
    Integrates between a and b until desired accuracy is reached Stores information to reduce function calls
    • real(dl) intent(in) :: a
    • real(dl) intent(in) :: b
    • real(dl) intent(in) :: y
    • real(dl) intent(in) :: z
    • integer intent(in) :: itype
    • HM_cosmology intent(in)  :: cosm
    • real(dl) intent(in) :: acc
    • integer intent(in) :: iorder
    • function  :: f(x, y, z, itype, cosm)
  • real(dl)
    inttab(n1, n2, x, y, iorder)
    Integrates tables y(x)dx
    • integer intent(in) :: n1
    • integer intent(in) :: n2
    • real(dl) intent(in) :: x(:)
    • real(dl) intent(in) :: y(:)
    • integer intent(in) :: iorder
  • real(dl)
    Lagrange_polynomial(n, x, xv, yv)
    Computes the result of the nth order Lagrange polynomial at point x, L(x)
    • integer intent(in) :: n
    • real(dl) intent(in) :: x
    • real(dl) intent(in) :: xv(n + 1)
    • real(dl) intent(in) :: yv(n + 1)
  • integer
    linear_table_integer(n, x, xtab)
    Assuming the table is exactly linear this gives you the integer position Clamped because the tabulated points can differ from the exactly equally spaced values by a last bit, which could otherwise put x in a neighbouring cell
    • integer intent(in) :: n
    • real(dl) intent(in) :: x
    • real(dl) intent(in) :: xtab(n)
  • real(dl)
    m_baryon(lut, cosm)
  • real(dl)
    mass(c)
    This calculates the (normalised) mass of a halo of concentration c The 'normalised' mass is that divided by the prefactor r_s^3 4*pi rho_n where rho_n is the profile normalisation [i.e, rho=rho_n/((r/r_s)*(1. + r/r_s)^2]
    • real(dl) intent(in) :: c
  • real(dl)
    neff(cosm, lut)
    Finds the effective spectral index at the collapse scale r_nl, where nu(r_nl) = 1.
  • real(dl)
    neff_integrand(t, R, z, itype, cosm)
    This is the integrand for the velocity dispersion integral
    • real(dl) intent(in) :: t
    • real(dl) intent(in) :: R
    • real(dl) intent(in) :: z
    • integer intent(in) :: itype
    • HM_cosmology intent(in)  :: cosm
  • real(dl)
    Omega_cold_hm(z, cosm)
    This calculates omega_cold variations with z (no neutrinos)
  • real(dl)
    Omega_m_hm(z, cosm)
    This calculates omega_m variations with z!
  • real(dl)
    p_dewiggle(k, z, p_linear, sigv, cosm)
    Call the dewiggled power spectrum, which is linear but with damped wiggles
    • real(dl) intent(in) :: k
    • real(dl) intent(in) :: z
    • real(dl) intent(in) :: p_linear
    • real(dl) intent(in) :: sigv
    • HM_cosmology intent(in)  :: cosm
  • real(dl)
    p_lin(itype, cosm)
    Looks up the value for the linear power spectrum
  • real(dl)
    Pk_nowiggle(k, s, alpha, cosm)
    Calculates the un-normalised no-wiggle power spectrum Comes from the Eisenstein & Hu approximation
    • real(dl) intent(in) :: k
    • real(dl) intent(in) :: s k-independent parameters from Tk_nw_init
    • real(dl) intent(in) :: alpha k-independent parameters from Tk_nw_init
    • HM_cosmology intent(in)  :: cosm
  • real(dl)
    r_nl(lut)
    Calculates R_nl, defined by nu(R_nl) = 1., nu = dc/sigma(R)
  • real(dl)
    radius_m(m, cosm)
    Calculates the co-moving radius that encloses a mass 'm' in the homogeneous Universe
  • real(dl)
    sigma_integral(r, z, itype, cosm)
    Gets sigma(R)
    • real(dl) intent(in) :: r
    • real(dl) intent(in) :: z
    • integer intent(in) :: itype
    • HM_cosmology intent(in)  :: cosm
  • real(dl)
    sigma_integrand(t, R, z, itype, cosm)
    The integrand for the sigma(R) integrals
    • real(dl) intent(in) :: t
    • real(dl) intent(in) :: R
    • real(dl) intent(in) :: z
    • integer intent(in) :: itype
    • HM_cosmology intent(in)  :: cosm
  • real(dl)
    sigma_lut(r, z, cosm)
    Finds sigma_cold(R) from look-up table
    • real(dl) intent(in) :: r
    • real(dl) intent(in) :: z
    • HM_cosmology intent(in)  :: cosm
  • real(dl)
    sigmaV(R, z, itype, cosm)
    • real(dl) intent(in) :: R
    • real(dl) intent(in) :: z
    • integer intent(in) :: itype
    • HM_cosmology intent(in)  :: cosm
  • real(dl)
    sigmaV_integrand(t, R, z, itype, cosm)
    This is the integrand for the velocity dispersion integral
    • real(dl) intent(in) :: t
    • real(dl) intent(in) :: R
    • real(dl) intent(in) :: z
    • integer intent(in) :: itype
    • HM_cosmology intent(in)  :: cosm
  • integer
    table_integer(n, x, xtab, imeth)
    Chooses between ways to find the integer location below some value in an array find() only calls this for x inside the table, and both methods clamp anyway
    • integer intent(in) :: n
    • real(dl) intent(in) :: x
    • real(dl) intent(in) :: xtab(n)
    • integer intent(in) :: imeth
  • real(dl)
    Tcb_Tcbnu_ratio(k, z, cosm)
    Calculates the ratio of T(k) for cold vs. all matter Uses approximations in Eisenstein & Hu (1999; arXiv 9710252) Note that this assumes that there are exactly 3 species of neutrinos with Nnu<=3 of these being massive, and with the mass split evenly between the number of massive species.
    • real(dl) intent(in) :: k
    • real(dl) intent(in) :: z
    • HM_cosmology intent(in)  :: cosm
  • real(dl)
    Tk_nw(k, s, alpha, cosm)
    No-wiggle transfer function from Eisenstein & Hu: astro-ph:9709112
    • real(dl) intent(in) :: k Wavenumber [h/Mpc]
    • real(dl) intent(in) :: s k-independent parameters from Tk_nw_init
    • real(dl) intent(in) :: alpha k-independent parameters from Tk_nw_init
    • HM_cosmology intent(in)  :: cosm
  • real(dl)
    w_de_hm(z, cosm)
    The dark energy w(a) function
  • real(dl)
    win(k, rv, c)
    Selects the halo window function (k-space halo profile)
    • real(dl) intent(in) :: k
    • real(dl) intent(in) :: rv
    • real(dl) intent(in) :: c
  • real(dl)
    winnfw(k, rv, c)
    The analytic Fourier Transform of the NFW profile; note W(k->0) = 1
    • real(dl) intent(in) :: k
    • real(dl) intent(in) :: rv
    • real(dl) intent(in) :: c
  • real(dl)
    wk_tophat(x)
    The normlaised Fourier Transform of a top-hat
    • real(dl) intent(in) :: x
  • real(dl)
    wk_tophat_deriv(x)
    The derivative of a normlaised Fourier Transform of a spherical top-hat
    • real(dl) intent(in) :: x
  • real(dl)
    X_de(z, cosm)
    The time evolution for dark energy: rho_de = rho_de,0 * X(a) X(a) = 1 for LCDM but changes for other models