ccl_massfunc.c File Reference

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

Include dependency graph for ccl_massfunc.c:

Go to the source code of this file.

Functions
double	dc_NakamuraSuto (ccl_cosmology *cosmo, double a, int *status)
double	Dv_BryanNorman (ccl_cosmology *cosmo, double a, int *status)
double	r_delta (ccl_cosmology *cosmo, double halomass, double a, double odelta, int *status)
static void	ccl_cosmology_compute_hmfparams (ccl_cosmology *cosmo, int *status)
static double	massfunc_f (ccl_cosmology *cosmo, double halomass, double a, double odelta, int *status)
static double	ccl_halo_b1 (ccl_cosmology *cosmo, double halomass, double a, double odelta, int *status)
void	ccl_cosmology_compute_sigma (ccl_cosmology *cosmo, int *status)
static double	ccl_dlninvsig_dlogm (ccl_cosmology *cosmo, double halomass, int *status)
double	ccl_massfunc (ccl_cosmology *cosmo, double halomass, double a, double odelta, int *status)
double	ccl_halo_bias (ccl_cosmology *cosmo, double halomass, double a, double odelta, int *status)
double	ccl_massfunc_m2r (ccl_cosmology *cosmo, double halomass, int *status)
double	ccl_sigmaM (ccl_cosmology *cosmo, double halomass, double a, int *status)

Function Documentation

static void ccl_cosmology_compute_hmfparams ( ccl_cosmology *  cosmo, int *  status  )

static

Definition at line 62 of file ccl_massfunc.c.

References ccl_data::accelerator_d, ccl_data::alphahmf, beta, ccl_data::betahmf, ccl_cosmology_set_status_message(), CCL_ERROR_SPLINE, ccl_tinker, ccl_tinker10, ccl_cosmology::computed_hmfparams, ccl_cosmology::config, D_SPLINE_TYPE, ccl_cosmology::data, ccl_data::etahmf, ccl_data::gammahmf, ccl_configuration::mass_function_method, and ccl_data::phihmf.

Referenced by massfunc_f().

{
  if(cosmo->computed_hmfparams)
    return;

  // declare parameter splines on case-by-case basis
  switch(cosmo->config.mass_function_method) {
  case ccl_tinker10:{
    double delta[9] = {200.0, 300.0, 400.0, 600.0, 800.0, 1200.0, 1600.0, 2400.0, 3200.0};
    double lgdelta[9];
    double alpha[9] = {0.368, 0.363, 0.385, 0.389, 0.393, 0.365, 0.379, 0.355, 0.327};
    double beta[9] = {0.589, 0.585, 0.544, 0.543, 0.564, 0.623, 0.637, 0.673, 0.702};
    double gamma[9] ={0.864, 0.922, 0.987, 1.09, 1.20, 1.34, 1.50, 1.68, 1.81};
    double phi[9] = {-0.729, -0.789, -0.910, -1.05, -1.20, -1.26, -1.45, -1.50, -1.49};
    double eta[9] = {-0.243, -0.261, -0.261, -0.273, -0.278, -0.301, -0.301, -0.319, -0.336};
    int nd = 9;
    int i;

    for(i=0; i<nd; i++) {
      lgdelta[i] = log10(delta[i]);
    }

    gsl_spline * alphahmf = gsl_spline_alloc(D_SPLINE_TYPE, nd);
    *status = gsl_spline_init(alphahmf, lgdelta, alpha, nd);
    if (*status) {
      gsl_spline_free(alphahmf);
      *status = CCL_ERROR_SPLINE ;
      ccl_cosmology_set_status_message(cosmo, "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error creating alpha(D) spline\n");
      return;
    }

    gsl_spline * betahmf  = gsl_spline_alloc(D_SPLINE_TYPE, nd);
    *status = gsl_spline_init(betahmf, lgdelta, beta, nd);
    if (*status) {
      gsl_spline_free(alphahmf);
      gsl_spline_free(betahmf);
      *status = CCL_ERROR_SPLINE ;
      ccl_cosmology_set_status_message(cosmo, "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error creating beta(D) spline\n");
      return;
    }

    gsl_spline * gammahmf = gsl_spline_alloc(D_SPLINE_TYPE, nd);
    *status = gsl_spline_init(gammahmf, lgdelta, gamma, nd);
    if (*status) {
      gsl_spline_free(alphahmf);
      gsl_spline_free(betahmf);
      gsl_spline_free(gammahmf);
      *status = CCL_ERROR_SPLINE ;
      ccl_cosmology_set_status_message(cosmo, "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error creating gamma(D) spline\n");
      return;
    }

    gsl_spline * phihmf   = gsl_spline_alloc(D_SPLINE_TYPE, nd);
    *status = gsl_spline_init(phihmf, lgdelta, phi, nd);
    if (*status) {
      gsl_spline_free(alphahmf);
      gsl_spline_free(betahmf);
      gsl_spline_free(gammahmf);
      gsl_spline_free(phihmf);
      *status = CCL_ERROR_SPLINE ;
      ccl_cosmology_set_status_message(cosmo, "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error creating phi(D) spline\n");
      return;
    }

    gsl_spline * etahmf   = gsl_spline_alloc(D_SPLINE_TYPE, nd);
    *status = gsl_spline_init(etahmf, lgdelta, eta, nd);
    if (*status) {
      gsl_spline_free(alphahmf);
      gsl_spline_free(betahmf);
      gsl_spline_free(gammahmf);
      gsl_spline_free(phihmf);
      gsl_spline_free(etahmf);
      *status = CCL_ERROR_SPLINE ;
      ccl_cosmology_set_status_message(cosmo, "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error creating eta(D) spline\n");
      return;
    }
    if(cosmo->data.accelerator_d==NULL)
      cosmo->data.accelerator_d=gsl_interp_accel_alloc();
    cosmo->data.alphahmf = alphahmf;
    cosmo->data.betahmf = betahmf;
    cosmo->data.gammahmf = gammahmf;
    cosmo->data.phihmf = phihmf;
    cosmo->data.etahmf = etahmf;
    cosmo->computed_hmfparams = true;

    break;
  }
  case ccl_tinker:{
    double delta[9] = {200.0, 300.0, 400.0, 600.0, 800.0, 1200.0, 1600.0, 2400.0, 3200.0};
    double lgdelta[9];
    double alpha[9] = {0.186, 0.200, 0.212, 0.218, 0.248, 0.255, 0.260, 0.260, 0.260};
    double beta[9] = {1.47, 1.52, 1.56, 1.61, 1.87, 2.13, 2.30, 2.53, 2.66};
    double gamma[9] ={2.57, 2.25, 2.05, 1.87, 1.59, 1.51, 1.46, 1.44, 1.41};
    double phi[9] = {1.19, 1.27, 1.34, 1.45, 1.58, 1.80, 1.97, 2.24, 2.44};
    int nd = 9;
    int i;

    for(i=0; i<nd; i++) {
      lgdelta[i] = log10(delta[i]);
    }

    gsl_spline * alphahmf = gsl_spline_alloc(D_SPLINE_TYPE, nd);
    *status = gsl_spline_init(alphahmf, lgdelta, alpha, nd);
    if (*status) {
      gsl_spline_free(alphahmf);
      *status = CCL_ERROR_SPLINE ;
      ccl_cosmology_set_status_message(cosmo, "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error creating alpha(D) spline\n");
      return;
    }

    gsl_spline * betahmf  = gsl_spline_alloc(D_SPLINE_TYPE, nd);
    *status = gsl_spline_init(betahmf, lgdelta, beta, nd);
    if (*status) {
      gsl_spline_free(alphahmf);
      gsl_spline_free(betahmf);
      *status = CCL_ERROR_SPLINE ;
      ccl_cosmology_set_status_message(cosmo, "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error creating beta(D) spline\n");
      return;
    }

    gsl_spline * gammahmf = gsl_spline_alloc(D_SPLINE_TYPE, nd);
    *status = gsl_spline_init(gammahmf, lgdelta, gamma, nd);
    if (*status) {
      gsl_spline_free(alphahmf);
      gsl_spline_free(betahmf);
      gsl_spline_free(gammahmf);
      *status = CCL_ERROR_SPLINE ;
      ccl_cosmology_set_status_message(cosmo, "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error creating gamma(D) spline\n");
      return;
    }

    gsl_spline * phihmf   = gsl_spline_alloc(D_SPLINE_TYPE, nd);
    *status = gsl_spline_init(phihmf, lgdelta, phi, nd);
    if (*status) {
      gsl_spline_free(alphahmf);
      gsl_spline_free(betahmf);
      gsl_spline_free(gammahmf);
      gsl_spline_free(phihmf);
      *status = CCL_ERROR_SPLINE ;
      ccl_cosmology_set_status_message(cosmo, "ccl_massfunc.c: ccl_cosmology_compute_hmfparams(): Error creating phi(D) spline\n");
      return;
    }

    if(cosmo->data.accelerator_d==NULL)
      cosmo->data.accelerator_d=gsl_interp_accel_alloc();
    cosmo->data.alphahmf = alphahmf;
    cosmo->data.betahmf = betahmf;
    cosmo->data.gammahmf = gammahmf;
    cosmo->data.phihmf = phihmf;
    cosmo->computed_hmfparams = true;

    break;
  }
  default:
    // Error message could go here if we decide to make this public facing. 
    // Currently not accessible from the API though.
    break;
  }
}

void ccl_cosmology_compute_sigma	(	ccl_cosmology *	cosmo,
		int *	status
	)

Definition at line 431 of file ccl_massfunc.c.

References ccl_data::accelerator_m, ccl_cosmology_set_status_message(), CCL_ERROR_LINSPACE, CCL_ERROR_SPLINE, ccl_linear_spacing(), ccl_massfunc_m2r(), ccl_raise_gsl_warning(), ccl_sigmaR(), ccl_splines, ccl_cosmology::computed_sigma, ccl_cosmology::data, ccl_data::dlnsigma_dlogm, ccl_spline_params::LOGM_SPLINE_DELTA, ccl_spline_params::LOGM_SPLINE_MAX, ccl_spline_params::LOGM_SPLINE_MIN, ccl_spline_params::LOGM_SPLINE_NM, ccl_data::logsigma, m, M_SPLINE_TYPE, and pow().

Referenced by ccl_dlninvsig_dlogm(), ccl_halo_bias(), and ccl_sigmaM().

{
  if(cosmo->computed_sigma)
    return;

  // create linearly-spaced values of the mass.
  int nm=ccl_splines->LOGM_SPLINE_NM;
  double * m = ccl_linear_spacing(ccl_splines->LOGM_SPLINE_MIN, ccl_splines->LOGM_SPLINE_MAX, nm);

  // create space for y, to be filled with sigma and dlnsigma_dlogm
  double * y = malloc(sizeof(double)*nm);
  double smooth_radius;
  double na, nb;

  // start up of GSL pointers
  int gslstatus = 0;
  gsl_spline *logsigma;
  gsl_spline *dlnsigma_dlogm;

  if (m==NULL ||
      (fabs(m[0]-ccl_splines->LOGM_SPLINE_MIN)>1e-5) ||
      (fabs(m[nm-1]-ccl_splines->LOGM_SPLINE_MAX)>1e-5) ||
      (m[nm-1]>10E17)
      ) {
    *status = CCL_ERROR_LINSPACE;
    ccl_cosmology_set_status_message(cosmo,"ccl_cosmology_compute_sigmas(): Error creating linear spacing in m\n");
  }

  // fill in sigma, if no errors have been triggered at this time.
  if (*status == 0) {
    for (int i=0; i<nm; i++) {
      smooth_radius = ccl_massfunc_m2r(cosmo, pow(10,m[i]), status);
      y[i] = log10(ccl_sigmaR(cosmo, smooth_radius, 1., status));
    }
    logsigma = gsl_spline_alloc(M_SPLINE_TYPE, nm);
    *status = gsl_spline_init(logsigma, m, y, nm);
  }

  if (*status !=0 ) {
    *status = CCL_ERROR_SPLINE ;
    ccl_cosmology_set_status_message(cosmo, "ccl_massfunc.c: ccl_cosmology_compute_sigma(): Error creating sigma(M) spline\n");
  }

  // again, making splines assuming nothing bad has happened to this point

  if (*status == 0 ) {
    for (int i=0; i<nm; i++) {
      if(i==0) {
        gslstatus |= gsl_spline_eval_e(logsigma, m[i], NULL,&na);
        gslstatus |= gsl_spline_eval_e(logsigma, m[i]+ccl_splines->LOGM_SPLINE_DELTA/2., NULL,&nb);
        y[i] = (na-nb)*log(10.);
        y[i] = 2.*y[i] / ccl_splines->LOGM_SPLINE_DELTA;
      }
      else if (i==nm-1) {
        gslstatus |= gsl_spline_eval_e(logsigma, m[i]-ccl_splines->LOGM_SPLINE_DELTA/2., NULL,&na);
        gslstatus |= gsl_spline_eval_e(logsigma, m[i], NULL,&nb);
        y[i] = (na-nb)*log(10.);
        y[i] = 2.*y[i] / ccl_splines->LOGM_SPLINE_DELTA;
      }
      else {
        gslstatus |= gsl_spline_eval_e(logsigma, m[i]-ccl_splines->LOGM_SPLINE_DELTA/2., NULL,&na);
        gslstatus |= gsl_spline_eval_e(logsigma, m[i]+ccl_splines->LOGM_SPLINE_DELTA/2., NULL,&nb);
        y[i] = (na-nb)*log(10.);
        y[i] = y[i] / ccl_splines->LOGM_SPLINE_DELTA;
      }
    }
  }

  if(gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_massfunc.c: ccl_cosmology_compute_sigma():");
    *status |= gslstatus;
    *status = CCL_ERROR_SPLINE ;
    ccl_cosmology_set_status_message(cosmo, "ccl_massfunc.c: ccl_cosmology_compute_sigma(): Error evaluating grid points for dlnsigma/dlogM spline\n");
  }

  if(*status==0) {
    dlnsigma_dlogm = gsl_spline_alloc(M_SPLINE_TYPE, nm);
    *status = gsl_spline_init(dlnsigma_dlogm, m, y, nm);
    if(cosmo->data.accelerator_m==NULL)
      cosmo->data.accelerator_m=gsl_interp_accel_alloc();
  }

  if(*status!=0) {
    *status = CCL_ERROR_SPLINE ;
    ccl_cosmology_set_status_message(cosmo, "ccl_massfunc.c: ccl_cosmology_compute_sigma(): Error creating dlnsigma/dlogM spline\n");
  }
  

  cosmo->data.logsigma = logsigma;
  cosmo->data.dlnsigma_dlogm = dlnsigma_dlogm;
  cosmo->computed_sigma = true;

  free(m);
  free(y);
  if(*status != 0) {
    gsl_spline_free(logsigma);
    gsl_spline_free(dlnsigma_dlogm);
  }
  return;
}

static double ccl_dlninvsig_dlogm ( ccl_cosmology *  cosmo, double  halomass, int *  status  )

static

Definition at line 537 of file ccl_massfunc.c.

References ccl_data::accelerator_m, ccl_check_status(), ccl_cosmology_compute_sigma(), ccl_raise_gsl_warning(), ccl_cosmology::computed_sigma, ccl_cosmology::data, and ccl_data::dlnsigma_dlogm.

Referenced by ccl_massfunc().

{
  if (!cosmo->computed_sigma) {
    ccl_cosmology_compute_sigma(cosmo, status);
    ccl_check_status(cosmo, status);
  }

  double val, logmass;

  logmass = log10(halomass);

  int gslstatus = gsl_spline_eval_e(cosmo->data.dlnsigma_dlogm, logmass, cosmo->data.accelerator_m,&val);
  if(gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_massfunc.c: ccl_massfunc():");
    *status |= gslstatus;
  }
  ccl_check_status(cosmo, status);

  return val;
}

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

static

Definition at line 374 of file ccl_massfunc.c.

References ccl_cosmology_set_status_message(), CCL_ERROR_HMF_DV, CCL_ERROR_MF, ccl_shethtormen, ccl_sigmaM(), ccl_tinker10, ccl_cosmology::config, dc_NakamuraSuto(), Dv_BryanNorman(), ccl_configuration::mass_function_method, and pow().

Referenced by ccl_halo_bias().

{
  double fit_A, fit_B, fit_C, fit_a, fit_b, fit_c, fit_p, overdensity_delta, y;
  double delta_c_Tinker, nu;
  double sigma=ccl_sigmaM(cosmo,halomass,a, status);
  switch(cosmo->config.mass_function_method) {

  // Equation (12) in  arXiv: 9901122
  // Derived using the peak-background split applied to the mass function in the same paper
  // Note that Sheth & Tormen (1999) use nu=(dc/sigma)^2 whereas we use nu=dc/sigma
  case ccl_shethtormen:

    // Check if Delta_v is the virial Delta_v
    if (odelta != Dv_BryanNorman(cosmo, a, status)) {
      *status = CCL_ERROR_HMF_DV;
      ccl_cosmology_set_status_message(cosmo, "ccl_massfunc.c: halo_b1(): Sheth-Tormen called with not virial Delta_v\n");
      return NAN;
    }

    // ST bias fitting parameters (which are the same as for the mass function)
    fit_p = 0.3;
    fit_a = 0.707;

    // Cosmology dependent delta_c and nu
    double delta_c = dc_NakamuraSuto(cosmo, a, status);
    nu = delta_c/ccl_sigmaM(cosmo, halomass, a, status);

    return 1.+(fit_a*pow(nu,2)-1.+2.*fit_p/(1.+pow(fit_a*pow(nu,2),fit_p)))/delta_c;

    //this version uses b(nu) parameterization, Eq. 6 in Tinker et al. 2010
    // use this for consistency with Tinker et al. 2010 fitting function for halo bias
  case ccl_tinker10:
    y = log10(odelta);
    //critical collapse overdensity assumed in this model
    delta_c_Tinker = 1.686;
    //peak height - note that this factorization is incorrect for e.g. massive neutrino cosmologies
    nu = delta_c_Tinker/(sigma);
    // Table 2 in https://arxiv.org/pdf/1001.3162.pdf
    fit_A = 1.0 + 0.24*y*exp(-pow(4./y,4.));
    fit_a = 0.44*y-0.88;
    fit_B = 0.183;
    fit_b = 1.5;
    fit_C = 0.019+0.107*y+0.19*exp(-pow(4./y,4.));
    fit_c = 2.4;

    return 1.-fit_A*pow(nu,fit_a)/(pow(nu,fit_a)+pow(delta_c_Tinker,fit_a))+fit_B*pow(nu,fit_b)+fit_C*pow(nu,fit_c);
    break;

  default:
    *status = CCL_ERROR_MF;
    ccl_cosmology_set_status_message(cosmo ,
        "ccl_massfunc.c: ccl_halo_b1(): No b(M) fitting function implemented for mass_function_method: %d \n",
      cosmo->config.mass_function_method);
    return 0;
  }
}

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

Definition at line 582 of file ccl_massfunc.c.

References ccl_check_status(), ccl_cosmology_compute_sigma(), ccl_cosmology_set_status_message(), CCL_ERROR_NOT_IMPLEMENTED, ccl_halo_b1(), ccl_cosmology::computed_sigma, ccl_parameters::N_nu_mass, and ccl_cosmology::params.

Referenced by main(), and two_halo_integrand().

{
  if (cosmo->params.N_nu_mass>0){
      *status = CCL_ERROR_NOT_IMPLEMENTED;
      ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): Support for the halo bias in cosmologies with massive neutrinos is not yet implemented.\n");
      return NAN;
  }


  if (!cosmo->computed_sigma) {
    ccl_cosmology_compute_sigma(cosmo, status);
    ccl_check_status(cosmo, status);
  }

  double f;
  f = ccl_halo_b1(cosmo,halomass,a,odelta, status);
  ccl_check_status(cosmo, status);
  return f;
}

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

Definition at line 562 of file ccl_massfunc.c.

References ccl_cosmology_set_status_message(), ccl_dlninvsig_dlogm(), CCL_ERROR_NOT_IMPLEMENTED, ccl_parameters::h, massfunc_f(), ccl_parameters::N_nu_mass, ccl_parameters::Omega_m, ccl_cosmology::params, and RHO_CRITICAL.

Referenced by compare_massfunc(), main(), one_halo_integrand(), and two_halo_integrand().

{
  if (cosmo->params.N_nu_mass>0){
      *status = CCL_ERROR_NOT_IMPLEMENTED;
      ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): Support for the halo mass function in cosmologies with massive neutrinos is not yet implemented.\n");
      return NAN;
  }

  double f, rho_m;

  rho_m = RHO_CRITICAL*cosmo->params.Omega_m*cosmo->params.h*cosmo->params.h;
  f=massfunc_f(cosmo,halomass,a,odelta,status);

  return f*rho_m*ccl_dlninvsig_dlogm(cosmo,halomass,status)/halomass;
}

double ccl_massfunc_m2r	(	ccl_cosmology *	cosmo,
		double	halomass,
		int *	status
	)

Definition at line 606 of file ccl_massfunc.c.

References ccl_rho_x(), ccl_species_m_label, M_PI, and pow().

Referenced by ccl_cosmology_compute_sigma().

{
  double rho_m, smooth_radius;

  // Comoving matter density
  //rho_m = RHO_CRITICAL*cosmo->params.Omega_m*cosmo->params.h*cosmo->params.h;
  rho_m = ccl_rho_x(cosmo, 1., ccl_species_m_label, 1, status);

  smooth_radius = pow((3.0*halomass) / (4*M_PI*rho_m), (1.0/3.0));

  return smooth_radius;
}

double ccl_sigmaM	(	ccl_cosmology *	cosmo,
		double	halomass,
		double	a,
		int *	status
	)

Definition at line 624 of file ccl_massfunc.c.

References ccl_data::accelerator_m, ccl_check_status(), ccl_cosmology_compute_sigma(), ccl_cosmology_set_status_message(), CCL_ERROR_NOT_IMPLEMENTED, ccl_growth_factor(), ccl_raise_gsl_warning(), ccl_cosmology::computed_sigma, ccl_cosmology::data, ccl_data::logsigma, ccl_parameters::N_nu_mass, ccl_cosmology::params, and pow().

Referenced by ccl_halo_b1(), ccl_halo_concentration(), compare_massfunc(), compare_sigmam(), and massfunc_f().

{
  if (cosmo->params.N_nu_mass>0){
      *status = CCL_ERROR_NOT_IMPLEMENTED;
      ccl_cosmology_set_status_message(cosmo, "ccl_background.c: ccl_cosmology_compute_growth(): Support for the sigma(M) function in cosmologies with massive neutrinos is not yet implemented.\n");
      return NAN;
  }

  double sigmaM;
  // Check if sigma has already been calculated
  if (!cosmo->computed_sigma) {
    ccl_cosmology_compute_sigma(cosmo, status);
    ccl_check_status(cosmo, status);
  }

  double lgsigmaM;

  int gslstatus = gsl_spline_eval_e(cosmo->data.logsigma, 
                                    log10(halomass), 
                                    cosmo->data.accelerator_m,&lgsigmaM);

  if(gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_massfunc.c: ccl_sigmaM():");
    *status |= gslstatus;
  }

  // Interpolate to get sigma
  sigmaM = pow(10,lgsigmaM)*ccl_growth_factor(cosmo, a, status);
  ccl_check_status(cosmo, status);
  return sigmaM;
}

double dc_NakamuraSuto	(	ccl_cosmology *	cosmo,
		double	a,
		int *	status
	)

Definition at line 21 of file ccl_massfunc.c.

References ccl_omega_x(), ccl_species_m_label, M_PI, and pow().

Referenced by ccl_halo_b1(), and massfunc_f().

                                                                   {

  double Om_mz = ccl_omega_x(cosmo, a, ccl_species_m_label, status);
  double dc0 = (3./20.)*pow(12.*M_PI,2./3.);
  double dc = dc0*(1.+0.012299*log10(Om_mz));

  return dc;

}

double Dv_BryanNorman	(	ccl_cosmology *	cosmo,
		double	a,
		int *	status
	)

Definition at line 37 of file ccl_massfunc.c.

References ccl_omega_x(), ccl_species_m_label, M_PI, pow(), and x.

Referenced by ccl_halo_b1(), ccl_halo_concentration(), ccl_twohalo_matter_power(), massfunc_f(), one_halo_integrand(), and two_halo_integrand().

                                                                  {

  double Om_mz = ccl_omega_x(cosmo, a, ccl_species_m_label, status);
  double x = Om_mz-1.;
  double Dv0 = 18.*pow(M_PI,2);
  double Dv = (Dv0+82.*x-39.*pow(x,2))/Om_mz;

  return Dv;

}

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

static

Definition at line 234 of file ccl_massfunc.c.

References ccl_data::accelerator_d, ccl_data::alphahmf, ccl_data::betahmf, ccl_angulo, ccl_check_status(), ccl_cosmology_compute_hmfparams(), ccl_cosmology_set_status_message(), CCL_ERROR_HMF_DV, CCL_ERROR_HMF_INTERP, CCL_ERROR_MF, ccl_omega_x(), ccl_raise_gsl_warning(), ccl_shethtormen, ccl_sigmaM(), ccl_species_m_label, ccl_tinker, ccl_tinker10, ccl_watson, ccl_cosmology::computed_hmfparams, ccl_cosmology::config, ccl_cosmology::data, dc_NakamuraSuto(), Dv_BryanNorman(), ccl_data::etahmf, ccl_data::gammahmf, ccl_configuration::mass_function_method, ccl_data::phihmf, and pow().

Referenced by ccl_massfunc().

{
  double fit_A, fit_a, fit_b, fit_c, fit_d, fit_p, overdensity_delta;
  double Omega_m_a;
  double delta_c_Tinker, nu;

  double sigma=ccl_sigmaM(cosmo, halomass, a, status);
  int gslstatus;

  switch(cosmo->config.mass_function_method) {

  // Equation (10) in arxiv: 9901122
  // Note that Sheth & Tormen (1999) use nu=(dc/sigma)^2 whereas we use nu=dc/sigma
  case ccl_shethtormen:

    // Check if odelta is outside the interpolated range
    if (odelta != Dv_BryanNorman(cosmo, a, status)) {
      *status = CCL_ERROR_HMF_DV;
      ccl_cosmology_set_status_message(cosmo, "ccl_massfunc.c: massfunc_f(): Sheth-Tormen called with not virial Delta_v\n");
      return NAN;
    }

    // ST mass function fitting parameters
    fit_A = 0.21616;
    fit_p = 0.3;
    fit_a = 0.707;

    // nu = delta_c(z) / sigma(M)
    nu = dc_NakamuraSuto(cosmo, a, status)/ccl_sigmaM(cosmo, halomass, a, status);

    return nu*fit_A*(1.+pow(fit_a*pow(nu,2),-fit_p))*exp(-fit_a*pow(nu,2)/2.);

  case ccl_tinker:

    // Check if odelta is outside the interpolated range
    if ((odelta < 200) || (odelta > 3200)) {
      *status = CCL_ERROR_HMF_INTERP;
      ccl_cosmology_set_status_message(cosmo, "ccl_massfunc.c: massfunc_f(): Tinker 2008 only supported in range of Delta = 200 to Delta = 3200.\n");
      return NAN;
    }

    // Compute HMF parameter (alpha, beta, gamma, phi) splines if they haven't
    // been computed already
    if (!cosmo->computed_hmfparams) {
      ccl_cosmology_compute_hmfparams(cosmo, status);
      ccl_check_status(cosmo, status);
    }
    gslstatus = gsl_spline_eval_e(cosmo->data.alphahmf, log10(odelta), cosmo->data.accelerator_d,&fit_A);
    gslstatus |= gsl_spline_eval_e(cosmo->data.betahmf, log10(odelta), cosmo->data.accelerator_d,&fit_a);
    gslstatus |= gsl_spline_eval_e(cosmo->data.gammahmf, log10(odelta), cosmo->data.accelerator_d,&fit_b);
    gslstatus |= gsl_spline_eval_e(cosmo->data.phihmf, log10(odelta), cosmo->data.accelerator_d,&fit_c);
    fit_d = pow(10, -1.0*pow(0.75 / log10(odelta / 75.0), 1.2));

    fit_A = fit_A*pow(a, 0.14);
    fit_a = fit_a*pow(a, 0.06);
    fit_b = fit_b*pow(a, fit_d);
    if(gslstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_massfunc.c: ccl_massfunc_f():");
      *status |= gslstatus;
      ccl_cosmology_set_status_message(cosmo, "ccl_massfunc.c: ccl_massfunc_f(): interpolation error for Tinker MF\n");
      return NAN;
    }
    return fit_A*(pow(sigma/fit_b,-fit_a)+1.0)*exp(-fit_c/sigma/sigma);
    break;
  case ccl_tinker10:
    // this version uses f(nu) parameterization from Eq. 8 in Tinker et al. 2010
    // use this for consistency with Tinker et al. 2010 fitting function for halo bias

    // Check if odelta is outside the interpolated range
    if ((odelta < 200) || (odelta > 3200)) {
      *status = CCL_ERROR_HMF_INTERP;
      ccl_cosmology_set_status_message(cosmo, "ccl_massfunc.c: massfunc_f(): Tinker 2010 only supported in range of Delta = 200 to Delta = 3200.\n");
      return 0;
    }

    if (!cosmo->computed_hmfparams) {
        ccl_cosmology_compute_hmfparams(cosmo, status);
        ccl_check_status(cosmo, status);
    }
    //critical collapse overdensity assumed in this model
    delta_c_Tinker = 1.686;
    nu = delta_c_Tinker/(sigma);

    gslstatus = gsl_spline_eval_e(cosmo->data.alphahmf, log10(odelta), cosmo->data.accelerator_d,&fit_A); //alpha in Eq. 8
    gslstatus |= gsl_spline_eval_e(cosmo->data.etahmf, log10(odelta), cosmo->data.accelerator_d,&fit_a); //eta in Eq. 8
    gslstatus |= gsl_spline_eval_e(cosmo->data.betahmf, log10(odelta), cosmo->data.accelerator_d,&fit_b); //beta in Eq. 8
    gslstatus |= gsl_spline_eval_e(cosmo->data.gammahmf, log10(odelta), cosmo->data.accelerator_d,&fit_c); //gamma in Eq. 8
    gslstatus |= gsl_spline_eval_e(cosmo->data.phihmf, log10(odelta), cosmo->data.accelerator_d,&fit_d); //phi in Eq. 8;

    fit_a *=pow(a, -0.27);
    fit_b *=pow(a, -0.20);
    fit_c *=pow(a, 0.01);
    fit_d *=pow(a, 0.08);
    if(gslstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_massfunc.c: ccl_massfunc_f():");
      *status |= gslstatus;
      ccl_cosmology_set_status_message(cosmo, "ccl_massfunc.c: ccl_massfunc_f(): interpolation error for Tinker 2010 MF\n");
      return NAN;
    }
    return nu*fit_A*(1.+pow(fit_b*nu,-2.*fit_d))*pow(nu, 2.*fit_a)*exp(-0.5*fit_c*nu*nu);
    break;

  case ccl_watson:
    if(odelta!=200.) {
      *status = CCL_ERROR_HMF_INTERP;
      ccl_cosmology_set_status_message(cosmo, "ccl_massfunc.c: ccl_massfunc_f(): Watson HMF only supported for Delta = 200.\n");
      return NAN;
    }
    // these parameters from: Angulo et al 2012 (arxiv 1203.3216 ) 
    Omega_m_a = ccl_omega_x(cosmo, a, ccl_species_m_label,status);
    fit_A = Omega_m_a*(0.990*pow(a,3.216)+0.074);
    fit_a = Omega_m_a*(5.907*pow(a,3.599)+2.344);
    fit_b = Omega_m_a*(3.136*pow(a,3.058)+2.349);
    fit_c = 1.318;

    return fit_A*(pow(sigma/fit_b,-fit_a)+1.0)*exp(-fit_c/sigma/sigma);

  case ccl_angulo:
    if(odelta!=200.) {
      *status = CCL_ERROR_HMF_INTERP;
      ccl_cosmology_set_status_message(cosmo, "ccl_massfunc.c: ccl_massfunc_f(): Angulo HMF only supported for Delta = 200.\n");
      return NAN;
    }
    // these parameters from: Watson et al 2012 (arxiv 1212.0095 )
    fit_A = 0.201;
    fit_a = 2.08;
    fit_b = 1.7;
    fit_c = 1.172;

    return fit_A*pow( (fit_a/sigma)+1.0, fit_b)*exp(-fit_c/sigma/sigma);

  default:
    *status = CCL_ERROR_MF;
    ccl_cosmology_set_status_message(cosmo ,
        "ccl_massfunc.c: ccl_massfunc(): Unknown or non-implemented mass function method: %d \n",
        cosmo->config.mass_function_method);
    return NAN;
  }
}

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

Definition at line 52 of file ccl_massfunc.c.

References ccl_rho_x(), M_PI, and pow().

Referenced by window_function().

                                                                                           {

  double rho_matter = ccl_rho_x(cosmo, 1., 1, 1, status);

  return pow(halomass*3.0/(4.0*M_PI*rho_matter*odelta),1.0/3.0);

}