ccl_power.h File Reference

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

Go to the source code of this file.

Functions
CCL_BEGIN_DECLS void	ccl_cosmology_write_power_class_z (char *filename, ccl_cosmology *cosmo, double z, int *status)
double	ccl_bcm_model_fka (ccl_cosmology *cosmo, double k, double a, int *status)
double	ccl_linear_matter_power (ccl_cosmology *cosmo, double k, double a, int *status)
double	ccl_nonlin_matter_power (ccl_cosmology *cosmo, double k, double a, int *status)
void	ccl_cosmology_compute_power (ccl_cosmology *cosmo, int *status)
double	ccl_sigmaR (ccl_cosmology *cosmo, double R, double a, int *status)
double	ccl_sigmaV (ccl_cosmology *cosmo, double R, double a, int *status)
double	ccl_sigma8 (ccl_cosmology *cosmo, int *status)

Function Documentation

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

Correction for the impact of baryonic physics on the matter power spectrum. Returns f(k,a) [dimensionless] for given cosmology, using the method specified for the baryonic transfer function. f(k,a) is the fractional change in the nonlinear matter power spectrum from the Baryon Correction Model (BCM) of Schenider & Teyssier (2015). The parameters of the model are passed as part of the cosmology class.

Parameters

cosmo	Cosmology parameters and configurations, including baryonic parameters.
k	Fourier mode, in [1/Mpc] units
a	scale factor, normalized to 1 for today
status	Status flag. 0 if there are no errors, nonzero otherwise. For specific cases see documentation for ccl_error.c

Returns

f(k,a).

Definition at line 574 of file ccl_power.c.

References ccl_parameters::bcm_etab, ccl_parameters::bcm_ks, ccl_parameters::bcm_log10Mc, ccl_parameters::h, ccl_cosmology::params, pow(), and z.

Referenced by ccl_nonlin_matter_power(), and compare_bcm().

                                                                                {

  double fka;
  double b0;
  double bfunc, bfunc4;
  double kg;
  double gf,scomp;
  double kh;
  double z;

  z=1./a-1.;
  kh = k/cosmo->params.h; //convert to h/Mpc
  b0 = 0.105*cosmo->params.bcm_log10Mc-1.27; //Eq. 4.4
  bfunc = b0/(1.+pow(z/2.3,2.5)); //Eq. 4.3
  bfunc4 = (1-bfunc)*(1-bfunc)*(1-bfunc)*(1-bfunc); //B^4(z)
  kg = 0.7*bfunc4*pow(cosmo->params.bcm_etab,-1.6); //Eq. 4.3
  gf = bfunc/(1+pow(kh/kg,3.))+1.-bfunc; //Eq. 4.2, k in h/Mpc
  scomp = 1+(kh/cosmo->params.bcm_ks)*(kh/cosmo->params.bcm_ks); //Eq 4.5, k in h/Mpc
  fka = gf*scomp;
  return fka;
}

void ccl_cosmology_compute_power	(	ccl_cosmology *	cosmo,
		int *	status
	)

Compute the power spectrum and create a 2d spline P(k,z) to be stored in the cosmology structure.

Parameters

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

Returns

void

Definition at line 1403 of file ccl_power.c.

References ccl_bbks, ccl_boltzmann_class, ccl_check_status(), ccl_cosmology_compute_power_bbks(), ccl_cosmology_compute_power_class(), ccl_cosmology_compute_power_eh(), ccl_cosmology_compute_power_emu(), ccl_cosmology_set_status_message(), ccl_eisenstein_hu, ccl_emulator, CCL_ERROR_INCONSISTENT, ccl_cosmology::computed_power, ccl_cosmology::config, and ccl_configuration::transfer_function_method.

Referenced by ccl_linear_matter_power(), and ccl_nonlin_matter_power().

{

  if (cosmo->computed_power) return;
    switch(cosmo->config.transfer_function_method){
        case ccl_bbks:
      ccl_cosmology_compute_power_bbks(cosmo,status);
      break;
        case ccl_eisenstein_hu:
      ccl_cosmology_compute_power_eh(cosmo,status);
      break;
        case ccl_boltzmann_class:
      ccl_cosmology_compute_power_class(cosmo,status);
      break;
        case ccl_emulator:
      ccl_cosmology_compute_power_emu(cosmo,status);
      break;
        default:
      *status = CCL_ERROR_INCONSISTENT;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power(): Unknown or non-implemented transfer function method: %d \n", cosmo->config.transfer_function_method);
    }

    ccl_check_status(cosmo,status);
    if (*status==0){
      cosmo->computed_power = true;
    }
  return;
}

CCL_BEGIN_DECLS void ccl_cosmology_write_power_class_z	(	char *	filename,
		ccl_cosmology *	cosmo,
		double	z,
		int *	status
	)

CLASS power spectrum without splines. Write k, P(k,z) [1/Mpc, Mpc^3] for given cosmology at the k values used within CLASS (spectra.ln_k[]), using the method specified in config.matter_power_spectrum_method.

Parameters

filename	File into which k, P(k,a) will be written
cosmo	Cosmology parameters and configurations
z	Redshift at which the power spectrum is evaluated
status	Status flag. 0 if there are no errors, nonzero otherwise.

Definition at line 596 of file ccl_power.c.

References ccl_cosmology_set_status_message(), CCL_ERROR_CLASS, ccl_fill_class_parameters(), ccl_free_class_structs(), ccl_run_class(), run_pk_param_space::nonlinear, run_pk_param_space::precision, and spectra().

Referenced by write_Pk_CLASS().

{
  struct precision pr;        // for precision parameters
  struct background ba;       // for cosmological background
  struct thermo th;           // for thermodynamics
  struct perturbs pt;         // for source functions
  struct transfers tr;        // for transfer functions
  struct primordial pm;       // for primordial spectra
  struct spectra sp;          // for output spectra
  struct nonlinear nl;        // for non-linear spectra
  struct lensing le;
  struct output op;
  struct file_content fc;

  ErrorMsg errmsg; // for error messages
  // generate file_content structure
  // CLASS configuration parameters will be passed through this structure,
  // to avoid writing and reading .ini files for every call
  int parser_length = 20;
  int init_arr[7]={0,0,0,0,0,0,0};
  if (parser_init(&fc,parser_length,"none",errmsg) == _FAILURE_) {
    *status = CCL_ERROR_CLASS;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: write_power_class_z(): parser init error:%s\n", errmsg);
    return;
  }

  ccl_fill_class_parameters(cosmo,&fc,parser_length, status);

  if (*status != CCL_ERROR_CLASS)
    ccl_run_class(cosmo, &fc,&pr,&ba,&th,&pt,&tr,&pm,&sp,&nl,&le,&op,init_arr,status);

  if (*status == CCL_ERROR_CLASS) {
    //printed error message while running CLASS
    ccl_free_class_structs(cosmo, &ba,&th,&pt,&tr,&pm,&sp,&nl,&le,init_arr,status);
    return;
  }
  if (parser_free(&fc)== _FAILURE_) {
    *status = CCL_ERROR_CLASS;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: write_power_class_z(): Error freeing CLASS parser\n");
    ccl_free_class_structs(cosmo, &ba,&th,&pt,&tr,&pm,&sp,&nl,&le,init_arr,status);
    return;
  }
  FILE *f;
  f = fopen(filename,"w");
  if (!f){
    *status = CCL_ERROR_CLASS;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: write_power_class_z(): Error opening output file\n");
    ccl_free_class_structs(cosmo, &ba,&th,&pt,&tr,&pm,&sp,&nl,&le,init_arr,status);
    fclose(f);
    return;
  }
  double psout_l,ic;
  int s=0;
  for (int i=0; i<sp.ln_k_size; i++) {
    s |= spectra_pk_at_k_and_z(&ba, &pm, &sp,exp(sp.ln_k[i]),z, &psout_l,&ic);
    fprintf(f,"%e %e\n",exp(sp.ln_k[i]),psout_l);
  }
  fclose(f);
  if(s) {
    *status = CCL_ERROR_CLASS;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: write_power_class_z(): Error writing CLASS power spectrum\n");
  }

  ccl_free_class_structs(cosmo, &ba,&th,&pt,&tr,&pm,&sp,&nl,&le,init_arr,status);
}

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

Linear matter power spectrum. Returns P_lin(k,a) [Mpc^3] for given cosmology, using the method specified in cosmo->config.transfer_function_method.

Parameters

cosmo	Cosmology parameters and configurations
k	Fourier mode, in [1/Mpc] units
a	scale factor, normalized to 1 for today
status	Status flag. 0 if there are no errors, nonzero otherwise. For specific cases see documentation for ccl_error.c

Returns

P_lin(k,a).

Definition at line 1506 of file ccl_power.c.

References A_MIN_EMU, ccl_spline_params::A_SPLINE_MINLOG_PK, ccl_cosmology_compute_power(), ccl_cosmology_set_status_message(), ccl_emulator, CCL_ERROR_INCONSISTENT, CCL_ERROR_SPLINE_EV, ccl_growth_factor(), ccl_power_extrapol_highk(), ccl_power_extrapol_lowk(), ccl_raise_gsl_warning(), ccl_splines, ccl_cosmology::computed_power, ccl_cosmology::config, ccl_cosmology::data, ccl_data::k_max_lin, ccl_data::k_min_lin, ccl_data::p_lin, and ccl_configuration::transfer_function_method.

Referenced by ccl_nonlin_matter_power(), ccl_twohalo_matter_power(), compare_bbks(), compare_eh(), compare_power_nu(), lin_pow_spec(), main(), sigmaR_integrand(), and sigmaV_integrand().

{
  if ((cosmo->config.transfer_function_method == ccl_emulator) && (a<A_MIN_EMU)){
    *status = CCL_ERROR_INCONSISTENT;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: the cosmic emulator cannot be used above a=%f\n",A_MIN_EMU);
    return NAN;
  }

  if (!cosmo->computed_power) ccl_cosmology_compute_power(cosmo, status);
  // Return if compilation failed
  //if (cosmo->data.p_lin == NULL) return NAN;
  if (!cosmo->computed_power) return NAN;

  double log_p_1;
  int gslstatus;

  if(a<ccl_splines->A_SPLINE_MINLOG_PK) {  //Extrapolate linearly at high redshift
    double pk0=ccl_linear_matter_power(cosmo,k,ccl_splines->A_SPLINE_MINLOG_PK,status);
    double gf=ccl_growth_factor(cosmo,a,status)/ccl_growth_factor(cosmo,ccl_splines->A_SPLINE_MINLOG_PK,status);

    return pk0*gf*gf;
  }
  if (*status!=CCL_ERROR_INCONSISTENT){
    if(k<=cosmo->data.k_min_lin) {
      log_p_1=ccl_power_extrapol_lowk(cosmo,k,a,cosmo->data.p_lin,cosmo->data.k_min_lin,status);

      return exp(log_p_1);
    }
    else if(k<cosmo->data.k_max_lin){
      gslstatus = gsl_spline2d_eval_e(cosmo->data.p_lin, log(k), a,NULL,NULL,&log_p_1);
      if(gslstatus != GSL_SUCCESS) {
        ccl_raise_gsl_warning(gslstatus, "ccl_power.c: ccl_linear_matter_power():");
        *status = CCL_ERROR_SPLINE_EV;
        ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_linear_matter_power(): Spline evaluation error\n");
        return NAN;
      }
      else {
        return exp(log_p_1);
      }
    }
    else { //Extrapolate using log derivative
      log_p_1 = ccl_power_extrapol_highk(cosmo,k,a,cosmo->data.p_lin,cosmo->data.k_max_lin,status);
      return exp(log_p_1);
    }
  }

  return exp(log_p_1);
}

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

Non-linear matter power spectrum. Returns P_NL(k,a) [Mpc^3] for given cosmology, using the method specified in cosmo->config.transfer_function_method and cosmo->config.matter_power_spectrum_method.

Parameters

cosmo	Cosmology parameters and configurations
k	Fourier mode, in [1/Mpc] units
a	scale factor, normalized to 1 for today
status	Status flag. 0 if there are no errors, nonzero otherwise. For specific cases see documentation for ccl_error.c

Returns

P_NL(k,a).

Definition at line 1562 of file ccl_power.c.

References A_MIN_EMU, ccl_spline_params::A_SPLINE_MINLOG_PK, ccl_configuration::baryons_power_spectrum_method, ccl_bcm, ccl_bcm_model_fka(), ccl_cosmology_compute_power(), ccl_cosmology_set_status_message(), ccl_emu, ccl_emulator, CCL_ERROR_EMULATOR_BOUND, CCL_ERROR_NOT_IMPLEMENTED, CCL_ERROR_SPLINE_EV, ccl_growth_factor(), ccl_halofit, ccl_linear, ccl_linear_matter_power(), ccl_power_extrapol_highk(), ccl_power_extrapol_lowk(), ccl_raise_gsl_warning(), ccl_raise_warning(), ccl_splines, ccl_cosmology::computed_power, ccl_cosmology::config, ccl_cosmology::data, ccl_data::k_max_nl, ccl_data::k_min_nl, ccl_configuration::matter_power_spectrum_method, ccl_data::p_nl, and ccl_configuration::transfer_function_method.

Referenced by ccl_correlation_3d(), cl_integrand(), compare_bcm(), compare_emu(), compare_emu_nu(), compare_power_nu_nl(), main(), non_lin_pow_spec(), and transfer_nc().

{
  double log_p_1, pk;

  switch(cosmo->config.matter_power_spectrum_method) {

  //If the matter PS specified was linear, then do the linear compuation
  case ccl_linear:
    return ccl_linear_matter_power(cosmo,k,a,status);

  case ccl_halofit:
    if (!cosmo->computed_power)
      ccl_cosmology_compute_power(cosmo, status);
    if (cosmo->data.p_nl == NULL) return NAN; // Return if computation failed

    if(a<ccl_splines->A_SPLINE_MINLOG_PK) { //Extrapolate linearly at high redshift
      double pk0=ccl_nonlin_matter_power(cosmo,k,ccl_splines->A_SPLINE_MINLOG_PK,status);
      double gf=ccl_growth_factor(cosmo,a,status)/ccl_growth_factor(cosmo,ccl_splines->A_SPLINE_MINLOG_PK,status);
      return pk0*gf*gf;
    }
        break;

  case ccl_emu:
    if ((cosmo->config.transfer_function_method == ccl_emulator) && (a<A_MIN_EMU)){
      *status = CCL_ERROR_EMULATOR_BOUND;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: the cosmic emulator cannot be used above z=2\n");
      return NAN;
    }

    // Compute power spectrum if needed; return if computation failed
    if (!cosmo->computed_power){
      ccl_cosmology_compute_power(cosmo,status);
    }
    if (cosmo->data.p_nl == NULL) return NAN;
        break;

  default:
    ccl_raise_warning(
            CCL_ERROR_NOT_IMPLEMENTED,
      "config.matter_power_spectrum_method = %d not yet supported "
            "continuing with linear power spectrum\n", cosmo->config.matter_power_spectrum_method);
    cosmo->config.matter_power_spectrum_method=ccl_linear;
    return ccl_linear_matter_power(cosmo,k,a,status);
  } // end switch

  // if we get here, try to evaluate the power spectrum
  // we need to account for bounds below and above
  if (k <= cosmo->data.k_min_nl) {
    // we assume no baryonic effects below k_min_nl
    log_p_1 = ccl_power_extrapol_lowk(cosmo, k, a, cosmo->data.p_nl, cosmo->data.k_min_nl, status);
    return exp(log_p_1);
  }

  if (k < cosmo->data.k_max_nl) {
    int gslstatus = gsl_spline2d_eval_e(cosmo->data.p_nl, log(k), a, NULL ,NULL, &log_p_1);
    if (gslstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_power.c: ccl_nonlin_matter_power():");
      *status = CCL_ERROR_SPLINE_EV;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_nonlin_matter_power(): Spline evaluation error\n");
      return NAN;
    } else {
      pk = exp(log_p_1);
    }
  } else {
    // Extrapolate NL regime using log derivative
    log_p_1 = ccl_power_extrapol_highk(cosmo, k, a, cosmo->data.p_nl, cosmo->data.k_max_nl, status);
    pk = exp(log_p_1);
  }

    // Add baryonic correction
  if (cosmo->config.baryons_power_spectrum_method == ccl_bcm) {
    int pwstatus = 0;
    double fbcm = ccl_bcm_model_fka(cosmo, k, a, &pwstatus);
    pk *= fbcm;
    if (pwstatus) {
      *status = CCL_ERROR_SPLINE_EV;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_nonlin_matter_power(): Error in BCM correction\n");
      return NAN;
    }
  }

  return pk;
}

double ccl_sigma8	(	ccl_cosmology *	cosmo,
		int *	status
	)

Computes sigma8, variance of the matter density field with (top-hat) smoothing scale R = 8 Mpc/h, from linear power spectrum. Returns sigma8 for specified cosmology.

Parameters

cosmo	Cosmology parameters and configurations
status	Status flag. 0 if there are no errors, nonzero otherwise. For specific cases see documentation for ccl_error.c

Returns

sigma8.

Definition at line 1777 of file ccl_power.c.

References ccl_sigmaR(), ccl_parameters::h, and ccl_cosmology::params.

Referenced by ccl_cosmology_compute_power_bbks(), ccl_cosmology_compute_power_eh(), main(), and norm_pwr().

{
  return ccl_sigmaR(cosmo,8/cosmo->params.h, 1.,status);
}

double ccl_sigmaR	(	ccl_cosmology *	cosmo,
		double	R,
		double	a,
		int *	status
	)

Variance of the matter density field with (top-hat) smoothing scale R [Mpc]. Returns sigma(R) for specified cosmology at a = 1.

Parameters

cosmo	Cosmology parameters and configurations
R	Smoothing scale, in [Mpc] units
a	scale factor
status	Status flag. 0 if there are no errors, nonzero otherwise. For specific cases see documentation for ccl_error.c

Returns

sigma(R).

Definition at line 1715 of file ccl_power.c.

References ccl_growth_factor(), ccl_gsl, ccl_raise_gsl_warning(), ccl_splines, cl_cmbl_bm::cosmo, SigmaR_pars::cosmo, ccl_gsl_params::INTEGRATION_SIGMAR_EPSREL, ccl_spline_params::K_MAX, ccl_spline_params::K_MIN, M_PI, ccl_gsl_params::N_ITERATION, SigmaR_pars::R, sigmaR_integrand(), sqrt(), and SigmaR_pars::status.

Referenced by ccl_cosmology_compute_sigma(), ccl_sigma8(), and main().

{
  SigmaR_pars par;
  par.status = status;

  par.cosmo=cosmo;
  par.R=R;
  gsl_integration_cquad_workspace *workspace=gsl_integration_cquad_workspace_alloc(ccl_gsl->N_ITERATION);
  gsl_function F;
  F.function=&sigmaR_integrand;
  F.params=∥
  double sigma_R;
  int gslstatus = gsl_integration_cquad(&F, log10(ccl_splines->K_MIN), log10(ccl_splines->K_MAX),
                                                0.0, ccl_gsl->INTEGRATION_SIGMAR_EPSREL,
                                        workspace,&sigma_R,NULL,NULL);
  if(gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_power.c: ccl_sigmaR():");
    *status |= gslstatus;
  }

  gsl_integration_cquad_workspace_free(workspace);

  return sqrt(sigma_R*M_LN10/(2*M_PI*M_PI))*ccl_growth_factor(cosmo, a, status);
}

double ccl_sigmaV	(	ccl_cosmology *	cosmo,
		double	R,
		double	a,
		int *	status
	)

Variance of the displacement field with (top-hat) smoothing scale R [Mpc] Returns sigma(V(R)) for specified cosmology at a = 1.

Parameters

cosmo	Cosmology parameters and configurations
R	smoothing scale, in [Mpc] units
a	scale factor
status	Status flag. 0 if there are no errors, nonzero otherwise. For specific cases see documentation for ccl_error.c

Returns

sigma(R).

Definition at line 1746 of file ccl_power.c.

References ccl_growth_factor(), ccl_gsl, ccl_raise_gsl_warning(), ccl_splines, cl_cmbl_bm::cosmo, SigmaV_pars::cosmo, ccl_gsl_params::INTEGRATION_SIGMAR_EPSREL, ccl_spline_params::K_MAX, ccl_spline_params::K_MIN, M_PI, ccl_gsl_params::N_ITERATION, SigmaV_pars::R, sigmaV_integrand(), sqrt(), and SigmaV_pars::status.

{
  SigmaV_pars par;
  par.status = status;

  par.cosmo=cosmo;
  par.R=R;
  gsl_integration_cquad_workspace *workspace=gsl_integration_cquad_workspace_alloc(ccl_gsl->N_ITERATION);
  gsl_function F;
  F.function=&sigmaV_integrand;
  F.params=∥
  double sigma_V;
    int gslstatus = gsl_integration_cquad(&F, log10(ccl_splines->K_MIN), log10(ccl_splines->K_MAX),
                                                                                0.0, ccl_gsl->INTEGRATION_SIGMAR_EPSREL,
                                                                                workspace,&sigma_V,NULL,NULL);

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

  gsl_integration_cquad_workspace_free(workspace);

  return sqrt(sigma_V*M_LN10/(2*M_PI*M_PI))*ccl_growth_factor(cosmo, a, status);
}