CLUMPY  v3.0.1
Developer documentation
cosmo.h File Reference

Cosmology base functions: distances, Omega parameters, and mass functions. More...

#include <vector>
#include <string>
Include dependency graph for cosmo.h:
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes

struct  Mmin_params_for_rootfinding
 
struct  pk_params_for_integration
 
struct  d_to_z_params_for_rootfinding
 

Functions

void compute_class_pk (double z)
 Calls CLASS to compute the matter power spectrum at redshift \(z\). More...
 
vector< double > compute_massfunction (const double z, const vector< double > &Mh_grid, const int card_mf, const vector< double > &linear_lnk_vec, const vector< double > &linear_lnp_vec, const int card_window, const int card_growth_mthd)
 Computes \(dn/d(\ln(M))\) (the number density per mass bin) for each z in vec_z and mass in vec_lnM. More...
 
double d_to_z (const double &d, const int switch_d)
 Computes the redshift from a given comoving distance. More...
 
double dh_c (const double &z)
 Computes the comoving distance by integration. More...
 
double dh_trans (const double &z)
 Computes the transverse comoving distance by integration. More...
 
double dh_a (const double &z)
 Angular diameter distance at redshift \(z\) in [ \(h^{-1} \; Mpc \)]. More...
 
double dh_l (const double &z)
 Luminosity distance at redshift \(z\) in [ \(h^{-1} \; Mpc \)]. More...
 
double dNdVhdlnMh_interpol (const double &z, const double &Mh)
 
void dNdVhdlnMh_integrand_Mh (double &M, double par_z[1], double &res)
 
void dNdVhdlnMh_integrand_logMh (double &logMh, double par_z[1], double &res)
 
void dNdVhdlnMh_sigmoid (double &Mh, double par[3], double &res)
 
double dNdVhdlnMh_literature (const double &z, const double &Mh, const int card_mf, vector< string > &label_literature)
 
double dNdOmega (double &zmin, double &zmax, double const &Mmin, double const &Mmax)
 
void dNdOmegadlnMh (double &lnMh, double par_z[2], double &res)
 
void dNdOmegadz (double &z, double par_Mh[2], double &res)
 
void dNdOmegadlnMhdz (double &z, double par_Mh[1], double &res)
 
double ds2dlnk (double lnk, void *p)
 Returns \( d(\sigma^2)/dk \propto k^3\;p(k)\; W(kr)\), where r is the scale of interest and \( W(kr)\) is the Fourier transform of the 3D top-hat function. More...
 
double dVhdzdOmega (double z)
 Differential comoving volume element in [ \(h^{-3} \; Mpc^3 \; sr^{-1}\)]. More...
 
string legend_for_delta ()
 
double lnsigma2 (double lnR, void *p)
 Returns \(\ln(\sigma^2)\) for the scale \( R \) passed as \(\ln(R)\). More...
 
double gamma_f (const int card_window)
 
void get_pk (const double &z, vector< double > &lnk_vec, vector< double > &lnp_vec, const bool is_nonlinear=false, const bool is_force_recompute=false)
 Loads the matter power spectrum from a file. Fills in vec_lnk with mode number \(\ln(k)\), \( k \) in [ \( h \;Mpc^{-1}\)]. Fills in linear_lnpkz with \(\ln(p_{\rm matter}(k,z))\) [ \( Mpc^3\;h^{-3}\)]. If the p_matter(k,z) file is not present, calls CLASS to compute it. More...
 
double growth_factor (const double &z, const int card_growth_mthd)
 Computes the linear growth rate factor. More...
 
double growth_integrand (const double z, void *params=NULL)
 
double mf (int card_mf, double const &nu, double const &z, double const &Delta_c)
 
double mf_bolshoi_planck (double const &nu, double const &z, double const &Delta_c)
 
double mf_magneticum (double const &nu, double const &z, double const &Delta_c, const bool is_hydro)
 
double mf_press_schechter (double const &nu, double const &z, double const &Delta_c)
 
double mf_sheth_tormen (double const &nu, double const &z, double const &Delta_c)
 
double mf_tinker (double const &nu, double const &z, double const &Delta_c)
 
double mf_tinker10 (double const &nu, double const &z, double const &Delta_c)
 
double mf_tinker_norm (double const &nu, double const &z, double const &Delta_c)
 
double mf_jenkins (double const &nu, double const &z, double const &Delta_c)
 
vector< vector< double > > nu_andderiv (const double z, const vector< double > &Mh_vec, const vector< double > &linear_lnk_vec, const vector< double > &linear_lnp_vec, const int card_window, const int card_growth_mthd)
 
double H2_over_H02 (const double &z)
 
double H_over_H0 (const double &z)
 \(H(z)/H_0\) where \(H(z)\) is the Hubble rate at redshift \(z\). More...
 
double H0_over_H (const double z, void *params=NULL)
 \(1/h(z) = H_0/H(z)\) where \(H(z)\) is the Hubble rate at redshift \(z\). More...
 
double Omega_r (const double &z)
 
double Omega_m (const double &z)
 
double Omega_cdm (const double &z)
 
double Omega_b (const double &z)
 
double Omega_lambda (const double &z)
 
double Omega_k (const double &z)
 
double rho_crit (const double &z)
 
vector< vector< double > > sigma2_andderiv (const double z, const vector< double > &Mh_vec, const vector< double > &linear_lnk_vec, const vector< double > &linear_lnp_vec, const int card_window, const int card_growth_mthd)
 
double sigma8 (vector< double > &lnk_vec, vector< double > &lnp_vec, const int card_window)
 
double solve_Mhmin (double Mhmin, void *p)
 
double solve_d_to_z (double d, void *p)
 
double window_fnct (const double &x, const int card_window)
 

Detailed Description

Cosmology base functions: distances, Omega parameters, and mass functions.

Definition in file cosmo.h.

Function Documentation

◆ compute_class_pk()

void compute_class_pk ( double  z)

Calls CLASS to compute the matter power spectrum at redshift \(z\).

Definition at line 32 of file cosmo.cc.

33 {
34  //--- Uses CLASS to compute the linear and nonlinear matter power spectrum p(k,z) at redshift z.
35  // Writes on disk file the output where [k]=h Mpc^-1 and [p(k,z)]=h^-3 Mpc^3.
36  // NOTE: We use CLASS for the computation up to kmax *h = 10000 and use the analytic
37  // result P(k) ~ k^(n_s-4) * ln^2(k) for larger k.
38  // z Redshift
39 
40  char char_tmp[512];
41 
42  bool is_extrapolate = false;
43  double kmax = gSIM_EXTRAGAL_KMAX_PRECOMP;
44  double kmax_class = 1e4;
45  if (z > 3.5) kmax_class = 1e4;
46 
47  if (kmax > kmax_class) {
48  sprintf(char_tmp, "Will use CLASS only to compute P(k) up to k = %d, "
49  "and extrapolate linear P(k) ~ k^(n_s-4) * ln^2(k) for larger k.", int(kmax_class));
50  print_warning("cosmo.cc", "compute_class_pk()", string(char_tmp));
51  is_extrapolate = true;
52  kmax = kmax_class;
53  }
54 
55  char const *gPATH_TO_CLASS_tmp = getenv("CLASS");
56 
57  string line = "";
58  ostringstream convert;
60  string h_str = convert.str();
61  convert.str("");
63  string OmegaM_str = convert.str();
64  convert.str("");
66  string OmegaB_str = convert.str();
67  convert.str("");
69  string OmegaL_str = convert.str();
70  convert.str("");
72  string OmegaCDM_str = convert.str();
73  convert.str("");
75  string OmegaK_str = convert.str();
76  convert.str("");
78  string ns_str = convert.str();
79  convert.str("");
81  string tau_str = convert.str();
82  convert.str("");
83  convert << kmax;
84  string kmax_str = convert.str();
85  convert.str("");
86 
87  sprintf(char_tmp, "%.3g", z);
88  if (z > 10) sprintf(char_tmp, "%.5g", z);
89  string zstr = string(char_tmp);
90 
91  string filename_base = "h" + h_str + "_OmegaB"
92  + OmegaB_str + "_OmegaM" + OmegaM_str + "_OmegaL" + OmegaL_str
93  + "_ns" + ns_str + "_tau" + tau_str + "_z" + zstr;
94 
95  string filename_template = "PREFIX_" + filename_base + "_XX.dat";
96 
97  if (gPATH_TO_CLASS_tmp == NULL) {
98  printf("\n====> ERROR: compute_class_pk() in cosmo.cc");
99  printf("\n Environmental variable CLASS not set.");
100  printf("\n You have two options:");
101  printf("\n 1.) Set CLASS variable and point it to an installation of CLASS, http://class-code.net/.");
102  printf("\n 2.) Provide manually an ASCII file with the name %s", filename_template.c_str());
103  printf("\n in the directory %sdata/pk_precomp/", gPATH_TO_CLUMPY_HOME.c_str());
104  printf("\n and PREFIX of your choice, and XX either 'lin' or 'nl', if the linear");
105  printf("\n or non-linear power spectrum is required, respectively.");
106  printf("\n The file must have two columns containing ");
107  printf("\n k in units of [h/Mpc] and the P(k) in units of [(Mpc/h)^3].");
108  printf("\n Make also sure that kmax >= %d in that case,", int(ceil(gSIM_EXTRAGAL_KMAX_PRECOMP)));
109  printf("\n or lower the hidden parameter gSIM_EXTRAGAL_KMAX_PRECOMP.\n");
110  printf("\n => abort()\n\n");
111  abort();
112  } else gPATH_TO_CLASS = string(gPATH_TO_CLASS_tmp);
113 
114  string infile = gPATH_TO_CLUMPY_HOME + "/data/pk_precomp/class_input_template.ini";
115  string inifile = gPATH_TO_CLASS + "/class_input_for_clumpy.ini";
116 
117  string command = "cd " + gPATH_TO_CLASS +
118  "; mkdir -p output_tmp/; mv output/* output_tmp/";
119 
120  sys_safe_execution(command);
121 
122  ifstream filein(infile.c_str());
123  ofstream myfile(inifile.c_str());
124 
125  while (getline(filein, line)) {
126  line.erase(remove_if(line.begin(), line.end(), isNotASCII), line.end());
127  if (line.substr(0, 1) == string("h")) line = "h = " + h_str;
128  else if (line.substr(0, 7) == string("Omega_b")) line = "Omega_b = " + OmegaB_str;
129  else if (line.substr(0, 9) == string("Omega_cdm")) line = "Omega_cdm = " + OmegaCDM_str;
130  else if (line.substr(0, 7) == string("Omega_k")) line = "Omega_k = " + OmegaK_str;
131  else if (line.substr(0, 3) == string("n_s")) line = "n_s = " + ns_str;
132  else if (line.substr(0, 8) == string("tau_reio")) line = "tau_reio = " + tau_str;
133  else if (line.substr(0, 13) == string("P_k_max_h/Mpc")) line = "P_k_max_h/Mpc = " + kmax_str;
134  else if (line[0] == 'z')line = "z_pk = " + zstr;
135  myfile << line << endl;
136  }
137  filein.close();
138  myfile.close();
139 
140  string filename_lin = "class_" + filename_base + "_lin.dat";
141 
142  string filename_nl = "class_" + filename_base + "_nl.dat";
143 
144  command = "cd " + gPATH_TO_CLASS +
145  "; ./class class_input_for_clumpy.ini; cp output/class_input_for_clumpy00_pk.dat "
146  + gPATH_TO_CLUMPY_HOME + "/data/pk_precomp/" + filename_lin + "; "
147  + "cp output/class_input_for_clumpy00_pk_nl.dat "
148  + gPATH_TO_CLUMPY_HOME + "/data/pk_precomp/" + filename_nl + "; "
149  + "rm -f output/*; mv output_tmp/* output; rm -rf output_tmp; "
150  ;//+ "rm class_input_for_clumpy.ini";
151 
152  sys_safe_execution(command);
153 
154  if (is_extrapolate) {
155  // we will append extrapolated P(k) ~ k^(n_s-4) * ln(k) to the file.
156  vector<double> linear_lnk_vec;
157  vector<double> linear_lnp_vec;
158  string file = gPATH_TO_CLUMPY_HOME + "/data/pk_precomp/" + filename_lin;
159  vector< vector<double> > table = read_ascii(file);
160  for (int i = 0; i < int(table.size()); ++i) {
161  for (int j = 0; j < int(table[i].size()); ++j) {
162  if (j == 0) linear_lnk_vec.push_back(log(table[i][j]));
163  else if (j == 1) linear_lnp_vec.push_back(log(table[i][j]));
164  }
165  }
166 
167  int n_k_old = int(linear_lnk_vec.size());
168  double delta_lnk = linear_lnk_vec[linear_lnk_vec.size() - 1] - linear_lnk_vec[linear_lnk_vec.size() - 2];
169 
170  // calculate normalisation:
171  double exponent = gCOSMO_N_S - 4.;
172  double lnnorm = linear_lnp_vec[linear_lnp_vec.size() - 1] - exponent * linear_lnk_vec[linear_lnk_vec.size() - 1] - 2. * log(linear_lnk_vec[linear_lnk_vec.size() - 1] + log(gCOSMO_HUBBLE));
173 
174 
175  while (linear_lnk_vec[linear_lnk_vec.size() - 1] < log(gSIM_EXTRAGAL_KMAX_PRECOMP)) {
176  linear_lnk_vec.push_back(linear_lnk_vec[linear_lnk_vec.size() - 1] + delta_lnk);
177  linear_lnp_vec.push_back(lnnorm + exponent * linear_lnk_vec[linear_lnk_vec.size() - 1] + 2. * log(linear_lnk_vec[linear_lnk_vec.size() - 1] + log(gCOSMO_HUBBLE)));
178  };
179  int n_k_new = int(linear_lnk_vec.size());
180 
181  FILE *fp;
182  fp = fopen(file.c_str(), "a");
183  for (int i = n_k_old; i < n_k_new; ++i) {
184  fprintf(fp, " %.12e %.12e \n", exp(linear_lnk_vec[i]), exp(linear_lnp_vec[i]));
185  }
186  fclose(fp);
187  }
188 
189  return;
190 }
string gPATH_TO_CLASS
Absolute path to CLASS executable.
Definition: params.cc:297
double gCOSMO_OMEGA0_LAMBDA
Present day dark energy density of the -CDM Universe (PDG)
Definition: params.cc:77
double gCOSMO_OMEGA0_K
See the user documentation
Definition: params.cc:65
double gCOSMO_HUBBLE
See the user documentation
Definition: params.cc:62
double gCOSMO_OMEGA0_B
See the user documentation
Definition: params.cc:64
void print_warning(string library, string function, string message)
Definition: misc.cc:986
def convert(infile='empty', i_ext=-1, i_col=-1, coordsys_out='-1')
vector< vector< double > > read_ascii(string const &filename, int const ncolumns=-1, bool const is_monotonous=false, int const column_start=0)
Definition: misc.cc:1003
void sys_safe_execution(string const &command)
Definition: misc.cc:1791
bool isNotASCII(const char s)
Returns true if a character is not ASCII, false otherwise.
Definition: misc.cc:593
string gPATH_TO_CLUMPY_HOME
Absolute path to CLUMPY home directory.
Definition: params.cc:299
double round(const double &x, const int digits)
Definition: misc.cc:1570
double gCOSMO_N_S
See the user documentation
Definition: params.cc:67
double gSIM_EXTRAGAL_KMAX_PRECOMP
See the user documentation
Definition: params.cc:232
double gCOSMO_TAU_REIO
See the user documentation
Definition: params.cc:68
double gCOSMO_OMEGA0_M
See the user documentation
Definition: params.cc:63

◆ compute_massfunction()

vector<double> compute_massfunction ( const double  z,
const vector< double > &  Mh_grid,
const int  card_mf,
const vector< double > &  linear_lnk_vec,
const vector< double > &  linear_lnp_vec,
const int  card_window,
const int  card_growth_mthd 
)

Computes \(dn/d(\ln(M))\) (the number density per mass bin) for each z in vec_z and mass in vec_lnM.

Result is in [ \(h^3 Mpc^{-3}\)] is stored in dndlnMh. The mass function parametrisation is selected by card_mf = kTINKER08, kJENKINS01, kSHETHTORMEN99. Warning: the mass is in units of \(\Omega_m^0\; h^{-1} \; M_\odot\)

Parameters
[in]zredshift \(z\)
[in]Mh_gridMass function grid, \([\Omega_{m0} h^{-1} M_{\odot}]\)
[in]card_mfmass function keyword
[in]linear_lnk_vec\(log(k/h)\) vector \([h \rm Mpc^{-1}]\)
[in]linear_lnp_vec\(log(P_{\rm lin})\) vector \([h^{-3} \rm Mpc^3 ]\)
[in]card_windowwindow function keyword
[in]card_growth_mthdselects the method to obtain \(P(k,z)\): either \(P(k,z)\) is directly given (kPKZ_FROMFILE),
or \(P(k,0)\) is given, and \(\sigma^2(z)\) is calculated from a growth factor model.

Definition at line 193 of file cosmo.cc.

196 {
197  // Computes dn/dln(Mh) (the number density per mass bin) at redshift z.
198  // The mass function parametrisation is selected by card_mf = kTINKER08, kJENKINS01, kSHETHTORMEN99....
199  // Note 1: Mh is in units of Omega_m0 h^-1 Msol
200  // Note 2: dn/dln(Mh) = Mh * dn/d(Mh) = M dn/dM, in units of halos/ (h^-1 Mpc)^3 (M in units of Msol).
201 
202  // Input:
203  // z redshift z
204  // Mh_grid mass values, length n_m, unit [Omega_m0 h^-1 Msol]
205  // card_mf selects the mass function to use (e.g. Tinker, Jenkins, etc.)
206  // linear_lnk_vec k values for p_k
207  // linear_lnp_vec lnp, must have same dimension as linear_lnk_vec
208  // card_window window function
209  // card_growth_mthd selects the method to obtain P(k,z): either P(k,z) is directly given (kPKZ_FROMFILE),
210  // or P(k,0) is given, and sigma^2(z) is calculated from a growth factor model.
211  // Output:
212  // dNdVhdlnMh_z mass function over tabulated masses at redshift z,
213  // length n_m, units [halos/ (h^-1 Mpc)^3].
214 
215  double rho_critperh2_MsolperMpc3 = RHO_CRITperh2_MSOLperKPC3 / KPC3_to_MPC3 ; // * in units of h^2 * Msol Mpc^-3
216  vector<double> dNdVhdlnMh_z;
217 
218  char char_tmp[80];
219  sprintf(char_tmp, "%.3g", z);
220  if (z > 10) sprintf(char_tmp, "%.5g", z);
221  string zstr = string(char_tmp);
222 
223  cout << " Computing mass function for z = " << zstr << " and given cosmology ...";
224 
225  // check sizes:
226  if (linear_lnk_vec.size() != linear_lnp_vec.size()) {
227  printf("\n====> ERROR: compute_massfunction() in cosmo.cc");
228  printf("\n Dimension of input arrays do not match:");
229  printf("\n len(linear_lnk_grid) = %d", int(linear_lnk_vec.size()));
230  printf("\n len(linear_lnplnkz[) = %d", int(linear_lnp_vec.size()));
231  printf("\n => abort()\n\n");
232  abort();
233  }
234 
235  // Determine now the reference m_delta for mf descriptions.
236  // We ca either (i) convert all masses via the mdelta1_to_mdelta2 trafo, or
237  // (ii) sometimes, different fitting formulae are provided in the literature for different Delta.
238  double Delta_c_mf;
239  vector<double> Mh_Delta_grid;
241  // use interpolations from literature, if available.
243  Mh_Delta_grid = Mh_grid;
244  } else {
245  printf("\n ... transform Delta from halo profile ...");
246  // calculate mf for Delta_crit = 200 (178, ...), and connect to input delta
247  // by mdelta1_to_mdelta2 conversion:
251  Delta_c_mf = 178; // w.r.t. the critical density
253  Delta_c_mf = delta_x_to_delta_crit(-1, kBRYANNORMAN98, z); // w.r.t. the critical density
254  } else {
255  Delta_c_mf = 200; // w.r.t. the critical density
256  }
257  double Delta_in = delta_x_to_delta_crit(gCOSMO_DELTA0, gCOSMO_FLAG_DELTA_REF, z); // w.r.t. the critical density
258  double par_prof[4] = {gHALO_SUBS_SHAPE_PARAMS[kEXTRAGAL][0],
262  };
264  print_warning("cosmo.cc", "compute_massfunction()",
265  "Mdelta conversion might fail for large masses and kPRADA12_200 cdelta-mdelta parametrization.");
266  }
267  for (int j = 0; j < (int)Mh_grid.size(); j++) {
268  // get the M_200,c (M_178,c) corresponding to user-chosen input mass_Delta_in,c:
269  double Mh_Delta_mf = gCOSMO_M_to_MH * mdelta1_to_mdelta2(Mh_grid[j] / gCOSMO_M_to_MH, Delta_in,
270  gHALO_SUBS_FLAG_CDELTAMDELTA[kEXTRAGAL], par_prof, z, Delta_c_mf);
271  // cout << Mh_Delta_mf << endl;
272  Mh_Delta_grid.push_back(Mh_Delta_mf);
273  }
274  printf(" done.\n");
275  }
276 
277  //Get nu (nu_deriv[0]) and d ln nu / dlnMh (nu_deriv[1]):
278  vector< vector<double> > nu_deriv = nu_andderiv(z, Mh_Delta_grid, linear_lnk_vec, linear_lnp_vec, card_window, card_growth_mthd);
279 
280  //Now, compute mass function!
281  for (int j = 0; j < (int)Mh_Delta_grid.size(); j++) {
282  double dNdVhdlnMh = nu_deriv[1][j] / Mh_Delta_grid[j] * rho_critperh2_MsolperMpc3 * mf(card_mf, nu_deriv[0][j], z, Delta_c_mf);
283  if (gDM_IS_IDM) {
284  double factor1 = pow(1 + gEXTRAGAL_IDM_MHALFMODE / Mh_Delta_grid[j] * gCOSMO_M_to_MH / gEXTRAGAL_IDM_BETA, gEXTRAGAL_IDM_ALPHA);
285  double factor2 = pow(1 + gEXTRAGAL_IDM_MHALFMODE / Mh_Delta_grid[j] * gCOSMO_M_to_MH / gEXTRAGAL_IDM_GAMMA, gEXTRAGAL_IDM_DELTA);
286  dNdVhdlnMh *= (factor1 * factor2);
287  }
288  if (dNdVhdlnMh < 1e-40) dNdVhdlnMh = 1e-40;
289 // if (Mh_grid[j] / gCOSMO_M_to_MH >= 1e10) {
290 // cout << "z = " << z << "\t M = " << Mh_grid[j] / gCOSMO_M_to_MH << "\t dNdVhdlnMh = " << dNdVhdlnMh * pow(gCOSMO_HUBBLE,3) << endl; abort();
291 // }
292  dNdVhdlnMh_z.push_back(dNdVhdlnMh);
293  }
294  cout << " done." << endl;
295 
296  return dNdVhdlnMh_z;
297 }
double mdelta1_to_mdelta2(double const &mdelta1, double const &Delta1, int const card_cdelta, double par_prof[4], double const &z, double Delta2=-1., double r=-1, double Rvir=0.)
Definition: clumps.cc:3029
double gEXTRAGAL_IDM_GAMMA
See the user documentation
Definition: params.cc:148
bool gSIM_EXTRAGAL_IS_DELTA_HMF_FROM_PROFILE
See the user documentation
Definition: params.cc:240
int gHALO_SUBS_FLAG_CDELTAMDELTA[gN_TYPEHALOES]
See the user documentation
Definition: params.cc:127
double gEXTRAGAL_IDM_BETA
See the user documentation
Definition: params.cc:147
double gCOSMO_DELTA0
See the user documentation
Definition: params.cc:71
int gHALO_SUBS_FLAG_PROFILE[gN_TYPEHALOES]
See the user documentation
Definition: params.cc:128
double gEXTRAGAL_IDM_DELTA
See the user documentation
Definition: params.cc:149
void print_warning(string library, string function, string message)
Definition: misc.cc:986
double gCOSMO_M_to_MH
conversion.
Definition: params.cc:79
int gCOSMO_FLAG_DELTA_REF
See the user documentation
Definition: params.cc:72
int gEXTRAGAL_FLAG_MASSFUNCTION
See the user documentation
Definition: params.cc:143
#define KPC3_to_MPC3
Definition: inlines.h:23
double gEXTRAGAL_IDM_ALPHA
See the user documentation
Definition: params.cc:146
See user documentation
Definition: params.h:121
See user documentation
Definition: params.h:123
bool gDM_IS_IDM
See the user documentation
Definition: params.cc:94
double delta_x_to_delta_crit(double const &Delta_x0, int const card_delta_ref, double const &z)
Definition: clumps.cc:102
double gEXTRAGAL_IDM_MHALFMODE
See the user documentation
Definition: params.cc:145
Generic cosmic halo object with same properties on all scales from cluster to mictrohalo size.
Definition: params.h:161
double mf(int card_mf, double const &nu, double const &z, double const &Delta_c)
Definition: cosmo.cc:1711
See user documentation
Definition: params.h:175
See user documentation
Definition: params.h:124
vector< vector< double > > nu_andderiv(const double z, const vector< double > &Mh_vec, const vector< double > &linear_lnk_vec, const vector< double > &linear_lnp_vec, const int card_window, const int card_growth_mthd)
Definition: cosmo.cc:2242
See user documentation
Definition: params.h:122
See user documentation
Definition: params.h:84
double gHALO_SUBS_SHAPE_PARAMS[gN_TYPEHALOES][gN_SHAPE_PARAMS]
See the user documentation
Definition: params.cc:129
#define RHO_CRITperh2_MSOLperKPC3
Critical density of the universe in units of .
Definition: inlines.h:7

◆ d_to_z()

double d_to_z ( const double &  d,
const int  switch_d 
)

Computes the redshift from a given comoving distance.

Definition at line 300 of file cosmo.cc.

301 {
302  // computes the redshift from a given distance.
303 
304  // Input:
305  // d distance [Mpc]
306  // switch_d which distance? 0: d_c, 1: d_trans, 2: d_l, 3: d_a
307  // Output:
308  // z redshift
309 
310  gsl_function F;
311  F.function = &solve_d_to_z;
312  d_to_z_params_for_rootfinding params = {d, switch_d};
313  F.params = &params;
314  double zmin = 0.;
315  double zmax = 1e6;
316  return rootsolver_gsl(gsl_root_fsolver_brent, F, zmin, zmax, gSIM_EPS);
317 }
double gSIM_EPS
See the user documentation
Definition: params.cc:196
double rootsolver_gsl(const gsl_root_fsolver_type *T, gsl_function Fx, double &x_lo, double &x_hi, double &precision)
Definition: misc.cc:1519
double solve_d_to_z(double z, void *p)
Definition: cosmo.cc:2481

◆ dh_a()

double dh_a ( const double &  z)

Angular diameter distance at redshift \(z\) in [ \(h^{-1} \; Mpc \)].

Definition at line 370 of file cosmo.cc.

371 {
372  // computes the angular diameter distance dh_a=dh_trans/(1+z) by integration
373  // for the cosmology defined by the global gCOSMO_OMEGA parameters.
374  // units = h^-1 Mpc
375 
376  // Input:
377  // z redshift
378  // Output:
379  // dh_a angular diameter distance [h^-1 Mpc]
380  return dh_trans(z) / (1 + z);
381 }
double dh_trans(const double &z)
Computes the transverse comoving distance by integration.
Definition: cosmo.cc:342

◆ dh_c()

double dh_c ( const double &  z)

Computes the comoving distance by integration.

Definition at line 320 of file cosmo.cc.

321 {
322  // computes the comoving line-of-sight distance by integration
323  // for the cosmology defined by the global gCOSMO_OMEGA parameters.
324  // units = h^-1 Mpc
325 
326  // Input:
327  // z redshift
328  // Output:
329  // dh_c comoving distance [h^-1 Mpc]
330 
331  gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
332  double eta, error;
333  gsl_function F;
334  F.function = &H0_over_H;
335 
336  gsl_integration_qags(&F, 0, z, 0, gSIM_EPS, 1000, w, &eta, &error);
337  gsl_integration_workspace_free(w);
338  return HUBBLE_LENGTHxh_Mpc * eta;
339 }
double gSIM_EPS
See the user documentation
Definition: params.cc:196
#define HUBBLE_LENGTHxh_Mpc
Hubble length in [ ].
Definition: inlines.h:5
double H0_over_H(const double z, void *params)
where is the Hubble rate at redshift .
Definition: cosmo.cc:2300

◆ dh_l()

double dh_l ( const double &  z)

Luminosity distance at redshift \(z\) in [ \(h^{-1} \; Mpc \)].

Definition at line 384 of file cosmo.cc.

385 {
386  // computes the luminosity distance dl=dc * (1+z) by integration
387  // for the cosmology defined by the global gCOSMO_OMEGA parameters.
388  // units = h^-1 Mpc
389 
390  // Input:
391  // z redshift
392  // Output:
393  // dh_l luminosity distance [h^-1 Mpc]
394  return dh_trans(z) * (1 + z);
395 }
double dh_trans(const double &z)
Computes the transverse comoving distance by integration.
Definition: cosmo.cc:342

◆ dh_trans()

double dh_trans ( const double &  z)

Computes the transverse comoving distance by integration.

Definition at line 342 of file cosmo.cc.

343 {
344  // computes the comoving transverse distance by integration
345  // for the cosmology defined by the global gCOSMO_OMEGA parameters.
346  // units = h^-1 Mpc
347 
348  // Input:
349  // z redshift
350  // Output:
351  // dh_trans comoving transverse distance [h^-1 Mpc]
352 
353  // do again computation to save a division in the result
354  gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
355  double eta, error;
356  gsl_function F;
357  F.function = &H0_over_H;
358 
359  gsl_integration_qags(&F, 0, z, 0, gSIM_EPS, 1000, w, &eta, &error);
360  gsl_integration_workspace_free(w);
361 
362  double result;
363  if ((abs(gCOSMO_OMEGA0_K)) < 1.e-5) result = HUBBLE_LENGTHxh_Mpc * eta; //ok0=gCOSMO_OMEGA0_K// //Hubble distance//
364  else if (gCOSMO_OMEGA0_K > 0) result = HUBBLE_LENGTHxh_Mpc * sinh(sqrt(gCOSMO_OMEGA0_K) * eta) / sqrt(gCOSMO_OMEGA0_K); //ok0=gCOSMO_OMEGA0_K// //Hubble distance//
365  else result = HUBBLE_LENGTHxh_Mpc * sin(sqrt(-gCOSMO_OMEGA0_K) * eta) / sqrt(-gCOSMO_OMEGA0_K); //ok0=gCOSMO_OMEGA0_K// //Hubble distance//
366  return result;
367 }
double gSIM_EPS
See the user documentation
Definition: params.cc:196
double gCOSMO_OMEGA0_K
See the user documentation
Definition: params.cc:65
#define HUBBLE_LENGTHxh_Mpc
Hubble length in [ ].
Definition: inlines.h:5
double H0_over_H(const double z, void *params)
where is the Hubble rate at redshift .
Definition: cosmo.cc:2300

◆ dNdOmega()

double dNdOmega ( double &  zmin,
double &  zmax,
double const &  Mmin,
double const &  Mmax 
)

Definition at line 1283 of file cosmo.cc.

1284 {
1285  //--- Returns the mean number of DM haloes per steradian, between zmin-zmax and Mmin-Mmax
1286  // INPUTS:
1287  // zmin min redshift of interval
1288  // zmax max redshift of interval
1289  // Mmin min mass of interval [Msol]
1290  // Mmax max mass of interval [Msol]
1291 
1292  // OUTPUTS:
1293  // dN/dOmega [sr^-1].
1294 
1295  double Mhmin = Mmin * gCOSMO_M_to_MH;
1296  double Mhmax = Mmax * gCOSMO_M_to_MH;
1297  double par[2] = {Mhmin, Mhmax};
1298  double res = 0.;
1299  simpson_lin_adapt(dNdOmegadz, zmin, zmax, par, res, gSIM_EPS);
1300  return res;
1301 }
void dNdOmegadz(double &z, double par_Mh[2], double &res)
Definition: cosmo.cc:1331
void simpson_lin_adapt(void fn(double &, double *, double &), double &xmin, double &xmax, double par[], double &s, double const &eps, bool is_verbose=false)
Returns the 1D integral of fn using Simpson rule in linear steps. Adaptative steps stop when the prec...
Definition: integr.cc:14
double gSIM_EPS
See the user documentation
Definition: params.cc:196
double gCOSMO_M_to_MH
conversion.
Definition: params.cc:79

◆ dNdOmegadlnMh()

void dNdOmegadlnMh ( double &  lnMh,
double  par_z[2],
double &  res 
)

Definition at line 1304 of file cosmo.cc.

1305 {
1306  //--- Returns the number of halos per lnMh and solid angle in the redshift
1307  // interval [zmin, zmax].
1308  // Can be integrated with simpson_lin/log_adapt over log mass.
1309 
1310  // DO NOT use for double integration first over z and then over lnMh.
1311  // Instead, it is much smarter (faster) to integrate first over lnMh
1312  // (with dNdOmegadz() over dNdVhdlnMh_integrand_logMh), and then over z
1313  // (dNdOmega() is already designed to be integrated over dNdOmegadz() after lnMh).
1314 
1315  // INPUTS:
1316  // lnMh log mass in units [Omega_m0 h^-1 Msol]
1317  // par_z[0] zmin
1318  // par_z[1] zmax
1319  // OUTPUTS:
1320  // dN/dz/dOmega [sr^-1].
1321 
1322  double par_Mh[1] = {exp(lnMh)};
1323  double zmin = par_z[0];
1324  double zmax = par_z[1];
1325  res = 0.;
1326  simpson_lin_adapt(dNdOmegadlnMhdz, zmin, zmax, par_Mh, res, gSIM_EPS);
1327  return;
1328 }
void dNdOmegadlnMhdz(double &z, double par_Mh[1], double &res)
Definition: cosmo.cc:1354
void simpson_lin_adapt(void fn(double &, double *, double &), double &xmin, double &xmax, double par[], double &s, double const &eps, bool is_verbose=false)
Returns the 1D integral of fn using Simpson rule in linear steps. Adaptative steps stop when the prec...
Definition: integr.cc:14
double gSIM_EPS
See the user documentation
Definition: params.cc:196

◆ dNdOmegadlnMhdz()

void dNdOmegadlnMhdz ( double &  z,
double  par_Mh[1],
double &  res 
)

Definition at line 1354 of file cosmo.cc.

1355 {
1356  //--- Returns the number of halos per steradian, log mass, and redshift.
1357  // Can be integrated with simpson_lin/log_adapt over z.
1358 
1359  // DO NOT use for double integration first over z and then over lnMh.
1360  // Instead, it is much faster to integrate first over lnMh
1361  // (with dNdOmegadz() over dNdVhdlnMh_integrand_logMh), and then over z
1362  // (dNdOmega() is already designed to be integrated over dNdOmegadz() after lnMh).
1363 
1364  // INPUTS:
1365  // z redshift
1366  // par_Mh[0] mass [Omega_m0 h^-1 Msol]
1367  // OUTPUTS:
1368  // dN/dln(Mh)/dz/dOmega = Mh * dN/d(Mh)/dz/dOmega = M dN/dM/dz/dOmega, in units of halos/ steradian.
1369 
1370  res = dNdVhdlnMh_interpol(z, par_Mh[0]); // h^3 Mpc^-3
1371  res *= dVhdzdOmega(z); // # of halos per sr per between z and z+dz
1372  return;
1373 }
double dVhdzdOmega(double z)
Differential comoving volume element in [ ].
Definition: cosmo.cc:1397
double dNdVhdlnMh_interpol(const double &z, const double &Mh)
Definition: cosmo.cc:398

◆ dNdOmegadz()

void dNdOmegadz ( double &  z,
double  par_Mh[2],
double &  res 
)

Definition at line 1331 of file cosmo.cc.

1332 {
1333  //--- Returns the number of halos per redshift and solid angle in the mass
1334  // interval [Mhmin, Mhmax].
1335  // Can be integrated with simpson_lin/log_adapt over z.
1336 
1337  // INPUTS:
1338  // z redshift
1339  // par_Mh[0] Mhmin [Omega_m0 h^-1 Msol]
1340  // par_Mh[1] Mhmax [Omega_m0 h^-1 Msol]
1341  // OUTPUTS:
1342  // dN/dz/dOmega [sr^-1].
1343 
1344  double par_z[1] = {z};
1345  double lnMhmin = log(par_Mh[0]);
1346  double lnMhmax = log(par_Mh[1]);
1347  res = 0.;
1348  simpson_lin_adapt(dNdVhdlnMh_integrand_logMh, lnMhmin, lnMhmax, par_z, res, gSIM_EPS);
1349  res *= dVhdzdOmega(z); // # of halos per sr per between z and z+dz
1350  return;
1351 }
void simpson_lin_adapt(void fn(double &, double *, double &), double &xmin, double &xmax, double par[], double &s, double const &eps, bool is_verbose=false)
Returns the 1D integral of fn using Simpson rule in linear steps. Adaptative steps stop when the prec...
Definition: integr.cc:14
double gSIM_EPS
See the user documentation
Definition: params.cc:196
double dVhdzdOmega(double z)
Differential comoving volume element in [ ].
Definition: cosmo.cc:1397
void dNdVhdlnMh_integrand_logMh(double &logMh, double par_z[1], double &res)
Definition: cosmo.cc:441

◆ dNdVhdlnMh_integrand_logMh()

void dNdVhdlnMh_integrand_logMh ( double &  logMh,
double  par_z[1],
double &  res 
)

Definition at line 441 of file cosmo.cc.

442 {
443  //--- Same as dNdVhdlnMh_interpol(), but written to be integrated with
444  // simpson_lin/log_adapt over logMh.
445 
446  // INPUTS:
447  // logMh log mass [Omega_m0 h^-1 Msol]
448  // par_z[0] z: redshift of the halo hosting the sub-clumps
449  // OUTPUTS:
450  // dn/dln(Mh) = Mh * dn/d(Mh) = M dn/dM, in units of halos/ (h^-1 Mpc)^3 (M in units of Msol).
451 
452  res = dNdVhdlnMh_interpol(par_z[0] /*redshift*/, exp(logMh)); // h^3 Mpc^-3
453  return;
454 }
double dNdVhdlnMh_interpol(const double &z, const double &Mh)
Definition: cosmo.cc:398

◆ dNdVhdlnMh_integrand_Mh()

void dNdVhdlnMh_integrand_Mh ( double &  M,
double  par_z[1],
double &  res 
)

Definition at line 425 of file cosmo.cc.

426 {
427  //--- Same as dNdVhdlnMh_interpol(), but written to be integrated with
428  // simpson_lin/log_adapt over Mh.
429 
430  // INPUTS:
431  // Mh mass [Omega_m0 h^-1 Msol]
432  // par_z[0] z: redshift of the halo hosting the sub-clumps
433  // OUTPUTS:
434  // dn/dln(Mh) = Mh * dn/d(Mh) = M dn/dM, in units of halos/ (h^-1 Mpc)^3 (M in units of Msol).
435 
436  res = dNdVhdlnMh_interpol(par_z[0] /*redshift*/, Mh); // h^3 Mpc^-3
437  return;
438 }
double dNdVhdlnMh_interpol(const double &z, const double &Mh)
Definition: cosmo.cc:398

◆ dNdVhdlnMh_interpol()

double dNdVhdlnMh_interpol ( const double &  z,
const double &  Mh 
)

Definition at line 398 of file cosmo.cc.

399 {
400  // Computes dn/dln(Mh) (the number density per mass bin) for redshift z and mass M.
401  // by interpolating the values stored in gCOSMO_DNDVHDLNMH_Z
402  // Note 1: Mh is in units of Omega_m0 h^-1 Msol
403  // Note 2: dn/dln(Mh) = Mh * dn/d(Mh) = M dn/dM, in units of halos/ (h^-1 Mpc)^3 (M in units of Msol).
404 
405  // INPUTS:
406  // z redshift z
407  // Mh mass in unit [Omega_m0 h^-1 Msol]
408  // OUTPUTS:
409  // dn/dln(Mh) = Mh * dn/d(Mh) = M dn/dM, in units of halos/ (h^-1 Mpc)^3 (M in units of Msol).
410 
411  // return error if base grids are empty (check only one of them; if this one is ok
412  // but the others not, an error will be drawn in the interpolation function):
413  if (gCOSMO_Z_GRID.empty()) {
414  printf("\n====> ERROR: dNdVhdlnMh_interpol() in cosmo.cc");
415  printf("\n gCOSMO_Z_GRID grid is empty. Please fill it before");
416  printf("\n via init_extragal().");
417  printf("\n => abort()\n\n");
418  abort();
419  }
420 
422 }
vector< double > gCOSMO_Z_GRID
Grid of redshift values for global tabulated values.
Definition: params.cc:80
x-log(y) linear interpolation
Definition: params.h:181
vector< vector< double > > gCOSMO_DNDVHDLNMH_Z
Definition: params.cc:85
vector< double > gCOSMO_MH_GRID
Definition: params.cc:81
double interp2D(double const &x_eval, double const &y_eval, const vector< double > &x_base, const vector< double > &y_base, const vector< vector< double > > &z_base, const int card_interp, bool no_error=false)
Definition: misc.cc:499

◆ dNdVhdlnMh_literature()

double dNdVhdlnMh_literature ( const double &  z,
const double &  Mh,
const int  card_mf,
vector< string > &  label_literature 
)

Definition at line 477 of file cosmo.cc.

478 {
479  // Computes dn/dln(Mh) (the number density per mass bin) for redshift z and mass M.
480  // from the curves in published papers.
481  // Note 1: Mh is in units of Omega_m0 h^-1 Msol
482  // Note 2: dn/dln(Mh) = Mh * dn/d(Mh) = M dn/dM, in units of halos/ (h^-1 Mpc)^3 (M in units of Msol).
483 
484  // Input:
485  // z redshift z
486  // Mh_var mass in unit [h^-1 Msol]
487  // card_mf mass function corresponding to a paper
488 
489  int n_points;
490  vector<double> log10_Mh_base; // [Msol/h]
491  vector<double> log10_dNdlnMh_base;
492  char Delta_c_char[50];
493 
494  switch (card_mf) {
496  string label_base = "Rodriguez-Puebla et al. (2016)";
497 
498  string Delta_str;
500  double tmp = delta_x_to_delta_crit(-1, kBRYANNORMAN98, z);
501  sprintf(Delta_c_char, "%d", int(round(delta_crit_to_delta_x(tmp, kRHO_MEAN, z))));
502  Delta_str = "#Delta_{m} = " + string(Delta_c_char);
503  } else {
504  sprintf(Delta_c_char, "%.1f", delta_x_to_delta_crit(-1, kBRYANNORMAN98, z));
505  Delta_str = "#Delta_{c} = " + string(Delta_c_char);
506  }
507  vector<vector<double> > rodriguez16_fig7;
508  if (z <= 0.75) {
509  n_points = 29;
510  double rodriguez16_fig7_arr_z0[29][2] = {
511  {6.69e9, 0.737},
512  {9.85e9, 0.522},
513  {1.41e10, 0.369},
514  {2.20e10, 0.245},
515  {3.47e10, 0.164},
516  {4.61e10, 0.126},
517  {6.48e10, 0.0944},
518  {1.74e11, 0.0392},
519  {2.56e11, 0.0273},
520  {4.66e11, 0.0163},
521  {1.24e12, 0.00695},
522  {2.59e12, 0.00360},
523  {4.03e12, 0.00241},
524  {6.35e12, 0.00161},
525  {1.01e13, 0.00106},
526  {1.62e13, 6.83e-4},
527  {2.59e13, 4.33e-4},
528  {4.11e13, 2.72e-4},
529  {6.62e13, 1.63e-4},
530  {9.95e13, 9.91e-5},
531  {1.58e14, 5.46e-5},
532  {2.57e14, 2.66e-5},
533  {3.24e14, 1.79e-5},
534  {4.16e14, 1.18e-5},
535  {5.16e14, 7.49e-6},
536  {6.35e14, 4.78e-6},
537  {7.72e14, 3.01e-6},
538  {9.21e14, 1.93e-6},
539  {1.15e15, 1.03e-6},
540  };
541  for (int i = 0; i < n_points; ++i) {
542  pushback_array2vector(rodriguez16_fig7, &rodriguez16_fig7_arr_z0[i][0], n_points);
543  log10_Mh_base.push_back(log10(rodriguez16_fig7_arr_z0[i][0]));
544  log10_dNdlnMh_base.push_back(log10(rodriguez16_fig7_arr_z0[i][1]));
545  }
546  static bool is_rodriguez0_added = false;
547  if (!is_rodriguez0_added) {
548  label_literature.push_back("#splitline{" + label_base + ", z = 0, " + Delta_str + "}"
549  "{(#Omega_{#Lambda},#Omega_{m},#Omega_{b},#sigma_{8},h,n_{s}) = (0.693,0.307,0.048,0.829,0.678,0.96)}");
550  is_rodriguez0_added = true;
551  }
552  } else if (z > 0.75 and z <= 1.5) {
553  n_points = 22;
554  double rodriguez16_fig7_arr_z1[22][2] = {
555  {1.72e10, 0.349},
556  {6.65e10, 0.104},
557  {1.03e11, 0.0672},
558  {1.47e11, 0.0489},
559  {2.28e11, 0.0325},
560  {4.18e11, 0.0186},
561  {7.49e11, 0.0107},
562  {1.47e12, 0.00564},
563  {2.65e12, 0.00309},
564  {4.22e12, 0.00189},
565  {6.58e12, 0.00117},
566  {1.08e13, 6.47e-4},
567  {1.61e13, 3.90e-4},
568  {2.55e13, 2.09e-4},
569  {4.16e13, 9.91e-5},
570  {6.55e13, 4.44e-5},
571  {9.35e13, 2.19e-5},
572  {1.27e14, 1.10e-5},
573  {1.66e14, 5.60e-6},
574  {2.09e14, 2.88e-6},
575  {2.43e14, 1.84e-6},
576  {2.88e14, 1.02e-6},
577  };
578  for (int i = 0; i < n_points; ++i) {
579  pushback_array2vector(rodriguez16_fig7, &rodriguez16_fig7_arr_z1[i][0], n_points);
580  log10_Mh_base.push_back(log10(rodriguez16_fig7_arr_z1[i][0]));
581  log10_dNdlnMh_base.push_back(log10(rodriguez16_fig7_arr_z1[i][1]));
582  }
583  static bool is_rodriguez1_added = false;
584  if (!is_rodriguez1_added) {
585  label_literature.push_back(label_base + ", z = 1, " + Delta_str);
586  is_rodriguez1_added = true;
587  }
588  } else if (z > 1.5 and z <= 2.5) {
589  n_points = 20;
590  double rodriguez16_fig7_arr_z2[20][2] = {
591  {6.51e9, 0.927},
592  {1.22e10, 0.500},
593  {2.57e10, 0.242},
594  {5.73e10, 0.113},
595  {1.24e11, 0.0521},
596  {2.43e11, 0.0262},
597  {5.14e11, 0.0116},
598  {1.01e12, 0.00545},
599  {2.01e12, 0.00234},
600  {4.05e12, 8.95e-4},
601  {6.55e12, 4.36e-4},
602  {1.01e13, 2.09e-4},
603  {1.40e13, 1.18e-4},
604  {2.08e13, 5.13e-5},
605  {3.36e13, 1.68e-5},
606  {4.16e13, 9.29e-6},
607  {5.06e13, 5.27e-6},
608  {6.09e13, 2.94e-6},
609  {7.26e13, 1.72e-6},
610  {8.35e13, 1.00e-6},
611  };
612  for (int i = 0; i < n_points; ++i) {
613  pushback_array2vector(rodriguez16_fig7, &rodriguez16_fig7_arr_z2[i][0], n_points);
614  log10_Mh_base.push_back(log10(rodriguez16_fig7_arr_z2[i][0]));
615  log10_dNdlnMh_base.push_back(log10(rodriguez16_fig7_arr_z2[i][1]));
616  }
617  static bool is_rodriguez2_added = false;
618  if (!is_rodriguez2_added) {
619  label_literature.push_back(label_base + ", z = 2, " + Delta_str);
620  is_rodriguez2_added = true;
621  }
622  } else if (z > 2.5 and z <= 3.5) {
623  n_points = 21;
624  double rodriguez16_fig7_arr_z3[21][2] = {
625  {6.55e9, 0.865},
626  {9.02e9, 0.616},
627  {1.51e10, 0.367},
628  {2.57e10, 0.209},
629  {4.29e10, 0.124},
630  {9.65e10, 0.0506},
631  {1.59e11, 0.0279},
632  {2.60e11, 0.0156},
633  {4.11e11, 0.00855},
634  {8.05e11, 0.00343},
635  {1.29e12, 0.00168},
636  {2.02e12, 8.12e-4},
637  {3.24e12, 3.57e-4},
638  {5.03e12, 1.45e-4},
639  {6.97e12, 7.21e-5},
640  {1.01e13, 2.91e-5},
641  {1.22e13, 1.70e-5},
642  {1.59e13, 8.14e-6},
643  {2.02e13, 3.91e-6},
644  {2.51e13, 1.81e-6},
645  {2.96e13, 9.93e-7},
646  };
647  for (int i = 0; i < n_points; ++i) {
648  pushback_array2vector(rodriguez16_fig7, &rodriguez16_fig7_arr_z3[i][0], n_points);
649  log10_Mh_base.push_back(log10(rodriguez16_fig7_arr_z3[i][0]));
650  log10_dNdlnMh_base.push_back(log10(rodriguez16_fig7_arr_z3[i][1]));
651  }
652  static bool is_rodriguez3_added = false;
653  if (!is_rodriguez3_added) {
654  label_literature.push_back(label_base + ", z = 3, " + Delta_str);
655  is_rodriguez3_added = true;
656  }
657  } else if (z > 3.5 and z <= 4.5) {
658  n_points = 17;
659  double rodriguez16_fig7_arr_z4[17][2] = {
660  {6.41e9, 0.763},
661  {1.85e10, 0.240},
662  {3.36e10, 0.117},
663  {6.22e10, 0.0554},
664  {1.22e11, 0.0228},
665  {2.46e11, 0.00873},
666  {4.18e11, 0.00388},
667  {8.01e11, 0.00131},
668  {1.31e12, 5.22e-4},
669  {2.14e12, 1.86e-4},
670  {3.08e12, 7.94e-5},
671  {4.05e12, 3.97e-5},
672  {5.35e12, 1.86e-5},
673  {6.97e12, 8.43e-6},
674  {8.83e12, 3.83e-6},
675  {1.09e13, 1.89e-6},
676  {1.28e13, 1.01e-6},
677  };
678  for (int i = 0; i < n_points; ++i) {
679  pushback_array2vector(rodriguez16_fig7, &rodriguez16_fig7_arr_z4[i][0], n_points);
680  log10_Mh_base.push_back(log10(rodriguez16_fig7_arr_z4[i][0]));
681  log10_dNdlnMh_base.push_back(log10(rodriguez16_fig7_arr_z4[i][1]));
682  }
683  static bool is_rodriguez4_added = false;
684  if (!is_rodriguez4_added) {
685  label_literature.push_back(label_base + ", z = 4, " + Delta_str);
686  is_rodriguez4_added = true;
687  }
688  } else if (z > 4.5 and z <= 5.5) {
689  n_points = 17;
690  double rodriguez16_fig7_arr_z5[17][2] = {
691  {6.41e9, 0.544},
692  {1.08e10, 0.286},
693  {2.01e10, 0.130},
694  {4.03e10, 0.0506},
695  {6.28e10, 0.0262},
696  {9.85e10, 0.0131},
697  {1.56e11, 0.00622},
698  {2.53e11, 0.00267},
699  {3.93e11, 0.00116},
700  {6.19e11, 4.48e-4},
701  {9.50e11, 1.69e-4},
702  {1.50e12, 5.28e-5},
703  {1.95e12, 2.50e-5},
704  {2.44e12, 1.26e-5},
705  {3.16e12, 5.57e-6},
706  {4.07e12, 2.29e-6},
707  {5.01e12, 1.03e-6},
708  };
709  for (int i = 0; i < n_points; ++i) {
710  pushback_array2vector(rodriguez16_fig7, &rodriguez16_fig7_arr_z5[i][0], n_points);
711  log10_Mh_base.push_back(log10(rodriguez16_fig7_arr_z5[i][0]));
712  log10_dNdlnMh_base.push_back(log10(rodriguez16_fig7_arr_z5[i][1]));
713  }
714  static bool is_rodriguez5_added = false;
715  if (!is_rodriguez5_added) {
716  label_literature.push_back(label_base + ", z = 5, " + Delta_str);
717  is_rodriguez5_added = true;
718  }
719  } else if (z > 5.5 and z <= 6.5) {
720  n_points = 13;
721  double rodriguez16_fig7_arr_z6[13][2] = {
722  {6.55e9, 0.333},
723  {1.60e10, 0.0918},
724  {4.05e10, 0.0219},
725  {6.32e10, 0.0101},
726  {9.85e10, 0.00455},
727  {1.58e11, 0.00178},
728  {3.00e11, 4.39e-4},
729  {4.71e11, 1.47e-4},
730  {7.53e11, 4.06e-5},
731  {1.02e12, 1.72e-5},
732  {1.33e12, 6.85e-6},
733  {1.90e12, 1.97e-6},
734  {2.25e12, 1.01e-6},
735 
736  };
737  for (int i = 0; i < n_points; ++i) {
738  pushback_array2vector(rodriguez16_fig7, &rodriguez16_fig7_arr_z6[i][0], n_points);
739  log10_Mh_base.push_back(log10(rodriguez16_fig7_arr_z6[i][0]));
740  log10_dNdlnMh_base.push_back(log10(rodriguez16_fig7_arr_z6[i][1]));
741  }
742  static bool is_rodriguez6_added = false;
743  if (!is_rodriguez6_added) {
744  label_literature.push_back(label_base + ", z = 6, " + Delta_str);
745  is_rodriguez6_added = true;
746  }
747  } else if (z > 6.5 and z <= 7.5) {
748  n_points = 12;
749  double rodriguez16_fig7_arr_z7[12][2] = {
750  {6.48e9, 0.191},
751  {1.42e10, 0.0554},
752  {2.64e10, 0.0193},
753  {4.24e10, 0.00814},
754  {7.97e10, 0.00231},
755  {1.53e11, 5.36e-4},
756  {2.68e11, 1.29e-4},
757  {3.89e11, 4.50e-5},
758  {5.35e11, 1.67e-5},
759  {7.11e11, 6.57e-6},
760  {9.40e11, 2.51e-6},
761  {1.19e12, 1.02e-6},
762  };
763  for (int i = 0; i < n_points; ++i) {
764  pushback_array2vector(rodriguez16_fig7, &rodriguez16_fig7_arr_z7[i][0], n_points);
765  log10_Mh_base.push_back(log10(rodriguez16_fig7_arr_z7[i][0]));
766  log10_dNdlnMh_base.push_back(log10(rodriguez16_fig7_arr_z7[i][1]));
767  }
768  static bool is_rodriguez7_added = false;
769  if (!is_rodriguez7_added) {
770  label_literature.push_back(label_base + ", z = 7, " + Delta_str);
771  is_rodriguez7_added = true;
772  }
773  } else if (z > 7.5 and z <= 8.5) {
774  n_points = 14;
775  double rodriguez16_fig7_arr_z8[14][2] = {
776  {6.45e9, 0.0893},
777  {1.24e10, 0.0277},
778  {2.02e10, 0.0111},
779  {3.31e10, 0.00406},
780  {4.96e10, 0.00168},
781  {6.28e10, 9.59e-4},
782  {9.55e10, 3.49e-4},
783  {1.53e11, 9.57e-5},
784  {1.96e11, 4.69e-5},
785  {2.55e11, 2.16e-5},
786  {3.20e11, 9.95e-6},
787  {4.22e11, 3.80e-6},
788  {5.22e11, 1.73e-6},
789  {6.00e11, 1.03e-6},
790  };
791  for (int i = 0; i < n_points; ++i) {
792  pushback_array2vector(rodriguez16_fig7, &rodriguez16_fig7_arr_z8[i][0], n_points);
793  log10_Mh_base.push_back(log10(rodriguez16_fig7_arr_z8[i][0]));
794  log10_dNdlnMh_base.push_back(log10(rodriguez16_fig7_arr_z8[i][1]));
795  }
796  static bool is_rodriguez8_added = false;
797  if (!is_rodriguez8_added) {
798  label_literature.push_back(label_base + ", z = 8, " + Delta_str);
799  is_rodriguez8_added = true;
800  }
801  } else if (z > 8.5 and z <= 9.5) {
802  n_points = 12;
803  double rodriguez16_fig7_arr_z9[12][2] = {
804  {6.45e9, 0.0301},
805  {1.03e10, 0.0118},
806  {1.72e10, 0.00402},
807  {2.59e10, 0.00154},
808  {3.97e10, 5.55e-4},
809  {6.09e10, 1.72e-4},
810  {9.95e10, 4.14e-5},
811  {1.30e11, 1.79e-5},
812  {1.60e11, 8.67e-6},
813  {1.99e11, 4.08e-6},
814  {2.48e11, 1.83e-6},
815  {2.90e11, 9.93e-7},
816  };
817  for (int i = 0; i < n_points; ++i) {
818  pushback_array2vector(rodriguez16_fig7, &rodriguez16_fig7_arr_z9[i][0], n_points);
819  log10_Mh_base.push_back(log10(rodriguez16_fig7_arr_z9[i][0]));
820  log10_dNdlnMh_base.push_back(log10(rodriguez16_fig7_arr_z9[i][1]));
821  }
822  static bool is_rodriguez9_added = false;
823  if (!is_rodriguez9_added) {
824  label_literature.push_back(label_base + ", z = 9, " + Delta_str);
825  is_rodriguez9_added = true;
826  }
827  } else return 0;
828  double log10_res = interp1D(log10(Mh_var), log10_Mh_base, log10_dNdlnMh_base, kSPLINE, true);
829  return pow(10, log10_res) / log(10);
830  }
831  case kTINKER08_N:
832  case kTINKER08: {
833  if (z > 0) {
834  return 0;
835  }
836 
838  double Delta_m0 = delta_crit_to_delta_x(Delta_c0, kRHO_MEAN, 0);
839  sprintf(Delta_c_char, "%d", int(round(delta_x_to_delta_crit(200, kRHO_MEAN, z))));
840  string Delta_m0_200_str = string(Delta_c_char);
841  sprintf(Delta_c_char, "%d", int(round(delta_x_to_delta_crit(800, kRHO_MEAN, z))));
842  string Delta_m0_800_str = string(Delta_c_char);
843  sprintf(Delta_c_char, "%d", int(round(delta_x_to_delta_crit(3200, kRHO_MEAN, z))));
844  string Delta_m0_3200_str = string(Delta_c_char);
845 
846  double tinker08_fig5_arr[22][2];
847  string Delta0_str;
848  n_points = 22;
849  if (Delta_m0 < 700) {
850  double tmp[22][2] = {
851  {10.01, -1.55 },
852  {10.51, -1.50 },
853  {11.05, -1.43 },
854  {11.59, -1.37 },
855  {12.11, -1.30 },
856  {12.63, -1.23 },
857  {13.08, -1.19 },
858  {13.46, -1.16 },
859  {13.79, -1.17 },
860  {14.08, -1.21 },
861  {14.36, -1.29 },
862  {14.55, -1.39 },
863  {14.77, -1.55 },
864  {14.93, -1.75 },
865  {15.06, -1.93 },
866  {15.18, -2.16 },
867  {15.28, -2.40 },
868  {15.36, -2.61 },
869  {15.43, -2.83 },
870  {15.50, -3.09 },
871  {15.56, -3.35 },
872  {15.61, -3.60 }
873  };
874  memcpy(tinker08_fig5_arr, tmp, sizeof(tinker08_fig5_arr));
875  if (gCOSMO_FLAG_DELTA_REF == kRHO_MEAN) Delta0_str = "#Delta_{m} = 200";
876  else Delta0_str = "#Delta_{c} = " + Delta_m0_200_str;
877  } else if (Delta_m0 < 2000) {
878  double tmp[22][2] = {
879  {10.02, -1.53 },
880  {10.29, -1.51 },
881  {10.63, -1.49 },
882  {10.96, -1.46 },
883  {11.44, -1.43 },
884  {11.78, -1.40 },
885  {12.22, -1.36 },
886  {12.72, -1.33 },
887  {13.07, -1.31 },
888  {13.49, -1.31 },
889  {13.93, -1.38 },
890  {14.22, -1.47 },
891  {14.44, -1.61 },
892  {14.66, -1.81 },
893  {14.80, -1.98 },
894  {14.94, -2.21 },
895  {15.07, -2.47 },
896  {15.16, -2.70 },
897  {15.23, -2.93 },
898  {15.30, -3.16 },
899  {15.36, -3.37 },
900  {15.41, -3.59 }
901  };
902  memcpy(tinker08_fig5_arr, tmp, sizeof(tinker08_fig5_arr));
903  if (gCOSMO_FLAG_DELTA_REF == kRHO_MEAN) Delta0_str = "#Delta_{m} = 800";
904  else Delta0_str = "#Delta_{c} = " + Delta_m0_800_str;
905  } else {
906  double tmp[22][2] = {
907  {10.01, -1.57 },
908  {10.55, -1.55 },
909  {10.98, -1.53 },
910  {11.51, -1.51 },
911  {11.85, -1.50 },
912  {12.30, -1.49 },
913  {12.72, -1.49 },
914  {13.01, -1.50 },
915  {13.35, -1.53 },
916  {13.64, -1.59 },
917  {13.85, -1.66 },
918  {14.07, -1.76 },
919  {14.24, -1.89 },
920  {14.40, -2.03 },
921  {14.53, -2.20 },
922  {14.64, -2.37 },
923  {14.71, -2.50 },
924  {14.79, -2.68 },
925  {14.87, -2.86 },
926  {14.95, -3.09 },
927  {15.03, -3.33 },
928  {15.10, -3.60 }
929  };
930  memcpy(tinker08_fig5_arr, tmp, sizeof(tinker08_fig5_arr));
931  if (gCOSMO_FLAG_DELTA_REF == kRHO_MEAN) Delta0_str = "#Delta_{m} = 3200";
932  else Delta0_str = "#Delta_{c} = " + Delta_m0_3200_str;
933  }
934  vector<vector<double> > tinker08_fig5;
935  for (int i = 0; i < n_points; ++i) {
936  pushback_array2vector(tinker08_fig5, &tinker08_fig5_arr[i][0], n_points);
937  log10_Mh_base.push_back(tinker08_fig5[i][0]);
938  log10_dNdlnMh_base.push_back(tinker08_fig5[i][1]);
939  }
940  double log10_res = interp1D(log10(Mh_var), log10_Mh_base, log10_dNdlnMh_base, kSPLINE, true);
941  static bool is_tinker_added = false;
942  if (!is_tinker_added) {
943  label_literature.push_back("#splitline{Tinker et al. (2008), WMAP1, z = 0, " + Delta0_str + "}"
944  "{(#Omega_{#Lambda},#Omega_{m},#Omega_{b},#sigma_{8},h,n_{s}) = (0.7,0.3,0.04,0.9,0.7,1)}");
945  is_tinker_added = true;
946  }
947  // rescale and return values:
948  return pow(10, log10_res) * RHO_CRITperh2_MSOLperKPC3 / KPC3_to_MPC3 / Mh_var * gCOSMO_OMEGA0_M;
949  }
950  case kTINKER10:
951  break;
952  case kBOCQUET16_DMONLY:
953  case kBOCQUET16_HYDRO: {
954  string label_base = "Bocquet et al. (2016, hydro)";
955  sprintf(Delta_c_char, "%d", int(round(delta_x_to_delta_crit(200, kRHO_MEAN, z))));
956  string Delta0_str;
957  if (gCOSMO_FLAG_DELTA_REF == kRHO_MEAN) Delta0_str = "#Delta_{m} = 200";
958  else Delta0_str = "#Delta_{c} = " + string(Delta_c_char);
959  vector<vector<double> > bocquet16_fig2;
960  if (z <= 0.3) {
961  n_points = 15;
962  double bocquet16_fig2_arr_z0[15][2] = {
963  {6.39e11, 0.00210},
964  {2.14e12, 7.37e-4},
965  {8.55e12, 2.22e-4},
966  {2.41e13, 8.11e-5},
967  {8.84e13, 2.19e-5},
968  {2.24e14, 7.03e-6},
969  {4.29e14, 2.57e-6},
970  {7.30e14, 9.22e-7},
971  {1.22e15, 2.77e-7},
972  {1.69e15, 1.04e-7},
973  {2.57e15, 2.35e-8},
974  {3.57e15, 5.56e-9},
975  {4.79e15, 1.18e-9},
976  {6.03e15, 2.92e-10},
977  {7.92e15, 3.74e-11},
978  };
979  for (int i = 0; i < n_points; ++i) {
980  pushback_array2vector(bocquet16_fig2, &bocquet16_fig2_arr_z0[i][0], n_points);
981  log10_Mh_base.push_back(log10(bocquet16_fig2_arr_z0[i][0] * gCOSMO_HUBBLE));
982  log10_dNdlnMh_base.push_back(log10(bocquet16_fig2_arr_z0[i][1]));
983  }
984  static bool is_bocqet1_added = false;
985  if (!is_bocqet1_added) {
986  label_literature.push_back("#splitline{" + label_base + ", z = 0, " + Delta0_str + "}"
987  "{(#Omega_{#Lambda},#Omega_{m},#Omega_{b},#sigma_{8},h,n_{s}) = (0.728,0.272,0.0456,0.809,0.704,0.963)}");
988  is_bocqet1_added = true;
989  }
990  } else if (z > 0.3 and z <= 0.8) {
991  n_points = 17;
992  double bocquet16_fig2_arr_z0_5[17][2] = {
993  {6.62e11, 0.00201},
994  {1.44e12, 9.78e-4},
995  {2.87e12, 5.31e-4},
996  {6.65e12, 2.37e-4},
997  {1.09e13, 1.46e-4},
998  {2.32e13, 6.60e-5},
999  {4.33e13, 3.10e-5},
1000  {9.22e13, 1.16e-5},
1001  {1.44e14, 5.65e-6},
1002  {3.11e14, 1.37e-6},
1003  {5.76e14, 2.90e-7},
1004  {8.94e14, 7.48e-8},
1005  {1.40e15, 1.22e-8},
1006  {1.89e15, 2.83e-9},
1007  {2.46e15, 6.69e-10},
1008  {3.13e15, 1.42e-10},
1009  {4.02e15, 2.03e-11},
1010  };
1011  for (int i = 0; i < n_points; ++i) {
1012  pushback_array2vector(bocquet16_fig2, &bocquet16_fig2_arr_z0_5[i][0], n_points);
1013  log10_Mh_base.push_back(log10(bocquet16_fig2_arr_z0_5[i][0] * gCOSMO_HUBBLE));
1014  log10_dNdlnMh_base.push_back(log10(bocquet16_fig2_arr_z0_5[i][1]));
1015  }
1016  static bool is_bocqet2_added = false;
1017  if (!is_bocqet2_added) {
1018  label_literature.push_back(label_base + ", z = 0.5, " + Delta0_str);
1019  is_bocqet2_added = true;
1020  }
1021  } else if (z > 0.8 and z <= 1.6) {
1022  n_points = 16;
1023  double bocquet16_fig2_arr_z1_2[16][2] = {
1024  {7.00e11, 0.00184},
1025  {3.08e12, 4.00e-4},
1026  {7.70e12, 1.40e-4},
1027  {1.81e13, 4.40e-5},
1028  {3.64e13, 1.54e-5},
1029  {6.64e13, 5.18e-6},
1030  {1.14e14, 1.59e-6},
1031  {1.87e14, 4.39e-7},
1032  {2.72e14, 1.44e-7},
1033  {3.60e14, 5.39e-8},
1034  {5.00e14, 1.36e-8},
1035  {6.81e14, 3.15e-9},
1036  {8.88e14, 7.30e-10},
1037  {1.08e15, 2.20e-10},
1038  {1.29e15, 6.32e-11},
1039  {1.50e15, 1.94e-11},
1040  };
1041  for (int i = 0; i < n_points; ++i) {
1042  pushback_array2vector(bocquet16_fig2, &bocquet16_fig2_arr_z1_2[i][0], n_points);
1043  log10_Mh_base.push_back(log10(bocquet16_fig2_arr_z1_2[i][0] * gCOSMO_HUBBLE));
1044  log10_dNdlnMh_base.push_back(log10(bocquet16_fig2_arr_z1_2[i][1]));
1045  }
1046  static bool is_bocqet3_added = false;
1047  if (!is_bocqet3_added) {
1048  label_literature.push_back(label_base + ", z = 1.2, " + Delta0_str);
1049  is_bocqet3_added = true;
1050  }
1051  } else if (z > 1.6 and z <= 2.5) {
1052  n_points = 14;
1053  double bocquet16_fig2_arr_z2[14][2] = {
1054  {6.90e11, 0.00142},
1055  {1.62e12, 5.19e-4},
1056  {3.13e12, 2.22e-4},
1057  {5.58e12, 9.66e-5},
1058  {9.24e12, 4.40e-5},
1059  {1.66e13, 1.58e-5},
1060  {2.91e13, 5.07e-6},
1061  {5.81e13, 8.82e-7},
1062  {1.01e14, 1.54e-7},
1063  {1.57e14, 2.86e-8},
1064  {2.35e14, 5.10e-9},
1065  {3.27e14, 8.50e-10},
1066  {4.41e14, 1.42e-10},
1067  {5.72e14, 2.03e-11},
1068  };
1069  for (int i = 0; i < n_points; ++i) {
1070  pushback_array2vector(bocquet16_fig2, &bocquet16_fig2_arr_z2[i][0], n_points);
1071  log10_Mh_base.push_back(log10(bocquet16_fig2_arr_z2[i][0] * gCOSMO_HUBBLE));
1072  log10_dNdlnMh_base.push_back(log10(bocquet16_fig2_arr_z2[i][1]));
1073  }
1074  static bool is_bocqet4_added = false;
1075  if (!is_bocqet4_added) {
1076  label_literature.push_back(label_base + ", z = 2, " + Delta0_str);
1077  is_bocqet4_added = true;
1078  }
1079  } else return 0;
1080 
1081  double log10_res = interp1D(log10(Mh_var), log10_Mh_base, log10_dNdlnMh_base, kSPLINE, true);
1082  return pow(10, log10_res) / pow(gCOSMO_HUBBLE, 3);
1083  }
1084  case kSHETHTORMEN99:
1085  break;
1086  case kPRESSSCHECHTER74:
1087  break;
1088  case kJENKINS01: {
1089  vector<vector<double> > springel05_fig2;
1090  string label_base = "Springel et al. (2005)";
1091  string Delta0_str;
1092  sprintf(Delta_c_char, "%d", int(round(delta_x_to_delta_crit(180, kRHO_MEAN, z))));
1093  if (gCOSMO_FLAG_DELTA_REF == kRHO_MEAN) Delta0_str = "#Delta_{m} = 180";
1094  else Delta0_str = "#Delta_{c} = " + string(Delta_c_char);
1095  if (z <= 1.) {
1096  n_points = 27;
1097  double springel05_fig2_arr_z0[27][2] = {
1098  {10.00, -1.64},
1099  {10.40, -1.55},
1100  {10.72, -1.49},
1101  {11.13, -1.43},
1102  {11.52, -1.38},
1103  {11.99, -1.34},
1104  {12.44, -1.30},
1105  {12.93, -1.27},
1106  {13.36, -1.23},
1107  {13.77, -1.23},
1108  {14.05, -1.25},
1109  {14.30, -1.31},
1110  {14.46, -1.39},
1111  {14.64, -1.49},
1112  {14.79, -1.62},
1113  {14.94, -1.79},
1114  {15.06, -1.98},
1115  {15.16, -2.16},
1116  {15.26, -2.39},
1117  {15.43, -2.85},
1118  {15.55, -3.25},
1119  {15.65, -3.65},
1120  {15.80, -4.28},
1121  {15.84, -4.57},
1122  {15.92, -5.04},
1123  {15.97, -5.37},
1124  {16.00, -5.52},
1125  };
1126  for (int i = 0; i < n_points; ++i) {
1127  pushback_array2vector(springel05_fig2, &springel05_fig2_arr_z0[i][0], n_points);
1128  log10_Mh_base.push_back(springel05_fig2_arr_z0[i][0]);
1129  log10_dNdlnMh_base.push_back(springel05_fig2_arr_z0[i][1]);
1130  }
1131  static bool is_springel1_added = false;
1132  if (!is_springel1_added) {
1133  label_literature.push_back("#splitline{" + label_base + ", z = 0, " + Delta0_str + "}"
1134  "{(#Omega_{#Lambda},#Omega_{m},#Omega_{b},#sigma_{8},h,n_{s}) = (0.75,0.25,0.045,0.9,0.73,1)}");
1135  is_springel1_added = true;
1136  }
1137  } else if (z > 1. and z <= 2.) {
1138  n_points = 24;
1139  double springel05_fig2_arr_z1_5[24][2] = {
1140  {10.00, -1.48},
1141  {10.33, -1.45},
1142  {10.76, -1.43},
1143  {11.11, -1.40},
1144  {11.56, -1.37},
1145  {11.97, -1.37},
1146  {12.32, -1.37},
1147  {12.68, -1.40},
1148  {12.97, -1.46},
1149  {13.20, -1.55},
1150  {13.45, -1.70},
1151  {13.62, -1.83},
1152  {13.79, -2.00},
1153  {13.93, -2.19},
1154  {14.15, -2.55},
1155  {14.29, -2.86},
1156  {14.43, -3.17},
1157  {14.57, -3.60},
1158  {14.69, -4.04},
1159  {14.77, -4.35},
1160  {14.88, -4.84},
1161  {14.96, -5.24},
1162  {15.02, -5.55},
1163  {15.05, -5.74},
1164  };
1165  for (int i = 0; i < n_points; ++i) {
1166  pushback_array2vector(springel05_fig2, &springel05_fig2_arr_z1_5[i][0], n_points);
1167  log10_Mh_base.push_back(springel05_fig2_arr_z1_5[i][0]);
1168  log10_dNdlnMh_base.push_back(springel05_fig2_arr_z1_5[i][1]);
1169  }
1170  static bool is_springel2_added = false;
1171  if (!is_springel2_added) {
1172  label_literature.push_back(label_base + ", z = 1.5, " + Delta0_str);
1173  is_springel2_added = true;
1174  }
1175  } else if (z > 2. and z <= 4.) {
1176  n_points = 26;
1177  double springel05_fig2_arr_z3_06[26][2] = {
1178  {10.00, -1.46},
1179  {10.44, -1.47},
1180  {10.79, -1.47},
1181  {11.13, -1.49},
1182  {11.49, -1.52},
1183  {11.71, -1.58},
1184  {11.93, -1.63},
1185  {12.08, -1.70},
1186  {12.28, -1.78},
1187  {12.43, -1.87},
1188  {12.59, -1.99},
1189  {12.76, -2.15},
1190  {12.90, -2.29},
1191  {13.07, -2.51},
1192  {13.15, -2.61},
1193  {13.28, -2.86},
1194  {13.37, -3.03},
1195  {13.47, -3.22},
1196  {13.58, -3.46},
1197  {13.66, -3.69},
1198  {13.77, -3.99},
1199  {13.84, -4.21},
1200  {13.93, -4.51},
1201  {14.03, -4.87},
1202  {14.11, -5.24},
1203  {14.22, -5.70},
1204  };
1205  for (int i = 0; i < n_points; ++i) {
1206  pushback_array2vector(springel05_fig2, &springel05_fig2_arr_z3_06[i][0], n_points);
1207  log10_Mh_base.push_back(springel05_fig2_arr_z3_06[i][0]);
1208  log10_dNdlnMh_base.push_back(springel05_fig2_arr_z3_06[i][1]);
1209  }
1210  static bool is_springel3_added = false;
1211  if (!is_springel3_added) {
1212  label_literature.push_back(label_base + ", z = 3.06, " + Delta0_str);
1213  is_springel3_added = true;
1214  }
1215  } else if (z > 4. and z <= 8.) {
1216  n_points = 15;
1217  double springel05_fig2_arr_z5_72[15][2] = {
1218  {10.00, -1.72},
1219  {10.20, -1.75},
1220  {10.44, -1.81},
1221  {10.72, -1.91},
1222  {10.98, -2.04},
1223  {11.25, -2.23},
1224  {11.52, -2.44},
1225  {11.76, -2.67},
1226  {11.90, -2.86},
1227  {12.08, -3.11},
1228  {12.23, -3.36},
1229  {12.45, -3.75},
1230  {12.70, -4.35},
1231  {12.93, -4.99},
1232  {13.12, -5.69},
1233  };
1234  for (int i = 0; i < n_points; ++i) {
1235  pushback_array2vector(springel05_fig2, &springel05_fig2_arr_z5_72[i][0], n_points);
1236  log10_Mh_base.push_back(springel05_fig2_arr_z5_72[i][0]);
1237  log10_dNdlnMh_base.push_back(springel05_fig2_arr_z5_72[i][1]);
1238  }
1239  static bool is_springel4_added = false;
1240  if (!is_springel4_added) {
1241  label_literature.push_back(label_base + ", z = 5.72, " + Delta0_str);
1242  is_springel4_added = true;
1243  }
1244  } else if (z > 8. and z <= 15.) {
1245  n_points = 10;
1246  double springel05_fig2_arr_z10_07[10][2] = {
1247  {10.00, -2.82},
1248  {10.19, -2.98},
1249  {10.41, -3.23},
1250  {10.63, -3.50},
1251  {10.80, -3.71},
1252  {10.95, -3.98},
1253  {11.17, -4.35},
1254  {11.34, -4.69},
1255  {11.51, -5.11},
1256  {11.79, -5.66},
1257  };
1258  for (int i = 0; i < n_points; ++i) {
1259  pushback_array2vector(springel05_fig2, &springel05_fig2_arr_z10_07[i][0], n_points);
1260  log10_Mh_base.push_back(springel05_fig2_arr_z10_07[i][0]);
1261  log10_dNdlnMh_base.push_back(springel05_fig2_arr_z10_07[i][1]);
1262  }
1263  static bool is_springel4_added = false;
1264  if (!is_springel4_added) {
1265  label_literature.push_back(label_base + ", z = 10.07, " + Delta0_str);
1266  is_springel4_added = true;
1267  }
1268  } else return 0;
1269 
1270  double log10_res = interp1D(log10(Mh_var), log10_Mh_base, log10_dNdlnMh_base, kSPLINE, true);
1271  return pow(10, log10_res) * RHO_CRITperh2_MSOLperKPC3 / KPC3_to_MPC3 / Mh_var * gCOSMO_OMEGA0_M;
1272  }
1273  default :
1274  printf("\n====> ERROR: dNdVhdlnMh_literature() in cosmo.cc");
1275  printf("\n card_mf = %d not a valid mass function!", card_mf);
1276  printf("\n => abort()\n\n");
1277  abort();
1278  }
1279  return 0;
1280 }
See user documentation
Definition: params.h:120
void pushback_array2vector(vector< vector< T > > &vec, const T *array, const int size)
Definition: misc.h:104
double gCOSMO_HUBBLE
See the user documentation
Definition: params.cc:62
See user documentation
Definition: params.h:116
double gCOSMO_DELTA0
See the user documentation
Definition: params.cc:71
int gCOSMO_FLAG_DELTA_REF
See the user documentation
Definition: params.cc:72
Spline interpolation.
Definition: params.h:184
#define KPC3_to_MPC3
Definition: inlines.h:23
See user documentation
Definition: params.h:119
double delta_crit_to_delta_x(double const &Delta_c, int const card_delta_ref, double const &z)
Definition: clumps.cc:182
See user documentation
Definition: params.h:121
See user documentation
Definition: params.h:123
See user documentation
Definition: params.h:117
double delta_x_to_delta_crit(double const &Delta_x0, int const card_delta_ref, double const &z)
Definition: clumps.cc:102
double round(const double &x, const int digits)
Definition: misc.cc:1570
See user documentation
Definition: params.h:175
See user documentation
Definition: params.h:174
See user documentation
Definition: params.h:124
See user documentation
Definition: params.h:118
See user documentation
Definition: params.h:122
double gCOSMO_OMEGA0_M
See the user documentation
Definition: params.cc:63
double interp1D(double const &x_eval, const vector< double > &x_base, const vector< double > &y_base, const int card_interp, bool no_error=false)
Definition: misc.cc:433
#define RHO_CRITperh2_MSOLperKPC3
Critical density of the universe in units of .
Definition: inlines.h:7

◆ dNdVhdlnMh_sigmoid()

void dNdVhdlnMh_sigmoid ( double &  Mh,
double  par[3],
double &  res 
)

Definition at line 457 of file cosmo.cc.

458 {
459  //--- Returns the mass function truncated at low masses with a sigmoid function.
460  // Can be integrated with simpson_lin/log_adapt over Mh.
461 
462  // INPUTS:
463  // M mass [Omega_m0 h^-1 Msol]
464  // par[0] z: redshift of the halo hosting the sub-clumps
465  // par[1] Mhmin [Omega_m0 h^-1 Msol], where the mass function is cut off.
466  // par[2] sigma [Omega_m0 h^-1 Msol], the sharpness of the cut-off.
467  // OUTPUTS:
468  // dn/dln(Mh) = Mh * dn/d(Mh) = M dn/dM, in units of halos/ (h^-1 Mpc)^3 (M in units of Msol).
469 
470  res = sigmoid_window(Mh, par[1] /*Mhmin*/, par[2] /*sigma*/) * dNdVhdlnMh_interpol(par[0] /*redshift*/, Mh); // h^3 Mpc^-3
471  return;
472 }
double dNdVhdlnMh_interpol(const double &z, const double &Mh)
Definition: cosmo.cc:398
double sigmoid_window(double const &x, double const &x0, double const &sigma)
Definition: misc.cc:1580

◆ ds2dlnk()

double ds2dlnk ( double  lnk,
void *  p 
)

Returns \( d(\sigma^2)/dk \propto k^3\;p(k)\; W(kr)\), where r is the scale of interest and \( W(kr)\) is the Fourier transform of the 3D top-hat function.

Definition at line 1376 of file cosmo.cc.

1377 {
1378  // returns d(sigma^2)/d(lnk) \propto k^3 P(k) W(kr) where r is the scale over which
1379 
1381  double rh = (params->rh);
1382  vector<double> vec_lnk = (params->vec_lnk);
1383  vector<double> linear_lnpkz = (params->linear_lnpkz);
1384  int card_window = params->card_window;
1385 
1386  double k = exp(lnk);
1387  double x = k * rh;
1388 
1389  // choose window function:
1390  double win = window_fnct(x, card_window);
1391  // Given tabulated linear P(k,z), interpolate the value of p(k) for in log(k):
1392  double linear_pk = exp(interp1D(log(k), vec_lnk, linear_lnpkz, kSPLINE));
1393  return pow(k, 3.) * (linear_pk * pow(win, 2));
1394 }
double window_fnct(const double &x, const int card_window)
Definition: cosmo.cc:2541
Spline interpolation.
Definition: params.h:184
vector< double > linear_lnpkz
Definition: cosmo.h:18
vector< double > vec_lnk
Definition: cosmo.h:17
double interp1D(double const &x_eval, const vector< double > &x_base, const vector< double > &y_base, const int card_interp, bool no_error=false)
Definition: misc.cc:433

◆ dVhdzdOmega()

double dVhdzdOmega ( double  z)

Differential comoving volume element in [ \(h^{-3} \; Mpc^3 \; sr^{-1}\)].

Definition at line 1397 of file cosmo.cc.

1398 {
1399  //--- Computes the differential comoving volume dV/dz element per steradian at
1400  // redshift z, in units if h^-3 Mpc^3 sr^-1.
1401  // INPUTS:
1402  // z redshift
1403  // OUTPUT:
1404  // dV/dz/dOmega co-moving volume element [h^-3 Mpc^3 sr^-1]
1405 
1406  return pow(dh_trans(z), 2.) * (HUBBLE_LENGTHxh_Mpc * H0_over_H(z)); //Hubble distance//
1407 }
double dh_trans(const double &z)
Computes the transverse comoving distance by integration.
Definition: cosmo.cc:342
#define HUBBLE_LENGTHxh_Mpc
Hubble length in [ ].
Definition: inlines.h:5
double H0_over_H(const double z, void *params)
where is the Hubble rate at redshift .
Definition: cosmo.cc:2300

◆ gamma_f()

double gamma_f ( const int  card_window)

Definition at line 1479 of file cosmo.cc.

1480 {
1481  //--- Sets res = gamma_f corresponding to the window function in k-space
1482  // Adapted from van den Bosch presentation slides,
1483  // http://www.astro.yale.edu/vdbosch/astro610_lecture9.pdf
1484  //
1485  // INPUTS:
1486  // card_window Choice of window function
1487  // OUTPUT:
1488  // res Value of normalisation gamma_f
1489 
1490  double res;
1491  // Calculate W(x)
1492  switch (card_window) {
1493  case kTOP_HAT:
1494  res = 4.* PI / 3.;
1495  break;
1496  case kGAUSS:
1497  res = pow(2. * PI, 1.5);
1498  break;
1499  case kSHARP_K:
1500  res = 6 * pow(PI, 2.);
1501  break;
1502  default :
1503  printf("\n====> ERROR: gamma_f() in cosmo.cc");
1504  printf("\n window card_window=%d does not correspond to any window function", card_window);
1505  printf("\n => abort()\n\n");
1506  abort();
1507  }
1508  return res;
1509 }
See user documentation
Definition: params.h:167
See user documentation
Definition: params.h:166
See user documentation
Definition: params.h:168
#define PI
What do you think it is?
Definition: inlines.h:34

◆ get_pk()

void get_pk ( const double &  z,
vector< double > &  lnk_vec,
vector< double > &  lnp_vec,
const bool  is_nonlinear = false,
const bool  is_force_recompute = false 
)

Loads the matter power spectrum from a file. Fills in vec_lnk with mode number \(\ln(k)\), \( k \) in [ \( h \;Mpc^{-1}\)]. Fills in linear_lnpkz with \(\ln(p_{\rm matter}(k,z))\) [ \( Mpc^3\;h^{-3}\)]. If the p_matter(k,z) file is not present, calls CLASS to compute it.

Definition at line 1512 of file cosmo.cc.

1513 {
1514  // - Computes the matter power spectrum at the redshifts z
1515  // --> Fills in vec_lnk with mode number ln(k) [k in h Mpc^-1]
1516  // --> Fills in linear_lnpkz with ln(p_matter(k,z)) [Mpc^3 h^-3]
1517  // If the p_matter(k,z) file is not present, calls CLASS to compute it
1518 
1519  // because we push back elements to a vector accessed by reference,
1520  // we want to make sure that it is empty at this point:
1521  vector<double>().swap(lnk_vec);
1522  vector<double>().swap(lnp_vec);
1523 
1524  ostringstream convert;
1525  convert << round(gCOSMO_HUBBLE, 3);
1526  string h_str = convert.str();
1527  convert.str("");
1528  convert << round(gCOSMO_OMEGA0_M, 3);
1529  string OmegaM_str = convert.str();
1530  convert.str("");
1531  convert << round(gCOSMO_OMEGA0_B, 4);
1532  string OmegaB_str = convert.str();
1533  convert.str("");
1535  string OmegaL_str = convert.str();
1536  convert.str("");
1537  convert << gCOSMO_N_S;
1538  string ns_str = convert.str();
1539  convert.str("");
1541  string tau_str = convert.str();
1542  convert.str("");
1543 
1544  double zz = z;
1545  if (z <= 1e-5) zz = 0.;
1546 
1547  char char_tmp[80];
1548  sprintf(char_tmp, "%.3g", zz);
1549  if (z > 10) sprintf(char_tmp, "%.5g", zz);
1550  string zstr = string(char_tmp);
1551 
1552  convert.str("");
1553 
1554  string file = "dummy";
1555  string path = gPATH_TO_CLUMPY_HOME + "/data/pk_precomp/";
1556 
1557  string filename_base = "h" + h_str + "_OmegaB"
1558  + OmegaB_str + "_OmegaM" + OmegaM_str + "_OmegaL" + OmegaL_str
1559  + "_ns" + ns_str + "_tau" + tau_str + "_z" + zstr;
1560 
1561  string file_pattern;
1562 
1563  if (is_nonlinear) {
1564  file_pattern = "_" + filename_base + "_nl.dat";
1565  } else {
1566  file_pattern = "_" + filename_base + "_lin.dat";
1567  }
1568 
1569  string command = "find " + path + "*" + file_pattern;
1570 
1571  string found_files = sys_command(command);
1572  string prefix;
1573  if (found_files != "") {
1574  prefix = found_files.substr(found_files.find_last_of("\\/") + 1, found_files.size());
1575  file = path + prefix;
1576  prefix = prefix.substr(0, prefix.find_first_of("_"));
1577  }
1578  file = removeblanks_from_startstop(file);
1579 
1580  ifstream ifile(file.c_str());
1581 
1582  if (!ifile or is_force_recompute) {
1583  ifile.close();
1584  if (is_force_recompute) cout << " => Force recomputation of P(k,z) file for" << endl;
1585  else cout << " => P(k,z) file does not exist for" << endl;
1586  cout << " (z,h,OmegaB,OmegaM,OmegaL,ns,tau) = ("
1587  + zstr + "," + h_str + "," + OmegaB_str + "," + OmegaM_str + "," + OmegaL_str
1588  + "," + ns_str + "," + tau_str + ")." << endl;
1589  cout << " => Running CLASS to compute it ..." << endl;
1590  compute_class_pk(zz);
1591  file = path + "class" + file_pattern;
1592  ifstream ifile_check(file.c_str());
1593  if (!ifile_check) {
1594  printf("\n====> ERROR: get_pk() in cosmo.cc");
1595  printf("\n Something went wrong in computing P(k,z) with CLASS (probably aborted.)");
1596  printf("\n => abort()\n\n");
1597  abort();
1598  }
1599  ifile.open(file.c_str());
1600  } else {
1601  if (is_nonlinear) cout << " Loading non-linear P(k) file '" << prefix << "' for z = " << zstr << " and given cosmology ";
1602  else cout << " Loading linear P(k) file '" << prefix << "' for z = " << zstr << " and given cosmology ";
1603  }
1604 
1605  vector< vector<double> > table = read_ascii(file);
1606 
1607 
1608  for (int i = 0; i < int(table.size()); ++i) {
1609  for (int j = 0; j < int(table[i].size()); ++j) {
1610  if (j == 0) lnk_vec.push_back(log(table[i][j]));
1611  else if (j == 1) {
1612  double pk = table[i][j];
1613  if (gDM_KMAX > 0) {
1614  pk *= (1 - sigmoid_window(table[i][0], gDM_KMAX, gDM_KMAX_SIGMA_CUTOFF));
1615  if (pk < 1e-40) pk = 1e-40;
1616  }
1617  lnp_vec.push_back(log(pk));
1618  } else {
1619  printf("\n====> ERROR: get_pk() in cosmo.cc");
1620  printf("\n File %s contains more than two columns", file.c_str());
1621  printf("\n => abort()\n\n");
1622  abort();
1623  }
1624  }
1625  }
1626 
1627  // check if loaded linear P(k) exists up to gSIM_EXTRAGAL_KMAX_PRECOMP:
1628  double kmax = exp(find_max_value(lnk_vec));
1629  static int please_break = 0; // don"t get lost in an infinite loop...
1630  if (!is_nonlinear and kmax < gSIM_EXTRAGAL_KMAX_PRECOMP) {
1631  please_break += 1;
1632  if (please_break < 2) {
1633  // load recursively:
1634  cout << "\n => P(k,z) file does not exist up to k_max = " << gSIM_EXTRAGAL_KMAX_PRECOMP << endl;
1635  cout << " => Running CLASS to compute it ..." << endl;
1636  compute_class_pk(zz);
1637  get_pk(z, lnk_vec, lnp_vec);
1638  } else {
1639  printf("\n====> ERROR: get_pk() in cosmo.cc");
1640  printf("\n Abort after one unsuccessful trial to recompute P(k) file.");
1641  printf("\n => abort()\n\n");
1642  abort();
1643  }
1644  } else {
1645  please_break = 0;
1646  cout << "...done." << endl;
1647  }
1648  return;
1649 }
double gDM_KMAX
See the user documentation
Definition: params.cc:95
void compute_class_pk(double z)
Calls CLASS to compute the matter power spectrum at redshift .
Definition: cosmo.cc:32
double gCOSMO_OMEGA0_LAMBDA
Present day dark energy density of the -CDM Universe (PDG)
Definition: params.cc:77
double gDM_KMAX_SIGMA_CUTOFF
See the user documentation
Definition: params.cc:96
double gCOSMO_HUBBLE
See the user documentation
Definition: params.cc:62
double gCOSMO_OMEGA0_B
See the user documentation
Definition: params.cc:64
string sys_command(string const &command)
Definition: misc.cc:1767
double find_max_value(vector< double > &array, int i_start=0, int i_end=-1)
Definition: misc.cc:115
def convert(infile='empty', i_ext=-1, i_col=-1, coordsys_out='-1')
string removeblanks_from_startstop(string const &str)
Definition: misc.cc:1085
vector< vector< double > > read_ascii(string const &filename, int const ncolumns=-1, bool const is_monotonous=false, int const column_start=0)
Definition: misc.cc:1003
void get_pk(const double &z, vector< double > &lnk_vec, vector< double > &lnp_vec, const bool is_nonlinear, const bool is_force_recompute)
Loads the matter power spectrum from a file. Fills in vec_lnk with mode number , in [ ]....
Definition: cosmo.cc:1512
string gPATH_TO_CLUMPY_HOME
Absolute path to CLUMPY home directory.
Definition: params.cc:299
double round(const double &x, const int digits)
Definition: misc.cc:1570
double gCOSMO_N_S
See the user documentation
Definition: params.cc:67
double sigmoid_window(double const &x, double const &x0, double const &sigma)
Definition: misc.cc:1580
double gSIM_EXTRAGAL_KMAX_PRECOMP
See the user documentation
Definition: params.cc:232
double gCOSMO_TAU_REIO
See the user documentation
Definition: params.cc:68
double gCOSMO_OMEGA0_M
See the user documentation
Definition: params.cc:63

◆ growth_factor()

double growth_factor ( const double &  z,
const int  card_growth_mthd 
)

Computes the linear growth rate factor.

Definition at line 1652 of file cosmo.cc.

1653 {
1654  //--- Computes the linear growth rate factor of linear density perturbations.
1655  //
1656  // INPUTS:
1657  // z redshift
1658  // card_growth_mthd Approximation/Model of the linear growth rate
1659  // OUTPUT:
1660  // res linear growth rate factor
1661 
1662  if (z == 0) return 1.;
1663 
1664  double res;
1665  // Calculate growth factor:
1666  switch (card_growth_mthd) {
1667  case kCARROLL92: {
1668  // approximation to integral solution from Heath (1977)
1669  double a = 1. / (1. + z);
1670  double omegaM0_var = 1. / (1 + gCOSMO_OMEGA0_LAMBDA / gCOSMO_OMEGA0_M);
1671  double omegaL0_var = 1 - omegaM0_var;
1672  double omegaM_var = 1. / (1 + gCOSMO_OMEGA0_LAMBDA / gCOSMO_OMEGA0_M * pow(a, 3));
1673  double omegaL_var = 1 - omegaM_var;
1674  // discard factor 2.5 in numerators for ratio:
1675  double g0 = omegaM0_var / (pow(omegaM0_var, 4. / 7.) - omegaL0_var + (1. + omegaM0_var / 2.) / (1. + omegaL0_var / 70.));
1676  double gz = omegaM_var * a / (pow(omegaM_var, 4. / 7.) - omegaL_var + (1. + omegaM_var / 2.) / (1. + omegaL_var / 70.));
1677  res = gz / g0;
1678  break;
1679  }
1680  case kHEATH77: {
1681  // exact analytical solution from Heath (1977)
1682  gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
1683  double res_int0, res_intz, error;
1684  gsl_function F;
1685  F.function = &growth_integrand;
1686  gsl_integration_qagiu(&F, 0, 0, gSIM_EPS, 1000, w, &res_int0, &error);
1687  gsl_integration_qagiu(&F, z, 0, gSIM_EPS, 1000, w, &res_intz, &error);
1688  gsl_integration_workspace_free(w);
1689  // discard factor 2.5 Omega_M in numerators for ratio:
1690  res = H_over_H0(z) * res_intz / res_int0;
1691  break;
1692  }
1693  default :
1694  printf("\n====> ERROR: growth_factor() in cosmo.cc");
1695  printf("\n growth factor card = %d does not correspond to any model", card_growth_mthd);
1696  printf("\n => abort()\n\n");
1697  abort();
1698  }
1699  return res;
1700 }
double gCOSMO_OMEGA0_LAMBDA
Present day dark energy density of the -CDM Universe (PDG)
Definition: params.cc:77
double gSIM_EPS
See the user documentation
Definition: params.cc:196
See user documentation
Definition: params.h:190
double growth_integrand(const double z, void *params)
Definition: cosmo.cc:1703
double H_over_H0(const double &z)
where is the Hubble rate at redshift .
Definition: cosmo.cc:2293
See user documentation
Definition: params.h:191
double gCOSMO_OMEGA0_M
See the user documentation
Definition: params.cc:63

◆ growth_integrand()

double growth_integrand ( const double  z,
void *  params = NULL 
)

Definition at line 1703 of file cosmo.cc.

1704 {
1705  // returns (1+z) * H_0^3/H(z)^3 to be integrated with gsl
1706 
1707  return (1 + z) * pow(H2_over_H02(z), -1.5);
1708 }
double H2_over_H02(const double &z)
Definition: cosmo.cc:2278

◆ H0_over_H()

double H0_over_H ( const double  z,
void *  params = NULL 
)

\(1/h(z) = H_0/H(z)\) where \(H(z)\) is the Hubble rate at redshift \(z\).

Definition at line 2300 of file cosmo.cc.

2301 {
2302  // returns H0_over_H(z) = H_0/H(z) approved by Celine. dfgfd
2303 
2304  return 1. / H_over_H0(z);
2305 }
double H_over_H0(const double &z)
where is the Hubble rate at redshift .
Definition: cosmo.cc:2293

◆ H2_over_H02()

double H2_over_H02 ( const double &  z)

Definition at line 2278 of file cosmo.cc.

2279 {
2280  // returns H(z)^2/H_0^2
2281  // x = 1/a = (1+z)
2282 
2283  double x = (1. + z);
2284  double fac_lambda;
2285  if (abs(gCOSMO_WDE + 1.) > 1e-5) fac_lambda = pow(x, (3. + 3.*gCOSMO_WDE));
2286  else fac_lambda = 1.;
2287 
2288  return gCOSMO_OMEGA0_R * pow(x, 4.) + gCOSMO_OMEGA0_M * pow(x, 3.)
2289  + gCOSMO_OMEGA0_K * pow(x, 2.) + gCOSMO_OMEGA0_LAMBDA * fac_lambda;
2290 }
double gCOSMO_OMEGA0_R
Present day radiation density of the Universe (PDG)
Definition: params.cc:76
double gCOSMO_OMEGA0_LAMBDA
Present day dark energy density of the -CDM Universe (PDG)
Definition: params.cc:77
double gCOSMO_OMEGA0_K
See the user documentation
Definition: params.cc:65
double gCOSMO_WDE
See the user documentation
Definition: params.cc:69
double gCOSMO_OMEGA0_M
See the user documentation
Definition: params.cc:63

◆ H_over_H0()

double H_over_H0 ( const double &  z)

\(H(z)/H_0\) where \(H(z)\) is the Hubble rate at redshift \(z\).

Definition at line 2293 of file cosmo.cc.

2294 {
2295  // returns H(z)/H_0
2296  return sqrt(H2_over_H02(z));
2297 }
double H2_over_H02(const double &z)
Definition: cosmo.cc:2278

◆ legend_for_delta()

string legend_for_delta ( )

Definition at line 1410 of file cosmo.cc.

1411 {
1412  //--- Returns a meaningful delta subscript according to the delta choice.
1413  //
1414  char Delta_char[50];
1415  sprintf(Delta_char, "%g", gCOSMO_DELTA0);
1416  string Delta_str = string(Delta_char);
1418  return Delta_str + ",c";
1419  } else if (gCOSMO_FLAG_DELTA_REF == kRHO_MEAN) {
1420  return Delta_str + ",m";
1421  } else if (gCOSMO_FLAG_DELTA_REF == kBRYANNORMAN98) {
1422  if (Delta_str == "-1") return "Bryan&Norman";
1423  else return Delta_str + ",Bryan&Norman";
1424  } else {
1425  return "unknown";
1426  }
1427 }
See user documentation
Definition: params.h:173
double gCOSMO_DELTA0
See the user documentation
Definition: params.cc:71
int gCOSMO_FLAG_DELTA_REF
See the user documentation
Definition: params.cc:72
See user documentation
Definition: params.h:175
See user documentation
Definition: params.h:174

◆ lnsigma2()

double lnsigma2 ( double  lnR,
void *  p 
)

Returns \(\ln(\sigma^2)\) for the scale \( R \) passed as \(\ln(R)\).

Definition at line 1430 of file cosmo.cc.

1431 {
1432  // returns ln(2pi^2/D(z)^2 * sigma^2) for the scale R (D(z): growth factor)
1433  // integration is performed between kmin=1.e-4 and kmax=20./R...
1434  // upper bound ensures enough modes are considered, without integrating over unnecessary region?
1435  //
1436  // INPUTS:
1437  // lnRh ln of scale radius in units of Mpc h^-1.
1438  // OUTPUT:
1439  // res ln(2pi^2 * sigma^2)
1440 
1442  params->rh = exp(lnRh);
1443 
1444  char char_tmp[512];
1445 
1446  double lnkmin_in = find_min_value(params->vec_lnk);
1447  double lnkmax_in = find_max_value(params->vec_lnk);
1448  double lnkmin = log(1.e-3 / params->rh);
1449  double lnkmax = log(1.e2 / params->rh);
1450  if (lnkmin_in > lnkmin) {
1451  sprintf(char_tmp, "Integration truncated by limits of input k grid: k_min = %g -> %g for scale %g",
1452  exp(lnkmin), exp(lnkmin_in), 2 * PI / params->rh);
1453  print_warning("cosmo.cc", "lnsigma2()", string(char_tmp));
1454  lnkmin = lnkmin_in;
1455  }
1456  if (lnkmax_in < lnkmax) {
1457  sprintf(char_tmp, "Integration truncated by limits of input k grid: k_max = %g -> %g for scale %g",
1458  exp(lnkmax), exp(lnkmax_in), 2 * PI / params->rh);
1459  print_warning("cosmo.cc", "lnsigma2()", string(char_tmp));
1460  lnkmax = lnkmax_in;
1461  }
1462 
1463  gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
1464  double res, error;
1465  gsl_function F;
1466  F.function = &ds2dlnk;
1467  F.params = params;
1468 
1469  double eps = gSIM_EPS;
1470  if (gSIM_EXTRAGAL_FLAG_WINDOWFUNC == kSHARP_K) eps *= 1e-2;
1471  // integration over ln(k):
1472  gsl_integration_qags(&F, lnkmin, lnkmax, 0, eps, 1000, w, &res, &error);
1473  gsl_integration_workspace_free(w);
1474 
1475  return log(res);
1476 }
double find_min_value(vector< double > &array, int i_start=0, int i_end=-1)
Definition: misc.cc:134
double gSIM_EPS
See the user documentation
Definition: params.cc:196
double eps
void print_warning(string library, string function, string message)
Definition: misc.cc:986
double find_max_value(vector< double > &array, int i_start=0, int i_end=-1)
Definition: misc.cc:115
double ds2dlnk(double lnk, void *p)
Returns , where r is the scale of interest and is the Fourier transform of the 3D top-hat function.
Definition: cosmo.cc:1376
See user documentation
Definition: params.h:168
vector< double > vec_lnk
Definition: cosmo.h:17
#define PI
What do you think it is?
Definition: inlines.h:34
int gSIM_EXTRAGAL_FLAG_WINDOWFUNC
See the user documentation
Definition: params.cc:230

◆ mf()

double mf ( int  card_mf,
double const &  nu,
double const &  z,
double const &  Delta_c 
)

Definition at line 1711 of file cosmo.cc.

1712 {
1713  // INPUTS:
1714  // card_mf mass function parametrization
1715  // lnnu
1716  // z redshift
1717  // Delta_c Delta w.r.t critical density
1718  // OUTPUT:
1719  //
1720 
1721  switch (card_mf) {
1723  return mf_bolshoi_planck(nu, z, Delta_c);
1724  case kBOCQUET16_HYDRO:
1725  return mf_magneticum(nu, z, Delta_c, true);
1726  case kBOCQUET16_DMONLY:
1727  return mf_magneticum(nu, z, Delta_c, false);
1728  case kTINKER08:
1729  return mf_tinker(nu, z, Delta_c);
1730  case kTINKER10:
1731  return mf_tinker10(nu, z, Delta_c);
1732  case kTINKER08_N:
1733  return mf_tinker_norm(nu, z, Delta_c);
1734  case kSHETHTORMEN99:
1735  return mf_sheth_tormen(nu, z, Delta_c);
1736  case kPRESSSCHECHTER74:
1737  return mf_press_schechter(nu, z, Delta_c);
1738  case kJENKINS01:
1739  return mf_jenkins(nu, z, Delta_c);
1740  default :
1741  printf("\n====> ERROR: mf() in cosmo.cc");
1742  printf("\n card_mf = %d not a valid mass function!", card_mf);
1743  printf("\n => abort()\n\n");
1744  abort();
1745  }
1746  return 0;
1747 }
See user documentation
Definition: params.h:120
double mf_tinker10(double const &nu, double const &z, double const &Delta_c)
Definition: cosmo.cc:2122
double mf_jenkins(double const &nu, double const &z, double const &Delta_c)
Definition: cosmo.cc:1800
See user documentation
Definition: params.h:116
double mf_tinker_norm(double const &nu, double const &z, double const &Delta_c)
Definition: cosmo.cc:2047
double mf_tinker(double const &nu, double const &z, double const &Delta_c)
Definition: cosmo.cc:1969
See user documentation
Definition: params.h:119
double mf_press_schechter(double const &nu, double const &z, double const &Delta_c)
Definition: cosmo.cc:1943
See user documentation
Definition: params.h:121
See user documentation
Definition: params.h:123
See user documentation
Definition: params.h:117
double mf_bolshoi_planck(double const &nu, double const &z, double const &Delta_c)
Definition: cosmo.cc:1750
double mf_magneticum(double const &nu, double const &z, double const &Delta_c, const bool is_hydro)
Definition: cosmo.cc:1830
See user documentation
Definition: params.h:124
See user documentation
Definition: params.h:118
See user documentation
Definition: params.h:122
double mf_sheth_tormen(double const &nu, double const &z, double const &Delta_c)
Definition: cosmo.cc:2211

◆ mf_bolshoi_planck()

double mf_bolshoi_planck ( double const &  nu,
double const &  z,
double const &  Delta_c 
)

Definition at line 1750 of file cosmo.cc.

1751 {
1752  // Halo multiplicity function per log(nu): 0.5*f(sigma),
1753  // where sigma=deltac/sqrt(nu).
1754  // \int_{-infty}^infty mf(lnnu,z) dlnnu = 1 is NOT satisfied!
1755  // Ref: Eq.(25) of Rodr�guez-Puebla (2016), after Tinker et al. (2008)
1756 
1757  double Delta_c_check = delta_x_to_delta_crit(-1, kBRYANNORMAN98, z);
1758 
1759  if (fabs(Delta_c - Delta_c_check) > 0.5) {
1760  printf("\n\n====> ERROR: mf_bolshoi_planck() in cosmo.cc");
1761  printf("\n Rodriguez-Puebla (2016) HMF parametrization is only");
1762  printf("\n given for values Delta_c(z=%.f) = %d", z, int(Delta_c_check));
1763  printf("\n (scaling according to Bryan & Norman, 1998)");
1764  printf("\n If you want to use a different Delta,");
1765  printf("\n you must choose gSIM_EXTRAGAL_IS_DELTA_HMF_FROM_PROFILE = 1.");
1766  printf("\n => abort()\n\n");
1767  abort();
1768  }
1769 
1770  double params[4][3] = {
1771  {0.144, -0.011, 0.003},
1772  {1.351, 0.068, 0.006},
1773  {3.113, -0.077, -0.013},
1774  {1.187, 0.009, 0.000},
1775  };
1776 
1777  double A, a, b, c;
1778 
1779  A = params[0][0] + params[0][1] * z + params[0][2] * pow(z, 2);
1780  a = params[1][0] + params[1][1] * z + params[1][2] * pow(z, 2);
1781  b = params[2][0] + params[2][1] * z + params[2][2] * pow(z, 2);
1782 
1783  if (b < 0) {
1784  static bool is_print_msg = false;
1785  if (!is_print_msg) {
1786  char char_tmp[512];
1787  sprintf(char_tmp, "b<0 obtained at z>=%g (causes negative root). Set b=0 resulting in extrapolation f(sigma)=A*exp(-c/sigma^2).", z);
1788  print_warning("cosmo.cc", "mf_bolshoi_planck()", string(char_tmp));
1789  is_print_msg = true;
1790  }
1791  b = 1e-40;
1792  }
1793  c = params[3][0] + params[3][1] * z + params[3][2] * pow(z, 2);
1794 
1795  double sigma = DELTA_COLLAPSE / sqrt(nu);
1796  return 0.5 * (A * (pow(sigma / b, -a) + 1.) * exp(-c / pow(sigma, 2.)));
1797 }
#define DELTA_COLLAPSE
Definition: inlines.h:8
void print_warning(string library, string function, string message)
Definition: misc.cc:986
double delta_x_to_delta_crit(double const &Delta_x0, int const card_delta_ref, double const &z)
Definition: clumps.cc:102
See user documentation
Definition: params.h:175

◆ mf_jenkins()

double mf_jenkins ( double const &  nu,
double const &  z,
double const &  Delta_c 
)

Definition at line 1800 of file cosmo.cc.

1801 {
1802 // Halo multiplicity function per log(nu): 0.5*f(sigma),
1803 // where sigma=deltac/sqrt(nu).
1804 // \int_{-infty}^infty mf(lnnu) dlnnu = 1 is NOT satisfied!
1805 // Ref: Eq.(B3) of Jenkins et al. (2001)
1806 // Delta_c = 178
1807 
1808  double Delta_m = delta_crit_to_delta_x(178, kRHO_MEAN, z);
1809  if (Delta_c < 177.5 or Delta_c > 178.5) {
1810  printf("\n\n====> ERROR: mf_jenkins() in cosmo.cc");
1811  printf("\n Jenkins et al. (2001) HMF parametrization is ");
1813  printf("\n given for Delta_m(z) = %d (Delta_c = 178)", int(round(Delta_m)));
1814  } else if (gCOSMO_FLAG_DELTA_REF == kBRYANNORMAN98) {
1815  double Delta_bn = delta_crit_to_delta_x(178, kBRYANNORMAN98, z);
1816  printf("\n given for Delta_Bryan() = %d (Delta_c = 178)", int(round(Delta_bn)));
1817  } else
1818  printf("\n given for values Delta_c = 178 (Delta_m(z) = %d)", int(round(Delta_m)));
1819  printf("\n If you want to use a different Delta,");
1820  printf("\n you must choose gSIM_EXTRAGAL_IS_DELTA_HMF_FROM_PROFILE = 1.");
1821  printf("\n => abort()\n\n");
1822  abort();
1823  }
1824 
1825  double sigma = DELTA_COLLAPSE / sqrt(nu);
1826  return 0.5 * 0.301 * exp(-pow(abs(log(1. / sigma) + 0.64), 3.82)); // Delta=180
1827 }
#define DELTA_COLLAPSE
Definition: inlines.h:8
int gCOSMO_FLAG_DELTA_REF
See the user documentation
Definition: params.cc:72
double delta_crit_to_delta_x(double const &Delta_c, int const card_delta_ref, double const &z)
Definition: clumps.cc:182
double round(const double &x, const int digits)
Definition: misc.cc:1570
See user documentation
Definition: params.h:175
See user documentation
Definition: params.h:174

◆ mf_magneticum()

double mf_magneticum ( double const &  nu,
double const &  z,
double const &  Delta_c,
const bool  is_hydro 
)

Definition at line 1830 of file cosmo.cc.

1831 {
1832 // Halo multiplicity function per log(nu): 0.5*f(sigma),
1833 // where sigma=deltac/sqrt(nu).
1834 // \int_{-infty}^infty mf(lnnu,z) dlnnu = 1 is NOT satisfied!
1835 // Coefficients correspond to Delta_mean=200, Hydro or DM only (Table 2 of Bocquet et al. (2016))
1836 // Ref: Eq.(1) of Bocquet et al. (2016)
1837 
1838  double Delta_c_of_Deltam_200 = delta_x_to_delta_crit(200, kRHO_MEAN, z);
1839 
1840  if (Delta_c < Delta_c_of_Deltam_200) {
1841  printf("\n\n====> ERROR: mf_magneticum() in cosmo.cc");
1842  printf("\n Bocquet et al. (2016) HMF parametrization is only");
1844  printf("\n given for values Delta_m >= 200");
1845  } else if (gCOSMO_FLAG_DELTA_REF == kBRYANNORMAN98) {
1846  double Delta_bn = delta_crit_to_delta_x(Delta_c_of_Deltam_200, kBRYANNORMAN98, z);
1847  printf("\n given for values Delta_Bryan >= %d", int(round(Delta_bn)));
1848  } else
1849  printf("\n given for values Delta_m >= 200.");
1850  printf("\n (This corresponds to Delta_c >= %d at z = %g)", int(round(Delta_c_of_Deltam_200)), z);
1851  printf("\n If you want to use a smaller Delta,");
1852  printf("\n you must choose gSIM_EXTRAGAL_IS_DELTA_HMF_FROM_PROFILE = 1.");
1853  printf("\n => abort()\n\n");
1854  abort();
1855  } else if (Delta_c > 500) {
1856  printf("\n\n====> ERROR: mf_magneticum() in cosmo.cc");
1857  printf("\n Bocquet et al. (2016) HMF parametrization is only");
1858  printf("\n given for values Delta_c <= 500");
1859  printf("\n If you want to use a larger Delta_c,");
1860  printf("\n you must choose gSIM_EXTRAGAL_IS_DELTA_HMF_FROM_PROFILE = 1.");
1861  printf("\n => abort()\n\n");
1862  abort();
1863  }
1864 
1865  // avoid problems if redshift-dependent Delta_c(Delta_m200) becomes larger than 200,
1866  // and also stop at 195 to avoid non-steady transition:
1867  if (Delta_c_of_Deltam_200 > 195) {
1868  Delta_c_of_Deltam_200 = 195;
1869  }
1870 
1871  double deltas[3] = {
1872  Delta_c_of_Deltam_200, 200, 500
1873  };
1874 
1875  vector<double> deltas_vec;
1876  deltas_vec.assign(deltas, deltas + 3);
1877  vector < vector<double> > params_vec;
1878 
1879  if (is_hydro) {
1880  double params[8][3] = {
1881  { 0.228, 0.202, 0.18 },
1882  { 2.15, 2.21, 2.29 },
1883  { 1.69, 2.00, 2.44 },
1884  { 1.30, 1.57, 1.97 },
1885  { 0.285, 1.147, 1.088},
1886  { -0.058, 0.375, 0.15 },
1887  { -0.366, -1.074, -1.008},
1888  { -0.045, -0.196, -0.322},
1889  };
1890  pushback_array2vector(params_vec, &params[0][0], sizeof params[0]);
1891  pushback_array2vector(params_vec, &params[1][0], sizeof params[1]);
1892  pushback_array2vector(params_vec, &params[2][0], sizeof params[2]);
1893  pushback_array2vector(params_vec, &params[3][0], sizeof params[3]);
1894  pushback_array2vector(params_vec, &params[4][0], sizeof params[4]);
1895  pushback_array2vector(params_vec, &params[5][0], sizeof params[5]);
1896  pushback_array2vector(params_vec, &params[6][0], sizeof params[6]);
1897  pushback_array2vector(params_vec, &params[7][0], sizeof params[7]);
1898  } else {
1899  double params[8][3] = {
1900  { 0.175, 0.222, 0.241},
1901  { 1.53, 1.71, 2.18 },
1902  { 2.55, 2.24, 2.35 },
1903  { 1.19, 1.46, 2.02 },
1904  { -0.012, 0.26, 0.370},
1905  { -0.040, 0.321, 0.251},
1906  { -0.194, -0.62, -0.698},
1907  { -0.021, -0.153, -0.310},
1908  };
1909  pushback_array2vector(params_vec, &params[0][0], sizeof params[0]);
1910  pushback_array2vector(params_vec, &params[1][0], sizeof params[1]);
1911  pushback_array2vector(params_vec, &params[2][0], sizeof params[2]);
1912  pushback_array2vector(params_vec, &params[3][0], sizeof params[3]);
1913  pushback_array2vector(params_vec, &params[4][0], sizeof params[4]);
1914  pushback_array2vector(params_vec, &params[5][0], sizeof params[5]);
1915  pushback_array2vector(params_vec, &params[6][0], sizeof params[6]);
1916  pushback_array2vector(params_vec, &params[7][0], sizeof params[7]);
1917  }
1918 
1919  double Ap0 = interp1D(Delta_c, deltas_vec, params_vec[0], kLINLIN);
1920  double a0 = interp1D(Delta_c, deltas_vec, params_vec[1], kLINLIN);
1921  double b0 = interp1D(Delta_c, deltas_vec, params_vec[2], kLINLIN);
1922  double c0 = interp1D(Delta_c, deltas_vec, params_vec[3], kLINLIN);
1923 
1924  double Apz = interp1D(Delta_c, deltas_vec, params_vec[4], kLINLIN);
1925  double az = interp1D(Delta_c, deltas_vec, params_vec[5], kLINLIN);
1926  double bz = interp1D(Delta_c, deltas_vec, params_vec[6], kLINLIN);
1927  double cz = interp1D(Delta_c, deltas_vec, params_vec[7], kLINLIN);
1928 
1929  double Ap, a, b, c;
1930 
1931  Ap = Ap0 * pow(1. + z, Apz);
1932  a = a0 * pow(1. + z, az);
1933  b = b0 * pow(1. + z, bz);
1934  c = c0 * pow(1. + z, cz);
1935 
1936  double sigma = DELTA_COLLAPSE / sqrt(nu);
1937  double mf = 0.5 * (Ap * (pow(sigma / b, -a) + 1.) * exp(-c / pow(sigma, 2.)));
1938  if (std::isnan(mf)) mf = 0.;
1939  return mf; // Delta_mean = 200
1940 }
#define DELTA_COLLAPSE
Definition: inlines.h:8
void pushback_array2vector(vector< vector< T > > &vec, const T *array, const int size)
Definition: misc.h:104
int gCOSMO_FLAG_DELTA_REF
See the user documentation
Definition: params.cc:72
double delta_crit_to_delta_x(double const &Delta_c, int const card_delta_ref, double const &z)
Definition: clumps.cc:182
double delta_x_to_delta_crit(double const &Delta_x0, int const card_delta_ref, double const &z)
Definition: clumps.cc:102
double round(const double &x, const int digits)
Definition: misc.cc:1570
double mf(int card_mf, double const &nu, double const &z, double const &Delta_c)
Definition: cosmo.cc:1711
See user documentation
Definition: params.h:175
See user documentation
Definition: params.h:174
x-y linear interpolation
Definition: params.h:180
double interp1D(double const &x_eval, const vector< double > &x_base, const vector< double > &y_base, const int card_interp, bool no_error=false)
Definition: misc.cc:433

◆ mf_press_schechter()

double mf_press_schechter ( double const &  nu,
double const &  z,
double const &  Delta_c 
)

Definition at line 1943 of file cosmo.cc.

1944 {
1945 // Halo multiplicity function per log(nu): nu*f(nu)
1946 // \int_{-infty}^infty mf(lnnu) dlnnu = 1
1947 
1948  double Delta_m = delta_crit_to_delta_x(178, kRHO_MEAN, z);
1949  if (Delta_c < 177.5 or Delta_c > 178.5) {
1950  printf("\n\n====> ERROR: mf_press_schechter() in cosmo.cc");
1951  printf("\n Press & Schechter HMF parametrization is ");
1953  printf("\n given for Delta_m(z) = %d (Delta_c = 178)", int(round(Delta_m)));
1954  } else if (gCOSMO_FLAG_DELTA_REF == kBRYANNORMAN98) {
1955  double Delta_bn = delta_crit_to_delta_x(178, kBRYANNORMAN98, z);
1956  printf("\n given for Delta_Bryan() = %d (Delta_c = 178)", int(round(Delta_bn)));
1957  } else
1958  printf("\n given for values Delta_c = 178 (Delta_m(z) = %d)", int(round(Delta_m)));
1959  printf("\n If you want to use a different Delta,");
1960  printf("\n you must choose gSIM_EXTRAGAL_IS_DELTA_HMF_FROM_PROFILE = 1.");
1961  printf("\n => abort()\n\n");
1962  abort();
1963  }
1964 
1965  return 1. / sqrt(2.*PI) * sqrt(nu) * exp(-nu / 2.);
1966 }
int gCOSMO_FLAG_DELTA_REF
See the user documentation
Definition: params.cc:72
double delta_crit_to_delta_x(double const &Delta_c, int const card_delta_ref, double const &z)
Definition: clumps.cc:182
double round(const double &x, const int digits)
Definition: misc.cc:1570
See user documentation
Definition: params.h:175
See user documentation
Definition: params.h:174
#define PI
What do you think it is?
Definition: inlines.h:34

◆ mf_sheth_tormen()

double mf_sheth_tormen ( double const &  nu,
double const &  z,
double const &  Delta_c 
)

Definition at line 2211 of file cosmo.cc.

2212 {
2213 // Halo multiplicity function per log(nu): nu*f(nu)
2214 // \int_{-infty}^infty mf(lnnu) dlnnu = 1
2215 // according to Sheth and Tormen (1999)
2216 
2217  double Delta_m = delta_crit_to_delta_x(178, kRHO_MEAN, z);
2218  if ((Delta_c < 177.5 or Delta_c > 178.5) and !gSIM_EXTRAGAL_IS_DELTA_HMF_FROM_PROFILE) {
2219  printf("\n\n====> ERROR: mf_sheth_tormen() in cosmo.cc");
2220  printf("\n Sheth & Tormen HMF parametrization is ");
2222  printf("\n given for Delta_m(z) = %d (Delta_c = 178)", int(round(Delta_m)));
2223  } else if (gCOSMO_FLAG_DELTA_REF == kBRYANNORMAN98) {
2224  double Delta_bn = delta_crit_to_delta_x(178, kBRYANNORMAN98, z);
2225  printf("\n given for Delta_Bryan() = %d (Delta_c = 178)", int(round(Delta_bn)));
2226  } else
2227  printf("\n given for values Delta_c = 178 (Delta_m(z) = %d)", int(round(Delta_m)));
2228  printf("\n If you want to use a different Delta,");
2229  printf("\n you must choose gSIM_EXTRAGAL_IS_DELTA_HMF_FROM_PROFILE = 1.");
2230  printf("\n => abort()\n\n");
2231  abort();
2232  }
2233 
2234  double q = 0.707;
2235  double Ap = 0.3222;
2236  double p = 0.3;
2237 
2238  return Ap / sqrt(2.*PI) * (1. + pow(q * nu, -p)) * sqrt(q * nu) * exp(-q * nu / 2.);
2239 }
bool gSIM_EXTRAGAL_IS_DELTA_HMF_FROM_PROFILE
See the user documentation
Definition: params.cc:240
int gCOSMO_FLAG_DELTA_REF
See the user documentation
Definition: params.cc:72
double delta_crit_to_delta_x(double const &Delta_c, int const card_delta_ref, double const &z)
Definition: clumps.cc:182
double round(const double &x, const int digits)
Definition: misc.cc:1570
See user documentation
Definition: params.h:175
See user documentation
Definition: params.h:174
#define PI
What do you think it is?
Definition: inlines.h:34

◆ mf_tinker()

double mf_tinker ( double const &  nu,
double const &  z,
double const &  Delta_c 
)

Definition at line 1969 of file cosmo.cc.

1970 {
1971  // Halo multiplicity function per log(nu): 0.5*f(sigma),
1972  // where sigma=deltac/sqrt(nu).
1973  // \int_{-infty}^infty mf(lnnu,z) dlnnu = 1 is NOT satisfied!
1974  // Ref: Eq.(3) of Tinker et al. (2008)
1975 
1976  double Delta_c_min = delta_x_to_delta_crit(200, kRHO_MEAN, z);
1977  double Delta_c_max = delta_x_to_delta_crit(3200, kRHO_MEAN, z);
1978 
1979  if (Delta_c < Delta_c_min) {
1980  printf("\n\n====> ERROR: mf_tinker() in cosmo.cc");
1981  printf("\n Tinker et al. (2008) HMF parametrization is only");
1983  printf("\n given for values Delta_c >= %d", int(round(Delta_c_min)));
1984  } else if (gCOSMO_FLAG_DELTA_REF == kRHO_MEAN) {
1985  double Delta_bn = delta_crit_to_delta_x(Delta_c_min, kBRYANNORMAN98, z);
1986  printf("\n given for values Delta_Bryan >= %d", int(round(Delta_bn)));
1987  } else
1988  printf("\n given for values Delta_m >= 200");
1989  printf("\n If you want to use a smaller Delta,");
1990  printf("\n you must choose gSIM_EXTRAGAL_IS_DELTA_HMF_FROM_PROFILE = 1.");
1991  printf("\n => abort()\n\n");
1992  abort();
1993  } else if (Delta_c > Delta_c_max) {
1994  printf("\n\n====> ERROR: mf_tinker() in cosmo.cc");
1995  printf("\n Tinker et al. (2008) HMF parametrization is only");
1996  printf("\n given for values Delta_m <= 3200");
1997  printf("\n If you want to use a larger Delta,");
1998  printf("\n you must choose gSIM_EXTRAGAL_IS_DELTA_HMF_FROM_PROFILE = 1.");
1999  printf("\n => abort()\n\n");
2000  abort();
2001  }
2002 
2003  double deltas[9] = {
2009  delta_x_to_delta_crit(1200, kRHO_MEAN, z),
2010  delta_x_to_delta_crit(1600, kRHO_MEAN, z),
2011  delta_x_to_delta_crit(2400, kRHO_MEAN, z),
2013  };
2014  double params[4][9] = {
2015  {0.186, 0.200, 0.212, 0.218, 0.248, 0.255, 0.260, 0.260, 0.260},
2016  {1.47, 1.52, 1.56, 1.61, 1.87, 2.13, 2.30, 2.53, 2.66},
2017  {2.57, 2.25, 2.05, 1.87, 1.59, 1.51, 1.46, 1.44, 1.41},
2018  {1.19, 1.27, 1.34, 1.45, 1.58, 1.80, 1.97, 2.24, 2.44},
2019  };
2020 
2021  vector<double> deltas_vec;
2022  deltas_vec.assign(deltas, deltas + 9);
2023  vector < vector<double> > params_vec;
2024  pushback_array2vector(params_vec, &params[0][0], sizeof params[0]);
2025  pushback_array2vector(params_vec, &params[1][0], sizeof params[1]);
2026  pushback_array2vector(params_vec, &params[2][0], sizeof params[2]);
2027  pushback_array2vector(params_vec, &params[3][0], sizeof params[3]);
2028 
2029  double Ap0 = interp1D(Delta_c, deltas_vec, params_vec[0], kLINLIN);
2030  double a0 = interp1D(Delta_c, deltas_vec, params_vec[1], kLINLIN);
2031  double b0 = interp1D(Delta_c, deltas_vec, params_vec[2], kLINLIN);
2032  double c = interp1D(Delta_c, deltas_vec, params_vec[3], kLINLIN);
2033 
2034  double Ap, a, b, alpha;
2035  double Delta_m = delta_crit_to_delta_x(Delta_c, kRHO_MEAN, z);
2036 
2037  Ap = Ap0 * pow(1. + z, -0.14);
2038  a = a0 * pow(1. + z, -0.06);
2039  alpha = exp(-pow(0.75 / log(Delta_m / 75.), 1.2));
2040  b = b0 * pow(1. + z, -alpha);
2041 
2042  double sigma = DELTA_COLLAPSE / sqrt(nu);
2043  return 0.5 * (Ap * (pow(sigma / b, -a) + 1.) * exp(-c / pow(sigma, 2.)));
2044 }
#define DELTA_COLLAPSE
Definition: inlines.h:8
void pushback_array2vector(vector< vector< T > > &vec, const T *array, const int size)
Definition: misc.h:104
See user documentation
Definition: params.h:173
int gCOSMO_FLAG_DELTA_REF
See the user documentation
Definition: params.cc:72
double delta_crit_to_delta_x(double const &Delta_c, int const card_delta_ref, double const &z)
Definition: clumps.cc:182
double delta_x_to_delta_crit(double const &Delta_x0, int const card_delta_ref, double const &z)
Definition: clumps.cc:102
double round(const double &x, const int digits)
Definition: misc.cc:1570
See user documentation
Definition: params.h:175
See user documentation
Definition: params.h:174
x-y linear interpolation
Definition: params.h:180
double interp1D(double const &x_eval, const vector< double > &x_base, const vector< double > &y_base, const int card_interp, bool no_error=false)
Definition: misc.cc:433

◆ mf_tinker10()

double mf_tinker10 ( double const &  nu,
double const &  z,
double const &  Delta_c 
)

Definition at line 2122 of file cosmo.cc.

2123 {
2124  // Halo multiplicity function per log(nu)
2125  // where sigma=deltac/sqrt(nu).
2126  // \int_{-infty}^infty mf(lnnu,z) dlnnu = 1 is satisfied!
2127  // Ref: Eq.(8) of Tinker et al. (2010)
2128  // Delta=500
2129 
2130  // NB: They define nu=deltac/sigma while we have nu=(deltac/sigma)^2
2131  // ==> what we call nu corresponds to their nu^2, hence the power of 2
2132  // missing here and there when comparing what we return to Eq.(8) of Tinker+ (2010)
2133  // Also, this implies f(nu_us)=0.5*f(nu_them)
2134 
2135  double zz = 0.;
2136  if (z < 3.) zz = z;
2137  else zz = 3.; // following recommendation in Tinker+ (2010)
2138 
2139  double Delta_c_min = delta_x_to_delta_crit(200, kRHO_MEAN, z);
2140  double Delta_c_max = delta_x_to_delta_crit(3200, kRHO_MEAN, z);
2141 
2142  if (Delta_c < Delta_c_min) {
2143  printf("\n\n====> ERROR: mf_tinker() in cosmo.cc");
2144  printf("\n Tinker et al. (2010) HMF parametrization is only");
2146  printf("\n given for values Delta_c >= %d", int(round(Delta_c_min)));
2147  } else if (gCOSMO_FLAG_DELTA_REF == kRHO_MEAN) {
2148  double Delta0_bn = delta_crit_to_delta_x(Delta_c_min, kBRYANNORMAN98, z);
2149  printf("\n given for values Delta_Bryan >= %d", int(round(Delta0_bn)));
2150  } else
2151  printf("\n given for values Delta_m >= 200");
2152  printf("\n If you want to use a smaller Delta,");
2153  printf("\n you must choose gSIM_EXTRAGAL_IS_DELTA_HMF_FROM_PROFILE = 1.");
2154  printf("\n => abort()\n\n");
2155  abort();
2156  } else if (Delta_c > Delta_c_max) {
2157  printf("\n\n====> ERROR: mf_tinker() in cosmo.cc");
2158  printf("\n Tinker et al. (2010) HMF parametrization is only");
2159  printf("\n given for values Delta_m <= 3200");
2160  printf("\n If you want to use a larger Delta,");
2161  printf("\n you must choose gSIM_EXTRAGAL_IS_DELTA_HMF_FROM_PROFILE = 1.");
2162  printf("\n => abort()\n\n");
2163  abort();
2164  }
2165 
2166  double deltas[9] = {
2172  delta_x_to_delta_crit(1200, kRHO_MEAN, z),
2173  delta_x_to_delta_crit(1600, kRHO_MEAN, z),
2174  delta_x_to_delta_crit(2400, kRHO_MEAN, z),
2176  };
2177  double params[5][9] = {
2178  {0.368, 0.363, 0.385, 0.389, 0.393, 0.365, 0.379, 0.355, 0.327 },
2179  {0.589, 0.585, 0.544, 0.543, 0.564, 0.623, 0.637, 0.673, 0.702 },
2180  {0.864, 0.922, 0.987, 1.090, 1.200, 1.340, 1.500, 1.680, 1.810 },
2181  { -0.729, -0.789, -0.910, -1.05, -1.20, -1.26, -1.45, -1.5, -1.49},
2182  { -0.243, -0.261, -0.261, -0.273, -0.278, -0.301, -0.301, -0.319, -0.336}
2183  };
2184 
2185  vector<double> deltas_vec;
2186  deltas_vec.assign(deltas, deltas + 9);
2187  vector < vector<double> > params_vec;
2188  pushback_array2vector(params_vec, &params[0][0], sizeof params[0]);
2189  pushback_array2vector(params_vec, &params[1][0], sizeof params[1]);
2190  pushback_array2vector(params_vec, &params[2][0], sizeof params[2]);
2191  pushback_array2vector(params_vec, &params[3][0], sizeof params[3]);
2192  pushback_array2vector(params_vec, &params[4][0], sizeof params[4]);
2193 
2194  double alpha = interp1D(Delta_c, deltas_vec, params_vec[0], kLINLIN);
2195  double beta0 = interp1D(Delta_c, deltas_vec, params_vec[1], kLINLIN);
2196  double gamma0 = interp1D(Delta_c, deltas_vec, params_vec[2], kLINLIN);
2197  double phi0 = interp1D(Delta_c, deltas_vec, params_vec[3], kLINLIN);
2198  double eta0 = interp1D(Delta_c, deltas_vec, params_vec[4], kLINLIN);
2199 
2200  double beta, gamma, phi, eta;
2201 
2202  beta = beta0 * pow(1. + zz, 0.2);
2203  gamma = gamma0 * pow(1. + zz, -0.01);
2204  phi = phi0 * pow(1. + zz, -0.08);
2205  eta = eta0 * pow(1. + zz, 0.27);
2206 
2207  return 0.5 * alpha * (1. + pow(beta, -2.*phi) * pow(nu, -phi)) * pow(nu, eta + 1. / 2.) * exp(-gamma * nu / 2.);
2208 }
void pushback_array2vector(vector< vector< T > > &vec, const T *array, const int size)
Definition: misc.h:104
See user documentation
Definition: params.h:173
int gCOSMO_FLAG_DELTA_REF
See the user documentation
Definition: params.cc:72
double delta_crit_to_delta_x(double const &Delta_c, int const card_delta_ref, double const &z)
Definition: clumps.cc:182
double delta_x_to_delta_crit(double const &Delta_x0, int const card_delta_ref, double const &z)
Definition: clumps.cc:102
double round(const double &x, const int digits)
Definition: misc.cc:1570
See user documentation
Definition: params.h:175
See user documentation
Definition: params.h:174
x-y linear interpolation
Definition: params.h:180
double interp1D(double const &x_eval, const vector< double > &x_base, const vector< double > &y_base, const int card_interp, bool no_error=false)
Definition: misc.cc:433

◆ mf_tinker_norm()

double mf_tinker_norm ( double const &  nu,
double const &  z,
double const &  Delta_c 
)

Definition at line 2047 of file cosmo.cc.

2048 {
2049  // Halo multiplicity function per log(nu): 0.5*f(sigma),
2050  // where sigma=deltac/sqrt(nu).
2051  // \int_{-infty}^infty mf(lnnu,z) dlnnu = 1
2052  // Coefficients correspond to Delta=200
2053  // Ref: Tinker et al. (2008) appendix C parametrisation. No z dependence.
2054  // This parametrisation IS normalised so that it does not diverge when going down to very low masses.
2055 
2056  double Delta_c_min = delta_x_to_delta_crit(200, kRHO_MEAN, z);
2057  double Delta_c_max = delta_x_to_delta_crit(3200, kRHO_MEAN, z);
2058 
2059  if (Delta_c < Delta_c_min) {
2060  printf("\n\n====> ERROR: mf_tinker() in cosmo.cc");
2061  printf("\n Tinker et al. (2008) HMF parametrization is only");
2063  printf("\n given for values Delta_c >= %d", int(round(Delta_c_min)));
2064  } else if (gCOSMO_FLAG_DELTA_REF == kRHO_MEAN) {
2065  double Delta_bn = delta_crit_to_delta_x(Delta_c_min, kBRYANNORMAN98, z);
2066  printf("\n given for values Delta_Bryan >= %d", int(round(Delta_bn)));
2067  } else
2068  printf("\n given for values Delta_m >= 200");
2069  printf("\n If you want to use a smaller Delta,");
2070  printf("\n you must choose gSIM_EXTRAGAL_IS_DELTA_HMF_FROM_PROFILE = 1.");
2071  printf("\n => abort()\n\n");
2072  abort();
2073  } else if (Delta_c > Delta_c_max) {
2074  printf("\n\n====> ERROR: mf_tinker() in cosmo.cc");
2075  printf("\n Tinker et al. (2008) HMF parametrization is only");
2076  printf("\n given for values Delta_m <= 3200");
2077  printf("\n If you want to use a larger Delta,");
2078  printf("\n you must choose gSIM_EXTRAGAL_IS_DELTA_HMF_FROM_PROFILE = 1.");
2079  printf("\n => abort()\n\n");
2080  abort();
2081  }
2082 
2083  double deltas[9] = {
2089  delta_x_to_delta_crit(1200, kRHO_MEAN, z),
2090  delta_x_to_delta_crit(1600, kRHO_MEAN, z),
2091  delta_x_to_delta_crit(2400, kRHO_MEAN, z),
2093  };
2094  double params[4][9] = {
2095  {1.97, 2.06, 2.30, 2.56, 2.83, 2.92, 3.29, 3.37, 3.30 },
2096  {1.00, 0.99, 0.93, 0.93, 0.96, 1.04, 1.07, 1.12, 1.16 },
2097  {0.51, 0.48, 0.48, 0.45, 0.44, 0.40, 0.40, 0.36, 0.33 },
2098  {1.228, 1.310, 1.403, 1.553, 1.702, 1.907, 2.138, 2.394, 2.572},
2099  };
2100 
2101  vector<double> deltas_vec;
2102  deltas_vec.assign(deltas, deltas + 9);
2103  vector < vector<double> > params_vec;
2104  pushback_array2vector(params_vec, &params[0][0], sizeof params[0]);
2105  pushback_array2vector(params_vec, &params[1][0], sizeof params[1]);
2106  pushback_array2vector(params_vec, &params[2][0], sizeof params[2]);
2107  pushback_array2vector(params_vec, &params[3][0], sizeof params[3]);
2108 
2109 
2110  double dd = interp1D(Delta_c, deltas_vec, params_vec[0], kSPLINE);
2111  double ee = interp1D(Delta_c, deltas_vec, params_vec[1], kSPLINE);
2112  double ff = interp1D(Delta_c, deltas_vec, params_vec[2], kSPLINE);
2113  double gg = interp1D(Delta_c, deltas_vec, params_vec[3], kSPLINE);
2114 
2115  double Bp = 2. / (pow(ee, dd) * pow(gg, -dd / 2.) * tgamma(dd / 2.) + pow(gg, -ff / 2.) * tgamma(ff / 2.));
2116 
2117  double sigma = DELTA_COLLAPSE / sqrt(nu);
2118  return 0.5 * Bp * (pow(sigma / ee, -dd) + pow(sigma, -ff)) * exp(-gg / (sigma * sigma));
2119 }
#define DELTA_COLLAPSE
Definition: inlines.h:8
void pushback_array2vector(vector< vector< T > > &vec, const T *array, const int size)
Definition: misc.h:104
See user documentation
Definition: params.h:173
int gCOSMO_FLAG_DELTA_REF
See the user documentation
Definition: params.cc:72
Spline interpolation.
Definition: params.h:184
double delta_crit_to_delta_x(double const &Delta_c, int const card_delta_ref, double const &z)
Definition: clumps.cc:182
double delta_x_to_delta_crit(double const &Delta_x0, int const card_delta_ref, double const &z)
Definition: clumps.cc:102
double round(const double &x, const int digits)
Definition: misc.cc:1570
See user documentation
Definition: params.h:175
See user documentation
Definition: params.h:174
double interp1D(double const &x_eval, const vector< double > &x_base, const vector< double > &y_base, const int card_interp, bool no_error=false)
Definition: misc.cc:433

◆ nu_andderiv()

vector< vector<double> > nu_andderiv ( const double  z,
const vector< double > &  Mh_vec,
const vector< double > &  linear_lnk_vec,
const vector< double > &  linear_lnp_vec,
const int  card_window,
const int  card_growth_mthd 
)

Definition at line 2242 of file cosmo.cc.

2245 {
2246  // Computes nu(Mh) and d ln(nu) / d ln(Mh) at redshift z for nu = nu=(delta_collapse/sigma)^2
2247  // Note 1: Mh is in units of Omega_m0 h^-1 Msol
2248  // Note 2: The masses Mh reflect all the masses contained in a collapsed region
2249 
2250  // Input:
2251  // z redshift z
2252  // Mh_vec mass values, length n_m, unit [Omega_m0 h^-1 Msol]
2253  // or P(k,0) is given, and sigma^2(z) is calculated from a growth factor model.
2254  // linear_lnk_grid k values for p_k
2255  // linear_lnplnk_z lnp/lnk, must have same dimension as linear_lnk_grid
2256  // card_window Real-space window for collapse region (kTOP_HAT, ..)
2257  // card_growth_mthd selects the method to obtain P(k,z): either P(k,z) is directly given (kPKZ_FROMFILE),
2258 
2259  // Output:
2260  // res vector of two vectors with each of length len(Mh_vec), containing:
2261  // res[0]: mu(Mh, z)
2262  // res[1]: d ln(nu) / d ln(Mh) (Mh, z)
2263 
2264  vector< vector<double> > res;
2265  vector< vector<double> > s2 = sigma2_andderiv(z, Mh_vec, linear_lnk_vec, linear_lnp_vec, card_window, card_growth_mthd);
2266  vector<double> nu, dlnnudlnMh;
2267  double delta2 = pow(DELTA_COLLAPSE, 2);
2268  for (int i = 0; i < (int)s2[0].size(); i++) {
2269  nu.push_back(delta2 / s2[0][i]); // nu
2270  dlnnudlnMh.push_back(-s2[1][i]); // d ln(nu) / d ln(Mh)
2271  }
2272  res.push_back(nu);
2273  res.push_back(dlnnudlnMh);
2274  return res;
2275 }
#define DELTA_COLLAPSE
Definition: inlines.h:8
vector< vector< double > > sigma2_andderiv(const double z, const vector< double > &Mh_vec, const vector< double > &linear_lnk_vec, const vector< double > &linear_lnp_vec, const int card_window, const int card_growth_mthd)
Definition: cosmo.cc:2355

◆ Omega_b()

double Omega_b ( const double &  z)

Definition at line 2326 of file cosmo.cc.

2327 {
2328  return gCOSMO_OMEGA0_B / H2_over_H02(z) * pow((1 + z), 3.);
2329 }
double gCOSMO_OMEGA0_B
See the user documentation
Definition: params.cc:64
double H2_over_H02(const double &z)
Definition: cosmo.cc:2278

◆ Omega_cdm()

double Omega_cdm ( const double &  z)

Definition at line 2320 of file cosmo.cc.

2321 {
2322  return (gCOSMO_OMEGA0_M - gCOSMO_OMEGA0_B) / H2_over_H02(z) * pow((1 + z), 3.);
2323 }
double gCOSMO_OMEGA0_B
See the user documentation
Definition: params.cc:64
double H2_over_H02(const double &z)
Definition: cosmo.cc:2278
double gCOSMO_OMEGA0_M
See the user documentation
Definition: params.cc:63

◆ Omega_k()

double Omega_k ( const double &  z)

Definition at line 2342 of file cosmo.cc.

2343 {
2344  return gCOSMO_OMEGA0_K / H2_over_H02(z) * pow((1 + z), 2.);
2345 }
double gCOSMO_OMEGA0_K
See the user documentation
Definition: params.cc:65
double H2_over_H02(const double &z)
Definition: cosmo.cc:2278

◆ Omega_lambda()

double Omega_lambda ( const double &  z)

Definition at line 2332 of file cosmo.cc.

2333 {
2334  double x = (1. + z);
2335  double fac_lambda;
2336  if (abs(gCOSMO_WDE + 1.) > 1e-5) fac_lambda = pow(x, (3. + 3.*gCOSMO_WDE));
2337  else fac_lambda = 1.;
2338  return gCOSMO_OMEGA0_LAMBDA / H2_over_H02(z) * fac_lambda;
2339 }
double gCOSMO_OMEGA0_LAMBDA
Present day dark energy density of the -CDM Universe (PDG)
Definition: params.cc:77
double gCOSMO_WDE
See the user documentation
Definition: params.cc:69
double H2_over_H02(const double &z)
Definition: cosmo.cc:2278

◆ Omega_m()

double Omega_m ( const double &  z)

Definition at line 2314 of file cosmo.cc.

2315 {
2316  return gCOSMO_OMEGA0_M / H2_over_H02(z) * pow((1 + z), 3.);
2317 }
double H2_over_H02(const double &z)
Definition: cosmo.cc:2278
double gCOSMO_OMEGA0_M
See the user documentation
Definition: params.cc:63

◆ Omega_r()

double Omega_r ( const double &  z)

Definition at line 2308 of file cosmo.cc.

2309 {
2310  return gCOSMO_OMEGA0_R / H2_over_H02(z) * pow((1 + z), 4.);
2311 }
double gCOSMO_OMEGA0_R
Present day radiation density of the Universe (PDG)
Definition: params.cc:76
double H2_over_H02(const double &z)
Definition: cosmo.cc:2278

◆ rho_crit()

double rho_crit ( const double &  z)

Definition at line 2348 of file cosmo.cc.

2349 {
2350  // returns rho_crit in units of Msol kpc^-3
2351  return gCOSMO_RHO0_C * H2_over_H02(z);
2352 }
double H2_over_H02(const double &z)
Definition: cosmo.cc:2278
double gCOSMO_RHO0_C
Present day critical density of the universe in unit of (PDG)
Definition: params.cc:75

◆ sigma2_andderiv()

vector< vector<double> > sigma2_andderiv ( const double  z,
const vector< double > &  Mh_vec,
const vector< double > &  linear_lnk_vec,
const vector< double > &  linear_lnp_vec,
const int  card_window,
const int  card_growth_mthd 
)

Definition at line 2355 of file cosmo.cc.

2358 {
2359  // Computes sigma^2(Mh) and d ln(sigma^2) / d ln(Mh) at redshift z.
2360  // Note 1: Mh is in units of Omega_m0 h^-1 Msol
2361  // Note 2: The masses Mh reflect all the masses contained in a collapsed region
2362 
2363  // Input:
2364  // z redshift z
2365  // Mh_vec mass values, length n_m, unit [Omega_m0 h^-1 Msol]
2366  // or P(k,0) is given, and sigma^2(z) is calculated from a growth factor model.
2367  // linear_lnk_grid k values for p_k
2368  // linear_lnplnk_z lnp/lnk, must have same dimension as linear_lnk_grid
2369  // card_window Real-space window for collapse region (kTOP_HAT, ..)
2370  // card_growth_mthd selects the method to obtain P(k,z): either P(k,z) is directly given (kPKZ_FROMFILE),
2371 
2372  // Output:
2373  // res vector of two vectors with each of length len(Mh_vec), containing:
2374  // res[0]: sigma^2(Mh, z)
2375  // res[1]: d ln(sigma^2) / d ln(Mh) (Mh, z)
2376 
2377  double Rh, lnRh, lns2, dlns2dlnRh, dlns2dlnMh;
2378  double rho_critperh2_MsolperMpc3 = RHO_CRITperh2_MSOLperKPC3 / KPC3_to_MPC3 ; // * in units of h^2 * Msol Mpc^-3
2379 
2380  vector<double> sigma2_vec;
2381  vector<double> dlns2dlnMh_vec;
2382  vector< vector<double> > res;
2383 
2384  bool is_rodriguez_sigmaapprox = false;
2385  if (gCOSMO_HUBBLE == 0.678 and
2386  gCOSMO_OMEGA0_M == 0.307 and
2387  gCOSMO_OMEGA0_B == 0.048 and
2388  gCOSMO_SIGMA8 == 0.829 and
2389  gCOSMO_N_S == 0.96 and
2390  gCOSMO_DELTA0 == -1 and
2392  card_growth_mthd != kPKZ_FROMFILE) {
2393  static bool is_rodriguez_sigmaapprox_msg = true; // have switched off analytic formula for now (chenge to false if you want to use it)
2394  is_rodriguez_sigmaapprox = false; // have switched off analytic formula for now (change to true if you want to use it)
2395  if (!is_rodriguez_sigmaapprox_msg) {
2396  printf("\n\n====> NOTE: compute_sigma2() in cosmo.cc");
2397  printf("\n Analytic approximation for sigma(Mvir)");
2398  printf("\n from Rodriguez-Puebla (2016) is used.\n");
2399  is_rodriguez_sigmaapprox_msg = true;
2400  }
2401  }
2402 
2403  // Following Komatsu, first compute the Chebyshev serie of lnsigma2 and its derivative.
2404  // This significantly speeds up the computation compared to performing the evaluation
2405  // and numerical derivative of lnsigma2 for each mass. We checked that the 2 approaches
2406  // give the same results.
2407 
2408  gsl_cheb_series *cs = gsl_cheb_alloc(30);
2409  gsl_cheb_series *deriv = gsl_cheb_alloc(30);
2410 
2411  if (!is_rodriguez_sigmaapprox) {
2412  pk_params_for_integration params = {0., linear_lnk_vec, linear_lnp_vec, card_window};
2413  gsl_function F;
2414  F.function = &lnsigma2;
2415  F.params = &params;
2416 
2417  double lnRhmin = log(pow(Mh_vec[0] / gamma_f(card_window) / rho_critperh2_MsolperMpc3, 1. / 3.));
2418  double lnRhmax = log(pow(Mh_vec[Mh_vec.size() - 1] / gamma_f(card_window) / rho_critperh2_MsolperMpc3, 1. / 3.));
2419  gsl_cheb_init(cs, &F, lnRhmin, lnRhmax); // cs contains the Chebyshev coefficients
2420  gsl_cheb_calc_deriv(deriv, cs); // deriv contains the derivatives coefficients
2421  }
2422 
2423  // Compute growth factor (only used when card_growth_mthd != kPKZ_FROMFILE)
2424  double log_growthfac = 0.;
2425  if (card_growth_mthd != kPKZ_FROMFILE) {
2426  log_growthfac = log(growth_factor(z, card_growth_mthd));
2427  }
2428 
2429  //Now, loop on all masses
2430  for (int j = 0; j < (int)Mh_vec.size(); j++) {
2431 
2432  if (is_rodriguez_sigmaapprox) {
2433  double Mh = Mh_vec[j] * gCOSMO_OMEGA0_M;
2434  double lnsigma = log(17.111 * pow(1e12 / Mh, 0.405) / (1 + 1.306 * pow(1e12 / Mh, 0.22) + 6.218 * pow(1e12 / Mh, 0.317)));
2435  lns2 = 2 * lnsigma + 2. * log_growthfac;
2436  dlns2dlnMh = -2. * (0.088 + (0.317 + 55.2987 * pow(Mh, -0.22)) /
2437  (1. + 570.09 * pow(Mh, -0.22) + 39595.9 * pow(Mh, -0.317)));
2438  } else {
2439  Rh = pow(Mh_vec[j] / gamma_f(card_window) / (rho_critperh2_MsolperMpc3), 1. / 3.); // in units of h^-1 Mpc
2440  lnRh = log(Rh);
2441  if (exp(linear_lnk_vec[linear_lnk_vec.size() - 1]) < 2.e1 / Rh) {
2442  cout << j << " Mass scale Mh = " << Mh_vec[j] << "<--> Rh = " << Rh <<
2443  " is too small for available k range --> abort!" << endl;
2444  cout << "Ask for larger P_k_max_h/Mpc when running CLASS" << endl;
2445  abort();
2446  }
2447 
2448  // evaluate lnsigma2 at lnRh using Chebyshev coefficients.
2449  // We had calculated ln(2pi^2/D(z)^2 sigma^2), so multiply with D(z)^2/(2pi^2) here
2450  lns2 = gsl_cheb_eval(cs, lnRh) + 2. * log_growthfac - log(2) - 2. * log(PI);
2451  // evaluates derivative at lnRh using Chebyshev coefficients.
2452  // Derivation in log space -> D(z)^2/(2pi^2) dissappears.
2453  dlns2dlnRh = gsl_cheb_eval(deriv, lnRh);
2454  dlns2dlnMh = dlns2dlnRh / 3.;
2455  }
2456  sigma2_vec.push_back(exp(lns2));
2457  dlns2dlnMh_vec.push_back(dlns2dlnMh);
2458  }
2459 
2460  if (!is_rodriguez_sigmaapprox) {
2461  gsl_cheb_free(cs);
2462  gsl_cheb_free(deriv);
2463  }
2464  res.push_back(sigma2_vec);
2465  res.push_back(dlns2dlnMh_vec);
2466  return res;
2467 }
double gCOSMO_SIGMA8
See the user documentation
Definition: params.cc:66
double gCOSMO_HUBBLE
See the user documentation
Definition: params.cc:62
double gCOSMO_OMEGA0_B
See the user documentation
Definition: params.cc:64
double gCOSMO_DELTA0
See the user documentation
Definition: params.cc:71
int gCOSMO_FLAG_DELTA_REF
See the user documentation
Definition: params.cc:72
#define KPC3_to_MPC3
Definition: inlines.h:23
See user documentation
Definition: params.h:175
double lnsigma2(double lnRh, void *p)
Returns for the scale passed as .
Definition: cosmo.cc:1430
double gCOSMO_N_S
See the user documentation
Definition: params.cc:67
See user documentation
Definition: params.h:189
#define PI
What do you think it is?
Definition: inlines.h:34
double gamma_f(const int card_window)
Definition: cosmo.cc:1479
double growth_factor(const double &z, const int card_growth_mthd)
Computes the linear growth rate factor.
Definition: cosmo.cc:1652
double gCOSMO_OMEGA0_M
See the user documentation
Definition: params.cc:63
#define RHO_CRITperh2_MSOLperKPC3
Critical density of the universe in units of .
Definition: inlines.h:7

◆ sigma8()

double sigma8 ( vector< double > &  lnk_vec,
vector< double > &  lnp_vec,
const int  card_window 
)

Definition at line 2470 of file cosmo.cc.

2471 {
2472  // - Returns sigma8 from given power spectrum
2473  gsl_function lnsigma2_func;
2474  lnsigma2_func.function = &lnsigma2;
2475  pk_params_for_integration params = { -999, lnk_vec, lnp_vec, card_window};
2476  lnsigma2_func.params = &params;
2477  return 1 / (sqrt(2) * PI) * exp(0.5 * GSL_FN_EVAL(&lnsigma2_func, log(8)));
2478 }
double lnsigma2(double lnRh, void *p)
Returns for the scale passed as .
Definition: cosmo.cc:1430
#define PI
What do you think it is?
Definition: inlines.h:34

◆ solve_d_to_z()

double solve_d_to_z ( double  d,
void *  p 
)

Definition at line 2481 of file cosmo.cc.

2482 {
2483  //--- GSL Function to solve the equation
2484  //
2485  // d(z0) = d0
2486  //
2487  // for z0.
2488 
2490  int switch_distance = params->switch_distance;
2491  double distance = params->distance;
2492  double d = 0.;
2493  if (switch_distance == 0)
2494  d = dh_c(z) / gCOSMO_HUBBLE;
2495  else if (switch_distance == 1)
2496  d = dh_trans(z) / gCOSMO_HUBBLE;
2497  else if (switch_distance == 2)
2498  d = dh_l(z) / gCOSMO_HUBBLE;
2499  else if (switch_distance == 3)
2500  d = dh_a(z) / gCOSMO_HUBBLE;
2501  else {
2502  printf("\n====> ERROR: solve_d_to_z() in cosmo.cc");
2503  printf("\n switch_distance < 4 required");
2504  printf("\n => abort()\n\n");
2505  abort();
2506  }
2507 
2508  return d - distance;
2509 
2510 }
double dh_trans(const double &z)
Computes the transverse comoving distance by integration.
Definition: cosmo.cc:342
double gCOSMO_HUBBLE
See the user documentation
Definition: params.cc:62
double dh_l(const double &z)
Luminosity distance at redshift in [ ].
Definition: cosmo.cc:384
double distance
input distance [Mpc]
Definition: cosmo.h:23
int switch_distance
which distance? 0: d_c, 1: d_trans, 2: d_l, 3: d_a
Definition: cosmo.h:24
double dh_a(const double &z)
Angular diameter distance at redshift in [ ].
Definition: cosmo.cc:370
double dh_c(const double &z)
Computes the comoving distance by integration.
Definition: cosmo.cc:320

◆ solve_Mhmin()

double solve_Mhmin ( double  Mhmin,
void *  p 
)

Definition at line 2513 of file cosmo.cc.

2514 {
2515  //--- GSL Function to solve the implicit equation
2516  //
2517  // Int_0^Infinity (1+ exp[ -(M - Mmin)/Delta_M])^-1 * M * dN/dMdV dM = rho_mass
2518  //
2519  // for Mhmin.
2520 
2522 
2523  double par[3] = {params->z, Mhmin, params->sigma_m};
2524 
2525  if (gCOSMO_MH_GRID.empty()) {
2526  printf("\n====> ERROR: solve_Mmin() in cosmo.cc");
2527  printf("\n gCOSMO_MH_GRID grid is empty. Please fill it before");
2528  printf("\n via init_extragal().");
2529  printf("\n => abort()\n\n");
2530  abort();
2531  }
2532 
2533  double rhoh_halos = 0.; // result will be in units of Omega_m0 h^-1 * h^3 Msol/Mpc^3
2534  double Mhmin_grid = gCOSMO_MH_GRID[0];
2535  double Mhmax_grid = gCOSMO_MH_GRID[gCOSMO_MH_GRID.size() - 1];
2536  simpson_log_adapt(dNdVhdlnMh_sigmoid, Mhmin_grid, Mhmax_grid, par, rhoh_halos, gSIM_EPS);
2538 }
double gSIM_EPS
See the user documentation
Definition: params.cc:196
double z
redshift
Definition: cosmo.h:11
double gDM_RHOHALOES_TO_RHOMEAN
See the user documentation
Definition: params.cc:93
#define KPC3_to_MPC3
Definition: inlines.h:23
vector< double > gCOSMO_MH_GRID
Definition: params.cc:81
void simpson_log_adapt(void fn(double &, double *, double &), double &xmin, double &xmax, double par[], double &s, double const &eps, bool is_verbose=false)
Returns 1D integral of fn using Simpson rule in logarithmic steps. Adaptative steps stop when the pre...
Definition: integr.cc:60
void dNdVhdlnMh_sigmoid(double &Mh, double par[3], double &res)
Definition: cosmo.cc:457
double sigma_m
Steepness of the sigmoid transition.
Definition: cosmo.h:12
#define RHO_CRITperh2_MSOLperKPC3
Critical density of the universe in units of .
Definition: inlines.h:7

◆ window_fnct()

double window_fnct ( const double &  x,
const int  card_window 
)

Definition at line 2541 of file cosmo.cc.

2542 {
2543  //--- Sets res = window function in k-space
2544  // Adapted from van den Bosch presentation slides,
2545  // http://www.astro.yale.edu/vdbosch/astro610_lecture9.pdf
2546  // See also Percival (2001), astro-ph/0107437
2547  //
2548  // INPUTS:
2549  // x k * r
2550  // card_window Choice of window function
2551  // OUTPUT:
2552  // res Value of W(x) in k space
2553 
2554  double res;
2555  // Calculate W(x)
2556  switch (card_window) {
2557  case kTOP_HAT:
2558  res = (3. / pow(x, 3)) * (sin(x) - x * cos(x));
2559  break;
2560  case kGAUSS:
2561  res = exp(- pow(x, 2) / 2.);
2562  break;
2563  case kSHARP_K:
2564  if (x <= 1) res = 1;
2565  else res = 0;
2566  break;
2567  default :
2568  printf("\n====> ERROR: window_fnct() in cosmo.cc");
2569  printf("\n window card_window=%d does not correspond to any window function", card_window);
2570  printf("\n => abort()\n\n");
2571  abort();
2572  }
2573  return res;
2574 }
See user documentation
Definition: params.h:167
See user documentation
Definition: params.h:166
See user documentation
Definition: params.h:168