ccl_neutrinos.c File Reference

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_const_mksa.h>
#include <gsl/gsl_roots.h>
#include "ccl.h"
#include "ccl_params.h"

Include dependency graph for ccl_neutrinos.c:

Go to the source code of this file.

Functions
static double	nu_integrand (double x, void *r)
gsl_spline *	calculate_nu_phasespace_spline (int *status)
double	nu_phasespace_intg (gsl_interp_accel *accel, double mnuOT, int *status)
double	ccl_Omeganuh2 (double a, int N_nu_mass, double *mnu, double T_CMB, gsl_interp_accel *accel, int *status)
double *	ccl_nu_masses (double OmNuh2, ccl_neutrino_mass_splits mass_split, double T_CMB, int *status)

Variables
gsl_spline *	nu_spline = ((void*)0)

Function Documentation

gsl_spline* calculate_nu_phasespace_spline

(

int *

status

)

Spline for the phasespace integral required for getting the fractional energy density of massive neutrinos. Returns a gsl spline for the phase space integral needed for massive neutrinos.

Parameters

status

Status flag. 0 if there are no errors, nonzero otherwise. For specific cases see documentation for ccl_error.c

Returns

spl, the gsl spline for the phasespace integral required for massive neutrino calculations.

Definition at line 33 of file ccl_neutrinos.c.

References A_SPLINE_TYPE, ccl_cosmology_read_config(), CCL_ERROR_MISSING_CONFIG_FILE, CCL_ERROR_NU_INT, ccl_gsl, ccl_linear_spacing(), CCL_NU_MNUT_MAX, CCL_NU_MNUT_MIN, CCL_NU_MNUT_N, ccl_raise_exception(), ccl_raise_gsl_warning(), ccl_splines, ccl_gsl_params::INTEGRATION_NU_EPSABS, ccl_gsl_params::INTEGRATION_NU_EPSREL, ccl_gsl_params::N_ITERATION, and nu_integrand().

Referenced by nu_phasespace_intg().

                                                        {
  int N=CCL_NU_MNUT_N;
  double *mnut = ccl_linear_spacing(log(CCL_NU_MNUT_MIN),log(CCL_NU_MNUT_MAX),N);
  double *y=malloc(sizeof(double)*CCL_NU_MNUT_N);
  if (y ==NULL) {
    // Not setting a status_message here because we can't easily pass a cosmology to this function - message printed in ccl_error.c.
    *status = CCL_ERROR_NU_INT;
  }
  
  // Check whether ccl_splines and ccl_gsl exist; exit gracefully if they 
  // can't be loaded
  if(ccl_splines==NULL || ccl_gsl==NULL) ccl_cosmology_read_config();
  if(ccl_splines==NULL || ccl_gsl==NULL) {
    ccl_raise_exception(CCL_ERROR_MISSING_CONFIG_FILE, 
                        "ccl_neutrinos.c: Failed to read config file.");
    *status = CCL_ERROR_MISSING_CONFIG_FILE;
    return NULL;
  }
  
  int stat=0, gslstatus;
  gsl_integration_cquad_workspace * workspace = \
                gsl_integration_cquad_workspace_alloc(ccl_gsl->N_ITERATION);
  gsl_function F;
  F.function = &nu_integrand;
  for (int i=0; i<CCL_NU_MNUT_N; i++) {
    double mnut_=exp(mnut[i]);
    F.params = &(mnut_);
    gslstatus = gsl_integration_cquad(&F, 0, 1000.0, 
                                      ccl_gsl->INTEGRATION_NU_EPSABS, 
                                      ccl_gsl->INTEGRATION_NU_EPSREL, 
                                      workspace,&y[i], NULL, NULL);
    if(gslstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_neutrinos.c: calculate_nu_phasespace_spline():");
      stat |= gslstatus;
    }
  }
  gsl_integration_cquad_workspace_free(workspace);
  double renorm=1./y[0];
  for (int i=0; i<CCL_NU_MNUT_N; i++)
    y[i]*=renorm;
  gsl_spline* spl = gsl_spline_alloc(A_SPLINE_TYPE, CCL_NU_MNUT_N);
  stat |= gsl_spline_init(spl, mnut, y, CCL_NU_MNUT_N);
  
  // Check for errors in creating the spline
  if (stat) {
    // Not setting a status_message here because we can't easily pass a cosmology to this function - message printed in ccl_error.c.  
    *status = CCL_ERROR_NU_INT;  
    free(mnut);
    free(y);
    gsl_spline_free(spl);
  }
  free(mnut);
  free(y);
  return spl;
}

double* ccl_nu_masses	(	double	OmNuh2,
		ccl_neutrino_mass_splits	mass_split,
		double	T_CMB,
		int *	status
	)

Returns mass of one neutrino species at a scale factor a.

Parameters

a	Scale factor
Neff	The effective number of species with neutrino mass mnu.
OmNuh2	Fractional energy density of neutrions with mass mnu, multiplied by h squared. (can be 0).
T_CMB	Temperature of the CMB
accel	- Interpolation accelerator to be used with phasespace spline. If not set yet, pass NULL.
status	Status flag. 0 if there are no errors, nonzero otherwise. For specific cases see documentation for ccl_error.c

Returns

Mnu Neutrino mass [eV].

Definition at line 168 of file ccl_neutrinos.c.

References ccl_check_status_nocosmo(), CCL_ERROR_MNU_UNPHYSICAL, ccl_nu_equal, ccl_nu_inverted, ccl_nu_normal, ccl_nu_sum, DELTAM12_sq, DELTAM13_sq_neg, DELTAM13_sq_pos, ccl_test_distances::mnu, and pow().

Referenced by __attribute__().

                                                                                                     {
  
  double sumnu;
  
  sumnu = 93.14 * OmNuh2;
  
  // Now split the sum up into three masses depending on the label given:
  
  if(mass_split==ccl_nu_normal){
      
     double *mnu;
     mnu = malloc(3*sizeof(double));
            
     // See CCL note for how we get these expressions for the neutrino masses in normal and inverted hierarchy. 
     mnu[0] = 2./3.* sumnu - 1./6. * pow(-6. * DELTAM12_sq + 12. * DELTAM13_sq_pos + 4. * sumnu*sumnu, 0.5) 
             - 0.25 * DELTAM12_sq / (2./3.* sumnu - 1./6. * pow(-6. * DELTAM12_sq + 12. * DELTAM13_sq_pos + 4. * sumnu*sumnu, 0.5));
     mnu[1] = 2./3.* sumnu - 1./6. * pow(-6. * DELTAM12_sq + 12. * DELTAM13_sq_pos + 4. * sumnu*sumnu, 0.5) 
             + 0.25 * DELTAM12_sq / (2./3.* sumnu - 1./6. * pow(-6. * DELTAM12_sq + 12. * DELTAM13_sq_pos + 4. * sumnu*sumnu, 0.5));
     mnu[2] = -1./3. * sumnu + 1./3 * pow(-6. * DELTAM12_sq + 12. * DELTAM13_sq_pos + 4. * sumnu*sumnu, 0.5); 
     if (mnu[0]<0 || mnu[1]<0 || mnu[2]<0){
        // The user has provided a sum that is below the physical limit.
        if (sumnu<1e-14){
            mnu[0] = 0.;
            mnu[1] = 0.;
            mnu[2] = 0.;
        }else{
         *status = CCL_ERROR_MNU_UNPHYSICAL;
     }
     }   
     ccl_check_status_nocosmo(status);
     return mnu; 
     
  } else if (mass_split==ccl_nu_inverted){ 

    double *mnu;
    mnu = malloc(3*sizeof(double));
    mnu[0] = 2./3.* sumnu - 1./6. * pow(-6. * DELTAM12_sq + 12. * DELTAM13_sq_neg + 4. * sumnu * sumnu, 0.5) 
             - 0.25 * DELTAM12_sq / (2./3.* sumnu - 1./6. * pow(-6. * DELTAM12_sq + 12. * DELTAM13_sq_neg + 4. * sumnu * sumnu, 0.5));
    mnu[1] = 2./3.* sumnu - 1./6. * pow(-6. * DELTAM12_sq + 12. * DELTAM13_sq_neg + 4. * sumnu*sumnu, 0.5) 
             + 0.25 * DELTAM12_sq / (2./3.* sumnu - 1./6. * pow(-6. * DELTAM12_sq + 12. * DELTAM13_sq_neg + 4. * sumnu*sumnu, 0.5));
    mnu[2] = -1./3. * sumnu + 1./3 * pow(-6. * DELTAM12_sq + 12. * DELTAM13_sq_neg + 4. * sumnu*sumnu, 0.5); 
    if(mnu[0]<0 || mnu[1]<0 || mnu[2]<0){
    // The user has provided a sum that is below the physical limit.
        if (sumnu<1e-14){
            mnu[0] = 0.;
            mnu[1] = 0.;
            mnu[2] = 0.;;
        }else{
        *status = CCL_ERROR_MNU_UNPHYSICAL;
        }
    }    
    
    ccl_check_status_nocosmo(status);
    return mnu; 
    
  } else if (mass_split==ccl_nu_equal){
      double *mnu;
      mnu = malloc(3*sizeof(double));
      mnu[0] = sumnu/3.;
      mnu[1] = sumnu/3.;
      mnu[2] = sumnu/3.;
      
      return mnu;
      
  } else if (mass_split == ccl_nu_sum){
      
      double *mnu;
      mnu = malloc(sizeof(double));
      
      mnu[0] = sumnu;
      
      return mnu;
      
  } else{
      
      printf("WARNING:  mass option = %d not yet supported\n continuing with normal hierarchy\n", mass_split);
      double *mnu;
      mnu = malloc(3*sizeof(double));
            
     if (sumnu>0.59){
         mnu[0] = (sumnu - 0.59) / 3.;
         mnu[1] = mnu[0] + 0.009;
         mnu[2] = mnu[0] + 0.05;
     }else{
         *status = CCL_ERROR_MNU_UNPHYSICAL;
     }
     
     ccl_check_status_nocosmo(status);
     return mnu; 
  } 
}     

double ccl_Omeganuh2	(	double	a,
		int	N_nu_mass,
		double *	mnu,
		double	T_CMB,
		gsl_interp_accel *	accel,
		int *	status
	)

Returns density of one neutrino species at a scale factor a. Users are encouraged to access this quantity via the function ccl_omega_x.

Parameters

a	Scale factor
Neff	The effective number of species with neutrino mass mnu.
mnu	Pointer to array containing neutrino mass (can be 0).
T_CMB	Temperature of the CMB
accel	- Interpolation accelerator to be used with phasespace spline. If not set yet, pass NULL.
status	Status flag. 0 if there are no errors, nonzero otherwise. For specific cases see documentation for ccl_error.c

Returns

OmNuh2 Fractional energy density of neutrions with mass mnu, multiplied by h squared.

Definition at line 125 of file ccl_neutrinos.c.

References EV_IN_J, KBOLTZ, NU_CONST, nu_phasespace_intg(), pow(), and TNCDM.

Referenced by ccl_cosmology_compute_power_emu(), ccl_omega_x(), ccl_parameters_fill_initial(), and h_over_h0().

{
  double Tnu, a4, prefix_massless, mnuone, OmNuh2;
  double Tnu_eff, mnuOT, intval, prefix_massive;
  double total_mass; // To check if this is the massless or massive case.
  
  // First check if N_nu_mass is 0
  if (N_nu_mass==0) return 0.0;  
  
  Tnu=T_CMB*pow(4./11.,1./3.);
  a4=a*a*a*a;  
  // Check if mnu=0. We assume that in the massless case mnu is a pointer to a single element and that element is 0. This should in principle never be called.
  if (mnu[0] < 0.00017) {  // Limit taken from Lesgourges et al. 2012
    prefix_massless = NU_CONST  * Tnu * Tnu * Tnu * Tnu; 
    return N_nu_mass*prefix_massless*7./8./a4;
  }
  
  // And the remaining massive case.
  // Tnu_eff is used in the massive case because CLASS uses an effective temperature of nonLCDM components to match to mnu / Omeganu =93.14eV. Tnu_eff = T_ncdm * T_CMB = 0.71611 * T_CMB
  Tnu_eff = Tnu * TNCDM / (pow(4./11.,1./3.));
  
  // Define the prefix using the effective temperature (to get mnu / Omega = 93.14 eV) for the massive case: 
  prefix_massive = NU_CONST * Tnu_eff * Tnu_eff * Tnu_eff * Tnu_eff;
  
  OmNuh2 = 0.; // Initialize to 0 - we add to this for each massive neutrino species.
  for(int i=0; i<N_nu_mass; i++){
    // Get mass over T (mass (eV) / ((kb eV/s/K) Tnu_eff (K)) 
    // This returns the density normalized so that we get nuh2 at a=0
    mnuOT = mnu[i] / (Tnu_eff/a) * (EV_IN_J / (KBOLTZ)); 
  
    // Get the value of the phase-space integral 
    intval=nu_phasespace_intg(accel,mnuOT, status);
    OmNuh2 = intval*prefix_massive/a4 + OmNuh2;
  }
  
  return OmNuh2;
}

static double nu_integrand ( double  x, void *  r  )

static

Definition at line 23 of file ccl_neutrinos.c.

References sqrt(), and x.

Referenced by calculate_nu_phasespace_spline().

{
  double rat=*((double*)(r));
  return sqrt(x*x+rat*rat)/(exp(x)+1.0)*x*x; 
}

double nu_phasespace_intg	(	gsl_interp_accel *	accel,
		double	mnuOT,
		int *	status
	)

Definition at line 94 of file ccl_neutrinos.c.

References calculate_nu_phasespace_spline(), ccl_check_status_nocosmo(), CCL_NU_MNUT_MAX, CCL_NU_MNUT_MIN, ccl_raise_gsl_warning(), and nu_spline.

Referenced by ccl_Omeganuh2().

{
  // Check if the global variable for the phasespace spline has been defined yet:
  if (nu_spline==NULL) nu_spline =calculate_nu_phasespace_spline(status);
  ccl_check_status_nocosmo(status);
  
  double integral_value =0.;
  
  // First check the cases where we are in the limits.
  if (mnuOT<CCL_NU_MNUT_MIN) {
    return 7./8.;
  }
  else if (mnuOT>CCL_NU_MNUT_MAX) {
    return 0.2776566337*mnuOT; 
  }
    
  // Evaluate the spline - this will use the accelerator if it has been defined.
  int gslstatus = gsl_spline_eval_e(nu_spline, log(mnuOT),accel, &integral_value);
  if(gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_neutrinos.c: nu_phasespace_intg():");
    *status |= gslstatus;
  }
  return integral_value*7./8.;
}

Variable Documentation

gsl_spline* nu_spline = ((void*)0)

Definition at line 17 of file ccl_neutrinos.c.

Referenced by nu_phasespace_intg().