ccl_halomod.c File Reference

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

Include dependency graph for ccl_halomod.c:

Go to the source code of this file.

Classes
struct	Int_one_halo_Par
struct	Int_two_halo_Par

Functions
static double	u_nfw_c (ccl_cosmology *cosmo, double rv, double c, double k, int *status)
double	ccl_halo_concentration (ccl_cosmology *cosmo, double halomass, double a, double odelta, int *status)
static double	window_function (ccl_cosmology *cosmo, double m, double k, double a, double odelta, ccl_win_label label, int *status)
static double	one_halo_integrand (double log10mass, void *params)
static double	one_halo_integral (ccl_cosmology *cosmo, double k, double a, int *status)
static double	two_halo_integrand (double log10mass, void *params)
static double	two_halo_integral (ccl_cosmology *cosmo, double k, double a, int *status)
double	ccl_twohalo_matter_power (ccl_cosmology *cosmo, double k, double a, int *status)
double	ccl_onehalo_matter_power (ccl_cosmology *cosmo, double k, double a, int *status)
double	ccl_halomodel_matter_power (ccl_cosmology *cosmo, double k, double a, int *status)

Function Documentation

double ccl_halo_concentration	(	ccl_cosmology *	cosmo,
		double	halomass,
		double	a,
		double	odelta,
		int *	status
	)

Computes the concentration of a halo of mass M. This is the ratio of virial raidus to scale radius for an NFW halo.

Parameters

cosmo	cosmology object containing parameters
halomass	halo mass in units of Msun
a	scale factor normalised to a=1 today
odelta	overdensity criteria (with respect to matter density) used for halo mass
status	Status flag: 0 if there are no errors, non-zero otherwise

Returns

halo_concentration: the halo concentration

Definition at line 52 of file ccl_halomod.c.

References ccl_bhattacharya2011, ccl_constant_concentration, ccl_cosmology_set_status_message(), ccl_duffy2008, CCL_ERROR_CONC_DV, CCL_ERROR_HALOCONC, ccl_growth_factor(), ccl_raise_exception(), ccl_sigmaM(), ccl_cosmology::config, Dv_BryanNorman(), ccl_parameters::h, ccl_configuration::halo_concentration_method, ccl_cosmology::params, and pow().

Referenced by window_function().

                                                                                                          {

  double gz, g0, nu, delta_c, a_form;
  double Mpiv, A, B, C;

  switch(cosmo->config.halo_concentration_method){

    // Bhattacharya et al. (2011; 1005.2239; Delta = 200rho_m; Table 2)
  case ccl_bhattacharya2011:

    if (odelta != 200.) {
      *status = CCL_ERROR_CONC_DV;
      ccl_cosmology_set_status_message(cosmo, "ccl_halomod.c: halo_concentration(): Bhattacharya (2011) concentration relation only valid for Delta_v = 200 \n");
      return NAN;
    }
    
    gz = ccl_growth_factor(cosmo,a,status);
    g0 = ccl_growth_factor(cosmo,1.0,status);
    delta_c = 1.686;
    nu = delta_c/ccl_sigmaM(cosmo, halomass, a, status);
    return 9.*pow(nu,-0.29)*pow(gz/g0,1.15);

    // Duffy et al. (2008; 0804.2486; Table 1)
  case ccl_duffy2008:

    Mpiv = 2e12/cosmo->params.h; // Pivot mass in Msun (note in the paper units are Msun/h)

    if (odelta == Dv_BryanNorman(cosmo, a, status)) {

      // Duffy et al. (2008) for virial density haloes (second section in Table 1)
      A = 7.85;
      B = -0.081;
      C = -0.71;
      return A*pow(halomass/Mpiv,B)*pow(a,-C);

    } else if (odelta == 200.) {

      // Duffy et al. (2008) for x200 mean-matter-density haloes (third section in Table 1)
      A = 10.14;
      B = -0.081;
      C = -1.01;
      return A*pow(halomass/Mpiv,B)*pow(a,-C);

    } else {

      *status = CCL_ERROR_CONC_DV;
      ccl_cosmology_set_status_message(cosmo, "ccl_halomod.c: halo_concentration(): Duffy (2008) virial concentration only valid for virial Delta_v or 200\n");
      return NAN;
      
    }

    // Constant concentration (good for tests)
  case ccl_constant_concentration:
    
    return 4.;

    // Something went wrong
  default:
    
    *status = CCL_ERROR_HALOCONC;
    ccl_raise_exception(*status, "ccl_halomod.c: concentration-mass relation specified incorrectly");
    return NAN;
          
  }
}

double ccl_halomodel_matter_power	(	ccl_cosmology *	cosmo,
		double	k,
		double	a,
		int *	status
	)

Computes the halo model density-density power spectrum as the sum of two- and one-halo terms.

Parameters

cosmo	cosmology object containing parameters
k	wavenumber in units of Mpc^{-1}
a	scale factor normalised to a=1 today
status	Status flag: 0 if there are no errors, non-zero otherwise

Returns

halomodel_matter_power: halo-model power spectrum, P(k), units of Mpc^{3}

Definition at line 314 of file ccl_halomod.c.

References ccl_onehalo_matter_power(), and ccl_twohalo_matter_power().

Referenced by compare_halomod().

                                                                                        {

  // Standard sum of two- and one-halo terms
  return ccl_twohalo_matter_power(cosmo, k, a, status)+ccl_onehalo_matter_power(cosmo, k, a, status);
  
}

double ccl_onehalo_matter_power	(	ccl_cosmology *	cosmo,
		double	k,
		double	a,
		int *	status
	)

Computes the halo model density-density power spectrum one-halo term.

Parameters

cosmo	cosmology object containing parameters
k	wavenumber in units of Mpc^{-1}
a	scale factor normalised to a=1 today
status	Status flag: 0 if there are no errors, non-zero otherwise

Returns

1halo_matter_power: halo-model one-halo matter power spectrum, P(k), units of Mpc^{3}

Definition at line 304 of file ccl_halomod.c.

References one_halo_integral().

Referenced by ccl_halomodel_matter_power().

                                                                                      {
  
  return one_halo_integral(cosmo, k, a, status);
  
}

double ccl_twohalo_matter_power	(	ccl_cosmology *	cosmo,
		double	k,
		double	a,
		int *	status
	)

Computes the halo model density-density power spectrum two-halo term.

Parameters

cosmo	cosmology object containing parameters
k	wavenumber in units of Mpc^{-1}
a	scale factor normalised to a=1 today
status	Status flag: 0 if there are no errors, non-zero otherwise

Returns

2halo_matter_power: halo-model two-halo matter power spectrum, P(k), units of Mpc^{3}

Definition at line 277 of file ccl_halomod.c.

References ccl_linear_matter_power(), ccl_nfw, Dv_BryanNorman(), HM_MMIN, two_halo_integral(), and window_function().

Referenced by ccl_halomodel_matter_power().

                                                                                      {
    
  // Get the integral
  double I2h = two_halo_integral(cosmo, k, a, status);
    
  // The addative correction is the missing part of the integral below the lower-mass limit
  double A = 1.-two_halo_integral(cosmo, 0., a, status);

  // Virial overdensity for haloes
  double odelta = Dv_BryanNorman(cosmo, a, status);

  // ...multiplied by the ratio of window functions
  double W1 = window_function(cosmo, HM_MMIN, k,  a, odelta, ccl_nfw, status);
  double W2 = window_function(cosmo, HM_MMIN, 0., a, odelta, ccl_nfw, status);
  A = A*W1/W2;    

  // Add the additive correction to the calculated integral
  I2h = I2h+A;
      
  return ccl_linear_matter_power(cosmo, k, a, status)*I2h*I2h;
  
}

static double one_halo_integral ( ccl_cosmology *  cosmo, double  k, double  a, int *  status  )

static

Definition at line 174 of file ccl_halomod.c.

References Int_one_halo_Par::a, ccl_cosmology_set_status_message(), CCL_ERROR_ONE_HALO_INT, ccl_raise_gsl_warning(), cl_cmbl_bm::cosmo, Int_one_halo_Par::cosmo, HM_EPSABS, HM_EPSREL, HM_INT_METHOD, HM_LIMIT, HM_MMAX, HM_MMIN, Int_one_halo_Par::k, one_halo_integrand(), Int_one_halo_Par::status, and w.

Referenced by ccl_onehalo_matter_power().

                                                                                      {

  int one_halo_integral_status = 0, qagstatus;
  double result = 0, eresult;
  double log10mmin = log10(HM_MMIN);
  double log10mmax = log10(HM_MMAX);
  Int_one_halo_Par ipar;
  gsl_function F;
  gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);

  // Structure required for the gsl integration
  ipar.cosmo = cosmo;
  ipar.k = k;
  ipar.a = a;
  ipar.status = &one_halo_integral_status;
  F.function = &one_halo_integrand;
  F.params = &ipar;

  // Actually does the integration
  qagstatus = gsl_integration_qag(&F, log10mmin, log10mmax, HM_EPSABS, HM_EPSREL, HM_LIMIT, HM_INT_METHOD, w, &result, &eresult);

  // Clean up
  gsl_integration_workspace_free(w);

  // Check for errors
  if (qagstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(qagstatus, "ccl_halomod.c: one_halo_integral():");
    *status = CCL_ERROR_ONE_HALO_INT;
    ccl_cosmology_set_status_message(cosmo, "ccl_halomod.c: one_halo_integral(): Integration failure\n");
    return NAN;      
  } else {
    return result;
  }
  
}

static double one_halo_integrand ( double  log10mass, void *  params  )

static

Definition at line 158 of file ccl_halomod.c.

References Int_one_halo_Par::a, ccl_massfunc(), ccl_nfw, Int_one_halo_Par::cosmo, Dv_BryanNorman(), Int_one_halo_Par::k, p, pow(), Int_one_halo_Par::status, and window_function().

Referenced by one_halo_integral().

                                                                {  
  
  Int_one_halo_Par *p = (Int_one_halo_Par *)params;;
  double halomass = pow(10,log10mass);
  double odelta = Dv_BryanNorman(p->cosmo, p->a, p->status); // Virial density for haloes

  // The normalised Fourier Transform of a halo density profile
  double wk = window_function(p->cosmo,halomass, p->k, p->a, odelta, ccl_nfw, p->status);

  // Fairly sure that there should be no ln(10) factor should be here since the integration is being specified in log10 range
  double dn_dlogM = ccl_massfunc(p->cosmo, halomass, p->a, odelta, p->status);
    
  return dn_dlogM*pow(wk,2);
}

static double two_halo_integral ( ccl_cosmology *  cosmo, double  k, double  a, int *  status  )

static

Definition at line 237 of file ccl_halomod.c.

References Int_two_halo_Par::a, ccl_cosmology_set_status_message(), CCL_ERROR_TWO_HALO_INT, ccl_raise_gsl_warning(), cl_cmbl_bm::cosmo, Int_two_halo_Par::cosmo, HM_EPSABS, HM_EPSREL, HM_INT_METHOD, HM_LIMIT, HM_MMAX, HM_MMIN, Int_two_halo_Par::k, Int_two_halo_Par::status, two_halo_integrand(), and w.

Referenced by ccl_twohalo_matter_power().

                                                                                      {

  int two_halo_integral_status = 0, qagstatus;
  double result = 0, eresult;
  double log10mmin = log10(HM_MMIN);
  double log10mmax = log10(HM_MMAX);
  Int_two_halo_Par ipar;
  gsl_function F;
  gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);

  // Structure required for the gsl integration
  ipar.cosmo = cosmo;
  ipar.k = k;
  ipar.a = a;
  ipar.status = &two_halo_integral_status;
  F.function = &two_halo_integrand;
  F.params = &ipar;

  // Actually does the integration
  qagstatus = gsl_integration_qag(&F, log10mmin, log10mmax, HM_EPSABS, HM_EPSREL, HM_LIMIT, HM_INT_METHOD, w, &result, &eresult);

  // Clean up
  gsl_integration_workspace_free(w);

  // Check for errors
  if (qagstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(qagstatus, "ccl_halomod.c: two_halo_integral():");
    *status = CCL_ERROR_TWO_HALO_INT;
    ccl_cosmology_set_status_message(cosmo, "ccl_halomod.c: two_halo_integral(): Integration failure\n");
    return NAN;      
  } else {
    return result;
  }
  
}

static double two_halo_integrand ( double  log10mass, void *  params  )

static

Definition at line 218 of file ccl_halomod.c.

References Int_two_halo_Par::a, ccl_halo_bias(), ccl_massfunc(), ccl_nfw, Int_two_halo_Par::cosmo, Dv_BryanNorman(), Int_two_halo_Par::k, p, pow(), Int_two_halo_Par::status, and window_function().

Referenced by two_halo_integral().

                                                                {  
  
  Int_two_halo_Par *p = (Int_two_halo_Par *)params;
  double halomass = pow(10,log10mass);
  double odelta = Dv_BryanNorman(p->cosmo, p->a, p->status); // Virial density for haloes

  // The normalised Fourier Transform of a halo density profile
  double wk = window_function(p->cosmo,halomass, p->k, p->a, odelta, ccl_nfw, p->status);

  // Fairly sure that there should be no ln(10) factor should be here since the integration is being specified in log10 range
  double dn_dlogM = ccl_massfunc(p->cosmo, halomass, p->a, odelta, p->status);

  // Halo bias
  double b = ccl_halo_bias(p->cosmo, halomass, p->a, odelta, p->status);
    
  return b*dn_dlogM*wk;
}

static double u_nfw_c ( ccl_cosmology *  cosmo, double  rv, double  c, double  k, int *  status  )

static

Definition at line 17 of file ccl_halomod.c.

References cl_cmbl_bm::c, and rs.

Referenced by window_function().

                                                                                       {
   
  double rs, ks;
  double f1, f2, f3, fc;

  // Special case to prevent numerical problems if k=0,
  // the result should be unity here because of the normalisation
  if (k==0.) {    
    return 1.;    
  }

  // The general k case
  else{
    
    // Scale radius for NFW (rs=rv/c)
    rs = rv/c;

    // Dimensionless wave-number variable
    ks = k*rs;

    // Various bits for summing together to get final result
    f1 = sin(ks)*(gsl_sf_Si(ks*(1.+c))-gsl_sf_Si(ks));
    f2 = cos(ks)*(gsl_sf_Ci(ks*(1.+c))-gsl_sf_Ci(ks));
    f3 = sin(c*ks)/(ks*(1.+c));
    fc = log(1.+c)-c/(1.+c);
  
    return (f1+f2-f3)/fc;
    
  }
}

static double window_function ( ccl_cosmology *  cosmo, double  m, double  k, double  a, double  odelta, ccl_win_label  label, int *  status  )

static

Definition at line 119 of file ccl_halomod.c.

References cl_cmbl_bm::c, CCL_ERROR_HALOWIN, ccl_halo_concentration(), ccl_nfw, ccl_raise_exception(), ccl_rho_x(), ccl_species_m_label, r_delta(), and u_nfw_c().

Referenced by ccl_twohalo_matter_power(), one_halo_integrand(), and two_halo_integrand().

                                                                                                                                  {

  double rho_matter, c, rv;
  
  switch(label){

  case ccl_nfw:

    // The mean background matter density in Msun/Mpc^3
    rho_matter = ccl_rho_x(cosmo, 1., ccl_species_m_label, 1, status);

    // The halo virial radius
    rv = r_delta(cosmo, m, a, odelta, status);

    // The halo concentration for this halo mass and at this scale factor  
    c = ccl_halo_concentration(cosmo, m, a, odelta, status);    

    // The function U is normalised to 1 for k<<1 so multiplying by M/rho turns units to overdensity
    return m*u_nfw_c(cosmo, rv, c, k, status)/rho_matter;

    // Something went wrong
  default:
    
    *status = CCL_ERROR_HALOWIN;
    ccl_raise_exception(*status, "ccl_halomod.c: Window function specified incorrectly");
    return NAN;
    
  }
  
}