omega_nu_single.c File Reference

#include "omega_nu_single.h"
#include <math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <string.h>
#include "physconst.h"
#include "utils/mymalloc.h"
#include "utils/endrun.h"

Include dependency graph for omega_nu_single.c:

Go to the source code of this file.

Macros
#define	HBAR   6.582119e-16 /*hbar in units of eV s*/
#define	STEFAN_BOLTZMANN   5.670373e-5
#define	NRHOTAB   200
#define	FLOAT_ACC   1e-6
#define	GSL_VAL   200
#define	NU_SW   100

Functions
void	init_omega_nu (_omega_nu *omnu, const double MNu[], const double a0, const double HubbleParam, const double tcmb0)
double	get_omega_nu (const _omega_nu *const omnu, const double a)
double	get_omega_nu_nopart (const _omega_nu *const omnu, const double a)
double	get_omegag (const _omega_nu *const omnu, const double a)
double	rho_nu_int (double q, void *params)
double	get_rho_nu_conversion ()
void	rho_nu_init (_rho_nu_single *const rho_nu_tab, double a0, const double mnu, const double kBtnu)
static double	non_rel_rho_nu (const double a, const double kT, const double amnu, const double kTamnu2)
static double	rel_rho_nu (const double a, const double kT)
double	rho_nu (const _rho_nu_single *rho_nu_tab, const double a, const double kT)
double	fermi_dirac_kernel (double x, void *params)
double	nufrac_low (const double qc)
void	init_hybrid_nu (_hybrid_nu *const hybnu, const double mnu[], const double vcrit, const double light, const double nu_crit_time, const double kBtnu)
double	particle_nu_fraction (const _hybrid_nu *const hybnu, const double a, const int i)
double	omega_nu_single (const _omega_nu *const omnu, const double a, int i)

Macro Definition Documentation

FLOAT_ACC

#define FLOAT_ACC   1e-6

Floating point accuracy

Definition at line 16 of file omega_nu_single.c.

GSL_VAL

#define GSL_VAL   200

Number of bins in integrations

Definition at line 18 of file omega_nu_single.c.

HBAR

#define HBAR   6.582119e-16 /*hbar in units of eV s*/

Definition at line 11 of file omega_nu_single.c.

NRHOTAB

#define NRHOTAB   200

Definition at line 14 of file omega_nu_single.c.

NU_SW

#define NU_SW   100

Value of kT/aM_nu on which to switch from the analytic expansion to the numerical integration

Definition at line 87 of file omega_nu_single.c.

STEFAN_BOLTZMANN

#define STEFAN_BOLTZMANN   5.670373e-5

Definition at line 12 of file omega_nu_single.c.

Function Documentation

fermi_dirac_kernel()

double fermi_dirac_kernel	(	double	x,
		void *	params
	)

Definition at line 212 of file omega_nu_single.c.

{
  return x * x / (exp(x) + 1);
}

Referenced by nufrac_low().

Here is the caller graph for this function:

get_omega_nu()

double get_omega_nu	(	const _omega_nu *const	omnu,
		const double	a
	)

Return the total matter density in neutrinos at scale factor a.

Definition at line 57 of file omega_nu_single.c.

{
        double rhonu=0;
        int mi;
        for(mi=0; mi<NUSPECIES; mi++) {
            if(omnu->nu_degeneracies[mi] > 0){
                 rhonu += omnu->nu_degeneracies[mi] * rho_nu(&omnu->RhoNuTab[mi], a, omnu->kBtnu);
            }
        }
        return rhonu/omnu->rhocrit;
}

References _omega_nu::kBtnu, _omega_nu::nu_degeneracies, NUSPECIES, rho_nu(), _omega_nu::rhocrit, and _omega_nu::RhoNuTab.

Referenced by check_units(), compute_mass(), delta_nu_from_power(), delta_tot_first_init(), get_long_range_timestep_dloga(), get_omega_nu_nopart(), growth(), hubble_function(), init_cosmology(), init_neutrinos_lra(), init_transfer_table(), main(), parse_transfer(), test_get_omega_nu(), and update_delta_tot().

Here is the call graph for this function:

Here is the caller graph for this function:

get_omega_nu_nopart()

double get_omega_nu_nopart	(	const _omega_nu *const	omnu,
		const double	a
	)

Return the total matter density in neutrinos at scale factor a , excluding active particles.

Definition at line 71 of file omega_nu_single.c.

{
    double omega_nu = get_omega_nu(omnu, a);
    double part_nu = get_omega_nu(omnu, 1) * particle_nu_fraction(&omnu->hybnu, a, 0) / (a*a*a);
    return omega_nu - part_nu;
}

References get_omega_nu(), _omega_nu::hybnu, and particle_nu_fraction().

Referenced by check_omega(), delta_nu_from_power(), delta_tot_first_init(), get_delta_nu_combined(), potential_transfer(), test_hybrid_neutrinos(), and update_delta_tot().

Here is the call graph for this function:

Here is the caller graph for this function:

get_omegag()

double get_omegag	(	const _omega_nu *const	omnu,
		const double	a
	)

Return the photon matter density at scale factor a

Definition at line 79 of file omega_nu_single.c.

{
    const double omegag = 4*STEFAN_BOLTZMANN/(LIGHTCGS*LIGHTCGS*LIGHTCGS)*pow(omnu->tcmb0,4)/omnu->rhocrit;
    return omegag/pow(a,4);
}

References LIGHTCGS, _omega_nu::rhocrit, STEFAN_BOLTZMANN, and _omega_nu::tcmb0.

Referenced by test_get_omegag().

Here is the caller graph for this function:

get_rho_nu_conversion()

double get_rho_nu_conversion

(

)

Definition at line 102 of file omega_nu_single.c.

{
        /*q has units of eV/c, so rho_nu now has units of (eV/c)^4*/
        double convert=4*M_PI*2; /* The factor of two is for antineutrinos*/
        /*rho_nu_val now has units of eV^4*/
        /*To get units of density, divide by (c*hbar)**3 in eV s and cm/s */
        const double chbar=1./(2*M_PI*LIGHTCGS*HBAR);
        convert*=(chbar*chbar*chbar);
        /*Now has units of (eV)/(cm^3)*/
        /* 1 eV = 1.60217646 × 10-12 g cm^2 s^(-2) */
        /* So 1eV/c^2 = 1.7826909604927859e-33 g*/
        /*So this is rho_nu_val in g /cm^3*/
        convert*=(1.60217646e-12/LIGHTCGS/LIGHTCGS);
        return convert;
}

References HBAR, and LIGHTCGS.

Referenced by do_exact_rho_nu_integration(), non_rel_rho_nu(), rel_rho_nu(), and rho_nu_init().

Here is the caller graph for this function:

init_hybrid_nu()

void init_hybrid_nu	(	_hybrid_nu *const	hybnu,
		const double	mnu[],
		const double	vcrit,
		const double	light,
		const double	nu_crit_time,
		const double	kBtnu
	)

Set up parameters for the hybrid neutrinos

Parameters

hybnu	initialised structure
mnu	array of neutrino masses in eV
vcrit	Critical velocity above which to treat neutrinos with particles. Note this is unperturbed velocity TODAY To get velocity at redshift z, multiply by (1+z)
light	speed of light in internal units
nu_crit_time	critical time to make neutrino particles live
kBtnu	Boltzmann constant times neutrino temperature. Dimensionful factor.

Definition at line 235 of file omega_nu_single.c.

{
    hybnu->enabled=1;
    int i;
    hybnu->nu_crit_time = nu_crit_time;
    hybnu->vcrit = vcrit / light;
    for(i=0; i< NUSPECIES; i++) {
        const double qc = mnu[i] * vcrit / light / kBtnu;
        hybnu->nufrac_low[i] = nufrac_low(qc);
    }
}

References _hybrid_nu::enabled, _hybrid_nu::nu_crit_time, nufrac_low(), _hybrid_nu::nufrac_low, NUSPECIES, and _hybrid_nu::vcrit.

Referenced by init_cosmology(), and test_hybrid_neutrinos().

Here is the call graph for this function:

Here is the caller graph for this function:

init_omega_nu()

void init_omega_nu	(	_omega_nu *const	omnu,
		const double	MNu[],
		const double	a0,
		const double	HubbleParam,
		const double	tcmb0
	)

Initialise the neutrino structure, do the time integration and allocate memory for the subclass rho_nu_single

Parameters

omnu	structure to initialise
MNu	array of neutrino masses in eV. Three entries.
a0	initial scale factor.
HubbleParam	Hubble parameter h0, eg, 0.7.
tcmb0	Redshift zero CMB temperature.

Definition at line 20 of file omega_nu_single.c.

{
    int mi;
    /*Explicitly disable hybrid neutrinos*/
    omnu->hybnu.enabled=0;
    /*CMB temperature*/
    omnu->tcmb0 = tcmb0;
    /*Neutrino temperature times k_B*/
    omnu->kBtnu = BOLEVK * TNUCMB * tcmb0;
    /*Store conversion between rho and omega*/
    omnu->rhocrit = (3 * HUBBLE * HubbleParam * HUBBLE * HubbleParam)/ (8 * M_PI * GRAVITY);
    /*First compute which neutrinos are degenerate with each other*/
    for(mi=0; mi<NUSPECIES; mi++){
        int mmi;
        omnu->nu_degeneracies[mi]=0;
        for(mmi=0; mmi<mi; mmi++){
            if(fabs(MNu[mi] -MNu[mmi]) < FLOAT_ACC){
                omnu->nu_degeneracies[mmi]+=1;
                break;
            }
        }
        if(mmi==mi) {
            omnu->nu_degeneracies[mi]=1;
        }
    }
    /*Now allocate a table for the species we want*/
    for(mi=0; mi<NUSPECIES; mi++){
        if(omnu->nu_degeneracies[mi]) {
            rho_nu_init(&omnu->RhoNuTab[mi], a0, MNu[mi], omnu->kBtnu);
        }
        else
            omnu->RhoNuTab[mi].loga = 0;
    }
}

References BOLEVK, _hybrid_nu::enabled, FLOAT_ACC, GRAVITY, HUBBLE, _omega_nu::hybnu, _omega_nu::kBtnu, _rho_nu_single::loga, _omega_nu::nu_degeneracies, NUSPECIES, rho_nu_init(), _omega_nu::rhocrit, _omega_nu::RhoNuTab, _omega_nu::tcmb0, and TNUCMB.

Referenced by init_cosmology(), test_allocate_delta_tot_table(), test_get_omega_nu(), test_get_omegag(), test_hybrid_neutrinos(), test_omega_nu_init_degenerate(), test_omega_nu_init_nondeg(), test_omega_nu_single(), and test_omega_nu_single_exact().

Here is the call graph for this function:

Here is the caller graph for this function:

non_rel_rho_nu()

static double non_rel_rho_nu ( const double  a, const double  kT, const double  amnu, const double  kTamnu2  )

inlinestatic

Definition at line 165 of file omega_nu_single.c.

{
    /*The constants are Riemann zetas: 3,5,7 respectively*/
    return amnu*(kT*kT*kT)/(a*a*a*a)*(1.5*1.202056903159594+kTamnu2*45./4.*1.0369277551433704+2835./32.*kTamnu2*kTamnu2*1.0083492773819229+80325/32.*kTamnu2*kTamnu2*kTamnu2*1.0020083928260826)*get_rho_nu_conversion();
}

References get_rho_nu_conversion().

Referenced by rho_nu().

Here is the call graph for this function:

Here is the caller graph for this function:

nufrac_low()

double nufrac_low

(

const double

qc

)

Integrate the fermi-dirac kernel between 0 and qc to find the fraction of neutrinos that are particles

Definition at line 219 of file omega_nu_single.c.

{
    /*These functions are so smooth that we don't need much space*/
    gsl_integration_workspace * w = gsl_integration_workspace_alloc (100);
    double abserr;
    gsl_function F;
    F.function = &fermi_dirac_kernel;
    F.params = NULL;
    double total_fd;
    gsl_integration_qag (&F, 0, qc, 0, 1e-6,100,GSL_INTEG_GAUSS61, w,&(total_fd), &abserr);
    /*divided by the total F-D probability (which is 3 Zeta(3)/2 ~ 1.8 if MAX_FERMI_DIRAC is large enough*/
    total_fd /= 1.5*1.202056903159594;
    gsl_integration_workspace_free (w);
    return total_fd;
}

References fermi_dirac_kernel().

Referenced by init_hybrid_nu(), Jfrac_high(), test_hybrid_neutrinos(), and test_nufrac_low().

Here is the call graph for this function:

Here is the caller graph for this function:

omega_nu_single()

double omega_nu_single	(	const _omega_nu *const	rho_nu_tab,
		const double	a,
		const int	i
	)

Return the matter density in a single neutrino species

Parameters

rho_nu_tab	structure containing pre-computed matter density values.
i	index of neutrino species we want
a	scale factor desired

Definition at line 267 of file omega_nu_single.c.

{
    /*Deal with case where we want a species degenerate with another one*/
    if(omnu->nu_degeneracies[i] == 0) {
        int j;
        for(j=i; j >=0; j--)
            if(omnu->nu_degeneracies[j]){
                i = j;
                break;
            }
    }
    double omega_nu = rho_nu(&omnu->RhoNuTab[i], a,omnu->kBtnu)/omnu->rhocrit;
    double omega_part = rho_nu(&omnu->RhoNuTab[i], 1,omnu->kBtnu)/omnu->rhocrit;
    omega_part *= particle_nu_fraction(&omnu->hybnu, a, i)/(a*a*a);
    omega_nu -= omega_part;
    return omega_nu;
 
}

References _omega_nu::hybnu, _omega_nu::kBtnu, _omega_nu::nu_degeneracies, particle_nu_fraction(), rho_nu(), _omega_nu::rhocrit, and _omega_nu::RhoNuTab.

Referenced by get_delta_nu_combined(), parse_transfer(), test_get_omega_nu(), test_hybrid_neutrinos(), test_omega_nu_single(), and test_omega_nu_single_exact().

Here is the call graph for this function:

Here is the caller graph for this function:

particle_nu_fraction()

double particle_nu_fraction	(	const _hybrid_nu *const	hybnu,
		const double	a,
		int	i
	)

Get fraction of neutrinos currently followed by particles.

Parameters

hybnu	Structure with hybrid neutrino parameters.
i	index of neutrino species to use.
a	redshift of interest.

Returns

the fraction of neutrinos currently traced by particles. 0 when neutrinos are fully analytic at early times.

Definition at line 251 of file omega_nu_single.c.

{
    /*Return zero if hybrid neutrinos not enabled*/
    if(!hybnu->enabled)
        return 0;
    /*Just use a redshift cut for now. Really we want something more sophisticated,
     * based on the shot noise and average overdensity.*/
    if (a > hybnu->nu_crit_time){
        return hybnu->nufrac_low[i];
    }
    return 0;
}

References _hybrid_nu::enabled, _hybrid_nu::nu_crit_time, and _hybrid_nu::nufrac_low.

Referenced by delta_nu_from_power(), delta_tot_first_init(), get_delta_nu(), get_delta_tot(), get_omega_nu_nopart(), gravpm_force(), omega_nu_single(), test_hybrid_neutrinos(), and update_delta_tot().

Here is the caller graph for this function:

rel_rho_nu()

static double rel_rho_nu ( const double  a, const double  kT  )

inlinestatic

Definition at line 173 of file omega_nu_single.c.

{
    return 7*pow(M_PI*kT/a,4)/120.*get_rho_nu_conversion();
}

References get_rho_nu_conversion().

Referenced by rho_nu().

Here is the call graph for this function:

Here is the caller graph for this function:

rho_nu()

double rho_nu	(	const _rho_nu_single *	rho_nu_tab,
		const double	a,
		const double	kT
	)

Computes the neutrino density for a single neutrino species at a given redshift, either by looking up in a table, or a simple calculation in the limits, or by direct integration.

Parameters

rho_nu_tab	Pre-computed table of values
a	Redshift desired.
kT	Boltzmann constant times neutrino temperature. Dimensionful factor.

Returns

neutrino density in g cm^-3

Definition at line 180 of file omega_nu_single.c.

{
        double rho_nu_val;
        double amnu=a*rho_nu_tab->mnu;
        const double kTamnu2=(kT*kT/amnu/amnu);
        /*Do it analytically if we are in a regime where we can
         * The next term is 141682 (kT/amnu)^8.
         * At kT/amnu = 8, higher terms are larger and the series stops converging.
         * Don't go lower than 50 here. */
        if(NU_SW*NU_SW*kTamnu2 < 1){
            /*Heavily non-relativistic*/
            rho_nu_val = non_rel_rho_nu(a, kT, amnu, kTamnu2);
        }
        else if(amnu < 1e-6*kT){
            /*Heavily relativistic*/
            rho_nu_val=rel_rho_nu(a, kT);
        }
        else{
            const double loga = log(a);
            /* Deal with early time case. In practice no need to be very accurate
             * so assume relativistic.*/
            if (!rho_nu_tab->loga || loga < rho_nu_tab->loga[0])
                rho_nu_val = rel_rho_nu(a,kT);
            else
                rho_nu_val=gsl_interp_eval(rho_nu_tab->interp,rho_nu_tab->loga,rho_nu_tab->rhonu,loga,rho_nu_tab->acc);
        }
        return rho_nu_val;
}

References _rho_nu_single::acc, _rho_nu_single::interp, _rho_nu_single::loga, _rho_nu_single::mnu, non_rel_rho_nu(), NU_SW, rel_rho_nu(), and _rho_nu_single::rhonu.

Referenced by get_omega_nu(), and omega_nu_single().

Here is the call graph for this function:

Here is the caller graph for this function:

rho_nu_init()

void rho_nu_init	(	_rho_nu_single *	rho_nu_tab,
		double	a0,
		const double	mnu,
		const double	kBtnu
	)

Initialise the tables for matter density in a single neutrino species, by doing numerical integration.

Parameters

rho_nu_tab	Structure to initialise
a0	First scale factor in rho_nu_tab. Do not try and evaluate rho_nu at higher redshift!
mnu	neutrino mass to compute neutrino matter density for in eV
kBtnu	Boltzmann constant times neutrino temperature. Dimensionful factor.

Definition at line 119 of file omega_nu_single.c.

{
     int i;
     double abserr;
     /* Tabulate to 1e-3, unless earlier requested.
      * Need early times for growth function.*/
     if(a0 > 1e-3)
         a0 = 1e-3;
     /* Do not need a table when relativistic*/
     if(a0 * mnu < 1e-6 * kBtnu)
        a0 = 1e-6 * kBtnu / mnu;
     /*Make the table over a slightly wider range than requested, in case there is roundoff error*/
     const double logA0=log(a0)-log(1.2);
     const double logaf=log(NU_SW*kBtnu/mnu)+log(1.2);
     gsl_function F;
     F.function = &rho_nu_int;
     /*Initialise constants*/
     rho_nu_tab->mnu = mnu;
     /*Shortcircuit if we don't need to do the integration*/
     if(mnu < 1e-6*kBtnu || logaf < logA0)
         return;
 
     /*Allocate memory for arrays*/
     rho_nu_tab->loga = (double *) mymalloc("rho_nu_table",2*NRHOTAB*sizeof(double));
     rho_nu_tab->rhonu = rho_nu_tab->loga+NRHOTAB;
     rho_nu_tab->acc = gsl_interp_accel_alloc();
     rho_nu_tab->interp=gsl_interp_alloc(gsl_interp_cspline,NRHOTAB);
     if(!rho_nu_tab->interp || !rho_nu_tab->acc || !rho_nu_tab->loga)
         endrun(2035,"Could not initialise tables for neutrino matter density\n");
 
     gsl_integration_workspace * w = gsl_integration_workspace_alloc (GSL_VAL);
     for(i=0; i< NRHOTAB; i++){
        double param[2];
        rho_nu_tab->loga[i]=logA0+i*(logaf-logA0)/(NRHOTAB-1);
        param[0]=mnu*exp(rho_nu_tab->loga[i]);
        param[1] = kBtnu;
        F.params = &param;
        gsl_integration_qag (&F, 0, 500*kBtnu,0 , 1e-9,GSL_VAL,6,w,&(rho_nu_tab->rhonu[i]), &abserr);
        rho_nu_tab->rhonu[i]=rho_nu_tab->rhonu[i]/pow(exp(rho_nu_tab->loga[i]),4)*get_rho_nu_conversion();
     }
     gsl_integration_workspace_free (w);
     gsl_interp_init(rho_nu_tab->interp,rho_nu_tab->loga,rho_nu_tab->rhonu,NRHOTAB);
     return;
}

References _rho_nu_single::acc, endrun(), get_rho_nu_conversion(), GSL_VAL, _rho_nu_single::interp, _rho_nu_single::loga, _rho_nu_single::mnu, mymalloc, NRHOTAB, NU_SW, rho_nu_int(), and _rho_nu_single::rhonu.

Referenced by init_omega_nu(), and test_rho_nu_init().

Here is the call graph for this function:

Here is the caller graph for this function:

rho_nu_int()

double rho_nu_int	(	double	q,
		void *	params
	)

Definition at line 91 of file omega_nu_single.c.

{
        double amnu = *((double *)params);
        double kT = *((double *)params+1);
        double epsilon = sqrt(q*q+amnu*amnu);
        double f0 = 1./(exp(q/kT)+1);
        return q*q*epsilon*f0;
}

Referenced by do_exact_rho_nu_integration(), and rho_nu_init().

Here is the caller graph for this function: