definitions of cosmological functions like power spectrum, growth, etc. More...

#include "core_power.h"
#include "class_data_vec.hpp"
#include "params.hpp"
#include <ccl_background.h>
#include <ccl_power.h>
#include <gsl/gsl_deriv.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_spline.h>

Include dependency graph for core_power.cpp:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	anonymous_namespace{core_power.cpp}::integrand_param< P >
	general integrand with callable function and one parameter More...

class	anonymous_namespace{core_power.cpp}::Integr_obj
	basic integration object, wrapper for GSL integration functions More...

class	anonymous_namespace{core_power.cpp}::Integr_obj_qag
	QAG adaptive integration. More...

class	anonymous_namespace{core_power.cpp}::Integr_obj_qagiu
	QAGI adaptive integration on infinite intervals. More...

class	anonymous_namespace{core_power.cpp}::Integr_obj_qawo
	QAWO adaptive integration for oscillatory functions. More...

class	anonymous_namespace{core_power.cpp}::Integr_obj_qawf
	QAWF adaptive integration for Fourier integrals. More...

class	anonymous_namespace{core_power.cpp}::ODE_Solver
	: Explicit embedded Runge-Kutta Prince-Dormand (8, 9) method. More...

Namespaces
	anonymous_namespace{core_power.cpp}

Functions
template<typename T >
T	anonymous_namespace{core_power.cpp}::hubble_param (T a, const Cosmo_Param &cosmo)

double	anonymous_namespace{core_power.cpp}::growth_factor_integrand (double a, void *params)

double	anonymous_namespace{core_power.cpp}::growth_factor (double a, void *params)

double	anonymous_namespace{core_power.cpp}::ln_growth_factor (double log_a, void *parameters)

template<typename T >
size_t	anonymous_namespace{core_power.cpp}::get_nearest (const T val, const std::vector< T > &vec)

template<class P >
double	anonymous_namespace{core_power.cpp}::xi_integrand_W (double k, void *params)
	integrand for correlation function when weight-function 'sin(kr)' is used in integration More...

template<class P >
double	anonymous_namespace{core_power.cpp}::xi_integrand_G (double k, void *params)
	integrand for correlation function when non-weighted integration is used More...

template<class P >
double	anonymous_namespace{core_power.cpp}::sigma_integrand_G (double k, void *params)
	integrand for amplitude of mass fluctuaton when non-weighted integration is used More...

template<class P , class T >
void	anonymous_namespace{core_power.cpp}::gen_fce_r_binned_gsl (const double x_min, const double x_max, const double lin_bin, const P &P_k, Data_Vec< double, 2 > &fce_binned, T &fce_r)
	compute given function in linear range and store results More...

template<class P , class T >
void	anonymous_namespace{core_power.cpp}::gen_corr_func_binned_gsl (const Sim_Param &sim, const P &P_k, Data_Vec< double, 2 > &corr_func_binned, T &xi_r)
	compute correlation function and store results More...

template<class P , class T >
void	anonymous_namespace{core_power.cpp}::gen_sigma_func_binned_gsl (const Sim_Param &sim, const P &P_k, Data_Vec< double, 2 > &sigma_binned, T &sigma_r)
	compute amplitude of density fluctuation and store results More...

void	norm_pwr (Cosmo_Param &cosmo)
	< end of anonymous namespace (private definitions) More...

double	norm_growth_factor (const Cosmo_Param &cosmo)
	when computing growth factor outside CCL range we need to normalize the growth factor; $D(a=1)\equiv1$ More...

double	growth_factor (double a, const Cosmo_Param &cosmo)

double	growth_rate (double a, const Cosmo_Param &cosmo)

double	growth_change (double a, const Cosmo_Param &cosmo)

double	Omega_lambda (double a, const Cosmo_Param &cosmo)

double	lin_pow_spec (double a, double k, const Cosmo_Param &cosmo)

double	non_lin_pow_spec (double a, double k, const Cosmo_Param &cosmo)

static double	truncation_fce (double k, double k2_G)

template<class P >
static double	get_gsl_error (const Cosmo_Param &cosmo, const P &P_k)

template<class P >
void	gen_corr_func_binned_gsl_qawf (const Sim_Param &sim, const P &P_k, Data_Vec< double, 2 > &corr_func_binned)
	compute correlation function and store results More...

void	gen_corr_func_binned_gsl_qawf_lin (const Sim_Param &sim, double a, Data_Vec< double, 2 > &corr_func_binned)
	compute linear correlation function and store results More...

void	gen_corr_func_binned_gsl_qawf_nl (const Sim_Param &sim, double a, Data_Vec< double, 2 > &corr_func_binned)
	compute non-linear correlation function and store results More...

template<class P >
void	gen_sigma_binned_gsl_qawf (const Sim_Param &sim, const P &P_k, Data_Vec< double, 2 > &sigma_binned)
	compute amplitude of density fluctuation and store results More...

void	gen_sigma_func_binned_gsl_qawf_lin (const Sim_Param &sim, double a, Data_Vec< double, 2 > &sigma_binned)
	compute linear amplitude of density fluctuation and store results More...

void	gen_sigma_func_binned_gsl_qawf_nl (const Sim_Param &sim, double a, Data_Vec< double, 2 > &sigma_binned)
	compute non-linear amplitude of density fluctuation and store results More...

template void	gen_corr_func_binned_gsl_qawf (const Sim_Param &, const Extrap_Pk< double, 2 > &, Data_Vec< double, 2 > &)

template void	gen_corr_func_binned_gsl_qawf (const Sim_Param &, const Extrap_Pk< double, 3 > &, Data_Vec< double, 2 > &)

template void	gen_corr_func_binned_gsl_qawf (const Sim_Param &, const Extrap_Pk_Nl< double, 2 > &, Data_Vec< double, 2 > &)

template void	gen_corr_func_binned_gsl_qawf (const Sim_Param &, const Extrap_Pk_Nl< double, 3 > &, Data_Vec< double, 2 > &)

template void	gen_sigma_binned_gsl_qawf (const Sim_Param &, const Extrap_Pk< double, 2 > &, Data_Vec< double, 2 > &)

template void	gen_sigma_binned_gsl_qawf (const Sim_Param &, const Extrap_Pk< double, 3 > &, Data_Vec< double, 2 > &)

template void	gen_sigma_binned_gsl_qawf (const Sim_Param &, const Extrap_Pk_Nl< double, 2 > &, Data_Vec< double, 2 > &)

template void	gen_sigma_binned_gsl_qawf (const Sim_Param &, const Extrap_Pk_Nl< double, 3 > &, Data_Vec< double, 2 > &)

Variables
constexpr double	anonymous_namespace{core_power.cpp}::GSL_EPSABS = 1e-7
	absolute error at z = 0, scales for higher redshifts More...

constexpr size_t	anonymous_namespace{core_power.cpp}::GSL_LIMIT = 1000
	max. number of subintervals for adaptive integration More...

constexpr size_t	anonymous_namespace{core_power.cpp}::GSL_N = 100
	max. number of bisections of integration interval (QAWO / QAWF) More...

Detailed Description

definitions of cosmological functions like power spectrum, growth, etc.

Author: Michal Vrastil

Date: 2018-06-23

Definition in file core_power.cpp.

Function Documentation

template<class P >

void gen_corr_func_binned_gsl_qawf	(	const Sim_Param &	sim,
		const P &	P_k,
		Data_Vec< double, 2 > &	corr_func_binned
	)

compute correlation function and store results

Template Parameters

P	callable object with 'operator()(double)'

Parameters

sim	simulation parameters
P_k	power spectrum (callable)
corr_func_binned	object to store binned correlation function

Definition at line 665 of file core_power.cpp.

References Sim_Param::cosmo, anonymous_namespace{core_power.cpp}::gen_corr_func_binned_gsl(), get_gsl_error(), anonymous_namespace{core_power.cpp}::GSL_LIMIT, and anonymous_namespace{core_power.cpp}::GSL_N.

Referenced by gen_corr_func_binned_gsl_qawf_lin(), gen_corr_func_binned_gsl_qawf_nl(), gen_sigma_func_binned_gsl_qawf_nl(), and App_Var< T >::Impl< T >::print_corr().

 {
     BOOST_LOG_TRIVIAL(debug) << "Computing correlation function via GSL integration QAWF...";
     const double gsp_epsabs = get_gsl_error(sim.cosmo, P_k);
     Integr_obj_qawf xi_r(&xi_integrand_W<P>, 0, gsp_epsabs,  GSL_LIMIT, GSL_N);
     gen_corr_func_binned_gsl(sim, P_k, corr_func_binned, xi_r);
 }

template void gen_corr_func_binned_gsl_qawf	(	const Sim_Param &	,
		const Extrap_Pk< double, 2 > &	,
		Data_Vec< double, 2 > &
	)

template void gen_corr_func_binned_gsl_qawf	(	const Sim_Param &	,
		const Extrap_Pk< double, 3 > &	,
		Data_Vec< double, 2 > &
	)

template void gen_corr_func_binned_gsl_qawf	(	const Sim_Param &	,
		const Extrap_Pk_Nl< double, 2 > &	,
		Data_Vec< double, 2 > &
	)

template void gen_corr_func_binned_gsl_qawf	(	const Sim_Param &	,
		const Extrap_Pk_Nl< double, 3 > &	,
		Data_Vec< double, 2 > &
	)

void gen_corr_func_binned_gsl_qawf_lin	(	const Sim_Param &	sim,
		double	a,
		Data_Vec< double, 2 > &	corr_func_binned
	)

compute linear correlation function and store results

Parameters

sim	simulation parameters
a	scale factor
corr_func_binned	object to store binned correlation function

Definition at line 673 of file core_power.cpp.

References Sim_Param::cosmo, gen_corr_func_binned_gsl_qawf(), and lin_pow_spec().

Referenced by App_Var< T >::Impl< T >::print_corr().

 {
     auto P_k = [&](FTYPE_t k){ return lin_pow_spec(a, k, sim.cosmo); };
     gen_corr_func_binned_gsl_qawf(sim, P_k, corr_func_binned);
 }

void gen_corr_func_binned_gsl_qawf_nl	(	const Sim_Param &	sim,
		double	a,
		Data_Vec< double, 2 > &	corr_func_binned
	)

compute non-linear correlation function and store results

Parameters

sim	simulation parameters
a	scale factor
corr_func_binned	object to store binned correlation function

Definition at line 679 of file core_power.cpp.

References Sim_Param::cosmo, gen_corr_func_binned_gsl_qawf(), and non_lin_pow_spec().

 {
     auto P_k = [&](FTYPE_t k){ return non_lin_pow_spec(a, k, sim.cosmo); };
     gen_corr_func_binned_gsl_qawf(sim, P_k, corr_func_binned);
 }

template<class P >

void gen_sigma_binned_gsl_qawf	(	const Sim_Param &	sim,
		const P &	P_k,
		Data_Vec< double, 2 > &	sigma_binned
	)

compute amplitude of density fluctuation and store results

Template Parameters

P	callable object with 'operator()(double)'

Parameters

sim	simulation parameters
P_k	power spectrum (callable)
sigma_binned	object to store binned correlation function

Definition at line 686 of file core_power.cpp.

References Sim_Param::cosmo, anonymous_namespace{core_power.cpp}::gen_sigma_func_binned_gsl(), get_gsl_error(), anonymous_namespace{core_power.cpp}::GSL_LIMIT, and anonymous_namespace{core_power.cpp}::GSL_N.

Referenced by gen_sigma_func_binned_gsl_qawf_lin(), and gen_sigma_func_binned_gsl_qawf_nl().

 {
     BOOST_LOG_TRIVIAL(debug) << "Computing mass fluctuations via GSL integration QAWF...";
     const double gsp_epsabs = get_gsl_error(sim.cosmo, P_k);
     Integr_obj_qagiu sigma_r(&sigma_integrand_G<P>, 0, gsp_epsabs,  GSL_LIMIT, GSL_N);
     gen_sigma_func_binned_gsl(sim, P_k, sigma_binned, sigma_r);
 }

template void gen_sigma_binned_gsl_qawf	(	const Sim_Param &	,
		const Extrap_Pk< double, 2 > &	,
		Data_Vec< double, 2 > &
	)

template void gen_sigma_binned_gsl_qawf	(	const Sim_Param &	,
		const Extrap_Pk< double, 3 > &	,
		Data_Vec< double, 2 > &
	)

template void gen_sigma_binned_gsl_qawf	(	const Sim_Param &	,
		const Extrap_Pk_Nl< double, 2 > &	,
		Data_Vec< double, 2 > &
	)

template void gen_sigma_binned_gsl_qawf	(	const Sim_Param &	,
		const Extrap_Pk_Nl< double, 3 > &	,
		Data_Vec< double, 2 > &
	)

void gen_sigma_func_binned_gsl_qawf_lin	(	const Sim_Param &	sim,
		double	a,
		Data_Vec< double, 2 > &	sigma_binned
	)

compute linear amplitude of density fluctuation and store results

Parameters

sim	simulation parameters
a	scale factor
sigma_binned	object to store binned correlation function

Definition at line 694 of file core_power.cpp.

References Sim_Param::cosmo, gen_sigma_binned_gsl_qawf(), and lin_pow_spec().

 {
     auto P_k = [&](FTYPE_t k){ return lin_pow_spec(a, k, sim.cosmo); };
     gen_sigma_binned_gsl_qawf(sim, P_k, sigma_binned);
 }

void gen_sigma_func_binned_gsl_qawf_nl	(	const Sim_Param &	sim,
		double	a,
		Data_Vec< double, 2 > &	sigma_binned
	)

compute non-linear amplitude of density fluctuation and store results

Parameters

sim	simulation parameters
a	scale factor
sigma_binned	object to store binned correlation function

Definition at line 700 of file core_power.cpp.

References Sim_Param::cosmo, gen_corr_func_binned_gsl_qawf(), gen_sigma_binned_gsl_qawf(), Interp_obj::init(), and non_lin_pow_spec().

 {
     auto P_k = [&](FTYPE_t k){ return non_lin_pow_spec(a, k, sim.cosmo); };
     gen_sigma_binned_gsl_qawf(sim, P_k, sigma_binned);
 }

template<class P >

static double get_gsl_error	(	const Cosmo_Param &	cosmo,
		const P &	P_k
	)

static

Definition at line 656 of file core_power.cpp.

References anonymous_namespace{core_power.cpp}::GSL_EPSABS, and lin_pow_spec().

Referenced by gen_corr_func_binned_gsl_qawf(), and gen_sigma_binned_gsl_qawf().

 {
     const FTYPE_t k_ref = 0.1;
     double ratio = P_k(k_ref) / lin_pow_spec(1, k_ref, cosmo);
     BOOST_LOG_TRIVIAL(debug) << "\tSet absolute error to " << GSL_EPSABS*ratio;
     return GSL_EPSABS*ratio;
 }

double growth_change	(	double	a,
		const Cosmo_Param &	cosmo
	)

Parameters

a
cosmo	cosmological parameters

Returns: FTYPE_t

Definition at line 447 of file core_power.cpp.

References cl_cmbl_bm::cosmo, growth_factor(), and growth_rate().

Referenced by kick_step_no_momentum(), App_Var_Chi::ChiImpl::kick_step_w_chi_no_momentum(), App_Var< T >::Impl< T >::print_vel_pwr(), set_pert_pos(), and TEST_CASE().

 {
     if (a){
         // compute with growth rate
         return growth_rate(a, cosmo)*growth_factor(a, cosmo)/a;
     }
     else{
         gsl_function F = {&growth_factor, static_cast<void*>(cosmo)};
         double dDda, error;
         int status = gsl_deriv_forward(&F, a, 1e-12, &dDda, &error);
         if (!status) return dDda;
         else throw std::runtime_error("GSL ODE error: " + std::string(gsl_strerror(status)));
     }
 }

double growth_factor	(	double	a,
		const Cosmo_Param &	cosmo
	)

Parameters

a
cosmo	cosmological parameters

Returns: FTYPE_t

Definition at line 414 of file core_power.cpp.

References ccl_growth_factor(), Cosmo_Param::cosmo, Cosmo_Param::D_norm, anonymous_namespace{core_power.cpp}::growth_factor_integrand(), and anonymous_namespace{core_power.cpp}::hubble_param().

Referenced by growth_change(), anonymous_namespace{core_power.cpp}::growth_factor(), lin_pow_spec(), and anonymous_namespace{core_power.cpp}::ln_growth_factor().

 {
     // D(0) == 0; D(1) == 1, do not check (not often, better performance)
     if (!a) return 0;
 
     // try ccl range
     int status = 0;
     FTYPE_t D_ccl = ccl_growth_factor(cosmo.cosmo, a, &status);
     if (!status) return D_ccl;
 
     // integrate outside range
     Integr_obj_qag D(&growth_factor_integrand, 0, a, 0, 1e-12, 1000, GSL_INTEG_GAUSS61);
     return hubble_param(a, cosmo)*D(static_cast<void*>(cosmo))/cosmo.D_norm;
 }

double growth_rate	(	double	a,
		const Cosmo_Param &	cosmo
	)

Parameters

a
cosmo	cosmological parameters

Returns: FTYPE_t

Definition at line 429 of file core_power.cpp.

References ccl_growth_rate(), cl_cmbl_bm::cosmo, Cosmo_Param::cosmo, and anonymous_namespace{core_power.cpp}::ln_growth_factor().

Referenced by growth_change(), and TEST_CASE().

 {
     // f(0) == 1
     if (!a) return 1;
 
     // try ccl range
     int status = 0;
     double f = ccl_growth_rate(cosmo.cosmo, a, &status);
     if (!status) return f;
     
     // logarithmic derivative outside range
     gsl_function F = {&ln_growth_factor, static_cast<void*>(cosmo)};
     double error;
     status = gsl_deriv_central(&F, log(a), 1e-12, &f, &error);
     if (!status) return f;
     else throw std::runtime_error("GSL ODE error: " + std::string(gsl_strerror(status)));
 }

double lin_pow_spec	(	double	a,
		double	k,
		const Cosmo_Param &	cosmo
	)

Parameters

a
k
cosmo	cosmological parameters

Returns: FTYPE_t

Definition at line 474 of file core_power.cpp.

References ccl_linear_matter_power(), Cosmo_Param::cosmo, growth_factor(), Cosmo_Param::h, pow(), and ccl_cosmology::status_message.

Referenced by Extrap_Pk< T, N >::fit_lin(), gen_corr_func_binned_gsl_qawf_lin(), gen_rho_w_pow_k(), gen_sigma_func_binned_gsl_qawf_lin(), get_gsl_error(), non_lin_pow_spec(), Extrap_Pk< T, N >::operator()(), and Extrap_Pk_Nl< T, N >::operator()().

 {
     int status = 0;
     FTYPE_t pk = ccl_linear_matter_power(cosmo.cosmo, k*cosmo.h, a, &status)*pow(cosmo.h, 3);
     if (!status) return pk;
     status = 0;
     pk = ccl_linear_matter_power(cosmo.cosmo, k*cosmo.h, 1, &status)*pow(cosmo.h, 3);
     FTYPE_t D = growth_factor(a, cosmo);
     if (!status) return D*D*pk;
     else throw std::runtime_error(cosmo.cosmo->status_message);
 }

double non_lin_pow_spec	(	double	a,
		double	k,
		const Cosmo_Param &	cosmo
	)

Parameters

a
k
cosmo	cosmological parameters

Returns: FTYPE_t

Definition at line 486 of file core_power.cpp.

References ccl_nonlin_matter_power(), Cosmo_Param::cosmo, Cosmo_Param::h, lin_pow_spec(), pow(), and ccl_cosmology::status_message.

Referenced by gen_corr_func_binned_gsl_qawf_nl(), gen_sigma_func_binned_gsl_qawf_nl(), and Extrap_Pk_Nl< T, N >::operator()().

 {
     if (a < 1e-3) return lin_pow_spec(a, k, cosmo);
     int status = 0;
     FTYPE_t pk = ccl_nonlin_matter_power(cosmo.cosmo, k*cosmo.h, a, &status)*pow(cosmo.h, 3);
     if (!status) return pk;
     else throw std::runtime_error(cosmo.cosmo->status_message);
 }

double norm_growth_factor ( const Cosmo_Param & cosmo )

when computing growth factor outside CCL range we need to normalize the growth factor; $D(a=1)\equiv1$

Parameters

cosmo cosmological parameters

Returns: FTYPE_t unnormalized growth factor at 'a = 1'

Definition at line 408 of file core_power.cpp.

References anonymous_namespace{core_power.cpp}::growth_factor_integrand().

Referenced by Cosmo_Param::init().

 {
     Integr_obj_qag D(&growth_factor_integrand, 0, 1, 0, 1e-12, 1000, GSL_INTEG_GAUSS61);
     return D(static_cast<void*>(cosmo));
 }

void norm_pwr ( Cosmo_Param & cosmo )

< end of anonymous namespace (private definitions)

declaration in class_data_vec.hpp

PUBLIC FUNCTIONS DEFINITIONS *

Definition at line 400 of file core_power.cpp.

References ccl_sigma8(), Cosmo_Param::cosmo, and ccl_cosmology::status_message.

Referenced by Cosmo_Param::init().

 {
     BOOST_LOG_TRIVIAL(debug) << "Initializing CCL power spectrum...";
     int status = 0;
     ccl_sigma8(cosmo.cosmo, &status);
     if (status) throw std::runtime_error(cosmo.cosmo->status_message);
 }

double Omega_lambda	(	double	a,
		const Cosmo_Param &	cosmo
	)

Parameters

a
cosmo	cosmological parameters

Returns: FTYPE_t

Definition at line 462 of file core_power.cpp.

References ccl_omega_x(), ccl_species_l_label, Cosmo_Param::cosmo, Cosmo_Param::Omega_L(), Cosmo_Param::Omega_m, and pow().

 {
     // try ccl range
     int status = 0;
     FTYPE_t OL = ccl_omega_x(cosmo.cosmo, a, ccl_species_l_label, &status);
     if(!status) return OL;
 
     // compute outside range
     OL = cosmo.Omega_L();
     return OL/(cosmo.Omega_m/pow(a, 3) + OL);
 }

static double truncation_fce	(	double	k,
		double	k2_G
	)

static

Definition at line 554 of file core_power.cpp.

Referenced by Extrap_Pk< T, N >::fit_lin().

 {
     return exp(-k*k/k2_G);
 }

Classes

Namespaces

Functions

Variables

Detailed Description

Function Documentation