CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
EisensteinHu.h
Go to the documentation of this file.
1 /********************************************************************
2  * Copyright (C) 2010 by Federico Marulli, Alfonso Veropalumbo *
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 
35 #ifndef __EH__
36 #define __EH__
37 
38 
39 namespace cbl {
40 
41  namespace cosmology {
42 
54  class EisensteinHu {
55 
56  protected:
57 
59  double alpha_gamma;
60 
62  double alpha_nu;
63 
65  double beta_c;
66 
68  double num_degen_hdm;
69 
71  double f_baryon;
72 
74  double f_bnu;
75 
77  double f_cb;
78 
80  double f_cdm;
81 
83  double f_hdm;
84 
86  double growth_k0;
87 
89  double growth_to_z0;
90 
92  double hhubble;
93 
95  double k_equality;
96 
98  double obhh;
99 
101  double omega_curv;
102 
104  double omega_lambda_z;
105 
107  double omega_matter_z;
108 
110  double omhh;
111 
113  double onhh;
114 
116  double p_c;
117 
119  double p_cb;
120 
123 
125  double theta_cmb;
126 
128  double y_drag;
129 
131  double z_drag;
132 
134  double z_equality;
135 
137  double gamma_eff;
138 
140  double growth_cb;
141 
143  double growth_cbnu;
144 
147 
149  double qq;
150 
152  double qq_eff;
153 
155  double qq_nu;
156 
158  double tf_master;
159 
161  double tf_sup;
162 
164  double y_freestream;
165 
167  double tf_cb;
168 
170  double tf_cbnu;
171 
174 
176  double ns;
177 
178  public:
179 
184 
188  EisensteinHu () = default;
189 
193  ~EisensteinHu () = default;
194 
196 
197 
231  int TFmdm_set_cosm (double omega_matter, double omega_baryon, double omega_hdm, int degen_hdm, double omega_lambda, double hubble, double redshift, double As=2.56e-9, double k_pivot=0.05, double n_spec=0.96)
232  {
233  double z_drag_b1, z_drag_b2, omega_denom;
234  int qwarn;
235  qwarn = 0;
236 
237  theta_cmb = 2.728/2.7; // assuming T_cmb = 2.728 K
238 
239  /* Look for strange input */
240 
241  if (omega_baryon<0.0) {
242  fprintf(stderr,
243  "TFmdm_set_cosm(): Negative omega_baryon set to trace amount.\n");
244  qwarn = 1;
245  }
246 
247  if (omega_hdm<0.0) {
248  fprintf(stderr,
249  "TFmdm_set_cosm(): Negative omega_hdm set to trace amount.\n");
250  qwarn = 1;
251  }
252 
253  if (hubble<=0.0) {
254  fprintf(stderr,"TFmdm_set_cosm(): Negative Hubble constant illegal.\n");
255  exit(1); /* Can't recover */
256  } else if (hubble>2.0) {
257  fprintf(stderr,"TFmdm_set_cosm(): Hubble constant should be in units of 100 km/s/Mpc.\n");
258  qwarn = 1;
259  }
260 
261  if (redshift<=-1.0) {
262  fprintf(stderr,"TFmdm_set_cosm(): Redshift < -1 is illegal.\n");
263  exit(1);
264  } else if (redshift>99.0) {
265  fprintf(stderr,
266  "TFmdm_set_cosm(): Large redshift entered. TF may be inaccurate.\n");
267  qwarn = 1;
268  }
269 
270  if (degen_hdm<1) degen_hdm = 1;
271  num_degen_hdm = (float) degen_hdm;
272  /* Have to save this for TFmdm_onek_mpc() */
273  /* This routine would crash if baryons or neutrinos were zero,
274  so don't allow that */
275  if (omega_baryon<=0) omega_baryon=1e-5;
276  if (omega_hdm<=0) omega_hdm=1e-5;
277 
278  omega_curv = 1.0-omega_matter-omega_lambda;
279  omhh = omega_matter*pow(hubble, 2);
280  obhh = omega_baryon*pow(hubble, 2);
281  onhh = omega_hdm*pow(hubble, 2);
282  f_baryon = omega_baryon/omega_matter;
283  f_hdm = omega_hdm/omega_matter;
284  f_cdm = 1.0-f_baryon-f_hdm;
285  f_cb = f_cdm+f_baryon;
286  f_bnu = f_baryon+f_hdm;
287 
288  /* Compute the equality scale. */
289  z_equality = 25000.0*omhh*pow(theta_cmb, -4); /* Actually 1+z_eq */
290  k_equality = 0.0746*omhh*pow(theta_cmb, -2);
291 
292  /* Compute the drag epoch and sound horizon */
293  z_drag_b1 = 0.313*pow(omhh,-0.419)*(1+0.607*pow(omhh,0.674));
294  z_drag_b2 = 0.238*pow(omhh,0.223);
295  z_drag = 1291*pow(omhh,0.251)/(1.0+0.659*pow(omhh,0.828))*(1.0+z_drag_b1*pow(obhh,z_drag_b2));
296  y_drag = z_equality/(1.0+z_drag);
297 
298  sound_horizon_fit = 44.5*log(9.83/omhh)/sqrt(1.0+10.0*pow(obhh,0.75));
299 
300  /* Set up for the free-streaming & infall growth function */
301  p_c = 0.25*(5.0-sqrt(1+24.0*f_cdm));
302  p_cb = 0.25*(5.0-sqrt(1+24.0*f_cb));
303 
304  omega_denom = omega_lambda+pow(1.0+redshift, 2)*(omega_curv+omega_matter*(1.0+redshift));
305  omega_lambda_z = omega_lambda/omega_denom;
306  omega_matter_z = omega_matter*pow(1.0+redshift,2)*(1.0+redshift)/omega_denom;
307 
308  growth_k0 = z_equality/(1.0+redshift)*2.5*omega_matter_z/(pow(omega_matter_z,4.0/7.0)-omega_lambda_z+(1.0+omega_matter_z/2.0)*(1.0+omega_lambda_z/70.0));
309  growth_to_z0 = z_equality*2.5*omega_matter/(pow(omega_matter,4.0/7.0)-omega_lambda + (1.0+omega_matter/2.0)*(1.0+omega_lambda/70.0));
310  double growth_z = 1./(1.0+redshift)*2.5*omega_matter_z/(pow(omega_matter_z,4.0/7.0)-omega_lambda_z + (1.0+omega_matter_z/2.0)*(1.0+omega_lambda_z/70.0));
312 
313  /*
314  growth_k0 = z_equality*5.*omega_matter_z/(2*(1+redshift))/(1./70.+209./140.*omega_matter_z-pow(omega_matter_z,2)/140.+pow(omega_matter_z,4./7.));
315  growth_to_z0 = z_equality*5.*omega_matter_z/(2*(1+redshift))/(1./70.+209./140.*omega_matter_z-pow(omega_matter_z,2)/140.+pow(omega_matter_z,4./7.));
316  growth_to_z0 = growth_k0/growth_to_z0;
317  */
318 
319  /* Compute small-scale suppression */
320  alpha_nu = f_cdm/f_cb*(5.0-2.*(p_c+p_cb))/(5.-4.*p_cb)*
321  pow(1+y_drag,p_cb-p_c)*
322  (1+f_bnu*(-0.553+0.126*f_bnu*f_bnu))/
323  (1-0.193*sqrt(f_hdm*num_degen_hdm)+0.169*f_hdm*pow(num_degen_hdm,0.2))*
324  (1+(p_c-p_cb)/2*(1+1/(3.-4.*p_c)/(7.-4.*p_cb))/(1+y_drag));
325  alpha_gamma = sqrt(alpha_nu);
326  beta_c = 1/(1-0.949*f_bnu);
327 
328  /* Done setting scalar variables */
329  hhubble = hubble; /* Need to pass Hubble constant to TFmdm_onek_hmpc() */
330 
331  ns = n_spec;
332 
333  //double norm = 2.*par::pi*par::pi*m_scalar_amp*pow(2.*pow(par::cc/100.,2)/(5*m_Omega_matter)*DN(redshift),2)*pow(k/m_scalar_pivot, m_n_spec-1)*k_ov_h*pow(fact, -3);
334  pk_normalization = 2.*par::pi*par::pi*As*pow(2.*pow(par::cc/100.,2)/(5*omega_matter)*growth_z,2)*pow(1./k_pivot, ns-1);
335 
336  return qwarn;
337  }
338 
351  double TFmdm_onek_mpc (double kk)
352  {
353  double tf_sup_L, tf_sup_C;
354  double temp1, temp2;
355 
356  qq = kk/omhh*pow(theta_cmb, 2);
357 
358  /* Compute the scale-dependent growth functions */
359  y_freestream = 17.2*f_hdm*(1+0.488*pow(f_hdm,-7.0/6.0))*pow(num_degen_hdm*qq/f_hdm, 2);
360  temp1 = pow(growth_k0, 1.0-p_cb);
361  temp2 = pow(growth_k0/(1+y_freestream),0.7);
362  growth_cb = pow(1.0+temp2, p_cb/0.7)*temp1;
363  growth_cbnu = pow(pow(f_cb,0.7/p_cb)+temp2, p_cb/0.7)*temp1;
364 
365  /* Compute the master function */
366  gamma_eff =omhh*(alpha_gamma+(1-alpha_gamma)/(1+pow(kk*sound_horizon_fit*0.43, 4)));
368 
369  tf_sup_L = log(2.71828+1.84*beta_c*alpha_gamma*qq_eff);
370  tf_sup_C = 14.4+325/(1+60.5*pow(qq_eff,1.11));
371  tf_sup = tf_sup_L/(tf_sup_L+tf_sup_C*pow(qq_eff, 2));
372 
373  qq_nu = 3.92*qq*sqrt(num_degen_hdm/f_hdm);
374  max_fs_correction = 1+1.2*pow(f_hdm,0.64)*pow(num_degen_hdm,0.3+0.6*f_hdm)/(pow(qq_nu,-1.6)+pow(qq_nu,0.8));
376 
377  /* Now compute the CDM+HDM+baryon transfer functions */
380 
381  return tf_cb;
382  }
383 
395  double TFmdm_onek_hmpc (double kk)
396  {
397  return TFmdm_onek_mpc(kk*hhubble);
398  }
399 
407  double Pk (double kk)
408  {
409  double k = kk*hhubble;
410 
411  double func = TFmdm_onek_mpc(k); // transfer function
412 
413  return pk_normalization*pow(func,2)*kk*pow(k, ns-1);
414  }
415 
424  std::vector<double> Pk(std::vector<double> kk) {
425 
426  std::vector<double> _Pk(kk.size(), 0);
427 
428  for (size_t i=0; i<kk.size(); i++)
429  _Pk[i] = Pk(kk[i]);
430 
431  return _Pk;
432  }
433 
434  };
435 
436  }
437 }
438 
439 #endif
The class EisensteinHu.
Definition: EisensteinHu.h:54
double f_cdm
CDM fraction.
Definition: EisensteinHu.h:80
double f_cb
baryon + CDM fraction
Definition: EisensteinHu.h:77
double growth_cbnu
growth factor for CDM+Baryon+Neutrino pert.
Definition: EisensteinHu.h:143
double omhh
Omega_matter * hubble^2.
Definition: EisensteinHu.h:110
double z_equality
redshift of matter-radiation equality
Definition: EisensteinHu.h:134
double num_degen_hdm
number of degenerate massive neutrino species
Definition: EisensteinHu.h:68
double pk_normalization
the power spectrum normalization
Definition: EisensteinHu.h:173
double omega_curv
= 1 - omega_matter - omega_lambda
Definition: EisensteinHu.h:101
double p_c
the correction to the exponent before drag epoch
Definition: EisensteinHu.h:116
double TFmdm_onek_mpc(double kk)
compute the transfer function
Definition: EisensteinHu.h:351
double z_drag
redshift of the drag epoch
Definition: EisensteinHu.h:131
double alpha_gamma
sqrt(alpha_nu)
Definition: EisensteinHu.h:59
double ns
the spectral index
Definition: EisensteinHu.h:176
double tf_sup
suppressed TF
Definition: EisensteinHu.h:161
double growth_cb
growth factor for CDM+Baryon perturbations
Definition: EisensteinHu.h:140
double y_drag
ratio of z_equality to z_drag
Definition: EisensteinHu.h:128
double k_equality
the comoving wave number of the horizon at equality
Definition: EisensteinHu.h:95
double growth_k0
D_1(z) – the growth function as k->0.
Definition: EisensteinHu.h:86
EisensteinHu()=default
default constructor
double omega_lambda_z
Omega_lambda at the given redshift.
Definition: EisensteinHu.h:104
double sound_horizon_fit
the sound horizon at the drag epoch
Definition: EisensteinHu.h:122
double omega_matter_z
Omega_matter at the given redshift.
Definition: EisensteinHu.h:107
double Pk(double kk)
get the power spectrum in unit of h^{-3}Mpc^{-3}
Definition: EisensteinHu.h:407
std::vector< double > Pk(std::vector< double > kk)
get the power spectrum in unit of h^{-3}Mpc^{-3}
Definition: EisensteinHu.h:424
double max_fs_correction
correction near maximal free streaming
Definition: EisensteinHu.h:146
double theta_cmb
the temperature of the CMB; in units of 2.7 K
Definition: EisensteinHu.h:125
double beta_c
the correction to the log in the small-scale
Definition: EisensteinHu.h:65
double f_hdm
massive Neutrino fraction
Definition: EisensteinHu.h:83
double qq
wavenumber rescaled by
Definition: EisensteinHu.h:149
double hhubble
need to pass Hubble constant to TFmdm_onek_hmpc()
Definition: EisensteinHu.h:92
~EisensteinHu()=default
default destructor
double f_baryon
baryon fraction
Definition: EisensteinHu.h:71
double y_freestream
the epoch of free-streaming for a given scale
Definition: EisensteinHu.h:164
double f_bnu
baryon + Massive Neutrino fraction
Definition: EisensteinHu.h:74
double qq_nu
wavenumber compared to maximal free streaming
Definition: EisensteinHu.h:155
double TFmdm_onek_hmpc(double kk)
compute the transfer function
Definition: EisensteinHu.h:395
int TFmdm_set_cosm(double omega_matter, double omega_baryon, double omega_hdm, int degen_hdm, double omega_lambda, double hubble, double redshift, double As=2.56e-9, double k_pivot=0.05, double n_spec=0.96)
set the cosmological parameters
Definition: EisensteinHu.h:231
double qq_eff
wavenumber rescaled by effective Gamma
Definition: EisensteinHu.h:152
double alpha_nu
the small-scale suppression
Definition: EisensteinHu.h:62
double onhh
Omega_hdm * hubble^2.
Definition: EisensteinHu.h:113
double p_cb
the correction to the exponent after drag epoch
Definition: EisensteinHu.h:119
double obhh
Omega_baryon * hubble^2.
Definition: EisensteinHu.h:98
double tf_cbnu
the transfer function for density-weighted CDM + Baryon + Massive Neutrino perturbations.
Definition: EisensteinHu.h:170
double tf_cb
the transfer function for density-weighted CDM + Baryon perturbations.
Definition: EisensteinHu.h:167
double growth_to_z0
D_1(z)/D_1(0) – the growth relative to z=0.
Definition: EisensteinHu.h:89
static const double pi
: the ratio of a circle's circumference to its diameter
Definition: Constants.h:199
static const double cc
: the speed of light in vacuum (the value is exact) [km/sec]
Definition: Constants.h:221
The global namespace of the CosmoBolognaLib
Definition: CAMB.h:38