#include "perturbations.h"

Include dependency graph for perturbations.c:

Functions
int	perturbations_sources_at_tau (struct perturbations ppt, int index_md, int index_ic, int index_tp, double tau, double psource_at_tau)

int	perturbations_sources_at_z (struct background pba, struct perturbations ppt, int index_md, int index_ic, int index_tp, double z, double *psource_at_z)

int	perturbations_sources_at_k_and_z (struct background pba, struct perturbations ppt, int index_md, int index_ic, int index_tp, double k, double z, double *psource_at_k_and_z)

int	perturbations_output_data_at_z (struct background pba, struct perturbations ppt, enum file_format output_format, double z, int number_of_titles, double *data)

int	perturbations_output_data_at_index_tau (struct background pba, struct perturbations ppt, enum file_format output_format, int index_tau, int number_of_titles, double *data)

int	perturbations_output_data (struct background pba, struct perturbations ppt, enum file_format output_format, double tkfull, int number_of_titles, double data)

int	perturbations_output_titles (struct background pba, struct perturbations ppt, enum file_format output_format, char titles[_MAXTITLESTRINGLENGTH_])

int	perturbations_output_firstline_and_ic_suffix (struct perturbations *ppt, int index_ic, char first_line[_LINE_LENGTH_MAX_], char ic_suffix[_SUFFIXNAMESIZE_])

int	perturbations_init (struct precision ppr, struct background pba, struct thermodynamics pth, struct perturbations ppt)

int	perturbations_free_input (struct perturbations *ppt)

int	perturbations_free (struct perturbations *ppt)

int	perturbations_indices (struct precision ppr, struct background pba, struct thermodynamics pth, struct perturbations ppt)

int	perturbations_timesampling_for_sources (struct precision ppr, struct background pba, struct thermodynamics pth, struct perturbations ppt)

int	perturbations_get_k_list (struct precision ppr, struct background pba, struct thermodynamics pth, struct perturbations ppt)

int	perturbations_workspace_init (struct precision ppr, struct background pba, struct thermodynamics pth, struct perturbations ppt, int index_md, struct perturbations_workspace *ppw)

int	perturbations_workspace_free (struct perturbations ppt, int index_md, struct perturbations_workspace ppw)

int	perturbations_solve (struct precision ppr, struct background pba, struct thermodynamics pth, struct perturbations ppt, int index_md, int index_ic, int index_k, struct perturbations_workspace *ppw)

int	perturbations_prepare_k_output (struct background pba, struct perturbations ppt)

int	perturbations_find_approximation_number (struct precision ppr, struct background pba, struct thermodynamics pth, struct perturbations ppt, int index_md, double k, struct perturbations_workspace ppw, double tau_ini, double tau_end, int interval_number, int *interval_number_of)

int	perturbations_find_approximation_switches (struct precision ppr, struct background pba, struct thermodynamics pth, struct perturbations ppt, int index_md, double k, struct perturbations_workspace ppw, double tau_ini, double tau_end, double precision, int interval_number, int interval_number_of, double interval_limit, int *interval_approx)

int	perturbations_vector_init (struct precision ppr, struct background pba, struct thermodynamics pth, struct perturbations ppt, int index_md, int index_ic, double k, double tau, struct perturbations_workspace ppw, int pa_old)

int	perturbations_vector_free (struct perturbations_vector *pv)

int	perturbations_initial_conditions (struct precision ppr, struct background pba, struct perturbations ppt, int index_md, int index_ic, double k, double tau, struct perturbations_workspace ppw)

int	perturbations_approximations (struct precision ppr, struct background pba, struct thermodynamics pth, struct perturbations ppt, int index_md, double k, double tau, struct perturbations_workspace *ppw)

int	perturbations_timescale (double tau, void parameters_and_workspace, double timescale, ErrorMsg error_message)

int	perturbations_einstein (struct precision ppr, struct background pba, struct thermodynamics pth, struct perturbations ppt, int index_md, double k, double tau, double y, struct perturbations_workspace ppw)

int	perturbations_total_stress_energy (struct precision ppr, struct background pba, struct thermodynamics pth, struct perturbations ppt, int index_md, double k, double y, struct perturbations_workspace ppw)

int	perturbations_sources (double tau, double y, double dy, int index_tau, void *parameters_and_workspace, ErrorMsg error_message)

int	perturbations_print_variables (double tau, double y, double dy, void *parameters_and_workspace, ErrorMsg error_message)

int	perturbations_derivs (double tau, double y, double dy, void *parameters_and_workspace, ErrorMsg error_message)

int	perturbations_tca_slip_and_shear (double y, void parameters_and_workspace, ErrorMsg error_message)

int	perturbations_rsa_delta_and_theta (struct precision ppr, struct background pba, struct thermodynamics pth, struct perturbations ppt, double k, double y, double a_prime_over_a, double pvecthermo, struct perturbations_workspace *ppw, ErrorMsg error_message)

int	perturbations_rsa_idr_delta_and_theta (struct precision ppr, struct background pba, struct thermodynamics pth, struct perturbations ppt, double k, double y, double a_prime_over_a, double pvecthermo, struct perturbations_workspace *ppw, ErrorMsg error_message)

Detailed Description

Documented perturbation module

Julien Lesgourgues, 23.09.2010

Deals with the perturbation evolution. This module has two purposes:

at the beginning; to initialize the perturbations, i.e. to integrate the perturbation equations, and store temporarily the terms contributing to the source functions as a function of conformal time. Then, to perform a few manipulations of these terms in order to infer the actual source functions $S^{X} (k, \tau)$ , and to store them as a function of conformal time inside an interpolation table.
at any time in the code; to evaluate the source functions at a given conformal time (by interpolating within the interpolation table).

Hence the following functions can be called from other modules:

perturbations_init() at the beginning (but after background_init() and thermodynamics_init())
perturbations_sources_at_tau() at any later time
perturbations_free() at the end, when no more calls to perturbations_sources_at_tau() are needed

Function Documentation

◆ perturbations_sources_at_tau()

int perturbations_sources_at_tau	(	struct perturbations *	ppt,
		int	index_md,
		int	index_ic,
		int	index_tp,
		double	tau,
		double *	psource_at_tau
	)

Source function $S^{X} (k, \tau)$ at a given conformal time tau.

Evaluate source functions at given conformal time tau by reading the pre-computed table and interpolating.

Parameters

ppt	Input: pointer to perturbation structure containing interpolation tables
index_md	Input: index of requested mode (for scalars, just pass ppt->index_md_scalars)
index_ic	Input: index of requested initial condition (for adiabatic, just pass ppt->index_ic_ad)
index_tp	Input: index of requested source function type
tau	Input: any value of conformal time
psource_at_tau	Output: vector (already allocated) of source function as a function of k, psource_at_tau[index_k]

Returns: the error status

Summary:

define local variables
If we have defined a z_max_pk > 0, then we have already an array of sources and of their second derivative with respect to time in the range 0 < z < z_max_pk, that can be used for an accurate spline interpolation at a given tau. Check whether we are in this situation and whether the value of tau is in the right range.
If yes, we do such a spline
otherwise, we just go for a quick linear interpolation. This is made available for developpers for completeness, but it is actually never used by the default version of CLASS

◆ perturbations_sources_at_z()

int perturbations_sources_at_z	(	struct background *	pba,
		struct perturbations *	ppt,
		int	index_md,
		int	index_ic,
		int	index_tp,
		double	z,
		double *	psource_at_z
	)

Source function $S^{X} (k, \tau)$ at a given redhsift z.

Evaluate source functions at given redhsift z by reading the pre-computed table and interpolating.

Parameters

pba	Input: pointer to background structure
ppt	Input: pointer to perturbation structure containing interpolation tables
index_md	Input: index of requested mode (for scalars, just pass ppt->index_md_scalars)
index_ic	Input: index of requested initial condition (for adiabatic, just pass ppt->index_ic_ad)
index_tp	Input: index of requested source function type
z	Input: any value of redshift
psource_at_z	Output: vector (already allocated) of source function as a function of k, psource_at_z[index_k]

Returns: the error status

◆ perturbations_sources_at_k_and_z()

int perturbations_sources_at_k_and_z	(	struct background *	pba,
		struct perturbations *	ppt,
		int	index_md,
		int	index_ic,
		int	index_tp,
		double	k,
		double	z,
		double *	psource_at_k_and_z
	)

Source function $S^{X} (k, \tau)$ at a given redhsift z and wavenumber k.

Evaluate source functions at given redhsift z and wavenumber k by reading the pre-computed table and interpolating.

Parameters

pba	Input: pointer to background structure
ppt	Input: pointer to perturbation structure containing interpolation tables
index_md	Input: index of requested mode (for scalars, just pass ppt->index_md_scalars)
index_ic	Input: index of requested initial condition (for adiabatic, just pass ppt->index_ic_ad)
index_tp	Input: index of requested source function type
k	Input: any value of wavenumber
z	Input: any value of redshift
psource_at_k_and_z	Output: pointer to the source function at (k,z) (so to just one number)

Returns: the error status

◆ perturbations_output_data_at_z()

int perturbations_output_data_at_z	(	struct background *	pba,
		struct perturbations *	ppt,
		enum file_format	output_format,
		double	z,
		int	number_of_titles,
		double *	data
	)

Function called by the output module or the wrappers, which returns the source functions $S^{X} (k, \tau)$ corresponding to densities ('dTk') and velocities ('vTk') at a given conformal time tau corresponding to the input redshift z.

Parameters

pba	Input: pointer to background structure
ppt	Input: pointer to perturbation structure
output_format	Input: choice of ordering and normalisation for the output quantities
z	Input: redshift
number_of_titles	Input: number of requested source functions (found in perturbations_output_titles)
data	Output: vector of all source functions for all k values and initial conditions (previously allocated with the right size)

Returns: the error status

allocate tkfull
compute for each k (if several ic's, compute it for each ic; if z_pk = 0, this is done by directly reading inside the pre-computed table; if not, this is done by interpolating the table at the correct value of tau.
store data
free tkfull

◆ perturbations_output_data_at_index_tau()

int perturbations_output_data_at_index_tau	(	struct background *	pba,
		struct perturbations *	ppt,
		enum file_format	output_format,
		int	index_tau,
		int	number_of_titles,
		double *	data
	)

Function called by the output module or the wrappers, which returns the source functions $S^{X} (k, \tau)$ corresponding to densities ('dTk') and velocities ('vTk') at a given conformal time tau corresponding to index index_tau in pre-computed table.

Parameters

pba	Input: pointer to background structure
ppt	Input: pointer to perturbation structure
output_format	Input: choice of ordering and normalisation for the output quantities
index_tau	Input: index pre-computed table ppt->ln_tau[index_tau]
number_of_titles	Input: number of requested source functions (found in perturbations_output_titles)
data	Output: vector of all source functions for all k values and initial conditions (previously allocated with the right size)

Returns: the error status

allocate and fill tkfull
store data
free tkfull

◆ perturbations_output_data()

int perturbations_output_data	(	struct background *	pba,
		struct perturbations *	ppt,
		enum file_format	output_format,
		double *	tkfull,
		int	number_of_titles,
		double *	data
	)

Function called by the output module or the wrappers, which returns the source functions $S^{X} (k, \tau)$ corresponding to densities ('dTk') and velocities ('vTk') in the correct order (matching that in perturbations_output_titles), given a vector containing the corresponding scalar source functions at a given time.

Parameters

pba	Input: pointer to background structure
ppt	Input: pointer to perturbation structure
output_format	Input: choice of ordering and normalisation for the output quantities
tkfull	Input: vector of scalar sources at given time, tk[(index_k * ppt->ic_size[index_md] + index_ic) * ppt->tp_size[index_md] + index_tp]
number_of_titles	Input: number of requested source functions (found in perturbations_output_titles)
data	Output: vector of all source functions for all k values and initial conditions (previously allocated with the right size)

Returns: the error status

store data

◆ perturbations_output_titles()

int perturbations_output_titles	(	struct background *	pba,
		struct perturbations *	ppt,
		enum file_format	output_format,
		char	titles[_MAXTITLESTRINGLENGTH_]
	)

Fill array of strings with the name of the requested 'mTk, vTk' functions (transfer functions as a function of wavenumber for fixed times).

Parameters

pba	Input: pointer to the background structure
ppt	Input: pointer to the perturbation structure
output_format	Input: flag for the format
titles	Output: name strings

Returns: the error status

◆ perturbations_output_firstline_and_ic_suffix()

int perturbations_output_firstline_and_ic_suffix	(	struct perturbations *	ppt,
		int	index_ic,
		char	first_line[_LINE_LENGTH_MAX_],
		char	ic_suffix[_SUFFIXNAMESIZE_]
	)

Fill strings that will be used when writing the transfer functions and the spectra in files (in the file names and in the comment at the beginning of each file).

Parameters

ppt	Input: pointer to the perturbation structure
index_ic	Input: index of the initial condition
first_line	Output: line of comment
ic_suffix	Output: suffix for the output file name

Returns: the error status

◆ perturbations_init()

int perturbations_init	(	struct precision *	ppr,
		struct background *	pba,
		struct thermodynamics *	pth,
		struct perturbations *	ppt
	)

Initialize the perturbations structure, and in particular the table of source functions.

Main steps:

given the values of the flags describing which kind of perturbations should be considered (modes: scalar/vector/tensor, initial conditions, type of source functions needed...), initialize indices and wavenumber list
define the time sampling for the output source functions
for each mode (scalar/vector/tensor): initialize the indices of relevant perturbations, integrate the differential system, compute and store the source functions.

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
pth	Input: pointer to thermodynamics structure
ppt	Output: Initialized perturbation structure

Returns: the error status

Summary:

define local variables
perform preliminary checks
initialize all indices and lists in perturbations structure using perturbations_indices()
define the common time sampling for all sources using perturbations_timesampling_for_sources()
if we want to store perturbations for given k values, write titles and allocate storage
create an array of workspaces in multi-thread case
loop over modes (scalar, tensors, etc). For each mode:
--> (a) create a workspace (one per thread in multi-thread case)
--> (b) initialize indices of vectors of perturbations with perturbations_indices_of_current_vectors()
--> (c) loop over initial conditions and wavenumbers; for each of them, evolve perturbations and compute source functions with perturbations_solve()
spline the source array with respect to the time variable

◆ perturbations_free_input()

int perturbations_free_input ( struct perturbations * ppt )

Free all memory space allocated by input.

Called by perturbations_free(), during shooting and if shooting failed

Parameters

ppt	Input: perturbation structure with input pointers to be freed

Returns: the error status

Here is the caller graph for this function:

◆ perturbations_free()

int perturbations_free ( struct perturbations * ppt )

Free all memory space allocated by perturbations_init().

To be called at the end of each run, only when no further calls to perturbations_sources_at_tau() are needed.

Parameters

ppt	Input: perturbation structure to be freed

Returns: the error status

Stuff related to perturbations output:

Free non-NULL pointers

Here is the call graph for this function:

◆ perturbations_indices()

int perturbations_indices	(	struct precision *	ppr,
		struct background *	pba,
		struct thermodynamics *	pth,
		struct perturbations *	ppt
	)

Initialize all indices and allocate most arrays in perturbations structure.

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
pth	Input: pointer to thermodynamics structure
ppt	Input/Output: Initialized perturbation structure

Returns: the error status

Summary:

define local variables
count modes (scalar, vector, tensor) and assign corresponding indices
allocate array of number of types for each mode, ppt->tp_size[index_md]
allocate array of number of initial conditions for each mode, ppt->ic_size[index_md]
allocate array of arrays of source functions for each mode, ppt->source[index_md]
initialize variables for the output of k values
initialization of all flags to false (will eventually be set to true later)
source flags and indices, for sources that all modes have in common (temperature, polarization, ...). For temperature, the term t2 is always non-zero, while other terms are non-zero only for scalars and vectors. For polarization, the term e is always non-zero, while the term b is only for vectors and tensors.
define k values with perturbations_get_k_list()
loop over modes. Initialize flags and indices which are specific to each mode.
(a) scalars
--> source flags and indices, for sources that are specific to scalars

gamma is not neccessary for converting output to Nbody gauge but is included anyway.

--> count scalar initial conditions (for scalars: ad, cdi, nid, niv; for tensors: only one) and assign corresponding indices
(b) vectors
--> source flags and indices, for sources that are specific to vectors
--> initial conditions for vectors
(c) tensors
--> source flags and indices, for sources that are specific to tensors
--> only one initial condition for tensors
(d) for each mode, allocate array of arrays of source functions for each initial conditions and wavenumber, (ppt->source[index_md])[index_ic][index_type]

◆ perturbations_timesampling_for_sources()

int perturbations_timesampling_for_sources	(	struct precision *	ppr,
		struct background *	pba,
		struct thermodynamics *	pth,
		struct perturbations *	ppt
	)

Define time sampling for source functions.

For each type, compute the list of values of tau at which sources will be sampled. Knowing the number of tau values, allocate all arrays of source functions.

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
pth	Input: pointer to thermodynamics structure
ppt	Input/Output: Initialized perturbation structure

Returns: the error status

Summary:

define local variables
allocate background/thermodynamics vectors
first, just count the number of sampling points in order to allocate the array containing all values
(a) if CMB requested, first sampling point = when the universe stops being opaque; otherwise, start sampling gravitational potential at recombination [however, if perturbed recombination is requested, we also need to start the system before recombination. Otherwise, the initial conditions for gas temperature and ionization fraction perturbations (delta_T = 1/3 delta_b, delta_x_e) are not valid].
(b) next sampling point = previous + ppr->perturbations_sampling_stepsize * timescale_source, where:
--> if CMB requested: timescale_source1 = $|g/\dot{g}| = |\dot{\kappa}-\ddot{\kappa}/\dot{\kappa}|^{-1}$ ; timescale_source2 = $|2\ddot{a}/a-(\dot{a}/a)^2|^{-1/2}$ (to sample correctly the late ISW effect; and timescale_source=1/(1/timescale_source1+1/timescale_source2); repeat till today.
--> if CMB not requested: timescale_source = 1/aH; repeat till today.
--> infer total number of time steps, ppt->tau_size
--> allocate array of time steps, ppt->tau_sampling[index_tau]
--> repeat the same steps, now filling the array with each tau value:
--> (b.1.) first sampling point = when the universe stops being opaque
--> (b.2.) next sampling point = previous + ppr->perturbations_sampling_stepsize * timescale_source, where timescale_source1 = $|g/\dot{g}| = |\dot{\kappa}-\ddot{\kappa}/\dot{\kappa}|^{-1}$ ; timescale_source2 = $|2\ddot{a}/a-(\dot{a}/a)^2|^{-1/2}$ (to sample correctly the late ISW effect; and timescale_source=1/(1/timescale_source1+1/timescale_source2); repeat till today. If CMB not requested: timescale_source = 1/aH; repeat till today.
last sampling point = exactly today
check the maximum redshift z_max_pk at which the Fourier transfer functions should be computable by interpolation. If it is equal to zero, only needs to be computed. If it is higher, we will store a table of log(tau) in the relevant time range, generously encompassing the range 0<z<z_max_pk, and used for the intepolation of sources
loop over modes, initial conditions and types. For each of them, allocate array of source functions.

◆ perturbations_get_k_list()

int perturbations_get_k_list	(	struct precision *	ppr,
		struct background *	pba,
		struct thermodynamics *	pth,
		struct perturbations *	ppt
	)

Define the number of comoving wavenumbers using the information passed in the precision structure.

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
pth	Input: pointer to thermodynamics structure
ppt	Input: pointer to perturbation structure

Returns: the error status

Summary:

allocate arrays related to k list for each mode
scalar modes
--> find k_max (as well as k_max_cmb[ppt->index_md_scalars], k_max_cl[ppt->index_md_scalars])
--> test that result for k_min, k_max make sense
vector modes
--> find k_max (as well as k_max_cmb[ppt->index_md_vectors], k_max_cl[ppt->index_md_vectors])
--> test that result for k_min, k_max make sense
tensor modes
--> find k_max (as well as k_max_cmb[ppt->index_md_tensors], k_max_cl[ppt->index_md_tensors])
--> test that result for k_min, k_max make sense
If user asked for k_output_values, add those to all k lists:
--> Find indices in ppt->k[index_md] corresponding to 'k_output_values'. We are assuming that ppt->k is sorted and growing, and we have made sure that ppt->k_output_values is also sorted and growing.
--> Decide if we should add k_output_value now. This has to be this complicated, since we can only compare the k-values when both indices are in range.
--> The two MIN statements are here because in a normal run, the cl and cmb arrays contain a single k value larger than their respective k_max. We are mimicking this behavior.
finally, find the global k_min and k_max for the ensemble of all modes 9scalars, vectors, tensors)

◆ perturbations_workspace_init()

int perturbations_workspace_init	(	struct precision *	ppr,
		struct background *	pba,
		struct thermodynamics *	pth,
		struct perturbations *	ppt,
		int	index_md,
		struct perturbations_workspace *	ppw
	)

Initialize a perturbations_workspace structure. All fields are allocated here, with the exception of the perturbations_vector '-->pv' field, which is allocated separately in perturbations_vector_init. We allocate one such perturbations_workspace structure per thread and per mode (scalar/../tensor). Then, for each thread, all initial conditions and wavenumbers will use the same workspace.

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
pth	Input: pointer to the thermodynamics structure
ppt	Input: pointer to the perturbation structure
index_md	Input: index of mode under consideration (scalar/.../tensor)
ppw	Input/Output: pointer to perturbations_workspace structure which fields are allocated or filled here

Returns: the error status

Summary:

define local variables
Compute maximum l_max for any multipole
Allocate [ ] array for freestreaming of multipoles (see arXiv:1305.3261) and initialize to 1.0, which is the K=0 value.
define indices of metric perturbations obeying constraint equations (this can be done once and for all, because the vector of metric perturbations is the same whatever the approximation scheme, unlike the vector of quantities to be integrated, which is allocated separately in perturbations_vector_init)
allocate some workspace in which we will store temporarily the values of background, thermodynamics, metric and source quantities at a given time
count number of approximations, initialize their indices, and allocate their flags
For definiteness, initialize approximation flags to arbitrary values (correct values are overwritten in pertub_find_approximation_switches)
allocate fields where some of the perturbations are stored

◆ perturbations_workspace_free()

int perturbations_workspace_free	(	struct perturbations *	ppt,
		int	index_md,
		struct perturbations_workspace *	ppw
	)

Free the perturbations_workspace structure (with the exception of the perturbations_vector '-->pv' field, which is freed separately in perturbations_vector_free).

Parameters

ppt	Input: pointer to the perturbation structure
index_md	Input: index of mode under consideration (scalar/.../tensor)
ppw	Input: pointer to perturbations_workspace structure to be freed

Returns: the error status

◆ perturbations_solve()

int perturbations_solve	(	struct precision *	ppr,
		struct background *	pba,
		struct thermodynamics *	pth,
		struct perturbations *	ppt,
		int	index_md,
		int	index_ic,
		int	index_k,
		struct perturbations_workspace *	ppw
	)

Solve the perturbation evolution for a given mode, initial condition and wavenumber, and compute the corresponding source functions.

For a given mode, initial condition and wavenumber, this function finds the time ranges over which the perturbations can be described within a given approximation. For each such range, it initializes (or redistributes) perturbations using perturbations_vector_init(), and integrates over time. Whenever a "source sampling time" is passed, the source terms are computed and stored in the source table using perturbations_sources().

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
pth	Input: pointer to the thermodynamics structure
ppt	Input/Output: pointer to the perturbation structure (output source functions S(k,tau) written here)
index_md	Input: index of mode under consideration (scalar/.../tensor)
index_ic	Input: index of initial condition under consideration (ad, iso...)
index_k	Input: index of wavenumber
ppw	Input: pointer to perturbations_workspace structure containing index values and workspaces

Returns: the error status

Summary:

define local variables
initialize indices relevant for back/thermo tables search
get wavenumber value
If non-zero curvature, update array of free-streaming coefficients ppw->s_l
maximum value of tau for which sources are calculated for this wavenumber
using bisection, compute minimum value of tau for which this wavenumber is integrated
find the number of intervals over which approximation scheme is constant
fill the structure containing all fixed parameters, indices and workspaces needed by perturbations_derivs
check whether we need to print perturbations to a file for this wavenumber
loop over intervals over which approximation scheme is uniform. For each interval:
--> (a) fix the approximation scheme
--> (b) get the previous approximation scheme. If the current interval starts from the initial time tau_ini, the previous approximation is set to be a NULL pointer, so that the function perturbations_vector_init() knows that perturbations must be initialized
--> (c) define the vector of perturbations to be integrated over. If the current interval starts from the initial time tau_ini, fill the vector with initial conditions for each mode. If it starts from an approximation switching point, redistribute correctly the perturbations from the previous to the new vector of perturbations.
--> (d) integrate the perturbations over the current interval.
if perturbations were printed in a file, close the file
fill the source terms array with zeros for all times between the last integrated time tau_max and tau_today.
free quantities allocated at the beginning of the routine

◆ perturbations_prepare_k_output()

int perturbations_prepare_k_output	(	struct background *	pba,
		struct perturbations *	ppt
	)

Fill array of strings with the name of the 'k_output_values' functions (transfer functions as a function of time, for fixed values of k).

Parameters

pba	Input: pointer to the background structure
ppt	Input/Output: pointer to the perturbation structure

Returns: the error status

Write titles for all perturbations that we would like to print/store.

Fluid

◆ perturbations_find_approximation_number()

int perturbations_find_approximation_number	(	struct precision *	ppr,
		struct background *	pba,
		struct thermodynamics *	pth,
		struct perturbations *	ppt,
		int	index_md,
		double	k,
		struct perturbations_workspace *	ppw,
		double	tau_ini,
		double	tau_end,
		int *	interval_number,
		int *	interval_number_of
	)

For a given mode and wavenumber, find the number of intervals of time between tau_ini and tau_end such that the approximation scheme (and the number of perturbation equations) is uniform.

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
pth	Input: pointer to the thermodynamics structure
ppt	Input: pointer to the perturbation structure
index_md	Input: index of mode under consideration (scalar/.../tensor)
k	Input: index of wavenumber
ppw	Input: pointer to perturbations_workspace structure containing index values and workspaces
tau_ini	Input: initial time of the perturbation integration
tau_end	Input: final time of the perturbation integration
interval_number	Output: total number of intervals
interval_number_of	Output: number of intervals with respect to each particular approximation

Returns: the error status

Summary:

fix default number of intervals to one (if no approximation switch)
loop over each approximation and add the number of approximation switching times

◆ perturbations_find_approximation_switches()

int perturbations_find_approximation_switches	(	struct precision *	ppr,
		struct background *	pba,
		struct thermodynamics *	pth,
		struct perturbations *	ppt,
		int	index_md,
		double	k,
		struct perturbations_workspace *	ppw,
		double	tau_ini,
		double	tau_end,
		double	precision,
		int	interval_number,
		int *	interval_number_of,
		double *	interval_limit,
		int **	interval_approx
	)

For a given mode and wavenumber, find the values of time at which the approximation changes.

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
pth	Input: pointer to the thermodynamics structure
ppt	Input: pointer to the perturbation structure
index_md	Input: index of mode under consideration (scalar/.../tensor)
k	Input: index of wavenumber
ppw	Input: pointer to perturbations_workspace structure containing index values and workspaces
tau_ini	Input: initial time of the perturbation integration
tau_end	Input: final time of the perturbation integration
precision	Input: tolerance on output values
interval_number	Input: total number of intervals
interval_number_of	Input: number of intervals with respect to each particular approximation
interval_limit	Output: value of time at the boundary of the intervals: tau_ini, tau_switch1, ..., tau_end
interval_approx	Output: value of approximations in each interval

Returns: the error status

Summary:

write in output arrays the initial time and approximation
if there are no approximation switches, just write final time and return
if there are switches, consider approximations one after each other. Find switching time by bisection. Store all switches in arbitrary order in array unsorted_tau_switch[ ]
now sort interval limits in correct order
store each approximation in chronological order

◆ perturbations_vector_init()

int perturbations_vector_init	(	struct precision *	ppr,
		struct background *	pba,
		struct thermodynamics *	pth,
		struct perturbations *	ppt,
		int	index_md,
		int	index_ic,
		double	k,
		double	tau,
		struct perturbations_workspace *	ppw,
		int *	pa_old
	)

Initialize the field '-->pv' of a perturbations_workspace structure, which is a perturbations_vector structure. This structure contains indices and values of all quantities which need to be integrated with respect to time (and only them: quantities fixed analytically or obeying constraint equations are NOT included in this vector). This routine distinguishes between two cases:

--> the input pa_old is set to the NULL pointer:

This happens when we start integrating over a new wavenumber and we want to set initial conditions for the perturbations. Then, it is assumed that ppw-->pv is not yet allocated. This routine allocates it, defines all indices, and then fills the vector ppw-->pv-->y with the initial conditions defined in perturbations_initial_conditions.

--> the input pa_old is not set to the NULL pointer and describes some set of approximations:

This happens when we need to change approximation scheme while integrating over a given wavenumber. The new approximation described by ppw-->pa is then different from pa_old. Then, this routine allocates a new vector with a new size and new index values; it fills this vector with initial conditions taken from the previous vector passed as an input in ppw-->pv, and eventually with some analytic approximations for the new variables appearing at this time; then the new vector comes in replacement of the old one, which is freed.

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
pth	Input: pointer to the thermodynamics structure
ppt	Input: pointer to the perturbation structure
index_md	Input: index of mode under consideration (scalar/.../tensor)
index_ic	Input: index of initial condition under consideration (ad, iso...)
k	Input: wavenumber
tau	Input: conformal time
ppw	Input/Output: workspace containing in input the approximation scheme, the background/thermodynamics/metric quantities, and eventually the previous vector y; and in output the new vector y.
pa_old	Input: NULL is we need to set y to initial conditions for a new wavenumber; points towards a perturbations_approximations if we want to switch of approximation.

Returns: the error status

Summary:

define local variables
allocate a new perturbations_vector structure to which ppw-->pv will point at the end of the routine
initialize pointers to NULL (they will be allocated later if needed), relevant for perturbations_vector_free()
define all indices in this new vector (depends on approximation scheme, described by the input structure ppw-->pa)
(a) metric perturbations V or depending on gauge
(b) metric perturbation h is a propagating degree of freedom, so h and hdot are included in the vector of ordinary perturbations, no in that of metric perturbations
allocate vectors for storing the values of all these quantities and their time-derivatives at a given time
specify which perturbations are needed in the evaluation of source terms
case of setting initial conditions for a new wavenumber
--> (a) check that current approximation scheme is consistent with initial conditions
--> (b) let ppw-->pv points towards the perturbations_vector structure that we just created
--> (c) fill the vector ppw-->pv-->y with appropriate initial conditions
case of switching approximation while a wavenumber is being integrated
--> (a) for the scalar mode:
—> (a.1.) check that the change of approximation scheme makes sense (note: before calling this routine there is already a check that we wish to change only one approximation flag at a time)
—> (a.2.) some variables (b, cdm, fld, ...) are not affected by any approximation. They need to be reconducted whatever the approximation switching is. We treat them here. Below we will treat other variables case by case.
--> (b) for the vector mode
—> (b.1.) check that the change of approximation scheme makes sense (note: before calling this routine there is already a check that we wish to change only one approximation flag at a time)
—> (b.2.) some variables (gw, gwdot, ...) are not affected by any approximation. They need to be reconducted whatever the approximation switching is. We treat them here. Below we will treat other variables case by case.
--> (c) for the tensor mode
—> (c.1.) check that the change of approximation scheme makes sense (note: before calling this routine there is already a check that we wish to change only one approximation flag at a time)
—> (c.2.) some variables (gw, gwdot, ...) are not affected by any approximation. They need to be reconducted whatever the approximation switching is. We treat them here. Below we will treat other variables case by case.
--> (d) free the previous vector of perturbations
--> (e) let ppw-->pv points towards the perturbations_vector structure that we just created

◆ perturbations_vector_free()

int perturbations_vector_free ( struct perturbations_vector * pv )

Free the perturbations_vector structure.

Parameters

pv	Input: pointer to perturbations_vector structure to be freed

Returns: the error status

◆ perturbations_initial_conditions()

int perturbations_initial_conditions	(	struct precision *	ppr,
		struct background *	pba,
		struct perturbations *	ppt,
		int	index_md,
		int	index_ic,
		double	k,
		double	tau,
		struct perturbations_workspace *	ppw
	)

For each mode, wavenumber and initial condition, this function initializes in the vector all values of perturbed variables (in a given gauge). It is assumed here that all values have previously been set to zero, only non-zero values are set here.

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
ppt	Input: pointer to the perturbation structure
index_md	Input: index of mode under consideration (scalar/.../tensor)
index_ic	Input: index of initial condition under consideration (ad, iso...)
k	Input: wavenumber
tau	Input: conformal time
ppw	Input/Output: workspace containing in input the approximation scheme, the background/thermodynamics/metric quantities, and eventually the previous vector y; and in output the new vector y.

Returns: the error status

Summary:

--> Declare local variables

--> For scalars

(a) compute relevant background quantities: compute rho_r, rho_m, rho_nu (= all relativistic except photons), and their ratio.
(b) starts by setting everything in synchronous gauge. If another gauge is needed, we will perform a gauge transformation below.
--> (b.1.) adiabatic
—> Canonical field (solving for the perturbations): initial perturbations set to zero, they should reach the attractor soon enough.
—> TODO: Incorporate the attractor IC from 1004.5509. delta_phi $= -(a/k)^2/\phi'(\rho + p)\theta$ , delta_phi_prime $= a^2/\phi'$ (delta_rho_phi + V'delta_phi), and assume theta, delta_rho as for perfect fluid with and w = 1/3 (ASSUMES radiation TRACKING)
--> (b.2.) Cold dark matter Isocurvature
--> (b.3.) Baryon Isocurvature
--> (b.4.) Neutrino density Isocurvature
--> (b.5.) Neutrino velocity Isocurvature
(c) If the needed gauge is really the synchronous gauge, we need to affect the previously computed value of eta to the actual variable eta
(d) If the needed gauge is the newtonian gauge, we must compute alpha and then perform a gauge transformation for each variable
(e) In any gauge, we should now implement the relativistic initial conditions in ur and ncdm variables

--> For tensors

tensor initial conditions take into account the fact that scalar (resp. tensor) 's are related to the real space power spectrum of curvature (resp. of the tensor part of metric perturbations)

$<R(x) R(x)> \ \ \sum_{ij} <h_{ij}(x) h^{ij}(x)>$

In momentum space it is conventional to use the modes R(k) and h(k) where the quantity h obeying to the equation of propagation:

$h'' + \frac{2a'}{a} h + [k2+2K] h = 12\pi Ga2 (\rho+p) \sigma = 8\pi Ga2 p \pi$

and the power spectra in real space and momentum space are related through:

$<R(x) R(x)> = \int \frac{dk}{k} \left[ \frac{k^3}{2\pi^2} <R(k)R(k)^*>\right] = \int \frac{dk}{k} \mathcal{P}_R(k)$

$\sum_{ij} <h_{ij}(x) h^{ij}(x)> = \frac{dk}{k} \left[ \frac{k^3}{2\pi^2} F\left(\frac{k^2}{K}\right) <h(k)h(k)^*>\right] = \int \frac{dk}{k} F\left(\frac{k^2}{K}\right) \mathcal{P}_h(k)$

where $\mathcal{P}_R$ and $\mathcal{P}_h$ are the dimensionless spectrum of curvature R, and F is a function of k2/K, where K is the curvature parameter. F is equal to one in flat space (K=0), and coming from the contraction of the laplacian eigentensor $Q_{ij}$ with itself. We will give F explicitly below.

Similarly the scalar (S) and tensor (T) 's are given by

$C_l^S = 4\pi \int \frac{dk}{k} [\Delta_l^S(q)]^2 \mathcal{P}_R(k)$

$C_l^T = 4\pi \int \frac{dk}{k} [\Delta_l^T(q)]^2 F\left(\frac{k^2}{K}\right) \mathcal{P}_h(k)$

The usual convention for the tensor-to-scalar ratio at pivot scale = 16 epsilon in single-field inflation is such that for constant $\mathcal{P}_R(k)$ and $\mathcal{P}_h(k)$ ,

$r = 6 \frac{\mathcal{P}_h(k)}{\mathcal{P}_R(k)}$

so

$\mathcal{P}_h(k) = \frac{\mathcal{P}_R(k) r}{6} = \frac{A_s r}{6} = \frac{A_t}{6}$

A priori it would make sense to say that for a power-law primordial spectrum there is an extra factor $(k/k_{pivot})^{n_t}$ (and eventually running and so on and so forth...)

However it has been shown that the minimal models of inflation in a negatively curved bubble lead to $\mathcal{P}_h(k)=\tanh(\pi*\nu/2)$ . In open models it is customary to define the tensor tilt in a non-flat universe as a deviation from this behavior rather than from true scale-invariance in the above sense.

Hence we should have

$\mathcal{P}_h(k) = \frac{A_t}{6} [ \tanh(\pi*\frac{\nu}{2})] (k/k_{pivot})^{(n_t+...)}$

where the brackets

$[...]$

mean "if K<0"

Then

$C_l^T = 4\pi \int \frac{dk}{k} [\Delta_l^T(q)]^2 F\left(\frac{k^2}{K}\right) \frac{A_t}{6} [\tanh(\pi*\frac{\nu}{2})] (k/k_{pivot})^{(n_t+...)}$

In the code, it is then a matter of choice to write:

In the primordial module: $\mathcal{P}_h(k) = \frac{A_t}{6} \tanh{(\pi*\frac{\nu}{2})} (k/k^*)^{n_T}$
In the perturbation initial conditions:
In the harmonic module: $C_l^T = \frac{4}{\pi} \int \frac{dk}{k} [\Delta_l^T(q)]^2 F\left(\frac{k^2}{K}\right) \mathcal{P}_h(k)$

or:

In the primordial module: $\mathcal{P}_h(k) = A_t (k/k^*)^{n_T}$
In the perturbation initial conditions: $h = \sqrt{[F\left(\frac{k^2}{K}\right) / 6] \tanh{(\pi*\frac{\nu}{2})}}$
In the harmonic module: $C_l^T = \frac{4}{\pi} \int \frac{dk}{k} [\Delta_l^T(q)]^2 \mathcal{P}_h(k)$

We choose this last option, such that the primordial and harmonic module differ minimally in flat and non-flat space. Then we must impose

$h = \sqrt{\left(\frac{F}{6}\right) \tanh{(\pi*\frac{\nu}{2})}}$

The factor F is found to be given by:

$\sum_{ij}<h_{ij}(x) h^{ij}(x)> = \int \frac{dk}{k} \frac{k2(k2-K)}{(k2+3K)(k2+2K)} \mathcal{P}_h(k)$

Introducing as usual and using qdq = kdk this gives

$\sum_{ij}<h_{ij}(x) h^{ij}(x)> = \int \frac{dk}{k} \frac{(q2-3K)(q2-4K)}{q2(q2-K)} \mathcal{P}_h(k)$

Using qdq = kdk this is equivalent to

$\sum_{ij}<h_{ij}(x) h^{ij}(x)> = \int \frac{dq}{q} \frac{q2-4K}{q2-K} \mathcal{P}_h(k(q))$

Finally, introducing $\nu=q/\sqrt{|K|}$ and sgnK=SIGN(k) $=\pm 1$ , this could also be written

$\sum_{ij}<h_{ij}(x) h^{ij}(x)> = \int \frac{d\nu}{\nu} \frac{(\nu2-4sgnK)}{(\nu2-sgnK)} \mathcal{P}_h(k(\nu))$

Equation (43,44) of Hu, Seljak, White, Zaldarriaga is equivalent to absorbing the above factor $(\nu2-4sgnK)/(\nu2-sgnK)$ in the definition of the primordial spectrum. Since the initial condition should be written in terms of k rather than nu, they should read

$h = \sqrt{ [k2(k2-K)]/[(k2+3K)(k2+2K)] / 6 * \tanh{(\pi*\frac{\nu}{2})} }$

We leave the freedom to multiply by an arbitrary number ppr->gw_ini. The standard convention corresponding to standard definitions of r, , is however ppr->gw_ini=1.

◆ perturbations_approximations()

int perturbations_approximations	(	struct precision *	ppr,
		struct background *	pba,
		struct thermodynamics *	pth,
		struct perturbations *	ppt,
		int	index_md,
		double	k,
		double	tau,
		struct perturbations_workspace *	ppw
	)

Evaluate background/thermodynamics at $\tau$ , infer useful flags / time scales for integrating perturbations.

Evaluate background quantities at $\tau$ , as well as thermodynamics for scalar mode; infer useful flags and time scales for integrating the perturbations:

check whether tight-coupling approximation is needed.
check whether radiation (photons, massless neutrinos...) perturbations are needed.
choose step of integration: step = ppr->perturbations_integration_stepsize * min_time_scale, where min_time_scale = smallest time scale involved in the equations. There are three time scales to compare:
1. that of recombination, $\tau_c = 1/\kappa'$
2. Hubble time scale, $\tau_h = a/a'$
3. Fourier mode, $\tau_k = 1/k$

So, in general, min_time_scale = $\min(\tau_c, \tau_b, \tau_h, \tau_k)$ .

However, if $\tau_c \ll \tau_h$ and $\tau_c \ll \tau_k$ , we can use the tight-coupling regime for photons and write equations in such way that the time scale $\tau_c$ becomes irrelevant (no effective mass term in $1/\tau_c$ ). Then, the smallest scale in the equations is only $\min(\tau_h, \tau_k)$ . In practise, it is sufficient to use only the condition $\tau_c \ll \tau_h$ .

Also, if $\rho_{matter} \gg \rho_{radiation}$ and $k \gg aH$ , we can switch off radiation perturbations (i.e. switch on the free-streaming approximation) and then the smallest scale is simply $\tau_h$ .

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
pth	Input: pointer to thermodynamics structure
ppt	Input: pointer to the perturbation structure
index_md	Input: index of mode under consideration (scalar/.../tensor)
k	Input: wavenumber
tau	Input: conformal time
ppw	Input/Output: in output contains the approximation to be used at this time

Returns: the error status

Summary:

define local variables
compute Fourier mode time scale = $\tau_k = 1/k$
evaluate background quantities with background_at_tau() and Hubble time scale $\tau_h = a/a'$
for scalar modes:
--> (a) evaluate thermodynamical quantities with thermodynamics_at_z()
—> (b.1.) if $\kappa'=0$ , recombination is finished; tight-coupling approximation must be off
—> (b.2.) if $\kappa' \neq 0$ , recombination is not finished: check tight-coupling approximation
-—> (b.2.a) compute recombination time scale for photons, $\tau_{\gamma} = 1/ \kappa'$
-—> (b.2.b) check whether tight-coupling approximation should be on
--> (c) free-streaming approximations
for tensor modes:
--> (a) evaluate thermodynamical quantities with thermodynamics_at_z()
—> (b.1.) if $\kappa'=0$ , recombination is finished; tight-coupling approximation must be off
—> (b.2.) if $\kappa' \neq 0$ , recombination is not finished: check tight-coupling approximation
-—> (b.2.a) compute recombination time scale for photons, $\tau_{\gamma} = 1/ \kappa'$
-—> (b.2.b) check whether tight-coupling approximation should be on

◆ perturbations_timescale()

int perturbations_timescale	(	double	tau,
		void *	parameters_and_workspace,
		double *	timescale,
		ErrorMsg	error_message
	)

Compute typical timescale over which the perturbation equations vary. Some integrators (e.g. Runge-Kunta) benefit from calling this routine at each step in order to adapt the next step.

This is one of the few functions in the code which is passed to the generic_integrator() routine. Since generic_integrator() should work with functions passed from various modules, the format of the arguments is a bit special:

fixed parameters and workspaces are passed through a generic pointer. generic_integrator() doesn't know the content of this pointer.
the error management is a bit special: errors are not written as usual to pth->error_message, but to a generic error_message passed in the list of arguments.

Parameters

tau	Input: conformal time
parameters_and_workspace	Input: fixed parameters (e.g. indices), workspace, approximation used, etc.
timescale	Output: perturbation variation timescale (given the approximation used)
error_message	Output: error message

Summary:

define local variables
extract the fields of the parameter_and_workspace input structure
compute Fourier mode time scale = $\tau_k = 1/k$
evaluate background quantities with background_at_tau() and Hubble time scale $\tau_h = a/a'$
for scalars modes:
--> compute recombination time scale for photons, $\tau_{\gamma} = 1/ \kappa'$
for vector modes:
--> compute recombination time scale for photons, $\tau_{\gamma} = 1/ \kappa'$
for tensor modes:
--> compute recombination time scale for photons, $\tau_{\gamma} = 1/ \kappa'$

◆ perturbations_einstein()

int perturbations_einstein	(	struct precision *	ppr,
		struct background *	pba,
		struct thermodynamics *	pth,
		struct perturbations *	ppt,
		int	index_md,
		double	k,
		double	tau,
		double *	y,
		struct perturbations_workspace *	ppw
	)

Compute metric perturbations (those not integrated over time) using Einstein equations

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
pth	Input: pointer to thermodynamics structure
ppt	Input: pointer to the perturbation structure
index_md	Input: index of mode under consideration (scalar/.../tensor)
k	Input: wavenumber
tau	Input: conformal time
y	Input: vector of perturbations (those integrated over time) (already allocated)
ppw	Input/Output: in output contains the updated metric perturbations

Returns: the error status

Summary:

define local variables
define wavenumber and scale factor related quantities
sum up perturbations from all species
for scalar modes:
--> infer metric perturbations from Einstein equations
for vector modes
for tensor modes

◆ perturbations_total_stress_energy()

int perturbations_total_stress_energy	(	struct precision *	ppr,
		struct background *	pba,
		struct thermodynamics *	pth,
		struct perturbations *	ppt,
		int	index_md,
		double	k,
		double *	y,
		struct perturbations_workspace *	ppw
	)

Summary:

define local variables

Variables used for FLD and PPF

wavenumber and scale factor related quantities
for scalar modes
--> (a) deal with approximation schemes
—> (a.1.) photons
-—> (a.1.1.) no approximation
-—> (a.1.2.) radiation streaming approximation
-—> (a.1.3.) tight coupling approximation
—> (a.2.) ur
—> (a.3.) baryon pressure perturbation
—> (a.4.) interacting dark radiation
--> (b) compute the total density, velocity and shear perturbations

We must gauge transform the pressure perturbation from the fluid rest-frame to the gauge we are working in

The equation is too stiff for Runge-Kutta when c_gamma_k_H_square is large. Use the asymptotic solution Gamma=Gamma'=0 in that case.

We must now check the stiffness criterion again and set Gamma_prime_fld accordingly.

Now construct the pressure perturbation, see 1903.xxxxx.

Construct energy density and pressure for DE (_fld) and the rest (_t). Also compute derivatives.

Compute background quantities X,Y,Z and their derivatives.

Construct theta_t and its derivative from the Euler equation

Analytic derivative of the equation for ppw->rho_plus_p_theta_fld above.

We can finally compute the pressure perturbation using the Euler equation for theta_fld

for vector modes
--> photon contribution to vector sources:
--> baryons
for tensor modes
--> photon contribution to gravitational wave source:
--> ur contribution to gravitational wave source:
--> ncdm contribution to gravitational wave source:

◆ perturbations_sources()

int perturbations_sources	(	double	tau,
		double *	y,
		double *	dy,
		int	index_tau,
		void *	parameters_and_workspace,
		ErrorMsg	error_message
	)

Compute the source functions (three terms for temperature, one for E or B modes, etc.)

This is one of the few functions in the code which is passed to the generic_integrator() routine. Since generic_integrator() should work with functions passed from various modules, the format of the arguments is a bit special:

fixed parameters and workspaces are passed through a generic pointer. generic_integrator() doesn't know the content of this pointer.
the error management is a bit special: errors are not written as usual to pth->error_message, but to a generic error_message passed in the list of arguments.

Parameters

tau	Input: conformal time
y	Input: vector of perturbations
dy	Input: vector of time derivative of perturbations
index_tau	Input: index in the array tau_sampling
parameters_and_workspace	Input/Output: in input, all parameters needed by perturbations_derivs, in output, source terms
error_message	Output: error message

Returns: the error status

Summary:

define local variables
rename structure fields (just to avoid heavy notations)
get background/thermo quantities in this point
for scalars
--> compute metric perturbations
--> compute quantities depending on approximation schemes
--> for each type, compute source terms

gamma in N-body gauge. Eq. A.2 in 1811.00904 gives k2gamma = (a'/a)H_T' + k2(phi-psi) - H_T''. The last term is cubersome to calculate (one would need finite derivatives) but usually small. Here we only compute an approximate k2gamma without this last term. If needed, the term could be restored: you can see how T. Tram did it in a previous commit beec79548877e1e43403d1f4de5ddee6741a3c16 (28.02.2019) - then it had to go to harmonic.c, now it could stay in this module. Later this feature was removed for simplicity. Note that to compute the transfer functions in the N-body gauge we do not need k2gamma anyway.

We follow the (debatable) CMBFAST/CAMB convention of not including rho_lambda in rho_tot

for tensors
--> compute quantities depending on approximation schemes

◆ perturbations_print_variables()

int perturbations_print_variables	(	double	tau,
		double *	y,
		double *	dy,
		void *	parameters_and_workspace,
		ErrorMsg	error_message
	)

When testing the code or a cosmological model, it can be useful to output perturbations at each step of integration (and not just the delta's at each source sampling point, which is achieved simply by asking for matter transfer functions). Then this function can be passed to the generic_evolver routine.

By default, instead of passing this function to generic_evolver, one passes a null pointer. Then this function is just not used.

Parameters

tau	Input: conformal time
y	Input: vector of perturbations
dy	Input: vector of its derivatives (already allocated)
parameters_and_workspace	Input: fixed parameters (e.g. indices)
error_message	Output: error message

Summary:

define local variables
ncdm sector begins
ncdm sector ends
rename structure fields (just to avoid heavy notations)
update background/thermo quantities in this point
update metric perturbations in this point
calculate perturbed recombination
for scalar modes
--> Get delta, deltaP/rho, theta, shear and store in array
--> TODO: gauge transformation of delta, deltaP/rho (?) and theta using -= 3aH(1+w_ncdm) alpha for delta.
--> Handle (re-)allocation

Fluid

for tensor modes:
--> Handle (re-)allocation

◆ perturbations_derivs()

int perturbations_derivs	(	double	tau,
		double *	y,
		double *	dy,
		void *	parameters_and_workspace,
		ErrorMsg	error_message
	)

Compute derivative of all perturbations to be integrated

For each mode (scalar/vector/tensor) and each wavenumber k, this function computes the derivative of all values in the vector of perturbed variables to be integrated.

This is one of the few functions in the code which is passed to the generic_integrator() routine. Since generic_integrator() should work with functions passed from various modules, the format of the arguments is a bit special:

fixed parameters and workspaces are passed through a generic pointer. generic_integrator() doesn't know what the content of this pointer is.
errors are not written as usual in pth->error_message, but in a generic error_message passed in the list of arguments.

Parameters

tau	Input: conformal time
y	Input: vector of perturbations
dy	Output: vector of its derivatives (already allocated)
parameters_and_workspace	Input/Output: in input, fixed parameters (e.g. indices); in output, background and thermo quantities evaluated at tau.
error_message	Output: error message

Summary:

define local variables
rename the fields of the input structure (just to avoid heavy notations)
get background/thermo quantities in this point
get metric perturbations with perturbations_einstein()
compute related background quantities
Compute 'generalised cotK function of argument $\sqrt{|K|}*\tau$ , for closing hierarchy. (see equation 2.34 in arXiv:1305.3261):
for scalar modes:
--> (a) define short-cut notations for the scalar perturbations
--> (b) perturbed recombination
--> (c) compute metric-related quantities (depending on gauge; additional gauges can be coded below)
Each continuity equation contains a term in (theta+metric_continuity) with metric_continuity = (h_prime/2) in synchronous gauge, (-3 phi_prime) in newtonian gauge
Each Euler equation contains a source term metric_euler with metric_euler = 0 in synchronous gauge, (k2 psi) in newtonian gauge
Each shear derivative equation contains a source term metric_shear equal to metric_shear = (h_prime+6eta_prime)/2 in synchronous gauge, 0 in newtonian gauge
metric_shear_prime is the derivative of metric_shear
In the ufa_class approximation, the leading-order source term is (h_prime/2) in synchronous gauge, (-3 (phi_prime+psi_prime)) in newtonian gauge: we approximate the later by (-6 phi_prime)
--> (d) if some approximation schemes are turned on, enforce a few y[] values computed in perturbations_einstein
--> (e) BEGINNING OF ACTUAL SYSTEM OF EQUATIONS OF EVOLUTION
—> photon temperature density
—> baryon density
—> baryon velocity (depends on tight-coupling approximation=tca)
-—> perturbed recombination has an impact
—> photon temperature higher momenta and photon polarization (depend on tight-coupling approximation)
-—> if photon tight-coupling is off
--—> define $\Pi = G_{\gamma 0} + G_{\gamma 2} + F_{\gamma 2}$
--—> photon temperature velocity
--—> photon temperature shear
--—> photon temperature l=3
--—> photon temperature l>3
--—> photon temperature lmax
--—> photon polarization l=0
--—> photon polarization l=1
--—> photon polarization l=2
--—> photon polarization l>2
--—> photon polarization lmax_pol
-—> if photon tight-coupling is on:
--—> in that case, only need photon velocity
—> cdm
-—> newtonian gauge: cdm density and velocity
-—> synchronous gauge: cdm density only (velocity set to zero by definition of the gauge)
—> interacting dark radiation
-—> idr velocity
-—> exact idr shear
-—> exact idr l=3
-—> exact idr l>3
-—> exact idr lmax_dr
—> dcdm and dr
-—> dcdm
—> dr
-—> dr F0
-—> dr F1
-—> exact dr F2
-—> exact dr l=3
-—> exact dr l>3
-—> exact dr lmax_dr
—> fluid (fld)
-—> factors w, w_prime, adiabatic sound speed ca2 (all three background-related), plus actual sound speed in the fluid rest frame cs2
-—> fluid density
-—> fluid velocity
—> scalar field (scf)
-—> field value
-—> Klein Gordon equation
—> ultra-relativistic neutrino/relics (ur)
-—> if radiation streaming approximation is off
--—> ur density
--—> ur velocity
--—> exact ur shear
--—> exact ur l=3
--—> exact ur l>3
--—> exact ur lmax_ur
--—> in fluid approximation (ufa): only ur shear needed
—> non-cold dark matter (ncdm): massive neutrinos, WDM, etc.
-—> first case: use a fluid approximation (ncdmfa)
--—> loop over species
--—> define intermediate quantitites
--—> exact continuity equation
--—> exact euler equation
--—> different ansatz for approximate shear derivative
--—> jump to next species
-—> second case: use exact equation (Boltzmann hierarchy on momentum grid)
--—> loop over species
--—> loop over momentum
--—> define intermediate quantities
--—> ncdm density for given momentum bin
--—> ncdm velocity for given momentum bin
--—> ncdm shear for given momentum bin
--—> ncdm l>3 for given momentum bin
--—> ncdm lmax for given momentum bin (truncation as in Ma and Bertschinger) but with curvature taken into account a la arXiv:1305.3261
--—> jump to next momentum bin or species
—> metric
—> eta of synchronous gauge
vector mode
--> baryon velocity
tensor modes:
--> non-cold dark matter (ncdm): massive neutrinos, WDM, etc.
—> loop over species
-—> loop over momentum
-—> define intermediate quantities
-—> ncdm density for given momentum bin
-—> ncdm l>0 for given momentum bin
-—> ncdm lmax for given momentum bin (truncation as in Ma and Bertschinger) but with curvature taken into account a la arXiv:1305.3261
-—> jump to next momentum bin or species
--> tensor metric perturbation h (gravitational waves)
--> its time-derivative

◆ perturbations_tca_slip_and_shear()

int perturbations_tca_slip_and_shear	(	double *	y,
		void *	parameters_and_workspace,
		ErrorMsg	error_message
	)

Compute the baryon-photon slip (theta_g - theta_b)' and the photon shear in the tight-coupling approximation

Parameters

y	Input: vector of perturbations
parameters_and_workspace	Input/Output: in input, fixed parameters (e.g. indices); in output, slip and shear
error_message	Output: error message

Summary:

define local variables
rename the fields of the input structure (just to avoid heavy notations)
compute related background quantities
--> (a) define short-cut notations for the scalar perturbations
--> (b) define short-cut notations used only in tight-coupling approximation
--> (c) compute metric-related quantities (depending on gauge; additional gauges can be coded below)
Each continuity equation contains a term in (theta+metric_continuity) with metric_continuity = (h_prime/2) in synchronous gauge, (-3 phi_prime) in newtonian gauge
Each Euler equation contains a source term metric_euler with metric_euler = 0 in synchronous gauge, (k2 psi) in newtonian gauge
Each shear derivative equation contains a source term metric_shear equal to metric_shear = (h_prime+6eta_prime)/2 in synchronous gauge, 0 in newtonian gauge
metric_shear_prime is the derivative of metric_shear
In the ufa_class approximation, the leading-order source term is (h_prime/2) in synchronous gauge, (-3 (phi_prime+psi_prime)) in newtonian gauge: we approximate the later by (-6 phi_prime)
--> (d) if some approximation schemes are turned on, enforce a few y[ ] values computed in perturbations_einstein
—> intermediate quantities for 2nd order tca (or 1st with idm_g): zero order for theta_b' = theta_g'
—> standard photon velocity derivative without tca (neglecting shear) - only needed for idm_g_dr behavior
—> like Ma & Bertschinger
—> relax assumption dkappa~a $^{-2}$ (like in CAMB)
—> also relax assumption cb2~a $^{-1}$
—> intermediate quantities for 2nd order tca: shear_g at first order in tight-coupling
—> intermediate quantities for 2nd order tca: shear_g_prime at first order in tight-coupling
—> 2nd order as in CRS
—> 2nd order like in CLASS paper
—> add only the most important 2nd order terms
—> store tight-coupling values of photon shear and its derivative

◆ perturbations_rsa_delta_and_theta()

int perturbations_rsa_delta_and_theta	(	struct precision *	ppr,
		struct background *	pba,
		struct thermodynamics *	pth,
		struct perturbations *	ppt,
		double	k,
		double *	y,
		double	a_prime_over_a,
		double *	pvecthermo,
		struct perturbations_workspace *	ppw,
		ErrorMsg	error_message
	)

Compute the density delta and velocity theta of photons and ultra-relativistic neutrinos in the radiation streaming approximation

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
pth	Input: pointer to thermodynamics structure
ppt	Input: pointer to perturbation structure
k	Input: wavenumber
y	Input: vector of perturbations
a_prime_over_a	Input: a'/a
pvecthermo	Input: vector of thermodynamics quantites
ppw	Input/Output: in input, fixed parameters (e.g. indices); in output, delta and theta
error_message	Output: error message

◆ perturbations_rsa_idr_delta_and_theta()

int perturbations_rsa_idr_delta_and_theta	(	struct precision *	ppr,
		struct background *	pba,
		struct thermodynamics *	pth,
		struct perturbations *	ppt,
		double	k,
		double *	y,
		double	a_prime_over_a,
		double *	pvecthermo,
		struct perturbations_workspace *	ppw,
		ErrorMsg	error_message
	)

Compute the density delta and velocity theta of interacting dark radiation in its streaming approximation

Parameters

ppr	Input: pointer to precision structure
pba	Input: pointer to background structure
pth	Input: pointer to thermodynamics structure
ppt	Input: pointer to perturbation structure
k	Input: wavenumber
y	Input: vector of perturbations
a_prime_over_a	Input: a'/a
pvecthermo	Input: vector of thermodynamics quantites
ppw	Input/Output: in input, fixed parameters (e.g. indices); in output, delta and theta
error_message	Output: error message

Functions

Detailed Description

Function Documentation

◆ perturbations_sources_at_tau()

◆ perturbations_sources_at_z()

◆ perturbations_sources_at_k_and_z()

◆ perturbations_output_data_at_z()

◆ perturbations_output_data_at_index_tau()

◆ perturbations_output_data()

◆ perturbations_output_titles()

◆ perturbations_output_firstline_and_ic_suffix()

◆ perturbations_init()

◆ perturbations_free_input()

◆ perturbations_free()

◆ perturbations_indices()

◆ perturbations_timesampling_for_sources()

◆ perturbations_get_k_list()

◆ perturbations_workspace_init()

◆ perturbations_workspace_free()

◆ perturbations_solve()

◆ perturbations_prepare_k_output()

◆ perturbations_find_approximation_number()

◆ perturbations_find_approximation_switches()

◆ perturbations_vector_init()

◆ perturbations_vector_free()

◆ perturbations_initial_conditions()

◆ perturbations_approximations()

◆ perturbations_timescale()

◆ perturbations_einstein()

◆ perturbations_total_stress_energy()

◆ perturbations_sources()

◆ perturbations_print_variables()

◆ perturbations_derivs()

◆ perturbations_tca_slip_and_shear()

◆ perturbations_rsa_delta_and_theta()

◆ perturbations_rsa_idr_delta_and_theta()