#include "fourier.h"

Include dependency graph for fourier.c:

Functions
int	fourier_pk_at_z (struct background pba, struct fourier pfo, enum linear_or_logarithmic mode, enum pk_outputs pk_output, double z, int index_pk, double out_pk, double out_pk_ic)

int	fourier_pk_at_k_and_z (struct background pba, struct primordial ppm, struct fourier pfo, enum pk_outputs pk_output, double k, double z, int index_pk, double out_pk, double *out_pk_ic)

int	fourier_pks_at_kvec_and_zvec (struct background pba, struct fourier pfo, enum pk_outputs pk_output, double kvec, int kvec_size, double zvec, int zvec_size, double out_pk, double out_pk_cb)

int	fourier_pk_tilt_at_k_and_z (struct background pba, struct primordial ppm, struct fourier pfo, enum pk_outputs pk_output, double k, double z, int index_pk, double pk_tilt)

int	fourier_sigmas_at_z (struct precision ppr, struct background pba, struct fourier pfo, double R, double z, int index_pk, enum out_sigmas sigma_output, double result)

int	fourier_k_nl_at_z (struct background pba, struct fourier pfo, double z, double k_nl, double k_nl_cb)

int	fourier_init (struct precision ppr, struct background pba, struct thermodynamics pth, struct perturbations ppt, struct primordial ppm, struct fourier pfo)

int	fourier_free (struct fourier *pfo)

int	fourier_indices (struct precision ppr, struct background pba, struct perturbations ppt, struct primordial ppm, struct fourier *pfo)

int	fourier_get_k_list (struct precision ppr, struct perturbations ppt, struct fourier *pfo)

int	fourier_get_tau_list (struct perturbations ppt, struct fourier pfo)

int	fourier_get_source (struct background pba, struct perturbations ppt, struct fourier pfo, int index_k, int index_ic, int index_tp, int index_tau, double sources, double source)

int	fourier_pk_linear (struct background pba, struct perturbations ppt, struct primordial ppm, struct fourier pfo, int index_pk, int index_tau, int k_size, double lnpk, double lnpk_ic)

int	fourier_sigmas (struct fourier pfo, double R, double lnpk_l, double ddlnpk_l, int k_size, double k_per_decade, enum out_sigmas sigma_output, double result)

int	fourier_sigma_at_z (struct background pba, struct fourier pfo, double R, double z, int index_pk, double k_per_decade, double *result)

int	fourier_halofit (struct precision ppr, struct background pba, struct perturbations ppt, struct primordial ppm, struct fourier pfo, int index_pk, double tau, double pk_nl, double lnpk_l, double ddlnpk_l, double k_nl, short nl_corr_not_computable_at_this_k)

int	fourier_halofit_integrate (struct fourier pfo, double integrand_array, int integrand_size, int ia_size, int index_ia_k, int index_ia_pk, int index_ia_sum, int index_ia_ddsum, double R, enum halofit_integral_type type, double *sum)

int	fourier_hmcode (struct precision ppr, struct background pba, struct perturbations ppt, struct primordial ppm, struct fourier pfo, int index_pk, int index_tau, double tau, double pk_nl, double lnpk_l, double ddlnpk_l, double k_nl, short nl_corr_not_computable_at_this_k, struct fourier_workspace *pnw)

int	fourier_hmcode_workspace_init (struct precision ppr, struct background pba, struct fourier pfo, struct fourier_workspace pnw)

int	fourier_hmcode_workspace_free (struct fourier pfo, struct fourier_workspace pnw)

int	fourier_hmcode_dark_energy_correction (struct precision ppr, struct background pba, struct fourier pfo, struct fourier_workspace pnw)

int	fourier_hmcode_baryonic_feedback (struct fourier *pfo)

int	fourier_hmcode_fill_sigtab (struct precision ppr, struct background pba, struct perturbations ppt, struct primordial ppm, struct fourier pfo, int index_tau, double lnpk_l, double ddlnpk_l, struct fourier_workspace pnw)

int	fourier_hmcode_fill_growtab (struct precision ppr, struct background pba, struct fourier pfo, struct fourier_workspace pnw)

int	fourier_hmcode_growint (struct precision ppr, struct background pba, struct fourier pfo, double a, double w0, double wa, double growth)

int	fourier_hmcode_window_nfw (struct fourier pfo, double k, double rv, double c, double window_nfw)

int	fourier_hmcode_halomassfunction (double nu, double *hmf)

int	fourier_hmcode_sigma8_at_z (struct background pba, struct fourier pfo, double z, double sigma_8, double sigma_8_cb, struct fourier_workspace *pnw)

int	fourier_hmcode_sigmadisp_at_z (struct background pba, struct fourier pfo, double z, double sigma_disp, double sigma_disp_cb, struct fourier_workspace *pnw)

int	fourier_hmcode_sigmadisp100_at_z (struct background pba, struct fourier pfo, double z, double sigma_disp_100, double sigma_disp_100_cb, struct fourier_workspace *pnw)

int	fourier_hmcode_sigmaprime_at_z (struct background pba, struct fourier pfo, double z, double sigma_prime, double sigma_prime_cb, struct fourier_workspace *pnw)

Detailed Description

Documented fourier module

Julien Lesgourgues, 6.03.2014

New module replacing an older one present up to version 2.0 The new module is located in a better place in the main, allowing it to compute non-linear correction to 's and not just . It will also be easier to generalize to new methods. The old implementation of one-loop calculations and TRG calculations has been dropped from this version, they can still be found in older versions.

Function Documentation

◆ fourier_pk_at_z()

int fourier_pk_at_z	(	struct background *	pba,
		struct fourier *	pfo,
		enum linear_or_logarithmic	mode,
		enum pk_outputs	pk_output,
		double	z,
		int	index_pk,
		double *	out_pk,
		double *	out_pk_ic
	)

Return the P(k,z) for a given redshift z and pk type (_m, _cb) (linear if pk_output = pk_linear, nonlinear if pk_output = pk_nonlinear)

In the linear case, if there are several initial conditions and the input pointer out_pk_ic is not set to NULL, the function also returns the decomposition into different IC contributions.

Hints on input index_pk:

a. if you want the total matter spectrum P_m(k,z), pass in input pfo->index_pk_total (this index is always defined)

b. if you want the power spectrum relevant for galaxy or halos, given by P_cb if there is non-cold-dark-matter (e.g. massive neutrinos) and to P_m otherwise, pass in input pfo->index_pk_cluster (this index is always defined)

c. there is another possible syntax (use it only if you know what you are doing): if pfo->has_pk_m == TRUE you may pass pfo->index_pk_m to get P_m if pfo->has_pk_cb == TRUE you may pass pfo->index_pk_cb to get P_cb

Output format:

if mode = logarithmic (most straightforward for the code): out_pk = ln(P(k)) out_pk_ic[diagonal] = ln(P_ic(k)) out_pk_ic[non-diagonal] = cos(correlation angle icxic)
if mode = linear (a conversion is done internally in this function) out_pk = P(k) out_pk_ic[diagonal] = P_ic(k) out_pk_ic[non-diagonal] = P_icxic(k)

Parameters

pba	Input: pointer to background structure
pfo	Input: pointer to fourier structure
mode	Input: linear or logarithmic
pk_output	Input: linear or nonlinear
z	Input: redshift
index_pk	Input: index of pk type (_m, _cb)
out_pk	Output: P(k) returned as out_pk_l[index_k]
out_pk_ic	Output: P_ic(k) returned as out_pk_ic[index_k * pfo->ic_ic_size + index_ic1_ic2]

Returns: the error status

check whether we need the decomposition into contributions from each initial condition
case z=0 requiring no interpolation in z
interpolation in z

--> get value of contormal time tau

-> check that tau is in pre-computed table

--> if ln(tau) much too small, raise an error

--> if ln(tau) too small but within tolerance, round it and get right values without interpolating

--> if ln(tau) much too large, raise an error

--> if ln(tau) too large but within tolerance, round it and get right values without interpolating

-> tau is in pre-computed table: interpolate

--> interpolate P_l(k) at tau from pre-computed array

--> interpolate P_ic_l(k) at tau from pre-computed array

--> interpolate P_nl(k) at tau from pre-computed array

so far, all output stored in logarithmic format. Eventually, convert to linear one.

--> loop over k

--> convert total spectrum

--> convert contribution of each ic (diagonal elements)

--> convert contribution of each ic (non-diagonal elements)

◆ fourier_pk_at_k_and_z()

int fourier_pk_at_k_and_z	(	struct background *	pba,
		struct primordial *	ppm,
		struct fourier *	pfo,
		enum pk_outputs	pk_output,
		double	k,
		double	z,
		int	index_pk,
		double *	out_pk,
		double *	out_pk_ic
	)

Return the P(k,z) for a given (k,z) and pk type (_m, _cb) (linear if pk_output = pk_linear, nonlinear if pk_output = pk_nonlinear)

In the linear case, if there are several initial conditions and the input pointer out_pk_ic is not set to NULL, the function also returns the decomposition into different IC contributions.

Hints on input index_pk:

a. if you want the total matter spectrum P_m(k,z), pass in input pfo->index_pk_total (this index is always defined)

b. if you want the power spectrum relevant for galaxy or halos, given by P_cb if there is non-cold-dark-matter (e.g. massive neutrinos) and to P_m otherwise, pass in input pfo->index_pk_cluster (this index is always defined)

c. there is another possible syntax (use it only if you know what you are doing): if pfo->has_pk_m == TRUE you may pass pfo->index_pk_m to get P_m if pfo->has_pk_cb == TRUE you may pass pfo->index_pk_cb to get P_cb

Output format:

out_pk = P(k)
out_pk_ic[diagonal] = P_ic(k)
out_pk_ic[non-diagonal] = P_icxic(k)

Parameters

pba	Input: pointer to background structure
ppm	Input: pointer to primordial structure
pfo	Input: pointer to fourier structure
pk_output	Input: linear or nonlinear
k	Input: wavenumber in 1/Mpc
z	Input: redshift
index_pk	Input: index of pk type (_m, _cb)
out_pk	Output: pointer to P
out_pk_ic	Ouput: P_ic returned as out_pk_ic_l[index_ic1_ic2]

Returns: the error status

preliminary: check whether we need the decomposition into contributions from each initial condition
first step: check that k is in valid range [0:kmax] (the test for z will be done when calling fourier_pk_linear_at_z())
deal with case k = 0 for which P(k) is set to zero (this non-physical result can be useful for interpolations)
deal with 0 < k <= kmax
deal with standard case kmin <= k <= kmax

--> First, get P(k) at the right z (in logarithmic format for more accurate interpolation, and then convert to linear format)

--> interpolate total spectrum

--> convert from logarithmic to linear format

--> interpolate each ic component

--> convert each ic component from logarithmic to linear format

--> deal with case 0 < k < kmin that requires extrapolation P(k) = [some number] * k * P_primordial(k) so P(k) = P(kmin) * (k P_primordial(k)) / (kmin P_primordial(kmin)) (note that the result is accurate only if kmin is such that [a0 kmin] << H0)

This is accurate for the synchronous gauge; TODO: write newtonian gauge case. Also, In presence of isocurvature modes, we assumes for simplicity that the mode with index_ic1_ic2=0 dominates at small k: exact treatment should be written if needed.

--> First, get P(k) at the right z (in linear format)

◆ fourier_pks_at_kvec_and_zvec()

int fourier_pks_at_kvec_and_zvec	(	struct background *	pba,
		struct fourier *	pfo,
		enum pk_outputs	pk_output,
		double *	kvec,
		int	kvec_size,
		double *	zvec,
		int	zvec_size,
		double *	out_pk,
		double *	out_pk_cb
	)

Return the P(k,z) for a grid of (k_i,z_j) passed in input, for all available pk types (_m, _cb), either linear or nonlinear depending on input.

If there are several initial conditions, this function is not designed to return individual contributions.

The main goal of this routine is speed. Unlike fourier_pk_at_k_and_z(), it performs no extrapolation when an input k_i falls outside the pre-computed range [kmin,kmax]: in that case, it just returns P(k,z)=0 for such a k_i

Parameters

pba	Input: pointer to background structure
pfo	Input: pointer to fourier structure
pk_output	Input: pk_linear or pk_nonlinear
kvec	Input: array of wavenumbers in ascending order (in 1/Mpc)
kvec_size	Input: size of array of wavenumbers
zvec	Input: array of redshifts in arbitrary order
zvec_size	Input: size of array of redshifts
out_pk	Output: P(k_i,z_j) for total matter (if available) in Mpc**3
out_pk_cb	Output: P_cb(k_i,z_j) for cdm+baryons (if available) in Mpc**3

Returns: the error status

Summary:

define local variables
Allocate arrays
Construct table of log(P(k_n,z_j)) for pre-computed wavenumbers but requested redshifts:
Spline it for interpolation along k
Construct ln(kvec):
Loop over first k values. If k<kmin, fill output with zeros. If not, go to next step.
Deal with case kmin<=k<=kmax. For better performance, do not loop through kvec, but through pre-computed k values.

--> Loop through k_i's that fall in interval [k_n,k_n+1]

--> for each of them, perform spine interpolation

Loop over possible remaining k values with k > kmax, to fill output with zeros.

◆ fourier_pk_tilt_at_k_and_z()

int fourier_pk_tilt_at_k_and_z	(	struct background *	pba,
		struct primordial *	ppm,
		struct fourier *	pfo,
		enum pk_outputs	pk_output,
		double	k,
		double	z,
		int	index_pk,
		double *	pk_tilt
	)

Return the logarithmic slope of P(k,z) for a given (k,z), a given pk type (_m, _cb) (computed with linear P_L if pk_output = pk_linear, nonlinear P_NL if pk_output = pk_nonlinear)

Parameters

pba	Input: pointer to background structure
ppm	Input: pointer to primordial structure
pfo	Input: pointer to fourier structure
pk_output	Input: linear or nonlinear
k	Input: wavenumber in 1/Mpc
z	Input: redshift
index_pk	Input: index of pk type (_m, _cb)
pk_tilt	Output: logarithmic slope of P(k,z)

Returns: the error status

◆ fourier_sigmas_at_z()

int fourier_sigmas_at_z	(	struct precision *	ppr,
		struct background *	pba,
		struct fourier *	pfo,
		double	R,
		double	z,
		int	index_pk,
		enum out_sigmas	sigma_output,
		double *	result
	)

This routine computes the variance of density fluctuations in a sphere of radius R at redshift z, sigma(R,z), or other similar derived quantitites, for one given pk type (_m, _cb).

The integral is performed until the maximum value of k_max defined in the perturbation module. Here there is not automatic checking that k_max is large enough for the result to be well converged. E.g. to get an accurate sigma8 at R = 8 Mpc/h, the user should pass at least about P_k_max_h/Mpc = 1.

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
pfo	Input: pointer to fourier structure
R	Input: radius in Mpc
z	Input: redshift
index_pk	Input: type of pk (_m, _cb)
sigma_output	Input: quantity to be computed (sigma, sigma', ...)
result	Output: result

Returns: the error status

allocate temporary array for P(k,z) as a function of k
get P(k,z) as a function of k, for the right z
spline it along k
calll the function computing the sigmas
free allocated arrays

◆ fourier_k_nl_at_z()

int fourier_k_nl_at_z	(	struct background *	pba,
		struct fourier *	pfo,
		double	z,
		double *	k_nl,
		double *	k_nl_cb
	)

Return the value of the non-linearity wavenumber k_nl for a given redshift z

Parameters

pba	Input: pointer to background structure
pfo	Input: pointer to fourier structure
z	Input: redshift
k_nl	Output: k_nl value
k_nl_cb	Ouput: k_nl value of the cdm+baryon part only, if there is ncdm

Returns: the error status

convert input redshift into a conformal time
interpolate the precomputed k_nl array at the needed valuetime
if needed, do the same for the baryon part only

◆ fourier_init()

int fourier_init	(	struct precision *	ppr,
		struct background *	pba,
		struct thermodynamics *	pth,
		struct perturbations *	ppt,
		struct primordial *	ppm,
		struct fourier *	pfo
	)

Initialize the fourier structure, and in particular the nl_corr_density and k_nl interpolation tables.

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
pth	Input: pointer to therodynamics structure
ppt	Input: pointer to perturbation structure
ppm	Input: pointer to primordial structure
pfo	Input/Output: pointer to initialized fourier structure

Returns: the error status

Do we want to compute P(k,z)? Propagate the flag has_pk_matter from the perturbations structure to the fourier structure
preliminary tests

--> This module only makes sense for dealing with scalar perturbations, so it should do nothing if there are no scalars

--> Nothing to be done if we don't want the matter power spectrum

--> check applicability of Halofit and HMcode

define indices in fourier structure (and allocate some arrays in the structure)
get the linear power spectrum at each time

--> loop over required pk types (_m, _cb)

--> get the linear power spectrum for this time and this type

--> if interpolation of $P(k,\tau)$ will be needed (as a function of tau), compute array of second derivatives in view of spline interpolation

compute and store sigma8 (variance of density fluctuations in spheres of radius 8/h Mpc at z=0, always computed by convention using the linear power spectrum)
get the non-linear power spectrum at each time

--> First deal with the case where non non-linear corrections requested

--> Then go through common preliminary steps to the HALOFIT and HMcode methods

--> allocate temporary arrays for spectra at each given time/redshift

--> Then go through preliminary steps specific to HMcode

--> Loop over decreasing time/growing redhsift. For each time/redshift, compute P_NL(k,z) using either Halofit or HMcode

--> fill the array of nonlinear power spectra (only if we are at a late time where P(k) and T(k) are supposed to be stored, i.e., such that z(tau < z_max_pk)

--> spline the array of nonlinear power spectrum

--> free the nonlinear workspace

if the nl_method could not be identified

◆ fourier_free()

int fourier_free ( struct fourier * pfo )

Free all memory space allocated by fourier_init().

Parameters

pfo	Input: pointer to fourier structure (to be freed)

Returns: the error status

◆ fourier_indices()

int fourier_indices	(	struct precision *	ppr,
		struct background *	pba,
		struct perturbations *	ppt,
		struct primordial *	ppm,
		struct fourier *	pfo
	)

Define indices in the fourier structure, and when possible, allocate arrays in this structure given the index sizes found here

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
ppt	Input: pointer to perturbation structure
ppm	Input: pointer to primordial structure
pfo	Input/Output: pointer to fourier structure

Returns: the error status

define indices for initial conditions (and allocate related arrays)
define flags indices for pk types (_m, _cb). Note: due to some dependencies in HMcode, when pfo->index_pk_cb exists, it must come first (e.g. the calculation of the non-linear P_m depends on sigma_cb so the cb-related quantitites must be evaluated first)
get list of k values
get list of tau values
given previous indices, we can allocate the array of linear power spectrum values
if interpolation of $P(k,\tau)$ will be needed (as a function of tau), compute also the array of second derivatives in view of spline interpolation
array of sigma8 values
if non-linear computations needed, allocate array of non-linear correction ratio R_nl(k,z), k_nl(z) and P_nl(k,z) for each P(k) type

◆ fourier_get_k_list()

int fourier_get_k_list	(	struct precision *	ppr,
		struct perturbations *	ppt,
		struct fourier *	pfo
	)

Copy list of k from perturbation module, and extended it if necessary to larger k for extrapolation (currently this extrapolation is required only by HMcode)

Parameters

ppr	Input: pointer to precision structure
ppt	Input: pointer to perturbation structure
pfo	Input/Output: pointer to fourier structure

Returns: the error status

if k extrapolation necessary, compute number of required extra values
otherwise, same number of values as in perturbation module
allocate array of k
fill array of k (not extrapolated)
fill additional values of k (extrapolated)

◆ fourier_get_tau_list()

int fourier_get_tau_list	(	struct perturbations *	ppt,
		struct fourier *	pfo
	)

Copy list of tau from perturbation module

Parameters

ppt	Input: pointer to perturbation structure
pfo	Input/Output: pointer to fourier structure

Returns: the error status

-> for linear calculations: only late times are considered, given the value z_max_pk inferred from the ionput

-> for non-linear calculations: we wills store a correction factor for all times

◆ fourier_get_source()

int fourier_get_source	(	struct background *	pba,
		struct perturbations *	ppt,
		struct fourier *	pfo,
		int	index_k,
		int	index_ic,
		int	index_tp,
		int	index_tau,
		double **	sources,
		double *	source
	)

Get sources for a given wavenumber (and for a given time, type, ic, mode...) either directly from precomputed valkues (computed ain perturbation module), or by analytic extrapolation

Parameters

pba	Input: pointer to background structure
ppt	Input: pointer to perturbation structure
pfo	Input: pointer to fourier structure
index_k	Input: index of required k value
index_ic	Input: index of required ic value
index_tp	Input: index of required tp value
index_tau	Input: index of required tau value
sources	Input: array containing the original sources
source	Output: desired value of source

Returns: the error status

use precomputed values
extrapolate

--> Get last source and k, which are used in (almost) all methods

--> Get previous source and k, which are used in best methods

--> Extrapolate by assuming the source to vanish Has terrible discontinuity

--> Extrapolate starting from the maximum value, assuming growth ~ ln(k) Has a terrible bend in log slope, discontinuity only in derivative

--> Extrapolate starting from the maximum value, assuming growth ~ ln(k) Here we use k in h/Mpc instead of 1/Mpc as it is done in the CAMB implementation of HMcode Has a terrible bend in log slope, discontinuity only in derivative

--> Extrapolate assuming source ~ ln(a*k) where a is obtained from the data at k_0 Mostly continuous derivative, quite good

--> Extrapolate assuming source ~ ln(e+a*k) where a is estimated like is done in original HMCode

--> If the user has a complicated model and wants to interpolate differently, they can define their interpolation here and switch to using it instead

◆ fourier_pk_linear()

int fourier_pk_linear	(	struct background *	pba,
		struct perturbations *	ppt,
		struct primordial *	ppm,
		struct fourier *	pfo,
		int	index_pk,
		int	index_tau,
		int	k_size,
		double *	lnpk,
		double *	lnpk_ic
	)

This routine computes all the components of the matter power spectrum P(k), given the source functions and the primordial spectra, at a given time within the pre-computed table of sources (= Fourier transfer functions) of the perturbation module, for a given type (total matter _m or baryon+CDM _cb), and for the same array of k values as in the pre-computed table.

If the input array of k values pfo->ln_k contains wavemumbers larger than those of the pre-computed table, the sources will be extrapolated analytically.

On the opther hand, if the primordial spectrum has sharp features and needs to be sampled on a finer grid than the sources, this function has to be modified to capture the features.

There are two output arrays, because we consider:

the total matter (_m) or CDM+baryon (_cb) power spectrum
in the quantitites labelled _ic, the splitting of one of these spectra in different modes for different initial conditions. If the pointer ln_pk_ic is NULL in input, the function will ignore this part; thus, to get the result, one should allocate the array before calling the function. Then the convention is the following:

– the index_ic1_ic2 labels ordered pairs (index_ic1, index_ic2) (since the primordial spectrum is symmetric in (index_ic1, index_ic2)).

– for diagonal elements (index_ic1 = index_ic2) this arrays contains ln[P(k)] where P(k) is positive by construction.

– for non-diagonal elements this arrays contains the k-dependent cosine of the correlation angle, namely P(k)_(index_ic1, index_ic2)/sqrt[P(k)_index_ic1 P(k)_index_ic2]. E.g. for fully correlated or anti-correlated initial conditions, this non-diagonal element is independent on k, and equal to +1 or -1.

Parameters

pba	Input: pointer to background structure
ppt	Input: pointer to perturbation structure
ppm	Input: pointer to primordial structure
pfo	Input: pointer to fourier structure
index_pk	Input: index of required P(k) type (_m, _cb)
index_tau	Input: index of time
k_size	Input: wavenumber array size
lnpk	Output: log of matter power spectrum for given type/time, for all wavenumbers
lnpk_ic	Output: log of matter power spectrum for given type/time, for all wavenumbers and initial conditions

Returns: the error status

allocate temporary vector where the primordial spectrum will be stored
loop over k values

--> get primordial spectrum

--> initialize a local variable for P_m(k) and P_cb(k) to zero

--> here we recall the relations relevant for the nomalization fo the power spectrum: For adiabatic modes, the curvature primordial spectrum thnat we just read was: P_R(k) = 1/(2pi^2) k^3 < R R > Thus the primordial curvature correlator is given by: < R R > = (2pi^2) k^-3 P_R(k) So the delta_m correlator reads: P(k) = < delta_m delta_m > = (source_m)^2 < R R > = (2pi^2) k^-3 (source_m)^2 P_R(k)

For isocurvature or cross adiabatic-isocurvature parts, one would just replace one or two 'R' by 'S_i's

--> get contributions to P(k) diagonal in the initial conditions

--> get contributions to P(k) non-diagonal in the initial conditions

◆ fourier_sigmas()

int fourier_sigmas	(	struct fourier *	pfo,
		double	R,
		double *	lnpk_l,
		double *	ddlnpk_l,
		int	k_size,
		double	k_per_decade,
		enum out_sigmas	sigma_output,
		double *	result
	)

Calculate intermediate quantities for hmcode (sigma, sigma', ...) for a given scale R and a given input P(k).

This function has several differences w.r.t. the standard external function non_linear_sigma (format of input, of output, integration stepsize, management of extrapolation at large k, ...) and is overall more precise for sigma(R).

Parameters

pfo	Input: pointer to fourier structure
R	Input: scale at which to compute sigma
lnpk_l	Input: array of ln(P(k))
ddlnpk_l	Input: its spline along k
k_size	Input: dimension of array lnpk_l, normally pfo->k_size, but inside hmcode it its increased by extrapolation to pfo->k_extra_size
k_per_decade	Input: logarithmic step for the integral (recommended: pass ppr->sigma_k_per_decade)
sigma_output	Input: quantity to be computed (sigma, sigma', ...)
result	Output: result

Returns: the error status

allocate temporary array for an integral over y(x)
fill the array with values of k and of the integrand
spline the integrand
integrate
preperly normalize the final result
free allocated array

◆ fourier_sigma_at_z()

int fourier_sigma_at_z	(	struct background *	pba,
		struct fourier *	pfo,
		double	R,
		double	z,
		int	index_pk,
		double	k_per_decade,
		double *	result
	)

This routine computes the variance of density fluctuations in a sphere of radius R at redshift z, sigma(R,z) for one given pk type (_m, _cb).

Try to use instead fourier_sigmas_at_z(). This function is just maintained for compatibility with the deprecated function harmonic_sigma()

The integral is performed until the maximum value of k_max defined in the perturbation module. Here there is not automatic checking that k_max is large enough for the result to be well converged. E.g. to get an accurate sigma8 at R = 8 Mpc/h, the user should pass at least about P_k_max_h/Mpc = 1.

Parameters

pba	Input: pointer to background structure
pfo	Input: pointer to fourier structure
R	Input: radius in Mpc
z	Input: redshift
index_pk	Input: type of pk (_m, _cb)
k_per_decade	Input: logarithmic step for the integral (recommended: pass ppr->sigma_k_per_decade)
result	Output: result

Returns: the error status

allocate temporary array for P(k,z) as a function of k
get P(k,z) as a function of k, for the right z
spline it along k
calll the function computing the sigmas
free allocated arrays

◆ fourier_halofit()

int fourier_halofit	(	struct precision *	ppr,
		struct background *	pba,
		struct perturbations *	ppt,
		struct primordial *	ppm,
		struct fourier *	pfo,
		int	index_pk,
		double	tau,
		double *	pk_nl,
		double *	lnpk_l,
		double *	ddlnpk_l,
		double *	k_nl,
		short *	nl_corr_not_computable_at_this_k
	)

Calculation of the nonlinear matter power spectrum with Halofit (includes Takahashi 2012 + Bird 2013 revisions).

At high redshift it is possible that the non-linear corrections are so small that they can be computed only by going to very large wavenumbers. Thius, for some combination of (z, k_max), the calculation is not possible. In this case a FALSE will be returned in the flag halofit_found_k_max.

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
ppt	Input: pointer to perturbation structure
ppm	Input: pointer to primordial structure
pfo	Input: pointer to fourier structure
index_pk	Input: index of component are we looking at (total matter or cdm+baryons?)
tau	Input: conformal time at which we want to do the calculation
pk_nl	Output: non linear spectrum at the relevant time
lnpk_l	Input: array of log(P(k)_linear)
ddlnpk_l	Input: array of second derivative of log(P(k)_linear) wrt k, for spline interpolation
k_nl	Output: non-linear wavenumber
nl_corr_not_computable_at_this_k	Ouput: flag concerning the status of the calculation (TRUE if not possible)

Returns: the error status

Determine non linear ratios (from pk)

◆ fourier_halofit_integrate()

int fourier_halofit_integrate	(	struct fourier *	pfo,
		double *	integrand_array,
		int	integrand_size,
		int	ia_size,
		int	index_ia_k,
		int	index_ia_pk,
		int	index_ia_sum,
		int	index_ia_ddsum,
		double	R,
		enum halofit_integral_type	type,
		double *	sum
	)

Internal routione of Halofit. In original Halofit, this is equivalent to the function wint(). It performs convolutions of the linear spectrum with two window functions.

Parameters

pfo	Input: pointer to non linear structure
integrand_array	Input: array with k, P_L(k) values
integrand_size	Input: one dimension of that array
ia_size	Input: other dimension of that array
index_ia_k	Input: index for k
index_ia_pk	Input: index for pk
index_ia_sum	Input: index for the result
index_ia_ddsum	Input: index for its spline
R	Input: radius
type	Input: which window function to use
sum	Output: result of the integral

Returns: the error status

◆ fourier_hmcode()

int fourier_hmcode	(	struct precision *	ppr,
		struct background *	pba,
		struct perturbations *	ppt,
		struct primordial *	ppm,
		struct fourier *	pfo,
		int	index_pk,
		int	index_tau,
		double	tau,
		double *	pk_nl,
		double **	lnpk_l,
		double **	ddlnpk_l,
		double *	k_nl,
		short *	nl_corr_not_computable_at_this_k,
		struct fourier_workspace *	pnw
	)

Computes the nonlinear correction on the linear power spectrum via the method presented in Mead et al. 1505.07833

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
ppt	Input: pointer to perturbation structure
ppm	Input: pointer to primordial structure
pfo	Input: pointer to fourier structure
index_pk	Input: index of the pk type, either index_m or index_cb
index_tau	Input: index of tau, at which to compute the nl correction
tau	Input: tau, at which to compute the nl correction
pk_nl	Output:nonlinear power spectrum
lnpk_l	Input: logarithm of the linear power spectrum for both index_m and index_cb
ddlnpk_l	Input: spline of the logarithm of the linear power spectrum for both index_m and index_cb
nl_corr_not_computable_at_this_k	Ouput: was the computation doable?
k_nl	Output: nonlinear scale for index_m and index_cb
pnw	Input/Output: pointer to nonlinear workspace

Returns: the error status

include precision parameters that control the number of entries in the growth and sigma tables

Compute background quantitites today

If index_pk_cb, choose Omega0_cb as the matter density parameter. If index_pk_m, choose Omega0_cbn as the matter density parameter.

Call all the relevant background parameters at this tau

Test whether pk_cb has to be taken into account (only if we have massive neutrinos)

Get sigma(R=8 Mpc/h), sigma_disp(R=0), sigma_disp(R=100 Mpc/h) and write them into pfo structure

Initialisation steps for the 1-Halo Power Integral

find nonlinear scales k_nl and r_nl and the effective spectral index n_eff

Calculate halo concentration-mass relation conc(mass) (Bullock et al. 2001)

Compute the nonlinear correction

◆ fourier_hmcode_workspace_init()

int fourier_hmcode_workspace_init	(	struct precision *	ppr,
		struct background *	pba,
		struct fourier *	pfo,
		struct fourier_workspace *	pnw
	)

allocate and fill arrays of nonlinear workspace (currently used only by HMcode)

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
pfo	Input: pointer to fourier structure
pnw	Output: pointer to nonlinear workspace

Returns: the error status

allocate arrays of the nonlinear workspace
fill table with scale independent growth factor

◆ fourier_hmcode_workspace_free()

int fourier_hmcode_workspace_free	(	struct fourier *	pfo,
		struct fourier_workspace *	pnw
	)

deallocate arrays in the nonlinear worksapce (currently used only by HMcode)

Parameters

pfo	Input: pointer to fourier structure
pnw	Input: pointer to nonlinear workspace

Returns: the error status

◆ fourier_hmcode_dark_energy_correction()

int fourier_hmcode_dark_energy_correction	(	struct precision *	ppr,
		struct background *	pba,
		struct fourier *	pfo,
		struct fourier_workspace *	pnw
	)

set the HMcode dark energy correction (if w is not -1)

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
pfo	Input: pointer to fourier structure
pnw	Output: pointer to nonlinear workspace

Returns: the error status

if there is dynamical Dark Energy (w is not -1) modeled as a fluid
otherwise, we assume no dynamical Dark Energy (w is -1)

◆ fourier_hmcode_baryonic_feedback()

int fourier_hmcode_baryonic_feedback ( struct fourier * pfo )

set the HMcode baryonic feedback parameters according to the chosen feedback model

Parameters

pfo	Output: pointer to fourier structure

Returns: the error status

◆ fourier_hmcode_fill_sigtab()

int fourier_hmcode_fill_sigtab	(	struct precision *	ppr,
		struct background *	pba,
		struct perturbations *	ppt,
		struct primordial *	ppm,
		struct fourier *	pfo,
		int	index_tau,
		double *	lnpk_l,
		double *	ddlnpk_l,
		struct fourier_workspace *	pnw
	)

Function that fills pnw->rtab, pnw->stab and pnw->ddstab with (r, sigma, ddsigma) logarithmically spaced in r. Called by fourier_init at for all tau to account for scale-dependant growth before fourier_hmcode is called

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
ppt	Input: pointer to perturbation structure
ppm	Input: pointer to primordial structure
pfo	Input: pointer to fourier structure
index_tau	Input: index of tau, at which to compute the nl correction
lnpk_l	Input: logarithm of the linear power spectrum for either index_m or index_cb
ddlnpk_l	Input: spline of the logarithm of the linear power spectrum for either index_m or index_cb
pnw	Output: pointer to nonlinear workspace

Returns: the error status

◆ fourier_hmcode_fill_growtab()

int fourier_hmcode_fill_growtab	(	struct precision *	ppr,
		struct background *	pba,
		struct fourier *	pfo,
		struct fourier_workspace *	pnw
	)

Function that fills pnw->tautable and pnw->growtable with (tau, D(tau)) linearly spaced in scalefactor a. Called by fourier_init at before the loop over tau

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure (will provide the scale independent growth factor)
pfo	Input/Output: pointer to fourier structure
pnw	Output: pointer to nonlinear workspace

Returns: the error status

◆ fourier_hmcode_growint()

int fourier_hmcode_growint	(	struct precision *	ppr,
		struct background *	pba,
		struct fourier *	pfo,
		double	a,
		double	w0,
		double	wa,
		double *	growth
	)

This function finds the scale independent growth factor by integrating the approximate relation d(lnD)/d(lna) = Omega_m(z)^gamma by Linder & Cahn 2007

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
pfo	Input: pointer to fourier structure
a	Input: scalefactor
w0	Input: dark energy equation of state today
wa	Input: dark energy equation of state varying with a: w=w0+(1-a)wa
growth	Output: scale independent growth factor at a

Returns: the error status

◆ fourier_hmcode_window_nfw()

int fourier_hmcode_window_nfw	(	struct fourier *	pfo,
		double	k,
		double	rv,
		double	c,
		double *	window_nfw
	)

This is the fourier transform of the NFW density profile.

Parameters

pfo	Input: pointer to fourier structure
k	Input: wave vector
rv	Input: virial radius
c	Input: concentration = rv/rs (with scale radius rs)
window_nfw	Output: Window Function of the NFW profile

Returns: the error status

◆ fourier_hmcode_halomassfunction()

int fourier_hmcode_halomassfunction	(	double	nu,
		double *	hmf
	)

This is the Sheth-Tormen halo mass function (1999, MNRAS, 308, 119)

Parameters

nu	Input: the $\nu$ parameter that depends on the halo mass via $\nu(M) = \delta_c/\sigma(M)$
hmf	Output: Value of the halo mass function at this $\nu$

Returns: the error status

◆ fourier_hmcode_sigma8_at_z()

int fourier_hmcode_sigma8_at_z	(	struct background *	pba,
		struct fourier *	pfo,
		double	z,
		double *	sigma_8,
		double *	sigma_8_cb,
		struct fourier_workspace *	pnw
	)

Compute sigma8(z)

Parameters

pba	Input: pointer to background structure
pfo	Input: pointer to fourier structure
z	Input: redshift
sigma_8	Output: sigma8(z)
sigma_8_cb	Output: sigma8_cb(z)
pnw	Output: pointer to nonlinear workspace

Returns: the error status

◆ fourier_hmcode_sigmadisp_at_z()

int fourier_hmcode_sigmadisp_at_z	(	struct background *	pba,
		struct fourier *	pfo,
		double	z,
		double *	sigma_disp,
		double *	sigma_disp_cb,
		struct fourier_workspace *	pnw
	)

Compute sigmadisp(z)

Parameters

pba	Input: pointer to background structure
pfo	Input: pointer to fourier structure
z	Input: redshift
sigma_disp	Output: sigmadisp(z)
sigma_disp_cb	Output: sigmadisp_cb(z)
pnw	Output: pointer to nonlinear workspace

Returns: the error status

◆ fourier_hmcode_sigmadisp100_at_z()

int fourier_hmcode_sigmadisp100_at_z	(	struct background *	pba,
		struct fourier *	pfo,
		double	z,
		double *	sigma_disp_100,
		double *	sigma_disp_100_cb,
		struct fourier_workspace *	pnw
	)

Compute sigmadisp100(z)

Parameters

pba	Input: pointer to background structure
pfo	Input: pointer to fourier structure
z	Input: redshift
sigma_disp_100	Output: sigmadisp100(z)
sigma_disp_100_cb	Output: sigmadisp100_cb(z)
pnw	Output: pointer to nonlinear workspace

Returns: the error status

◆ fourier_hmcode_sigmaprime_at_z()

int fourier_hmcode_sigmaprime_at_z	(	struct background *	pba,
		struct fourier *	pfo,
		double	z,
		double *	sigma_prime,
		double *	sigma_prime_cb,
		struct fourier_workspace *	pnw
	)

Compute sigma'(z)

Parameters

pba	Input: pointer to background structure
pfo	Input: pointer to fourier structure
z	Input: redshift
sigma_prime	Output: sigma'(z)
sigma_prime_cb	Output: sigma'_cb(z)
pnw	Output: pointer to nonlinear workspace

Returns: the error status

Functions

Detailed Description

Function Documentation

◆ fourier_pk_at_z()

◆ fourier_pk_at_k_and_z()

◆ fourier_pks_at_kvec_and_zvec()

◆ fourier_pk_tilt_at_k_and_z()

◆ fourier_sigmas_at_z()

◆ fourier_k_nl_at_z()

◆ fourier_init()

◆ fourier_free()

◆ fourier_indices()

◆ fourier_get_k_list()

◆ fourier_get_tau_list()

◆ fourier_get_source()

◆ fourier_pk_linear()

◆ fourier_sigmas()

◆ fourier_sigma_at_z()

◆ fourier_halofit()

◆ fourier_halofit_integrate()

◆ fourier_hmcode()

◆ fourier_hmcode_workspace_init()

◆ fourier_hmcode_workspace_free()

◆ fourier_hmcode_dark_energy_correction()

◆ fourier_hmcode_baryonic_feedback()

◆ fourier_hmcode_fill_sigtab()

◆ fourier_hmcode_fill_growtab()

◆ fourier_hmcode_growint()

◆ fourier_hmcode_window_nfw()

◆ fourier_hmcode_halomassfunction()

◆ fourier_hmcode_sigma8_at_z()

◆ fourier_hmcode_sigmadisp_at_z()

◆ fourier_hmcode_sigmadisp100_at_z()

◆ fourier_hmcode_sigmaprime_at_z()