#include "fourier.h"
#include "halofit.h"
#include "hmcode.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 primordial ppm, 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_pk_analytic_nowiggle (struct precision ppr, struct background pba, struct primordial ppm, struct fourier pfo)

int	fourier_wnw_split (struct precision ppr, struct background pba, struct primordial ppm, struct fourier pfo)

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)

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 $ C_l$ 's and not just $ P(k)$ . 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, nowiggle linear spectrum if pk_output = pk_numerical_nowiggle, analytic approximation to nowiggle linear spectrum if pk_output = pk_analytic_nowiggle)

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.

In the pk_analytic_nowiggle case, the overall normalisation of the spectrum is currently arbitrary and independent of redhsift.

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, nonlinear, nowiggle...
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. The pk_analytic_nowiggle can also be computed here because, in the current implementation, it is only computed at z=0.
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

--> we requested P_nl(k) at tau where P_l is computed but the NL correction is NOT -> Return P_l

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

--> interpolate P_l_nw(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, nowiggle linear spectrum if pk_output = pk_numerical_nowiggle, analytic approximation to linear nowiggle spectrum if pk_output = pk_analytic_nowiggle)

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, nonlinear, nowiggle...
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, pk_nonlinear, nowiggle...
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, nowiggle linear spectrum if pk_output = pk_numerical_nowiggle, analytic approximation to linear nowiggle spectrum if pk_output = pk_analytic_nowiggle)

Parameters

pba	Input: pointer to background structure
ppm	Input: pointer to primordial structure
pfo	Input: pointer to fourier structure
pk_output	Input: linear, nonlinear, nowiggle...
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

--> one more call to get the linear power spectrum extrapolated up to very large k, for this time and type, but ignoring the case of multiple initial conditions. Result stored in ln_pk_l_extra (different from non-extrapolated ln_pk_l)

if interpolation of $P(k,\tau)$ will be needed (as a function of tau), compute array of second derivatives in view of spline interpolation
get the analytic_nowiggle power spectrum
get the dewiggled power spectrum at each time in ln_tau
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

--> spline the array of nonlinear power spectrum (only above the first index where it could be computed)

--> 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
do we we want to compute and store the analytic_nowiggle power spectrum? This flag was already set by the user in input.c. If we need HMcode, overwrite the user request and compute it. Will soon do the same for Oneloop.
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
allocate arrays for P_nw(k,z) and its splines in log(tau)

◆ fourier_get_k_list()

int fourier_get_k_list	(	struct precision *	ppr,
		struct primordial *	ppm,
		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
ppm	Input: pointer to primordial 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_pk_analytic_nowiggle()

int fourier_pk_analytic_nowiggle	(	struct precision *	ppr,
		struct background *	pba,
		struct primordial *	ppm,
		struct fourier *	pfo
	)

Compute a smooth analytic approximation to the matter power spectrum today. Store the results in the fourier structure.

Parameters

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

Returns: the error status

◆ fourier_wnw_split()

int fourier_wnw_split	(	struct precision *	ppr,
		struct background *	pba,
		struct primordial *	ppm,
		struct fourier *	pfo
	)

Compute the decomposition of the linear power spectrum into a wiggly and a non-wiggly part. Store the results in the fourier structure.

Parameters

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

Returns: the error status

spline the nowiggle spectrum with respect to time

◆ 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

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_pk_analytic_nowiggle()

◆ fourier_wnw_split()

◆ fourier_sigmas()

◆ fourier_sigma_at_z()