CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
Cosmology.cpp
Go to the documentation of this file.
1 /********************************************************************
2  * Copyright (C) 2010 by Federico Marulli *
3  * federico.marulli3@unibo.it *
4  * *
5  * This program is free software; you can redistribute it and/or *
6  * modify it under the terms of the GNU General Public License as *
7  * published by the Free Software Foundation; either version 2 of *
8  * the License, or (at your option) any later version. *
9  * *
10  * This program is distributed in the hope that it will be useful, *
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13  * GNU General Public License for more details. *
14  * *
15  * You should have received a copy of the GNU General Public *
16  * License along with this program; if not, write to the Free *
17  * Software Foundation, Inc., *
18  * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
19  ********************************************************************/
20 
34 #include "Cosmology.h"
35 
36 using namespace std;
37 
38 using namespace cbl;
39 using namespace glob;
40 
41 
42 // =====================================================================================
43 
44 
45 std::string cbl::cosmology::CosmologicalParameter_name (const CosmologicalParameter parameter)
46 {
47  string name;
48 
49  switch (parameter) {
50 
51  case (CosmologicalParameter::_Omega_matter_LCDM_):
52  name = "Omega_matter_LCDM";
53  break;
54 
55  case (CosmologicalParameter::_Omega_matter_):
56  name = "Omega_matter";
57  break;
58 
59  case (CosmologicalParameter::_Omega_baryon_):
60  name = "Omega_baryon";
61  break;
62 
63  case (CosmologicalParameter::_Omega_baryon_h2_):
64  name = "Omega_baryon_h2";
65  break;
66 
67  case (CosmologicalParameter::_Omega_neutrinos_):
68  name = "Omega_matter";
69  break;
70 
71  case (CosmologicalParameter::_massless_neutrinos_):
72  name = "massless_neutrinos";
73  break;
74 
75  case (CosmologicalParameter::_massive_neutrinos_):
76  name = "massive_neutrinos";
77  break;
78 
79  case (CosmologicalParameter::_neutrino_mass_):
80  name = "neutrino_mass";
81  break;
82 
83  case (CosmologicalParameter::_Omega_DE_):
84  name = "Omega_DE";
85  break;
86 
87  case (CosmologicalParameter::_Omega_radiation_):
88  name = "Omega_radiation";
89  break;
90 
91  case (CosmologicalParameter::_H0_):
92  name = "H0";
93  break;
94 
95  case (CosmologicalParameter::_hh_):
96  name = "hh";
97  break;
98 
99  case (CosmologicalParameter::_scalar_amp_):
100  name = "scalar_amp";
101  break;
102 
103  case (CosmologicalParameter::_ln_scalar_amp_):
104  name = "ln_scalar_amp";
105  break;
106 
107  case (CosmologicalParameter::_scalar_pivot_):
108  name = "scalar_pivot";
109  break;
110 
111  case (CosmologicalParameter::_n_spec_):
112  name = "n_spec";
113  break;
114 
115  case (CosmologicalParameter::_w0_):
116  name = "w0";
117  break;
118 
119  case (CosmologicalParameter::_wa_):
120  name = "wa";
121  break;
122 
123  case (CosmologicalParameter::_fNL_):
124  name = "fNL";
125  break;
126 
127  case (CosmologicalParameter::_sigma8_):
128  name = "sigma8";
129  break;
130 
131  case (CosmologicalParameter::_tau_):
132  name = "tau";
133  break;
134 
135  case (CosmologicalParameter::_rs_):
136  name = "rs";
137  break;
138 
139  default:
140  ErrorCBL("no such a variable in the list!", "CosmologicalParameter_name", "Cosmology.cpp");
141  }
142 
143  return name;
144 }
145 
146 // =====================================================================================
147 
148 
149 void cbl::cosmology::Cosmology::set_default ()
150 {
151  if (m_Omega_matter==0) ErrorCBL("Omega_matter=0!", "set_default", "Cosmology.cpp");
152 
153  m_Omega_k = 1.-m_Omega_matter-m_Omega_radiation-m_Omega_DE;
154  m_Omega_CDM = m_Omega_matter-m_Omega_baryon-m_Omega_neutrinos;
155  m_H0 = (m_unit) ? 100. : 100.*m_hh;
156  m_t_H = 1./m_H0;
157  m_D_H = par::cc*m_t_H;
158  m_RhoZero = rho_m(0.);
159  m_Pk0_EH = 1., m_Pk0_CAMB = 1., m_Pk0_MPTbreeze = 1., m_Pk0_CLASS = 1.;
160 }
161 
162 // =====================================================================================
163 
164 
165 cbl::cosmology::Cosmology::Cosmology (const double Omega_matter, const double Omega_baryon, const double Omega_neutrinos, const double massless_neutrinos, const int massive_neutrinos, const double Omega_DE, const double Omega_radiation, const double hh, const double scalar_amp, const double scalar_pivot, const double n_spec, const double w0, const double wa, const double fNL, const int type_NG, const double tau, const std::string model, const bool unit)
166  : m_Omega_matter(Omega_matter), m_Omega_baryon(Omega_baryon), m_Omega_neutrinos(Omega_neutrinos), m_massless_neutrinos(massless_neutrinos), m_massive_neutrinos(massive_neutrinos), m_Omega_DE(Omega_DE), m_Omega_radiation(Omega_radiation), m_hh(hh), m_sigma8(-1.), m_scalar_amp(scalar_amp), m_scalar_pivot(scalar_pivot), m_n_spec(n_spec), m_w0(w0), m_wa(wa), m_fNL(fNL), m_type_NG(type_NG), m_tau(tau), m_model(model), m_unit(unit)
167 { set_default(); }
168 
169 
170 // =====================================================================================
171 
172 
173 cbl::cosmology::Cosmology::Cosmology (const CosmologicalModel cosmoModel, const std::string model, const bool unit)
174  : m_Omega_neutrinos(0.), m_massless_neutrinos(3.04), m_massive_neutrinos(0.), m_sigma8(-1.), m_n_spec(0.96), m_w0(-1.), m_wa(0.), m_fNL(0.), m_type_NG(1.), m_tau(0.09), m_model(model), m_unit(unit)
175 {
176  switch(cosmoModel) {
177 
178  // Komatsu et al. 2009: Table 1, WMAP 5 Year Mean
179  case(CosmologicalModel::_WMAP5_):
180  m_Omega_radiation = 2.469e-5*(1.+0.2271*m_massless_neutrinos); // Eq.8
181  m_Omega_DE = 0.742; // Omega_DE = 0.742 ± 0.030
182  m_Omega_matter = 1.-m_Omega_radiation-m_Omega_DE; // flat LCDM (e.g. Eq. 11)
183  m_Omega_baryon = 0.0441; // Omega_b = 0.0441 ± 0.0030
184  m_hh = 0.719; // h = 0.719 ± 0.027 Km/s/Mpc
185  m_scalar_pivot = 0.002; // baseline
186  m_scalar_amp = 2.41e-9; // scalar amplitude = (2.41 ± 0.11)e-9 -> sigma8 = 0.796 ± 0.036
187  m_n_spec = 0.963; // n = 0.963 ± 0.015
188  m_tau = 0.087; // tau = 0.087 ± 0.017
189  break;
190 
191  // Komatsu et al. 2011: Table 1, WMAP Seven-year Mean
192  case(CosmologicalModel::_WMAP7_):
193  m_Omega_radiation = 2.469e-5*(1.+0.2271*m_massless_neutrinos);
194  m_Omega_DE = 0.727; // Omega_DE = 0.727 ± 0.030
195  m_Omega_matter = 1.-m_Omega_radiation-m_Omega_DE; // flat LCDM
196  m_Omega_baryon = 0.0455; // Omega_b = 0.0455 ± 0.0028
197  m_hh = 0.704; // h = 0.704 ± 0.025 Km/s/Mpc
198  m_scalar_amp = 2.43e-9; // scalar amplitude = (2.43 ± 0.11)e-9 -> sigma8 = 0.811 ± 0.031
199  m_scalar_pivot = 0.002; // baseline
200  m_n_spec = 0.968; // n = 0.968 ± 0.012
201  m_tau = 0.088; // tau = 0.088 ± 0.015
202  break;
203 
204  // Hinshaw et al. 2013: Table 3, WMAP-only Nine-year
205  case(CosmologicalModel::_WMAP9_):
206  m_Omega_baryon = 0.0463; // Omega_b = 0.0463 ± 0.0024
207  m_Omega_DE = 0.721; // Omega_DE = 0.721 ± 0.025
208  m_Omega_matter = 0.233+m_Omega_baryon; // Omega_CDM = 0.233 ± 0.023
209  m_Omega_radiation = OmegaR_zeq(3273.); // Table 2
210  m_hh = 0.70; // h = 0.700 ± 0.022 Km/s/Mpc
211  m_scalar_amp = 2.41e-9; // scalar amplitude = (2.41 ± 0.10)e-9 -> sigma8 = 0.821 ± 0.023
212  m_scalar_pivot = 0.002; // baseline
213  m_n_spec = 0.972; // n = 0.972 ± 0.013
214  m_tau = 0.089; // tau = 0.089 ± 0.014
215  break;
216 
217  // Planck Collab 2013, Paper XVI: Table 2, Planck+WP
218  case(CosmologicalModel::_Planck13_):
219  m_Omega_matter = 0.315; // Omega_M = 0.315 ± 0.018
220  m_Omega_baryon = 0.0487; // Omega_b*h^2 = 0.02205 ± 0.00028
221  m_massless_neutrinos = 2.04; // baseline (see Table 1)
222  m_massive_neutrinos = 1; // baseline (see Table 1)
223  m_Omega_radiation = OmegaR_zeq(3391.); // z_eq = 3391 ± 60
224  m_Omega_DE = 1.-m_Omega_matter-m_Omega_radiation; // assuming Omega_k = 0 (in Table 3: Omega_DE = 0.685 ± 0.018)
225  m_hh = 0.673; // h = 0.673 ± 0.012 Km/s/Mpc
226  m_scalar_amp = 2.196e-9; // scalar amplitude = (2.196 ± 0.06)e-9 -> sigma8 = 0.829 ± 0.012
227  m_scalar_pivot = 0.05; // baseline
228  m_n_spec = 0.9603; // n = 0.9603 ± 0.0073
229  m_tau = 0.089; // tau = 0.089 ± 0.014
230  m_Omega_neutrinos = Omega_neutrinos(0.06); // baseline (see Table 1)
231  break;
232 
233  // Planck Collab 2015, Paper XIII: Table 4, TT,TE,EE+lowP+lensing
234  case(CosmologicalModel::_Planck15_):
235  m_Omega_matter = 0.3121; // Omega_M = 0.3121 ± 0.0087
236  m_Omega_baryon = 0.0488; // Omega_b*h^2 = 0.02226 ± 0.00016
237  m_massless_neutrinos = 2.04; // baseline
238  m_massive_neutrinos = 1; // baseline
239  m_Omega_radiation = OmegaR_zeq(3382.); // z_eq = 3382 ± 32
240  m_Omega_DE = 1.-m_Omega_matter-m_Omega_radiation; // assuming Omega_k = 0 (in Table 4: Omega_DE = 0.6879 ± 0.0087)
241  m_hh = 0.6751; // h = 0.6751 ± 0.0064 Km/s/Mpc
242  m_scalar_amp = 2.13e-9; // scalar amplitude = (2.130 ± 0.053)e-9 -> sigma8 = 0.8150 ± 0.0087
243  m_scalar_pivot = 0.05; // baseline
244  m_n_spec = 0.9653; // n = 0.9653 ± 0.0048
245  m_tau = 0.063; // tau = 0.063 ± 0.014
246  m_Omega_neutrinos = Omega_neutrinos(0.06); // baseline
247  break;
248 
249  // Planck Collab 2015, Paper XIII: Table 4, TT+lowP+lensing
250  case(CosmologicalModel::_Planck15_TT_):
251  m_Omega_matter = 0.308; // Omega_M = 0.308 ± 0.012
252  m_Omega_baryon = 0.0484; // Omega_b*h^2 = 0.02226 ± 0.0023
253  m_massless_neutrinos = 2.04; // baseline
254  m_massive_neutrinos = 1; // baseline
255  m_Omega_radiation = OmegaR_zeq(3365.); // z_eq = 3365 ± 44
256  m_Omega_DE = 1.-m_Omega_matter-m_Omega_radiation; // assuming Omega_k = 0 (in Table 4: Omega_DE = 0.692 ± 0.012)
257  m_hh = 0.6781; // h = 0.6781 ± 0.092 Km/s/Mpc
258  m_scalar_amp = 2.139e-9; // scalar amplitude = (2.139 ± 0.063)e-9 -> sigma8 = 0.8149 ± 0.0093
259  m_scalar_pivot = 0.05; // baseline
260  m_n_spec = 0.9677; // n = 0.9677 ± 0.060
261  m_tau = 0.066; // tau = 0.066 ± 0.016
262  m_Omega_neutrinos = Omega_neutrinos(0.06); // baseline
263  break;
264 
265  // Planck Collab 2018, Paper VI: Table 2, TT,TE,EE+lowE+lensing
266  case(CosmologicalModel::_Planck18_):
267  m_Omega_matter = 0.3153; // Omega_M = 0.3153 ± 0.0073
268  m_Omega_baryon = 0.0486; // Omega_b*h^2 = 0.02237 ± 0.00015
269  m_massless_neutrinos = 2.04; // baseline
270  m_massive_neutrinos = 1; // baseline
271  m_Omega_radiation = OmegaR_zeq(3402.); // z_eq = 3402 ± 26
272  m_Omega_DE = 1.-m_Omega_matter-m_Omega_radiation; // assuming Omega_k = 0 (in Table 4: Omega_DE = 0.692 ± 0.012)
273  m_hh = 0.6736; // h = 0.6736 ± 0.054 Km/s/Mpc
274  m_scalar_amp = 2.100e-9; // scalar amplitude = (2.100 ± 0.030)e-9 -> sigma8 = 0.8111 ± 0.0060
275  m_scalar_pivot = 0.05; // baseline
276  m_n_spec = 0.9649; // n = 0.9649 ± 0.0042
277  m_tau = 0.0544; // tau = 0.0544 ± 0.0073
278  m_Omega_neutrinos = Omega_neutrinos(0.06); // baseline
279  break;
280 
281  default:
282  ErrorCBL("the chosen built-in cosmological model is not implemented", "Cosmology", "Cosmology.cpp");
283 
284  }
285 
286  set_default();
287 }
288 
289 
290 // =====================================================================================
291 
292 
293 void cbl::cosmology::Cosmology::print_parameters () const
294 {
295  coutCBL << "Omega_matter = " << m_Omega_matter << std::endl;
296  coutCBL << "Omega_baryon = " << m_Omega_baryon << std::endl;
297  coutCBL << "Omega_neutrinos = " << m_Omega_neutrinos << std::endl;
298  coutCBL << "massless_neutrinos = " << m_massless_neutrinos << std::endl;
299  coutCBL << "massive_neutrinos = " << m_massive_neutrinos << std::endl;
300  coutCBL << "Omega_DE = " << m_Omega_DE << std::endl;
301  coutCBL << "Omega_radiation = " << m_Omega_radiation << std::endl;
302  coutCBL << "Omega_k = " << m_Omega_k << std::endl;
303  coutCBL << "Omega_CDM = " << m_Omega_CDM << std::endl;
304  coutCBL << "h = " << m_hh << std::endl;
305  coutCBL << "scalar_amp = " << m_scalar_amp << std::endl;
306  coutCBL << "scalar_pivot = " << m_scalar_pivot << std::endl;
307  coutCBL << "n_spec = " << m_n_spec << std::endl;
308  coutCBL << "w0 = " << m_w0 << std::endl;
309  coutCBL << "wa = " << m_wa << std::endl;
310  coutCBL << "fNL = " << m_fNL << std::endl;
311  coutCBL << "type_NG = " << m_type_NG << std::endl;
312  coutCBL << "tau = " << m_tau << std::endl;
313 }
314 
315 
316 // =====================================================================================
317 
318 
319 double cbl::cosmology::Cosmology::value (const CosmologicalParameter parameter) const
320 {
321  double param_value;
322 
323  switch (parameter) {
324 
325  case (CosmologicalParameter::_Omega_matter_LCDM_):
326  param_value = m_Omega_matter;
327  break;
328 
329  case (CosmologicalParameter::_Omega_matter_):
330  return m_Omega_matter;
331  break;
332 
333  case (CosmologicalParameter::_Omega_baryon_):
334  param_value = m_Omega_baryon;
335  break;
336 
337  case (CosmologicalParameter::_Omega_baryon_h2_):
338  param_value = m_Omega_baryon*m_hh*m_hh;
339  break;
340 
341  case (CosmologicalParameter::_Omega_neutrinos_):
342  param_value = m_Omega_neutrinos;
343  break;
344 
345  case (CosmologicalParameter::_massless_neutrinos_):
346  param_value = m_massless_neutrinos;
347  break;
348 
349  case (CosmologicalParameter::_massive_neutrinos_):
350  param_value = m_massive_neutrinos;
351  break;
352 
353  case (CosmologicalParameter::_neutrino_mass_):
354  param_value = neutrino_mass();
355  break;
356 
357  case (CosmologicalParameter::_Omega_DE_):
358  param_value = m_Omega_DE;
359  break;
360 
361  case (CosmologicalParameter::_Omega_radiation_):
362  param_value = m_Omega_radiation;
363  break;
364 
365  case (CosmologicalParameter::_H0_):
366  param_value = m_H0;
367  break;
368 
369  case (CosmologicalParameter::_hh_):
370  param_value = m_hh;
371  break;
372 
373  case (CosmologicalParameter::_scalar_amp_):
374  param_value = m_scalar_amp;
375  break;
376 
377  case (CosmologicalParameter::_ln_scalar_amp_):
378  param_value = log(1.e10*m_scalar_amp);
379  break;
380 
381  case (CosmologicalParameter::_scalar_pivot_):
382  param_value = m_scalar_pivot;
383  break;
384 
385  case (CosmologicalParameter::_n_spec_):
386  param_value = m_n_spec;
387  break;
388 
389  case (CosmologicalParameter::_w0_):
390  param_value = m_w0;
391  break;
392 
393  case (CosmologicalParameter::_wa_):
394  param_value = m_wa;
395  break;
396 
397  case (CosmologicalParameter::_fNL_):
398  param_value = m_fNL;
399  break;
400 
401  case (CosmologicalParameter::_sigma8_):
402  param_value = m_sigma8;
403  break;
404 
405  case (CosmologicalParameter::_tau_):
406  param_value = m_tau;
407  break;
408 
409  case (CosmologicalParameter::_rs_):
410  param_value = m_rs;
411  break;
412 
413  default:
414  ErrorCBL("no such a variable in the list!", "value", "Cosmology.cpp");
415  }
416 
417  return param_value;
418 }
419 
420 
421 // =====================================================================================
422 
423 
424 void cbl::cosmology::Cosmology::set_parameter (const CosmologicalParameter parameter, const double value)
425 {
426  switch (parameter) {
427 
428  case (CosmologicalParameter::_Omega_matter_LCDM_):
429  set_Omega(value);
430  break;
431 
432  case (CosmologicalParameter::_Omega_matter_):
433  set_OmegaM(value);
434  break;
435 
436  case (CosmologicalParameter::_Omega_baryon_):
437  set_OmegaB(value);
438  break;
439 
440  case (CosmologicalParameter::_Omega_baryon_h2_):
441  set_OmegaB_h2(value);
442  break;
443 
444  case (CosmologicalParameter::_Omega_neutrinos_):
445  set_OmegaNu(value, m_massless_neutrinos, m_massive_neutrinos);
446  break;
447 
448  case (CosmologicalParameter::_massless_neutrinos_):
449  set_OmegaNu(m_Omega_neutrinos, value, m_massive_neutrinos);
450  break;
451 
452  case (CosmologicalParameter::_massive_neutrinos_):
453  set_OmegaNu(m_Omega_neutrinos, m_massless_neutrinos, int(value));
454  break;
455 
456  case (CosmologicalParameter::_neutrino_mass_):
457  set_OmegaNu(Omega_neutrinos(value), m_massless_neutrinos, m_massive_neutrinos);
458  break;
459 
460  case (CosmologicalParameter::_Omega_DE_):
461  set_OmegaDE(value);
462  break;
463 
464  case (CosmologicalParameter::_Omega_radiation_):
465  set_Omega_radiation(value);
466  break;
467 
468  case (CosmologicalParameter::_H0_):
469  set_H0(value);
470  break;
471 
472  case (CosmologicalParameter::_hh_):
473  set_hh(value, false);
474  break;
475 
476  case (CosmologicalParameter::_scalar_amp_):
477  set_scalar_amp(value);
478  break;
479 
480  case (CosmologicalParameter::_ln_scalar_amp_):
481  set_scalar_amp(exp(value)*1.e-10);
482  break;
483 
484  case (CosmologicalParameter::_scalar_pivot_):
485  set_scalar_pivot(value);
486  break;
487 
488  case (CosmologicalParameter::_n_spec_):
489  set_n_spec(value);
490  break;
491 
492  case (CosmologicalParameter::_w0_):
493  set_w0(value);
494  break;
495 
496  case (CosmologicalParameter::_wa_):
497  set_wa(value);
498  break;
499 
500  case (CosmologicalParameter::_fNL_):
501  set_fNL(value);
502  break;
503 
504  case (CosmologicalParameter::_sigma8_):
505  set_sigma8(value);
506  break;
507 
508  case (CosmologicalParameter::_tau_):
509  set_tau(value);
510  break;
511 
512  case (CosmologicalParameter::_rs_):
513  set_rs(value);
514  break;
515 
516  default:
517  ErrorCBL("no such a variable in the list!", "set_parameter", "Cosmology.cpp");
518  }
519 }
520 
521 
522 // =====================================================================================
523 
524 
525 void cbl::cosmology::Cosmology::set_parameters (const std::vector<CosmologicalParameter> parameter, const std::vector<double> value)
526 {
527  for (size_t i=0; i<parameter.size(); i++)
528  set_parameter(parameter[i], value[i]);
529 }
530 
531 
532 // =====================================================================================
533 
534 
535 double cbl::cosmology::Cosmology::w_CPL (const double redshift) const
536 {
537  return m_w0+m_wa*redshift/(1.+redshift);
538 }
539 
540 
541 // =====================================================================================
542 
543 
544 double cbl::cosmology::Cosmology::f_DE (const double redshift) const
545 {
546  // CPL parameterisation (see e.g. Bassett & Hlozek 2010)
547 
548  return pow((1.+redshift),3.*(1.+m_w0+m_wa))*exp(-3.*m_wa*redshift/(1.+redshift));
549 
550  /*
551  // direct calculation, useful for future generalizations:
552  func_fDE Func(m_w0,m_wa);
553  return exp(3.*qromb(Func,0.,redshift));
554  */
555 }
556 
557 
558 // =====================================================================================
559 
560 
561 double cbl::cosmology::Cosmology::EE (const double redshift) const
562 {
563  return sqrt(m_Omega_matter*pow((1.+redshift), 3)+m_Omega_DE*f_DE(redshift)+m_Omega_k*pow((1.+redshift), 2)+m_Omega_radiation*pow(1.+redshift, 4));
564 }
565 
566 
567 // =====================================================================================
568 
569 
570 double cbl::cosmology::Cosmology::HH (const double redshift) const
571 {
572  return m_H0*EE(redshift);
573 }
574 
575 
576 // =====================================================================================
577 
578 
579 double cbl::cosmology::Cosmology::OmegaM (const double redshift) const
580 {
581  return m_Omega_matter/EE2(redshift)*(1.+redshift)*(1.+redshift)*(1.+redshift);
582 }
583 
584 
585 // =====================================================================================
586 
587 
588 double cbl::cosmology::Cosmology::OmegaDE (const double redshift) const
589 {
590  return m_Omega_DE/EE2(redshift)*f_DE(redshift);
591 }
592 
593 
594 // =====================================================================================
595 
596 
597 double cbl::cosmology::Cosmology::OmegaR (const double redshift) const
598 {
599  return m_Omega_radiation/EE2(redshift)*pow(1.+redshift,4.);
600 }
601 
602 
603 // =====================================================================================
604 
605 
606 double cbl::cosmology::Cosmology::OmegaR_zeq (const double z_eq) const
607 {
608  return m_Omega_matter/(z_eq+1.);
609 }
610 
611 
612 // =====================================================================================
613 
614 
615 double cbl::cosmology::Cosmology::OmegaK (const double redshift) const
616 {
617  return m_Omega_k/EE2(redshift)*(1.+redshift)*(1.+redshift);
618 }
619 
620 
621 // =====================================================================================
622 
623 
624 double cbl::cosmology::Cosmology::OmegaNu (const double redshift) const
625 {
626  return m_Omega_neutrinos/EE2(redshift)*(1.+redshift)*(1.+redshift)*(1.+redshift);
627 }
628 
629 
630 // =====================================================================================
631 
632 
633 double cbl::cosmology::Cosmology::Omega (const double redshift) const
634 {
635  return (m_Omega_radiation*pow(1.+redshift,4.)+m_Omega_matter*pow(1.+redshift,3.)+m_Omega_DE*f_DE(redshift)*pow(1.+redshift,4.))/EE2(redshift);
636 }
637 
638 
639 // =====================================================================================
640 
641 
642 double cbl::cosmology::Cosmology::Omega_neutrinos (const double Mnu) const
643 {
644  if (m_hh<1.e-33) ErrorCBL("m_hh should be >0", "Omega_neutrinos", "Cosmology.cpp");
645  return Mnu/(93.14*pow(m_hh, 2));
646 }
647 
648 
649 // =====================================================================================
650 
651 
652 double cbl::cosmology::Cosmology::neutrino_mass () const
653 {
654  if (m_hh<1.e-33) ErrorCBL("m_hh should be >0", "neutrino_mass", "Cosmology.cpp");
655  return m_Omega_neutrinos*93.14*pow(m_hh, 2);
656 }
657 
658 
659 // =====================================================================================
660 
661 
662 double cbl::cosmology::Cosmology::linear_growth_rate (const double redshift, const double prec) const
663 {
664  if (m_w0==-1. and m_wa==0. and m_Omega_neutrinos==0.)
665  return pow(OmegaM(redshift), 0.545);
666 
667  else {
668  // Dodelson eq. (7.59)
669  const double a_in = 1.e-8;
670  double ff = (m_Omega_radiation==0.) ? 1. : (a_in/m_Omega_radiation*m_Omega_matter)/(a_in/m_Omega_radiation*m_Omega_matter+2./3.);
671 
672  auto func = [&] (const double y, double &dyda, const double ln_aa)
673  {
674  const double zz = 1./exp(ln_aa)-1.;
675  dyda = -y*y-y*(1.-0.5*(OmegaM(zz)+OmegaR(zz)+(1.+3.*m_w0+3.*m_wa*zz/(1.+zz))*OmegaDE(zz)))+1.5*OmegaM(zz);
676  };
677 
678 
679  typedef boost::numeric::odeint::runge_kutta_dopri5<double> stepper_type;
680 
681  boost::numeric::odeint::integrate_adaptive(make_controlled(1e-8, 1e-8, stepper_type() ),
682  func, ff, log(a_in), log(1./(1.+redshift)), prec);
683  return ff;
684  }
685 }
686 
687 
688 // =====================================================================================
689 
690 
691 double cbl::cosmology::Cosmology::DN (const double redshift, const double redshift_norm, const double prec) const
692 {
693  auto func = [prec, this] (const double aa) { return linear_growth_rate(1./aa-1., prec)/aa; };
694 
695  return exp(wrapper::gsl::GSL_integrate_qag(func, 1./(1.+redshift_norm), 1./(1.+redshift)));
696 }
697 
698 
699 // =====================================================================================
700 
701 
702 double cbl::cosmology::Cosmology::DD (const double redshift) const
703 {
704  if (m_w0==-1. and m_wa==0. and m_Omega_neutrinos==0.) {
705  const double aa = 1./(1.+redshift);
706 
707  auto func = [&] (const double aa) { return pow(aa*EE(1./aa-1.),-3.); };
708 
709  return cbl::wrapper::gsl::GSL_integrate_qag(func, 0., aa)*2.5*Omega_matter()*EE(redshift);
710  }
711 
712  else
713  return ErrorCBL("the current implementation is valid only for LCDM cosmologies", "DD", "Cosmology.cpp", ExitCode::_workInProgress_);
714 }
715 
716 
717 // =====================================================================================
718 
719 
720 double cbl::cosmology::Cosmology::gg (const double redshift) const
721 {
722  return DD(redshift)*(1.+redshift);
723 }
724 
725 
726 // =====================================================================================
727 
728 
729 double cbl::cosmology::Cosmology::sigma8 (const double redshift) const
730 {
731  if (m_sigma8<0)
732  ErrorCBL("sigma8 at z=0 is not set!", "sigma8", "Cosmology.cpp");
733 
734  return m_sigma8*DN(redshift);
735 }
736 
737 
738 // =====================================================================================
739 
740 
741 double cbl::cosmology::Cosmology::D_C (const double redshift) const
742 {
743  if (redshift<0) ErrorCBL("redshift have to be >=0!", "D_C", "Cosmology.cpp");
744 
745  double Dc;
746 
747  if (m_model=="LCDM") {
748  function<double(double)> integrand = bind(&Cosmology::EE_inv, this, std::placeholders::_1);
749  Dc = wrapper::gsl::GSL_integrate_qag(integrand, 0, redshift);
750  }
751 
752  else {
753  cbl::Path path;
754  string dir = path.DirCosmo()+"Cosmology/Tables/dc_cDE/";
755  string file_in;
756  if (m_model=="LCDM_Baldi_wmap7") file_in = dir+"LCDM-wmap7-comovingdist.dat";
757  else if (m_model=="EXP005_Baldi_wmap7") file_in = dir+"EXP005-wmap7-comovingdist.dat";
758  else if (m_model=="EXP010e2_Baldi_wmap7") file_in = dir+"EXP010e2-wmap7-comovingdist.dat";
759  else if (m_model=="LCDM_Baldi_CoDECS") file_in = dir+"LCDM_CoDECS-comovingdist.dat";
760  else if (m_model=="EXP001_Baldi_CoDECS") file_in = dir+"EXP001_CoDECS-comovingdist.dat";
761  else if (m_model=="EXP002_Baldi_CoDECS") file_in = dir+"EXP002_CoDECS-comovingdist.dat";
762  else if (m_model=="EXP003_Baldi_CoDECS") file_in = dir+"EXP003_CoDECS-comovingdist.dat";
763  else if (m_model=="EXP008e3_Baldi_CoDECS") file_in = dir+"EXP008e3_CoDECS-comovingdist.dat";
764  else if (m_model=="EXP010e2_Baldi_CoDECS") file_in = dir+"EXP010e2_CoDECS-comovingdist.dat";
765  else if (m_model=="SUGRA003_Baldi_CoDECS") file_in = dir+"SUGRA003_CoDECS-comovingdist.dat";
766  else { ErrorCBL("model = " + m_model + "!", "D_C", "Cosmology.cpp"); }
767 
768  ifstream fin(file_in.c_str()); checkIO(fin, file_in);
769 
770  double Red, DC;
771  vector<double> Redshift, dc;
772  while (fin >>Red>>DC) {
773  Redshift.push_back(Red);
774  dc.push_back(DC);
775  }
776  fin.clear(); fin.close();
777 
778  double err = -1.;
779  Dc = interpolated(redshift, Redshift, dc, "Rat");
780 
781  if (err/Dc>0.1) {
782  string Err = "Error in cbl::cosmology::Cosmology::D_C of Cosmology.cpp: " + conv(redshift, par::fDP3) + " " + conv(Redshift.size(), par::fINT) + " " + conv(dc.size(), par::fINT);
783  ErrorCBL(conv(redshift, par::fDP3) + " " + conv(Redshift.size(), par::fINT) + " " + conv(dc.size(), par::fINT), "D_C", "Cosmology.cpp");
784  }
785  }
786 
787  return m_D_H*Dc;
788 }
789 
790 
791 // =====================================================================================
792 
793 
794 void cbl::cosmology::Cosmology::D_C_table (const std::string file_table, const double z_min, const double z_max, const int step, std::vector<double> &Redshift, std::vector<double> &dc) const
795 {
796  cbl::Path path;
797  string File_table = path.DirCosmo()+"Cosmology/Tables/dc/"+file_table;
798 
799  ifstream fin;
800  fin.open (File_table.c_str());
801  if (!fin) {
802 
803  ofstream fout(File_table.c_str()); checkIO(fout, File_table);
804 
805  double delta_z = (z_max-z_min)/step;
806  double z1 = z_min;
807  double z2 = z_min+delta_z;
808 
809  for (int i=0; i<step; i++) {
810  double zmean = (z1+z2)*0.5;
811  fout <<zmean<<" "<<D_C(zmean)<<endl;
812  z1 = z2; z2 += delta_z;
813  }
814 
815  fout.clear(); fout.close(); coutCBL <<"I wrote the file: "<<File_table<<endl;
816  }
817  fin.clear(); fin.close();
818 
819  fin.open(File_table.c_str());
820  double Red, DC;
821  while (fin >>Red>>DC) {
822  Redshift.push_back(Red);
823  dc.push_back(DC);
824  }
825  fin.clear(); fin.close();
826 }
827 
828 
829 // =====================================================================================
830 
831 
832 double cbl::cosmology::Cosmology::D_M (const double redshift) const
833 {
834  if (m_Omega_k>1.e-10)
835  return m_D_H/sqrt(m_Omega_k)*sinh(sqrt(m_Omega_k)*D_C(redshift)/m_D_H);
836 
837  else if (fabs(m_Omega_k)<1.e-10)
838  return D_C(redshift);
839 
840  else
841  return m_D_H/sqrt(-m_Omega_k)*sin(sqrt(-m_Omega_k)*D_C(redshift)/m_D_H);
842 }
843 
844 
845 // =====================================================================================
846 
847 
848 double cbl::cosmology::Cosmology::D_A (const double redshift) const
849 {
850  return D_M(redshift)/(1.+redshift);
851 }
852 
853 
854 // =====================================================================================
855 
856 
857 double cbl::cosmology::Cosmology::D_A (const double z1, const double z2) const
858 {
859  if (m_Omega_k<0)
860  ErrorCBL("the implemented formula is not correct for Omega_k<0!", "D_A", "Cosmology.cpp");
861 
862  const double _z1 = (z1<z2) ? z1 : z2;
863  const double _z2 = (z1<z2) ? z2 : z1;
864 
865  const double D1 = D_M(_z1);
866  const double D2 = D_M(_z2);
867 
868  return 1./(1.+_z2)*(D2*sqrt(1.+m_Omega_k*pow(D1/m_D_H, 2))-D1*sqrt(1.+m_Omega_k*pow(D2/m_D_H, 2)));
869 }
870 
871 
872 // =====================================================================================
873 
874 
875 double cbl::cosmology::Cosmology::D_L (const double redshift) const
876 {
877  return (1.+redshift)*D_M(redshift);
878 }
879 
880 
881 // =====================================================================================
882 
883 
884 double cbl::cosmology::Cosmology::D_V (const double redshift) const
885 {
886  return pow(pow(D_M(redshift),2)*par::cc*redshift/HH(redshift),1./3.);
887 }
888 
889 
890 // =====================================================================================
891 
892 
893 double cbl::cosmology::Cosmology::F_AP (const double redshift) const
894 {
895  return D_M(redshift)*HH(redshift)/par::cc;
896 }
897 
898 
899 // =====================================================================================
900 
901 
902 double cbl::cosmology::Cosmology::Distance (const double redshift, const std::string distance_type) const
903 {
904  if (distance_type=="DC")
905  return D_C(redshift);
906 
907  else if (distance_type=="DL")
908  return D_L(redshift);
909 
910  else if (distance_type=="DA")
911  return D_A(redshift);
912 
913  else if (distance_type=="Dv")
914  return D_V(redshift);
915 
916  else if (distance_type=="Dvrs")
917  return D_V(redshift)/rs_CAMB();
918 
919  else if (distance_type=="rsDv")
920  return rs_CAMB()/D_V(redshift);
921 
922  else {
923  ErrorCBL("no such a distance type!", "Distance", "Cosmology.cpp");
924  return -1;
925  }
926 }
927 
928 
929 // =====================================================================================
930 
931 
932 double cbl::cosmology::Cosmology::lookback_time (const double redshift) const
933 {
934  function<double(double)> integrand = bind(&Cosmology::EE_inv2, this, std::placeholders::_1);
935  double tt = wrapper::gsl::GSL_integrate_qag(integrand,0, redshift);
936 
937  double Mpc = par::mega*par::pc*1.e-3; // in Km;
938  double Gyr = par::giga*par::yr; // in sec
939 
940  return 1./(m_hh*100.)*tt*Mpc/Gyr;
941 }
942 
943 
944 // =====================================================================================
945 
946 
947 double cbl::cosmology::Cosmology::cosmic_time (const double redshift) const
948 {
949  function<double(double)> integrand = bind(&Cosmology::EE_inv3, this, std::placeholders::_1);
950 
951  double aa = 1./(1.+redshift);
952  double tt = wrapper::gsl::GSL_integrate_qag(integrand,0, aa);
953 
954  double Mpc = par::mega*par::pc*1.e-3; // in Km;
955  double Gyr = par::giga*par::yr; // in sec
956 
957  return 1./(m_hh*100.)*tt*Mpc/Gyr;
958 }
959 
960 
961 // =====================================================================================
962 
963 
964 double cbl::cosmology::Cosmology::EE2 (const double redshift) const
965 {
966  return m_Omega_matter*pow((1.+redshift),3)+m_Omega_DE*f_DE(redshift)+m_Omega_k*pow((1.+redshift),2)+m_Omega_radiation*pow(1.+redshift,4);
967 }
968 
969 
970 // =====================================================================================
971 
972 
973 double cbl::cosmology::Cosmology::qq (const double redshift) const
974 {
975  if (m_wa!=0) ErrorCBL("w_a!=0", "qq", "Cosmology.cpp", ExitCode::_workInProgress_);
976 
977  double aa = 1./(1.+redshift);
978  return (m_Omega_matter*aa+2.*m_Omega_radiation+(1.+3.*m_w0)*m_Omega_DE*pow(aa,1.-3.*m_w0))/(2.*EE2(redshift));
979 }
980 
981 
982 // =====================================================================================
983 
984 
985 double cbl::cosmology::Cosmology::Hdot (const double redshift) const
986 {
987  return -pow(HH(redshift),2)*(1.+qq(redshift));
988 }
989 
990 
991 // =====================================================================================
992 
993 
994 double cbl::cosmology::Cosmology::z_acc () const
995 {
996  if (m_wa!=0) ErrorCBL("w_a!=0", "z_acc", "Cosmology.cpp", ExitCode::_workInProgress_);
997 
998  double z_acc = pow(-(1.+3.*m_w0)*m_Omega_DE/m_Omega_matter,-1./(3.*m_w0))-1.;
999  if (std::isnan(z_acc)) ErrorCBL("", "z_acc", "Cosmology.cpp");
1000 
1001  return z_acc;
1002 }
1003 
1004 
1005 // =====================================================================================
1006 
1007 
1008 double cbl::cosmology::Cosmology::z_eq () const
1009 {
1010  if (m_wa!=0) ErrorCBL("w_a!=0", "z_eq", "Cosmology.cpp", ExitCode::_workInProgress_);
1011  return pow(m_Omega_DE/m_Omega_matter,-1./(3.*m_w0))-1.;
1012 }
1013 
1014 
1015 // =====================================================================================
1016 
1017 
1018 double cbl::cosmology::Cosmology::z_eq_rad (const double T_CMB) const
1019 {
1020  return 2.5e+4*m_Omega_matter*m_hh*m_hh*pow((T_CMB/2.7),-4.);
1021 }
1022 
1023 
1024 // =====================================================================================
1025 
1026 
1027 double cbl::cosmology::Cosmology::Mag_Volume_limited (const double z_max, const double mag_lim) const
1028 {
1029  return (mag_lim-5.*log10(D_L(z_max)))-25.;
1030 }
1031 
1032 
1033 // =====================================================================================
1034 
1035 
1036 double cbl::cosmology::Cosmology::Lum_bol (const double redshift, const double flux) const
1037 {
1038  return 4.*par::pi*flux*pow(D_L(redshift),2);
1039 }
1040 
1041 
1042 // =====================================================================================
1043 
1044 
1045 double cbl::cosmology::Cosmology::Redshift (const double d_c, const double z1_guess, const double z2_guess, const double prec) const
1046 {
1047  double Prec = prec;
1048  if (Prec<1.e-5) {
1049  WarningMsgCBL("prec has been set to 1.e-5", "Redshift", "Cosmology.h");
1050  Prec = 1.e-5;
1051  }
1052 
1053  double redshift = -1.;
1054 
1055  if (m_model=="LCDM") {
1056  function<double(double)> func = bind(&Cosmology::D_C, this, std::placeholders::_1);
1057  redshift = wrapper::gsl::GSL_root_brent(func, d_c, z1_guess, z2_guess, prec);
1058  }
1059 
1060  else {
1061  WarningMsgCBL("the quantity prec is not used", "Redshift", "Cosmology.h");
1062 
1063  cbl::Path path;
1064  string dir = path.DirCosmo()+"Cosmology/Tables/dc_cDE/";
1065  string file_in;
1066 
1067  if (m_model=="LCDM_Baldi_wmap7") file_in = dir+"LCDM-wmap7-comovingdist.dat";
1068  else if (m_model=="EXP005_Baldi_wmap7") file_in = dir+"EXP005-wmap7-comovingdist.dat";
1069  else if (m_model=="EXP010e2_Baldi_wmap7") file_in = dir+"EXP010e2-wmap7-comovingdist.dat";
1070  else if (m_model=="LCDM_Baldi_CoDECS") file_in = dir+"LCDM_CoDECS-comovingdist.dat";
1071  else if (m_model=="EXP001_Baldi_CoDECS") file_in = dir+"EXP001_CoDECS-comovingdist.dat";
1072  else if (m_model=="EXP002_Baldi_CoDECS") file_in = dir+"EXP002_CoDECS-comovingdist.dat";
1073  else if (m_model=="EXP003_Baldi_CoDECS") file_in = dir+"EXP003_CoDECS-comovingdist.dat";
1074  else if (m_model=="EXP008e3_Baldi_CoDECS") file_in = dir+"EXP003_CoDECS-comovingdist.dat";
1075  else if (m_model=="EXP010e2_Baldi_CoDECS") file_in = dir+"EXP010e2_CoDECS-comovingdist.dat";
1076  else if (m_model=="SUGRA003_Baldi_CoDECS") file_in = dir+"SUGRA003_CoDECS-comovingdist.dat";
1077  else { ErrorCBL("model = " + m_model + "!", "Redshift", "Cosmology.cpp"); }
1078 
1079  ifstream fin(file_in.c_str()); checkIO(fin, file_in);
1080 
1081  double Red, DC;
1082  vector<double> Redshift, dc;
1083  while (fin >> Red >> DC) {
1084  Redshift.push_back(Red);
1085  dc.push_back(m_D_H*DC);
1086  }
1087  fin.clear(); fin.close();
1088 
1089  double err = -1;
1090  redshift = interpolated(d_c, dc, Redshift, "Rat");
1091  if (err/redshift>0.1) ErrorCBL("", "Redshift", "Cosmology.cpp");
1092  }
1093 
1094  return redshift;
1095 }
1096 
1097 
1098 // =====================================================================================
1099 
1100 
1101 double cbl::cosmology::Cosmology::Redshift_LCDM (const double d_c, const double z1_guess, const double z2_guess, const bool go_fast, const double prec) const
1102 {
1103  if (m_model!="LCDM")
1104  ErrorCBL("this method works only for a LambdaCDM universe", "Redshift_LCDM", "Cosmology.cpp");
1105 
1106  double Prec = prec;
1107  if (Prec<1.e-5) {
1108  WarningMsgCBL("prec has been set to 1.e-5", "Redshift_LCDM", "Cosmology.h");
1109  Prec = 1.e-5;
1110  }
1111 
1112  double fact = 1./(par::cc/100.); // factor to convert [Mpc/h] into [c/H0]
1113 
1114  double D_c = d_c * fact; // [Mpc/h] -> [c/H0]
1115 
1116 
1117  // first interval definition: based on user definition and cosmology
1118 
1119  double z0 = max(max(z1_guess, D_c), 4./pow(sqrt(m_Omega_matter)*D_c-2.,2)-1.); // z is always > than these two numbers
1120 
1121  double z1 = (D_c<2.) ? min(z2_guess, 4./pow(D_c-2.,2)-1.) : z2_guess; // z is always < than this number, but valid only when d_c<2
1122 
1123  z1 = max(z1, z0);
1124 
1125 
1126  // choice of the method
1127 
1128  function<double(double)> DC = (go_fast) ? bind(&Cosmology::D_C_LCDM, this, std::placeholders::_1) : bind(&Cosmology::D_C, this, std::placeholders::_1);
1129 
1130 
1131  // interval check: needs z0<=z and z1>=z
1132 
1133  double d0 = DC(z0)*fact;
1134 
1135  if (d0>D_c) {
1136  z0 = D_c/d0*z0/EE(z0);
1137  d0 = DC(z0)*fact;
1138  }
1139 
1140  double d1 = DC(z1)*fact;
1141  if (d1<D_c) {
1142  z1 = EE(z1)*D_c/d1*z1;
1143  d1 = DC(z1)*fact;
1144  }
1145 
1146 
1147  // calculation
1148 
1149  double diffz = z1-z0;
1150  double diffd = d1-d0;
1151  double prec2 = 2.*Prec; // prec is the tolerance on the result, prec2 on the interval length
1152 
1153  while (diffz>prec2 && diffd>0.) {
1154  z1 = z0+(D_c-d0)*(diffz/diffd); // given the shape this estimate is always > z
1155  d1 = DC(z1)*fact;
1156  z0 = z0+EE(z0)*(D_c-d0); // given the shape this estimate is always < z
1157  d0 = DC(z0)*fact;
1158  diffz = z1-z0;
1159  diffd = d1-d0;
1160  }
1161 
1162  return 0.5*(z0+z1);
1163 }
1164 
1165 
1166 // =====================================================================================
1167 
1168 
1169 double cbl::cosmology::Cosmology::Redshift_time (const double time, const double z1_guess, const double z2_guess) const
1170 {
1171  if (m_model!="LCDM") ErrorCBL("model!=LCDM", "Redshift_time", "Cosmology.cpp", ExitCode::_workInProgress_);
1172 
1173  double prec = 0.0001;
1174 
1175  function<double(double)> func = bind(&Cosmology::cosmic_time, this, std::placeholders::_1);
1176  return wrapper::gsl::GSL_root_brent(func, time, z1_guess, z2_guess, prec);
1177 }
1178 
1179 
1180 // =====================================================================================
1181 
1182 
1183 double cbl::cosmology::Cosmology::Volume (const double z1, const double z2, const double Area) const
1184 {
1185  const double Area_steradians = Area*pow(par::pi/180., 2);
1186 
1187  return 4./3.*par::pi*fabs(pow(D_C(z1), 3)-pow(D_C(z2), 3))*Area_steradians/(4.*par::pi);
1188 }
1189 
1190 
1191 // =====================================================================================
1192 
1193 
1194 double cbl::cosmology::Cosmology::Volume (const double z1, const double z2, const double RA_min, const double RA_max, const double Dec_min, const double Dec_max) const
1195 {
1196  const double Area = degrees(haversine_distance(RA_min, RA_max, Dec_min, Dec_min))*degrees(haversine_distance(RA_min, RA_min, Dec_min, Dec_max));
1197 
1198  return Volume(z1, z2, Area);
1199 }
1200 
1201 
1202 // =====================================================================================
1203 
1204 
1205 double cbl::cosmology::Cosmology::Volume (const double zz) const
1206 {
1207  const double DDMM = D_M(zz);
1208 
1209  if (m_Omega_k>1.e-10)
1210  return (4*par::pi*pow(m_D_H, 3)/(2*m_Omega_k))*(DDMM/m_D_H*sqrt(1+m_Omega_k*pow(DDMM/m_D_H, 2))-pow(fabs(m_Omega_k), -0.5)*asinh(sqrt(fabs(m_Omega_k))*DDMM/m_D_H));
1211 
1212  else if (fabs(m_Omega_k)<1.e-10)
1213  return 4*par::pi*pow(DDMM, 3)/3;
1214 
1215  else
1216  return (4*par::pi*pow(m_D_H, 3)/(2*m_Omega_k))*(DDMM/m_D_H*sqrt(1+m_Omega_k*pow(DDMM/m_D_H, 2))-pow(fabs(m_Omega_k), -0.5)*asin(sqrt(fabs(m_Omega_k))*DDMM/m_D_H));
1217 }
1218 
1219 
1220 // =====================================================================================
1221 
1222 
1223 double cbl::cosmology::Cosmology::max_redshift (const double Volume, const double Area, const double z_min) const
1224 {
1225  double Area_steradians = Area*pow(par::pi/180.,2);
1226  double dcz2 = pow(3*Volume/Area_steradians+pow(D_C(z_min),3),1./3);
1227 
1228  function<double(double)> func = bind(&Cosmology::D_C, this, std::placeholders::_1);
1229  return wrapper::gsl::GSL_root_brent(func, dcz2, z_min, 10, 1.e-9);
1230 }
1231 
1232 
1233 // =====================================================================================
1234 
1235 
1236 double cbl::cosmology::Cosmology::dV_dZdOmega (const double redshift, const bool angle_rad) const
1237 {
1238  // angle_rad: true -> Omega in steradians; false -> Omega in square degrees
1239  double conv = (angle_rad) ? 1. : 3282.80635;
1240 
1241  return m_D_H*pow((1.+redshift)*D_A(redshift),2)/EE(redshift)/conv;
1242 }
1243 
1244 
1245 // =====================================================================================
1246 
1247 
1248 double cbl::cosmology::Cosmology::deltac (const double redshift) const
1249 {
1250  return 1.686*(1.+0.012299*log10(OmegaM(redshift)));
1251 }
1252 
1253 
1254 // =====================================================================================
1255 
1256 
1257 double cbl::cosmology::Cosmology::rho_crit (const double redshift, const bool unit1) const
1258 {
1259  // km sec^-1 Mpc^-1 -> sec^-1
1260  double HHc = HH(redshift)*pow(cbl::par::kilo*cbl::par::pc, -1);
1261 
1262  // force cosmological units
1263  if (!m_unit && unit1) HHc /= m_hh;
1264 
1265  // m^3 Kg^-1 sec^-2 -> Mpc^3 Msun^-1 sec^-2
1266  const double GNc = cbl::par::GN*(cbl::par::Msol*pow(cbl::par::mega*cbl::par::pc, -3));
1267 
1268  return 3.*pow(HHc, 2)/(8.*par::pi*GNc);
1269 }
1270 
1271 // =====================================================================================
1272 
1273 
1274 double cbl::cosmology::Cosmology::rho_m (const double redshift, const bool unit1, const bool nu) const
1275 {
1276  return rho_crit(redshift, unit1) * ((!nu) ? OmegaM(redshift) : OmegaM(redshift)-OmegaNu(redshift));
1277 }
1278 
1279 
1280 // =====================================================================================
1281 
1282 
1283 double cbl::cosmology::Cosmology::Delta_c (const double redshift, const std::string author) const
1284 {
1285  if (author=="BryanNorman") {
1286  const double xx = OmegaM(redshift)-1.;
1287 
1288  if (fabs(m_Omega_DE)<1.e-30)
1289  return 18.*par::pi*par::pi+60.*xx-32.*xx*xx;
1290 
1291  else if (fabs(m_Omega_k)<1.e-30)
1292  return 18.*par::pi*par::pi+82.*xx-39.*xx*xx;
1293 
1294  else return ErrorCBL("cosmological parameters not allowed for the current implementation", "Delta_c", "Cosmology.cpp");
1295  }
1296 
1297  else if (author=="Eke") {
1298  if (fabs(m_Omega_DE)<1.e-30)
1299  return 178.*pow(OmegaM(redshift), 0.30);
1300 
1301  else if (fabs(m_Omega_k)<1.e-30)
1302  return 178.*pow(OmegaM(redshift), 0.45);
1303 
1304  else return ErrorCBL("cosmological parameters not allowed for the current implementation", "Delta_c", "Cosmology.cpp");
1305  }
1306 
1307  else if (author=="NakamuraSuto") {
1308  if (fabs(m_Omega_DE)<1.e-30) {
1309  const double x = pow(1./m_Omega_matter-1., 1./3.)/(1.+redshift);
1310  return 18.*pow(par::pi, 2)*(1+0.4093*pow(x, 2.7152))*OmegaM(redshift);
1311  }
1312  else return ErrorCBL("cosmological parameters not allowed for the current implementation", "Delta_c", "Cosmology.cpp");
1313  }
1314 
1315  else return ErrorCBL("author not allowed", "Delta_c", "Cosmology.cpp");
1316 }
1317 
1318 
1319 // =====================================================================================
1320 
1321 
1322 double cbl::cosmology::Cosmology::M_vir (const double r_vir, const double redshift, const std::string author, const bool unit1) const
1323 {
1324  return 4./3.*par::pi*pow(r_vir, 3)*Delta_c(redshift, author)*rho_crit(redshift, unit1);
1325 }
1326 
1327 
1328 // =====================================================================================
1329 
1330 
1331 double cbl::cosmology::Cosmology::r_vir (const double M_vir, const double redshift, const std::string author, const bool unit1) const
1332 {
1333  return pow(3*M_vir/(4.*par::pi*Delta_c(redshift, author)*rho_crit(redshift, unit1)), 1./3.);
1334 }
1335 
1336 
1337 // =====================================================================================
1338 
1339 
1340 double cbl::cosmology::Cosmology::c_vir (const double c_200, const double redshift, const std::string author) const
1341 {
1342  const double a = -1.119*log10(Delta_c(redshift, author))+3.537;
1343  const double b = -0.967*log10(Delta_c(redshift, author))+2.181;
1344  return a*c_200+b;
1345 }
1346 
1347 
1348 // =====================================================================================
1349 
1350 
1351 double cbl::cosmology::Cosmology::D_C_LCDM (const double redshift) const
1352 {
1353  if (m_model!="LCDM" || fabs(m_Omega_k)>1.e-10)
1354  ErrorCBL("this method works only for a flat LambdaCDM universe; it does not work for non-standard dark energy or non-flat models", "D_C_LCDM", "Cosmology.h");
1355 
1356  if (redshift>10)
1357  WarningMsgCBL("this function works only at low enough redshifts where &Omega<SUB>r</SUB> is negligible", "D_C_LCDM","Cosmology.h");
1358 
1359  double CC = 1./pow(3.,0.25)/(pow(m_Omega_matter,1./3.)*pow(1.-m_Omega_matter,1./6.));
1360  double f_m = (1.-sqrt(3.))*pow((1.-m_Omega_matter)/m_Omega_matter,1./3.);
1361  double f_p = (1.+sqrt(3.))*pow((1.-m_Omega_matter)/m_Omega_matter,1./3.);
1362  double phi0 = acos((1.+f_m)/(1.+f_p));
1363  double F_phi0 = m_elf_dz(phi0);
1364 
1365  double aa = 1./(1.+redshift);
1366  double phi1 = acos((1.+f_m*aa)/(1.+f_p*aa));
1367 
1368  return CC*(F_phi0-m_elf_dz(phi1))*2997.9199*m_t_H*100.;
1369 }
1370 
1371 
1372 // =====================================================================================
1373 
1374 
1375 double cbl::cosmology::Cosmology::m_elf_dz (const double phi) const
1376 {
1377  int jj = round(phi/M_PI);
1378  double phi0 = phi-jj*M_PI;
1379 
1380  // then, it has to be reduced again to [0,pi/2], taking the sign into account
1381  double ss = phi0/abs(phi0);
1382  phi0 = ss*phi0;
1383 
1384  // optimisation parameters
1385  double phiS = 1.249;
1386  double yS = 0.9;
1387  double phic = 1.5707963-phi0; // pi/2 - phi0
1388 
1389  if (phi0 < phiS)
1390  return ss*m_asn_dz(sin(phi0))+jj*5.5361264;
1391  else {
1392  double cc = sin(phic);
1393  double xx = cc*cc;
1394  double d2 = 0.066987298 + 0.93301270 * xx;
1395  if (xx < yS*d2)
1396  return ss*(2.7680632-m_asn_dz(cc/sqrt(d2)))+jj*5.5361264;
1397  else {
1398  double vv = 0.066987298 * (1.-xx);
1399  if (vv < xx*d2) return ss*m_acn_dz(cc);
1400  else return ss*m_acn_dz(sqrt(vv/d2)) + jj*5.5361264;
1401  }
1402  }
1403 }
1404 
1405 
1406 // =====================================================================================
1407 
1408 
1409 double cbl::cosmology::Cosmology::m_acn_dz (const double cc) const
1410 {
1411  double pp = 1.;
1412  double xx = cc*cc;
1413  for (int j=1; j<10; j++) {
1414  if (xx > 0.5) return pp*m_asn_dz(sqrt(1.-xx));
1415  double dd = sqrt(0.066987298+0.93301270*xx);
1416  xx = (sqrt(xx)+dd)/(1.+dd);
1417  pp *= 2.;
1418  }
1419 
1420  ErrorCBL("too many half argument transformations of cn", "m_acn_dz", "Cosmology.cpp");
1421  return 0;
1422 }
1423 
1424 
1425 // =====================================================================================
1426 
1427 
1428 double cbl::cosmology::Cosmology::m_asn_dz (const double ss) const
1429 {
1430  //double yA = 0.034856757;
1431  double yA = 0.153532;
1432 
1433  double yy = ss*ss;
1434  if (yy < yA)
1435  return ss*m_serf_dz(yy);
1436 
1437  double pp = 1.;
1438  for (int j=1; j<10; j++) {
1439  yy /= ((1.+sqrt(1.-yy))*(1.+sqrt(1.-0.93301270*yy)));
1440  pp *= 2.;
1441  if (yy < yA)
1442  return pp*sqrt(yy)*m_serf_dz(yy);
1443  }
1444 
1445  return ErrorCBL("too many half argument transformations of sn", "m_asn_dz", "Cosmology.cpp");
1446 }
1447 
1448 
1449 // =====================================================================================
1450 
1451 
1452 double cbl::cosmology::Cosmology::m_serf_dz (const double yy) const
1453 {
1454  return (1.+yy*(0.32216878+yy*(0.18693909+yy*(0.12921048+yy*(0.097305732+yy*(0.077131543+yy*0.063267775+yy*(0.053185339+yy*(0.045545557+yy*0.039573617))))))));
1455 }
1456 
1457 
1458 // =====================================================================================
1459 
1460 
1461 double cbl::cosmology::Cosmology::concentration_NFW_Duffy (const double Mass, const double redshift, const string halo_def) const
1462 {
1463  double AA, BB, CC;
1464 
1465  if (redshift>2) ErrorCBL("The concentration-mass relation by Duffy et al. has been tested only at z<2", "concentration_NFW_Duffy", "Cosmology.cpp");
1466 
1467  if (halo_def=="200") {
1468  AA = 5.71;
1469  BB = -0.084;
1470  CC = -0.47;
1471  }
1472 
1473  else if (halo_def=="vir") {
1474  AA = 7.85;
1475  BB = -0.081;
1476  CC = -0.71;
1477  }
1478 
1479  else if (halo_def=="mean") {
1480  AA = 10.14;
1481  BB = -0.081;
1482  CC = -1.01;
1483  }
1484  else return ErrorCBL("halo_def not allowed!", "concentration_NFW_Duffy", "Cosmology.cpp");
1485 
1486  const double Mpivot = 2.e12; // in [Msun/h]
1487 
1488  return AA*pow(Mass/Mpivot, BB)*pow(1.+redshift, CC);
1489 }
1490 
1491 
1492 // ============================================================================================
Cosmology.h
The class Cosmology.
coutCBL
#define coutCBL
CBL print message.
Definition: Kernel.h:734
cbl::Path
The class Path.
Definition: Path.h:145
cbl::Path::DirCosmo
std::string DirCosmo()
get the directory where the CosmoBolognaLbi are stored
Definition: Path.h:178
cbl::cosmology::Cosmology::rho_crit
double rho_crit(const double redshift, const bool unit1=false) const
the critical cosmic density
Definition: Cosmology.cpp:1257
cbl::cosmology::Cosmology::EE_inv3
double EE_inv3(const double aa) const
inverse of the auxiliary function used to compute the Hubble function integrand of the cosmic time
Definition: Cosmology.h:1928
cbl::cosmology::Cosmology::DN
double DN(const double redshift, const double redshift_norm=0., const double prec=1.e-4) const
the normalised amplitude of the growing mode at a given redshift,
Definition: Cosmology.cpp:691
cbl::cosmology::Cosmology::m_hh
double m_hh
: the Hubble parameter,
Definition: Cosmology.h:312
cbl::cosmology::Cosmology::set_parameters
void set_parameters(const std::vector< CosmologicalParameter > parameter, const std::vector< double > value)
set the value of some cosmological paramters
Definition: Cosmology.cpp:525
cbl::cosmology::Cosmology::c_vir
double c_vir(const double c200, const double redshift, const std::string author="BryanNorman") const
virial halo concentration given
Definition: Cosmology.cpp:1340
cbl::cosmology::Cosmology::neutrino_mass
double neutrino_mass() const
the total neutrino mass
Definition: Cosmology.cpp:652
cbl::cosmology::Cosmology::m_scalar_pivot
double m_scalar_pivot
the scalar pivot k in
Definition: Cosmology.h:327
cbl::cosmology::Cosmology::rho_m
double rho_m(const double redshift=0., const bool unit1=false, const bool nu=false) const
the mean cosmic background density
Definition: Cosmology.cpp:1274
cbl::cosmology::Cosmology::Mag_Volume_limited
double Mag_Volume_limited(const double z_max=1., const double mag_lim=-20.) const
maximum absolute magnitude to have a volume-limited catalogue
Definition: Cosmology.cpp:1027
cbl::cosmology::Cosmology::Omega
double Omega(const double redshift=0.) const
the cosmic density at a given redshift
Definition: Cosmology.cpp:633
cbl::cosmology::Cosmology::m_Omega_radiation
double m_Omega_radiation
: the radiation density at z=0
Definition: Cosmology.h:300
cbl::cosmology::Cosmology::Delta_c
double Delta_c(const double redshift, const std::string author="BryanNorman") const
the critical overdensity
Definition: Cosmology.cpp:1283
cbl::cosmology::Cosmology::m_n_spec
double m_n_spec
: the primordial spectral index
Definition: Cosmology.h:330
cbl::cosmology::Cosmology::z_acc
double z_acc() const
redshift at which the Universe begins to accelerate
Definition: Cosmology.cpp:994
cbl::cosmology::Cosmology::gg
double gg(const double redshift=0.) const
the linear growth factor at a given redshift,
Definition: Cosmology.cpp:720
cbl::cosmology::Cosmology::D_C_LCDM
double D_C_LCDM(const double redshift) const
the comoving line-of-sight distance at a given redshift
Definition: Cosmology.cpp:1351
cbl::cosmology::Cosmology::m_Omega_baryon
double m_Omega_baryon
: the baryon density at z=0
Definition: Cosmology.h:285
cbl::cosmology::Cosmology::set_parameter
void set_parameter(const CosmologicalParameter parameter, const double value)
set the value of one cosmological paramter
Definition: Cosmology.cpp:424
cbl::cosmology::Cosmology::w_CPL
double w_CPL(const double redshift=0.) const
the DE equation of state in the CPL parameterisation, as a function of redshift
Definition: Cosmology.cpp:535
cbl::cosmology::Cosmology::m_elf_dz
double m_elf_dz(const double phi) const
the incomplete elliptic integral
Definition: Cosmology.cpp:1375
cbl::cosmology::Cosmology::sigma8
double sigma8() const
get the private member Cosmology::m_sigma8
Definition: Cosmology.h:1212
cbl::cosmology::Cosmology::max_redshift
double max_redshift(const double Volume, const double Area, const double z_min) const
maximum redshift for a given volume, sky area and minimum redshift
Definition: Cosmology.cpp:1223
cbl::cosmology::Cosmology::Volume
double Volume(const double z1, const double z2, const double Area) const
comoving volume for a given redshift range and sky area
Definition: Cosmology.cpp:1183
cbl::cosmology::Cosmology::EE
double EE(const double redshift=0.) const
auxiliary function used to compute the Hubble function
Definition: Cosmology.cpp:561
cbl::cosmology::Cosmology::F_AP
double F_AP(const double redshift) const
F_AP, the ALCOCK-PACZYNSKI distortion parameter.
Definition: Cosmology.cpp:893
cbl::cosmology::Cosmology::D_V
double D_V(const double redshift) const
the average distance at a given redshift, used to rescale the correlation function
Definition: Cosmology.cpp:884
cbl::cosmology::Cosmology::z_eq_rad
double z_eq_rad(const double T_CMB=2.7255) const
redshift of matter-radiation equality
Definition: Cosmology.cpp:1018
cbl::cosmology::Cosmology::qq
double qq(const double redshift=0.) const
the deceleration parameter at a given redshift
Definition: Cosmology.cpp:973
cbl::cosmology::Cosmology::m_massless_neutrinos
double m_massless_neutrinos
: the effective number (for QED + non-instantaneous decoupling)
Definition: Cosmology.h:291
cbl::cosmology::Cosmology::deltac
double deltac(const double redshift) const
spherical collapse density threshold at a given redshift
Definition: Cosmology.cpp:1248
cbl::cosmology::Cosmology::OmegaR
double OmegaR(const double redshift=0.) const
the radiation density at a given redshift
Definition: Cosmology.cpp:597
cbl::cosmology::Cosmology::Cosmology
Cosmology(const double Omega_matter=0.27, const double Omega_baryon=0.046, const double Omega_neutrinos=0., const double massless_neutrinos=3.04, const int massive_neutrinos=0, const double Omega_DE=0.73, const double Omega_radiation=0., const double hh=0.7, const double scalar_amp=2.46e-9, const double scalar_pivot=0.05, const double n_spec=0.96, const double w0=-1., const double wa=0., const double fNL=0., const int type_NG=1, const double tau=0.09, const std::string model="LCDM", const bool unit=true)
constructor
Definition: Cosmology.cpp:165
cbl::cosmology::Cosmology::Omega_neutrinos
double Omega_neutrinos() const
get the private member Cosmology::m_Omega_neutrinos
Definition: Cosmology.h:1134
cbl::cosmology::Cosmology::dV_dZdOmega
double dV_dZdOmega(const double redshift, const bool angle_rad) const
the derivative of the comoving volume, d2V/(dz*dΩ) at a given redshift
Definition: Cosmology.cpp:1236
cbl::cosmology::Cosmology::m_asn_dz
double m_asn_dz(const double ss) const
the inverse sine amplitude of the Jacobian elliptic function
Definition: Cosmology.cpp:1428
cbl::cosmology::Cosmology::m_Omega_neutrinos
double m_Omega_neutrinos
: the density of massive neutrinos at z=0
Definition: Cosmology.h:288
cbl::cosmology::Cosmology::HH
double HH(const double redshift=0.) const
the Hubble function
Definition: Cosmology.cpp:570
cbl::cosmology::Cosmology::EE2
double EE2(const double redshift=0.) const
auxiliary function used to compute the deceleration parameter
Definition: Cosmology.cpp:964
cbl::cosmology::Cosmology::m_acn_dz
double m_acn_dz(const double cc) const
the inverse cosine amplitude of the Jacobian elliptic function
Definition: Cosmology.cpp:1409
cbl::cosmology::Cosmology::D_A
double D_A(const double redshift) const
the angular diameter distance at a given redshift
Definition: Cosmology.cpp:848
cbl::cosmology::Cosmology::Redshift
double Redshift(const double d_c=1., const double z1_guess=0., const double z2_guess=10., const double prec=0.0001) const
redshift at a given comoving distance
Definition: Cosmology.cpp:1045
cbl::cosmology::Cosmology::print_parameters
void print_parameters() const
print the values of the private members on the screen
Definition: Cosmology.cpp:293
cbl::cosmology::Cosmology::EE_inv2
double EE_inv2(const double redshift=0.) const
inverse of the auxiliary function used to compute the Hubble function, integrand of the lookback time
Definition: Cosmology.h:1919
cbl::cosmology::Cosmology::m_serf_dz
double m_serf_dz(const double yy) const
the inverse truncated series necessary to compute sn-1(s|m) in ASN_DZ
Definition: Cosmology.cpp:1452
cbl::cosmology::Cosmology::f_DE
double f_DE(const double redshift=0.) const
auxiliary function used to compute the Hubble function
Definition: Cosmology.cpp:544
cbl::cosmology::Cosmology::cosmic_time
double cosmic_time(const double redshift=0.) const
cosmic time at a given redshift
Definition: Cosmology.cpp:947
cbl::cosmology::Cosmology::m_Omega_DE
double m_Omega_DE
: the dark energy density at z=0
Definition: Cosmology.h:297
cbl::cosmology::Cosmology::lookback_time
double lookback_time(const double redshift=0.) const
lookback time at a given redshift
Definition: Cosmology.cpp:932
cbl::cosmology::Cosmology::Redshift_LCDM
double Redshift_LCDM(const double d_c=1., const double z1_guess=0., const double z2_guess=10., const bool go_fast=1, const double prec=0.0001) const
redshift at a given comoving distance
Definition: Cosmology.cpp:1101
cbl::cosmology::Cosmology::set_default
void set_default()
internal function to set default values
Definition: Cosmology.cpp:149
cbl::cosmology::Cosmology::Redshift_time
double Redshift_time(const double time, const double z1_guess, const double z2_guess) const
redshift at a given cosmic time
Definition: Cosmology.cpp:1169
cbl::cosmology::Cosmology::EE_inv
double EE_inv(const double redshift=0.) const
inverse of the auxiliary function used to compute the Hubble function integrand of the comoving dista...
Definition: Cosmology.h:1910
cbl::cosmology::Cosmology::Lum_bol
double Lum_bol(const double redshift=0., const double flux=1.) const
bolometric luminosity
Definition: Cosmology.cpp:1036
cbl::cosmology::Cosmology::OmegaR_zeq
double OmegaR_zeq(const double z_eq=3395.) const
the radiation density, as a function of the redshift of radiation-matter equality
Definition: Cosmology.cpp:606
cbl::cosmology::Cosmology::Hdot
double Hdot(const double redshift=0.) const
derivative of the Hubble function at a given redshift
Definition: Cosmology.cpp:985
cbl::cosmology::Cosmology::m_scalar_amp
double m_scalar_amp
: the initial scalar amplitude of the power spectrum
Definition: Cosmology.h:324
cbl::cosmology::Cosmology::m_Omega_matter
double m_Omega_matter
: the density of baryons, cold dark matter and massive neutrinos (in units of the critical density) a...
Definition: Cosmology.h:282
cbl::cosmology::Cosmology::M_vir
double M_vir(const double r_vir, const double redshift, const std::string author="BryanNorman", const bool unit1=false) const
the virial mass, given the virial radius and the redshift
Definition: Cosmology.cpp:1322
cbl::cosmology::Cosmology::D_C
double D_C(const double redshift) const
the comoving line-of-sight distance at a given redshift
Definition: Cosmology.cpp:741
cbl::cosmology::Cosmology::D_C_table
void D_C_table(const std::string file_table, const double z_min, const double z_max, const int step, std::vector< double > &Redshift, std::vector< double > &dc) const
create a table of [redshift, comoving line-of-sight distance]
Definition: Cosmology.cpp:794
cbl::cosmology::Cosmology::D_L
double D_L(const double redshift) const
the luminosity distance at a given redshift
Definition: Cosmology.cpp:875
cbl::cosmology::Cosmology::OmegaK
double OmegaK(const double redshift=0.) const
the density of curvature energy at a given redshift
Definition: Cosmology.cpp:615
cbl::cosmology::Cosmology::linear_growth_rate
double linear_growth_rate(const double redshift, const double prec=1.e-4) const
the linear growth rate at a given redshift,
Definition: Cosmology.cpp:662
cbl::cosmology::Cosmology::r_vir
double r_vir(const double M_vir, const double redshift, const std::string author="BryanNorman", const bool unit1=false) const
the virial radius, given the virial mass and the redshift
Definition: Cosmology.cpp:1331
cbl::cosmology::Cosmology::Distance
double Distance(const double redshift, const std::string distance_type) const
the distance at a given redshift
Definition: Cosmology.cpp:902
cbl::cosmology::Cosmology::OmegaNu
double OmegaNu(const double redshift=0.) const
the neutrino density at a given redshift
Definition: Cosmology.cpp:624
cbl::cosmology::Cosmology::m_massive_neutrinos
int m_massive_neutrinos
the number of degenerate massive neutrino species
Definition: Cosmology.h:294
cbl::cosmology::Cosmology::OmegaM
double OmegaM(const double redshift=0.) const
the matter density at a given redshift
Definition: Cosmology.cpp:579
cbl::cosmology::Cosmology::D_M
double D_M(const double redshift) const
the comoving transverse distance at a given redshift
Definition: Cosmology.cpp:832
cbl::cosmology::Cosmology::m_tau
double m_tau
: Thomson scattering optical depth due to reionization
Definition: Cosmology.h:348
cbl::cosmology::Cosmology::OmegaDE
double OmegaDE(const double redshift=0.) const
the dark energy density at a given redshift
Definition: Cosmology.cpp:588
cbl::cosmology::Cosmology::DD
double DD(const double redshift) const
the amplitude of the growing mode at a given redshift,
Definition: Cosmology.cpp:702
cbl::cosmology::Cosmology::value
double value(const CosmologicalParameter parameter) const
get the private member specified by the enum CosmologicalParameter
Definition: Cosmology.cpp:319
cbl::cosmology::Cosmology::z_eq
double z_eq() const
redshift of matter-dark energy equality
Definition: Cosmology.cpp:1008
cbl::cosmology::Cosmology::concentration_NFW_Duffy
double concentration_NFW_Duffy(const double Mass, const double redshift, const std::string halo_def="vir") const
the halo concentration-mass relation for NFW prfile and Duffy model
Definition: Cosmology.cpp:1461
cbl::par::mega
static const double mega
conversion factor
Definition: Constants.h:78
cbl::par::fDP3
static const char fDP3[]
conversion symbol for: double -> std::string
Definition: Constants.h:136
cbl::par::giga
static const double giga
conversion factor
Definition: Constants.h:75
cbl::par::fINT
static const char fINT[]
conversion symbol for: int -> std::string
Definition: Constants.h:121
cbl::par::kilo
static const double kilo
conversion factor
Definition: Constants.h:81
cbl::par::pi
static const double pi
: the ratio of a circle's circumference to its diameter
Definition: Constants.h:199
cbl::par::pc
static const double pc
: the parsec, defined as 1au/1arc sec [m]
Definition: Constants.h:272
cbl::par::GN
static const double GN
: the Newtonian constant of gravitation [m3 Kg-1 s-2]
Definition: Constants.h:245
cbl::par::cc
static const double cc
: the speed of light in vacuum (the value is exact) [km/sec]
Definition: Constants.h:221
cbl::par::yr
static const double yr
: one year [s]
Definition: Constants.h:278
cbl::par::Msol
static const double Msol
: the solar mass [kg]
Definition: Constants.h:257
cbl::cosmology::CosmologicalParameter
CosmologicalParameter
the cosmological parameters
Definition: Cosmology.h:134
cbl::cosmology::CosmologicalParameter::_sigma8_
@ _sigma8_
: the power spectrum normalisation
cbl::cosmology::CosmologicalParameter::_massive_neutrinos_
@ _massive_neutrinos_
the number of degenerate massive neutrino species
cbl::cosmology::CosmologicalParameter::_hh_
@ _hh_
: the Hubble constant at z=0 divided by 100
cbl::cosmology::CosmologicalParameter::_H0_
@ _H0_
: the Hubble constant at z=0 [km/sec/Mpc]
cbl::cosmology::CosmologicalParameter::_Omega_baryon_h2_
@ _Omega_baryon_h2_
: the baryon density times at z=0
cbl::cosmology::CosmologicalParameter::_ln_scalar_amp_
@ _ln_scalar_amp_
: the logarithm of 1e10 times the initial scalar amplitude of the power spectrum
cbl::cosmology::CosmologicalParameter::_Omega_radiation_
@ _Omega_radiation_
: the radiation density at z=0
cbl::cosmology::CosmologicalParameter::_massless_neutrinos_
@ _massless_neutrinos_
: the effective number (for QED + non-instantaneous decoupling)
cbl::cosmology::CosmologicalParameter::_Omega_matter_LCDM_
@ _Omega_matter_LCDM_
: the density of baryons, cold dark matter and massive neutrinos (in units of the critical density) a...
cbl::cosmology::CosmologicalParameter::_fNL_
@ _fNL_
: the non-Gaussian amplitude
cbl::cosmology::CosmologicalParameter::_neutrino_mass_
@ _neutrino_mass_
the total neutrino mass
cbl::cosmology::CosmologicalParameter::_scalar_pivot_
@ _scalar_pivot_
the scalar pivot k in
cbl::cosmology::CosmologicalParameter::_Omega_baryon_
@ _Omega_baryon_
: the baryon density at z=0
cbl::cosmology::CosmologicalParameter::_scalar_amp_
@ _scalar_amp_
: the initial scalar amplitude of the power spectrum
cbl::cosmology::CosmologicalParameter::_wa_
@ _wa_
: the parameter of the dark energy equation of state (CPL parameterisation)
cbl::cosmology::CosmologicalParameter::_Omega_matter_
@ _Omega_matter_
: the density of baryons, cold dark matter and massive neutrinos (in units of the critical density) a...
cbl::cosmology::CosmologicalParameter::_w0_
@ _w0_
: the parameter of the dark energy equation of state (CPL parameterisation)
cbl::cosmology::CosmologicalParameter::_tau_
@ _tau_
: Thomson scattering optical depth due to reionization
cbl::cosmology::CosmologicalParameter::_Omega_DE_
@ _Omega_DE_
: the dark energy density at z=0
cbl::cosmology::CosmologicalParameter::_Omega_neutrinos_
@ _Omega_neutrinos_
: the density of massive neutrinos at z=0
cbl::cosmology::CosmologicalParameter::_n_spec_
@ _n_spec_
: the primordial spectral index
cbl::cosmology::CosmologicalModel
CosmologicalModel
built-in cosmological models
Definition: Cosmology.h:62
cbl::cosmology::CosmologicalModel::_Planck15_TT_
@ _Planck15_TT_
Planck collaboration 2015, paper XIII: Table 4, TT+lowP+lensing.
cbl::cosmology::CosmologicalModel::_WMAP7_
@ _WMAP7_
Komatsu et al. 2011: Table 1, WMAP Seven-year Mean.
cbl::cosmology::CosmologicalModel::_Planck13_
@ _Planck13_
Planck collaboration 2013, paper XVI: Table 3, Planck+WP.
cbl::cosmology::CosmologicalModel::_WMAP5_
@ _WMAP5_
Komatsu et al. 2009: Table 1, WMAP 5 Year Mean.
cbl::cosmology::CosmologicalModel::_Planck18_
@ _Planck18_
Planck collaboration 2018, Paper VI: Table 2, TT,TE,EE+lowE+lensing.
cbl::cosmology::CosmologicalModel::_Planck15_
@ _Planck15_
Planck collaboration 2015, paper XIII: Table 4, TT,TE,EE+lowP+lensing.
cbl::cosmology::CosmologicalModel::_WMAP9_
@ _WMAP9_
Hinshaw et al. 2013: Table 3, WMAP-only Nine-year.
cbl::cosmology::CosmologicalParameter_name
std::string CosmologicalParameter_name(const CosmologicalParameter parameter)
name of the cosmological parameter
Definition: Cosmology.cpp:45
cbl::wrapper::gsl::GSL_root_brent
double GSL_root_brent(gsl_function Func, const double low_guess, const double up_guess, const double rel_err=1.e-3, const double abs_err=0)
function to find roots using GSL qag method
Definition: GSLwrapper.cpp:430
cbl::wrapper::gsl::GSL_integrate_qag
double GSL_integrate_qag(gsl_function Func, const double a, const double b, const double rel_err=1.e-3, const double abs_err=0, const int limit_size=1000, const int rule=6)
integral, computed using the GSL qag method
Definition: GSLwrapper.cpp:180
cbl
The global namespace of the CosmoBolognaLib
Definition: CAMB.h:38
cbl::conv
std::string conv(const T val, const char *fact)
convert a number to a std::string
Definition: Kernel.h:903
cbl::Mass
T Mass(const T RR, const T Rho)
the mass of a sphere of a given radius and density
Definition: Func.h:1243
cbl::haversine_distance
double haversine_distance(const double ra1, const double ra2, const double dec1, const double dec2)
the haversine angular separation
Definition: Func.cpp:286
cbl::checkIO
void checkIO(const std::ifstream &fin, const std::string file="NULL")
check if an input file can be opened
Definition: Kernel.cpp:160
cbl::ErrorCBL
int ErrorCBL(const std::string msg, const std::string functionCBL, const std::string fileCBL, const cbl::glob::ExitCode exitCode=cbl::glob::ExitCode::_error_)
throw an exception: it is used for handling exceptions inside the CosmoBolognaLib
Definition: Kernel.h:780
cbl::interpolated
double interpolated(const double _xx, const std::vector< double > xx, const std::vector< double > yy, const std::string type)
1D interpolation
Definition: Func.cpp:445
cbl::degrees
double degrees(const double angle, const CoordinateUnits inputUnits=CoordinateUnits::_radians_)
conversion to degrees
Definition: Func.cpp:104
cbl::WarningMsgCBL
void WarningMsgCBL(const std::string msg, const std::string functionCBL, const std::string fileCBL)
internal CBL warning message
Definition: Kernel.h:747