cosmology.h File Reference

#include <stddef.h>
#include "omega_nu_single.h"
#include "utils/unitsystem.h"

Include dependency graph for cosmology.h:

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

Go to the source code of this file.

Classes
struct	Cosmology
struct	FunctionOfK

Functions
void	function_of_k_normalize_sigma (FunctionOfK *fk, double R, double sigma)
double	function_of_k_eval (FunctionOfK *fk, double k)
double	function_of_k_tophat_sigma (FunctionOfK *fk, double R)
double	hubble_function (const Cosmology *CP, double a)
double	GrowthFactor (Cosmology *CP, double astart, double aend)
double	F_Omega (Cosmology *CP, double a)
int	hybrid_nu_tracer (const Cosmology *CP, double atime)
void	init_cosmology (Cosmology *CP, double TimeBegin, const struct UnitSystem units)
void	check_units (const Cosmology *CP, const struct UnitSystem units)

Function Documentation

check_units()

void check_units	(	const Cosmology *	CP,
		const struct UnitSystem	units
	)

Check and print properties of the cosmological unit system.

Definition at line 260 of file cosmology.c.

{
    /* Detect cosmologies that are likely to be typos in the parameter files*/
    if(CP->HubbleParam < 0.1 || CP->HubbleParam > 10 ||
        CP->OmegaLambda < 0 || CP->OmegaBaryon < 0 || CP->OmegaG < 0 || CP->OmegaCDM < 0)
        endrun(5, "Bad cosmology: H0 = %g OL = %g Ob = %g Og = %g Ocdm = %g\n",
               CP->HubbleParam, CP->OmegaLambda, CP->OmegaBaryon, CP->OmegaCDM);
 
    message(0, "Hubble (internal units) = %g\n", CP->Hubble);
    message(0, "G (internal units) = %g\n", CP->GravInternal);
    message(0, "UnitLength_in_cm = %g \n", units.UnitLength_in_cm);
    message(0, "UnitMass_in_g = %g \n", units.UnitMass_in_g);
    message(0, "UnitTime_in_s = %g \n", units.UnitTime_in_s);
    message(0, "UnitVelocity_in_cm_per_s = %g \n", units.UnitVelocity_in_cm_per_s);
    message(0, "UnitDensity_in_cgs = %g \n", units.UnitDensity_in_cgs);
    message(0, "UnitEnergy_in_cgs = %g \n", units.UnitEnergy_in_cgs);
    message(0, "Dark energy model: OmegaL = %g OmegaFLD = %g\n",CP->OmegaLambda, CP->Omega_fld);
    message(0, "Photon density OmegaG = %g\n",CP->OmegaG);
    if(!CP->MassiveNuLinRespOn)
        message(0, "Massless Neutrino density OmegaNu0 = %g\n",get_omega_nu(&CP->ONu, 1));
    message(0, "Curvature density OmegaK = %g\n",CP->OmegaK);
    if(CP->RadiationOn) {
        /* note that this value is inaccurate if there is a massive neutrino. */
        double OmegaTot = CP->OmegaG + CP->OmegaK + CP->Omega0 + CP->OmegaLambda;
        if(!CP->MassiveNuLinRespOn)
            OmegaTot += get_omega_nu(&CP->ONu, 1);
        message(0, "Radiation is enabled in Hubble(a). "
               "Following CAMB convention: Omega_Tot - 1 = %g\n", OmegaTot - 1);
    }
    message(0, "\n");
}

References CP, endrun(), get_omega_nu(), Cosmology::GravInternal, Cosmology::Hubble, Cosmology::HubbleParam, Cosmology::MassiveNuLinRespOn, message(), Cosmology::Omega0, Cosmology::Omega_fld, Cosmology::OmegaBaryon, Cosmology::OmegaCDM, Cosmology::OmegaG, Cosmology::OmegaK, Cosmology::OmegaLambda, Cosmology::ONu, Cosmology::RadiationOn, UnitSystem::UnitDensity_in_cgs, UnitSystem::UnitEnergy_in_cgs, UnitSystem::UnitLength_in_cm, UnitSystem::UnitMass_in_g, UnitSystem::UnitTime_in_s, and UnitSystem::UnitVelocity_in_cm_per_s.

Referenced by begrun().

Here is the call graph for this function:

Here is the caller graph for this function:

F_Omega()

double F_Omega	(	Cosmology *	CP,
		double	a
	)

Definition at line 148 of file cosmology.c.

{
    double dD1da=0;
    double D1 = growth(CP, a, &dD1da);
    return a / D1 * dD1da;
}

References CP, and growth().

Referenced by displacement_fields(), test_cosmology(), and test_growth_numerical().

Here is the call graph for this function:

Here is the caller graph for this function:

function_of_k_eval()

double function_of_k_eval	(	FunctionOfK *	fk,
		double	k
	)

Definition at line 191 of file cosmology.c.

{
    /* ignore the 0 mode */
 
    if(k == 0) return 1;
 
    int l = 0;
    int r = fk->size - 1;
 
    while(r - l > 1) {
        int m = (r + l) / 2;
        if(k < fk->table[m].k)
            r = m;
        else
            l = m;
    }
    double k2 = fk->table[r].k,
           k1 = fk->table[l].k;
    double p2 = fk->table[r].Pk,
           p1 = fk->table[l].Pk;
 
    if(l == r) {
        return fk->table[l].Pk;
    }
 
    if(p1 == 0 || p2 == 0 || k1 == 0 || k2 == 0) {
        /* if any of the p is zero, use linear interpolation */
        double p = (k - k1) * p2 + (k2 - k) * p1;
        p /= (k2 - k1);
        return p;
    } else {
        k = log(k);
        p1 = log(p1);
        p2 = log(p2);
        k1 = log(k1);
        k2 = log(k2);
        double p = (k - k1) * p2 + (k2 - k) * p1;
        p /= (k2 - k1);
        return exp(p);
    }
}

References sigma2_params::fk, FunctionOfK::k, FunctionOfK::Pk, FunctionOfK::size, and FunctionOfK::table.

Referenced by sigma2_int().

Here is the caller graph for this function:

function_of_k_normalize_sigma()

void function_of_k_normalize_sigma	(	FunctionOfK *	fk,
		double	R,
		double	sigma
	)

Definition at line 249 of file cosmology.c.

                                                                             {
    double old = function_of_k_tophat_sigma(fk, R);
    size_t i;
    for(i = 0; i < fk->size; i ++) {
        fk->table[i].Pk *= sigma / old;
    };
}

References sigma2_params::fk, function_of_k_tophat_sigma(), FunctionOfK::Pk, sigma2_params::R, sigma, FunctionOfK::size, and FunctionOfK::table.

Here is the call graph for this function:

function_of_k_tophat_sigma()

double function_of_k_tophat_sigma	(	FunctionOfK *	fk,
		double	R
	)

Definition at line 233 of file cosmology.c.

{
    gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);
    struct sigma2_params params = {fk, R};
    double result,abserr;
    gsl_function F;
    F.function = &sigma2_int;
    F.params = &params;
 
    /* note: 500/R is here chosen as integration boundary (infinity) */
    gsl_integration_qags (&F, 0, 500. / R, 0, 1e-4,1000,w,&result, &abserr);
    //   printf("gsl_integration_qng in TopHatSigma2. Result %g, error: %g, intervals: %lu\n",result, abserr,w->size);
    gsl_integration_workspace_free (w);
    return sqrt(result);
}

References sigma2_params::fk, sigma2_params::R, and sigma2_int().

Referenced by function_of_k_normalize_sigma().

Here is the call graph for this function:

Here is the caller graph for this function:

GrowthFactor()

double GrowthFactor	(	Cosmology *	CP,
		double	astart,
		double	aend
	)

Definition at line 82 of file cosmology.c.

{
    return growth(CP, astart, NULL) / growth(CP, aend, NULL);
}

References CP, and growth().

Referenced by gravpm_force(), init_powerspectrum(), print_spec(), and test_cosmology().

Here is the call graph for this function:

Here is the caller graph for this function:

hubble_function()

double hubble_function	(	const Cosmology *	CP,
		double	a
	)

Definition at line 58 of file cosmology.c.

{
 
    double hubble_a;
 
    /* first do the terms in SQRT */
    hubble_a = CP->OmegaLambda;
 
    hubble_a += OmegaFLD(CP, a);
    hubble_a += CP->OmegaK / (a * a);
    hubble_a += (CP->OmegaCDM + CP->OmegaBaryon) / (a * a * a);
 
    if(CP->RadiationOn) {
        hubble_a += CP->OmegaG / (a * a * a * a);
        hubble_a += get_omega_nu(&CP->ONu, a);
    }
    hubble_a += CP->Omega_ur/(a*a*a*a);
    /* Now finish it up. */
    hubble_a = CP->Hubble * sqrt(hubble_a);
    return (hubble_a);
}

References CP, get_omega_nu(), Cosmology::Hubble, Cosmology::Omega0, Cosmology::Omega_ur, Cosmology::OmegaBaryon, Cosmology::OmegaCDM, OmegaFLD(), Cosmology::OmegaG, Cosmology::OmegaK, Cosmology::OmegaLambda, Cosmology::ONu, and Cosmology::RadiationOn.

Referenced by atime_integ(), blackhole(), cooling_and_starformation(), displacement_fields(), drift_integ(), fac_integ(), find_timesteps(), fof_save_particles(), fof_write_header(), fslength_int(), get_delta_nu(), get_delta_nu_int(), get_long_range_timestep_dloga(), gravkick_integ(), growth(), growth_ode(), hydro_force(), hydrokick_integ(), init_transfer_table(), kernel(), petaio_read_snapshot(), petaio_save_snapshot(), petaio_write_header(), and test_cosmology().

Here is the call graph for this function:

Here is the caller graph for this function:

hybrid_nu_tracer()

int hybrid_nu_tracer	(	const Cosmology *	CP,
		double	atime
	)

Definition at line 53 of file cosmology.c.

{
    return CP->HybridNeutrinosOn && (atime <= CP->HybridNuPartTime);
}

References CP, Cosmology::HybridNeutrinosOn, and Cosmology::HybridNuPartTime.

Referenced by run(), runfof(), and runpower().

Here is the caller graph for this function:

init_cosmology()

void init_cosmology	(	Cosmology *	CP,
		double	TimeBegin,
		const struct UnitSystem	units
	)

Definition at line 15 of file cosmology.c.

{
    CP->Hubble = HUBBLE * units.UnitTime_in_s;
    CP->UnitTime_in_s = units.UnitTime_in_s;
    CP->GravInternal = GRAVITY / pow(units.UnitLength_in_cm, 3) * units.UnitMass_in_g * pow(units.UnitTime_in_s, 2);
 
    /*With slightly relativistic massive neutrinos, for consistency we need to include radiation.
     * A note on normalisation (as of 08/02/2012):
     * CAMB appears to set Omega_Lambda + Omega_Matter+Omega_K = 1,
     * calculating Omega_K in the code and specifying Omega_Lambda and Omega_Matter in the paramfile.
     * This means that Omega_tot = 1+ Omega_r + Omega_g, effectively
     * making h0 (very) slightly larger than specified, and the Universe is no longer flat!
     */
    CP->OmegaCDM = CP->Omega0 - CP->OmegaBaryon;
    CP->OmegaK = 1.0 - CP->Omega0 - CP->OmegaLambda;
 
    /* Omega_g = 4 \sigma_B T_{CMB}^4 8 \pi G / (3 c^3 H^2) */
 
    CP->OmegaG = 4 * STEFAN_BOLTZMANN
                  * pow(CP->CMBTemperature, 4)
                  * (8 * M_PI * GRAVITY)
                  / (3*pow(LIGHTCGS, 3)*HUBBLE*HUBBLE)
                  / (CP->HubbleParam*CP->HubbleParam);
 
    init_omega_nu(&CP->ONu, CP->MNu, TimeBegin, CP->HubbleParam, CP->CMBTemperature);
    /*Initialise the hybrid neutrinos, after Omega_nu*/
    if(CP->HybridNeutrinosOn)
        init_hybrid_nu(&CP->ONu.hybnu, CP->MNu, CP->HybridVcrit, LIGHTCGS/1e5, CP->HybridNuPartTime, CP->ONu.kBtnu);
 
    /* Neutrinos will be included in Omega0, if massive.
     * This ensures that OmegaCDM contains only non-relativistic species.*/
    if(CP->MNu[0] + CP->MNu[1] + CP->MNu[2] > 0) {
        CP->OmegaCDM -= get_omega_nu(&CP->ONu, 1);
    }
}

References Cosmology::CMBTemperature, CP, get_omega_nu(), Cosmology::GravInternal, GRAVITY, Cosmology::Hubble, HUBBLE, Cosmology::HubbleParam, _omega_nu::hybnu, Cosmology::HybridNeutrinosOn, Cosmology::HybridNuPartTime, Cosmology::HybridVcrit, init_hybrid_nu(), init_omega_nu(), _omega_nu::kBtnu, LIGHTCGS, Cosmology::MNu, Cosmology::Omega0, Cosmology::OmegaBaryon, Cosmology::OmegaCDM, Cosmology::OmegaG, Cosmology::OmegaK, Cosmology::OmegaLambda, Cosmology::ONu, STEFAN_BOLTZMANN, UnitSystem::UnitLength_in_cm, UnitSystem::UnitMass_in_g, Cosmology::UnitTime_in_s, and UnitSystem::UnitTime_in_s.

Referenced by begrun(), do_density_test(), do_force_test(), main(), setup(), and setup_cosmology().

Here is the call graph for this function:

Here is the caller graph for this function: