| [3315] | 1 | /* Fitting Formulae for CDM + Baryon + Massive Neutrino (MDM) cosmologies. */
|
|---|
| 2 | /* Daniel J. Eisenstein & Wayne Hu, Institute for Advanced Study */
|
|---|
| 3 |
|
|---|
| 4 | /* There are two primary routines here, one to set the cosmology, the
|
|---|
| 5 | other to construct the transfer function for a single wavenumber k.
|
|---|
| 6 | You should call the former once (per cosmology) and the latter as
|
|---|
| 7 | many times as you want. */
|
|---|
| 8 |
|
|---|
| 9 | /* TFmdm_set_cosm() -- User passes all the cosmological parameters as
|
|---|
| 10 | arguments; the routine sets up all of the scalar quantites needed
|
|---|
| 11 | computation of the fitting formula. The input parameters are:
|
|---|
| 12 | 1) omega_matter -- Density of CDM, baryons, and massive neutrinos,
|
|---|
| 13 | in units of the critical density.
|
|---|
| 14 | 2) omega_baryon -- Density of baryons, in units of critical.
|
|---|
| 15 | 3) omega_hdm -- Density of massive neutrinos, in units of critical
|
|---|
| 16 | 4) degen_hdm -- (Int) Number of degenerate massive neutrino species
|
|---|
| 17 | 5) omega_lambda -- Cosmological constant
|
|---|
| 18 | 6) hubble -- Hubble constant, in units of 100 km/s/Mpc
|
|---|
| 19 | 7) redshift -- The redshift at which to evaluate */
|
|---|
| 20 |
|
|---|
| 21 | /* TFmdm_onek_mpc() -- User passes a single wavenumber, in units of Mpc^-1.
|
|---|
| 22 | Routine returns the transfer function from the Eisenstein & Hu
|
|---|
| 23 | fitting formula, based on the cosmology currently held in the
|
|---|
| 24 | internal variables. The routine returns T_cb (the CDM+Baryon
|
|---|
| 25 | density-weighted transfer function), although T_cbn (the CDM+
|
|---|
| 26 | Baryon+Neutrino density-weighted transfer function) is stored
|
|---|
| 27 | in the global variable tf_cbnu. */
|
|---|
| 28 |
|
|---|
| 29 | /* We also supply TFmdm_onek_hmpc(), which is identical to the previous
|
|---|
| 30 | routine, but takes the wavenumber in units of h Mpc^-1. */
|
|---|
| 31 |
|
|---|
| 32 | /* We hold the internal scalar quantities in global variables, so that
|
|---|
| 33 | the user may access them in an external program, via "extern" declarations. */
|
|---|
| 34 |
|
|---|
| 35 | /* Please note that all internal length scales are in Mpc, not h^-1 Mpc! */
|
|---|
| 36 |
|
|---|
| 37 | /* -------------------------- Prototypes ----------------------------- */
|
|---|
| 38 |
|
|---|
| 39 | int TFmdm_set_cosm(double omega_matter, double omega_baryon, double omega_hdm,
|
|---|
| 40 | int degen_hdm, double omega_lambda, double hubble, double redshift);
|
|---|
| 41 | double TFmdm_onek_mpc(double kk);
|
|---|
| 42 | double TFmdm_onek_hmpc(double kk);
|
|---|
| 43 |
|
|---|
| 44 | #include <math.h>
|
|---|
| 45 | #include <stdio.h>
|
|---|
| 46 | #include <stdlib.h>
|
|---|
| 47 |
|
|---|
| 48 | /* Convenience from Numerical Recipes in C, 2nd edition
|
|---|
| 49 | ------------ COMPLETEMENT DEBILE (cmv) ????????? ------------
|
|---|
| 50 | static double sqrarg;
|
|---|
| 51 | #define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
|
|---|
| 52 | */
|
|---|
| 53 | #define SQR(a) (a*a)
|
|---|
| 54 |
|
|---|
| 55 | /* ------------------------- Global Variables ------------------------ */
|
|---|
| 56 |
|
|---|
| 57 | /* The following are set in TFmdm_set_cosm() */
|
|---|
| 58 | double alpha_gamma, /* sqrt(alpha_nu) */
|
|---|
| 59 | alpha_nu, /* The small-scale suppression */
|
|---|
| 60 | beta_c, /* The correction to the log in the small-scale */
|
|---|
| 61 | num_degen_hdm, /* Number of degenerate massive neutrino species */
|
|---|
| 62 | f_baryon, /* Baryon fraction */
|
|---|
| 63 | f_bnu, /* Baryon + Massive Neutrino fraction */
|
|---|
| 64 | f_cb, /* Baryon + CDM fraction */
|
|---|
| 65 | f_cdm, /* CDM fraction */
|
|---|
| 66 | f_hdm, /* Massive Neutrino fraction */
|
|---|
| 67 | growth_k0, /* D_1(z) -- the growth function as k->0 */
|
|---|
| 68 | growth_to_z0, /* D_1(z)/D_1(0) -- the growth relative to z=0 */
|
|---|
| 69 | hhubble, /* Need to pass Hubble constant to TFmdm_onek_hmpc() */
|
|---|
| 70 | k_equality, /* The comoving wave number of the horizon at equality*/
|
|---|
| 71 | obhh, /* Omega_baryon * hubble^2 */
|
|---|
| 72 | omega_curv, /* = 1 - omega_matter - omega_lambda */
|
|---|
| 73 | omega_lambda_z, /* Omega_lambda at the given redshift */
|
|---|
| 74 | omega_matter_z, /* Omega_matter at the given redshift */
|
|---|
| 75 | omhh, /* Omega_matter * hubble^2 */
|
|---|
| 76 | onhh, /* Omega_hdm * hubble^2 */
|
|---|
| 77 | p_c, /* The correction to the exponent before drag epoch */
|
|---|
| 78 | p_cb, /* The correction to the exponent after drag epoch */
|
|---|
| 79 | sound_horizon_fit, /* The sound horizon at the drag epoch */
|
|---|
| 80 | theta_cmb, /* The temperature of the CMB, in units of 2.7 K */
|
|---|
| 81 | y_drag, /* Ratio of z_equality to z_drag */
|
|---|
| 82 | z_drag, /* Redshift of the drag epoch */
|
|---|
| 83 | z_equality; /* Redshift of matter-radiation equality */
|
|---|
| 84 |
|
|---|
| 85 | /* The following are set in TFmdm_onek_mpc() */
|
|---|
| 86 | double gamma_eff, /* Effective \Gamma */
|
|---|
| 87 | growth_cb, /* Growth factor for CDM+Baryon perturbations */
|
|---|
| 88 | growth_cbnu, /* Growth factor for CDM+Baryon+Neutrino pert. */
|
|---|
| 89 | max_fs_correction, /* Correction near maximal free streaming */
|
|---|
| 90 | qq, /* Wavenumber rescaled by \Gamma */
|
|---|
| 91 | qq_eff, /* Wavenumber rescaled by effective Gamma */
|
|---|
| 92 | qq_nu, /* Wavenumber compared to maximal free streaming */
|
|---|
| 93 | tf_master, /* Master TF */
|
|---|
| 94 | tf_sup, /* Suppressed TF */
|
|---|
| 95 | y_freestream; /* The epoch of free-streaming for a given scale */
|
|---|
| 96 |
|
|---|
| 97 | /* Finally, TFmdm_onek_mpc() and TFmdm_onek_hmpc() give their answers as */
|
|---|
| 98 | double tf_cb, /* The transfer function for density-weighted
|
|---|
| 99 | CDM + Baryon perturbations. */
|
|---|
| 100 | tf_cbnu; /* The transfer function for density-weighted
|
|---|
| 101 | CDM + Baryon + Massive Neutrino perturbations. */
|
|---|
| 102 |
|
|---|
| 103 | /* By default, these functions return tf_cb */
|
|---|
| 104 |
|
|---|
| 105 | /* ------------------------- TFmdm_set_cosm() ------------------------ */
|
|---|
| 106 | int TFmdm_set_cosm(double omega_matter, double omega_baryon, double omega_hdm,
|
|---|
| 107 | int degen_hdm, double omega_lambda, double hubble, double redshift)
|
|---|
| 108 | /* This routine takes cosmological parameters and a redshift and sets up
|
|---|
| 109 | all the internal scalar quantities needed to compute the transfer function. */
|
|---|
| 110 | /* INPUT: omega_matter -- Density of CDM, baryons, and massive neutrinos,
|
|---|
| 111 | in units of the critical density. */
|
|---|
| 112 | /* omega_baryon -- Density of baryons, in units of critical. */
|
|---|
| 113 | /* omega_hdm -- Density of massive neutrinos, in units of critical */
|
|---|
| 114 | /* degen_hdm -- (Int) Number of degenerate massive neutrino species */
|
|---|
| 115 | /* omega_lambda -- Cosmological constant */
|
|---|
| 116 | /* hubble -- Hubble constant, in units of 100 km/s/Mpc */
|
|---|
| 117 | /* redshift -- The redshift at which to evaluate */
|
|---|
| 118 | /* OUTPUT: Returns 0 if all is well, 1 if a warning was issued. Otherwise,
|
|---|
| 119 | sets many global variables for use in TFmdm_onek_mpc() */
|
|---|
| 120 | {
|
|---|
| 121 | double z_drag_b1, z_drag_b2, omega_denom;
|
|---|
| 122 | int qwarn;
|
|---|
| 123 | qwarn = 0;
|
|---|
| 124 |
|
|---|
| 125 | theta_cmb = 2.728/2.7; /* Assuming T_cmb = 2.728 K */
|
|---|
| 126 |
|
|---|
| 127 | /* Look for strange input */
|
|---|
| 128 | if (omega_baryon<0.0) {
|
|---|
| 129 | fprintf(stderr,
|
|---|
| 130 | "TFmdm_set_cosm(): Negative omega_baryon set to trace amount.\n");
|
|---|
| 131 | qwarn = 1;
|
|---|
| 132 | }
|
|---|
| 133 | if (omega_hdm<0.0) {
|
|---|
| 134 | fprintf(stderr,
|
|---|
| 135 | "TFmdm_set_cosm(): Negative omega_hdm set to trace amount.\n");
|
|---|
| 136 | qwarn = 1;
|
|---|
| 137 | }
|
|---|
| 138 | if (hubble<=0.0) {
|
|---|
| 139 | fprintf(stderr,"TFmdm_set_cosm(): Negative Hubble constant illegal.\n");
|
|---|
| 140 | exit(1); /* Can't recover */
|
|---|
| 141 | } else if (hubble>2.0) {
|
|---|
| 142 | fprintf(stderr,"TFmdm_set_cosm(): Hubble constant should be in units of 100 km/s/Mpc.\n");
|
|---|
| 143 | qwarn = 1;
|
|---|
| 144 | }
|
|---|
| 145 | if (redshift<=-1.0) {
|
|---|
| 146 | fprintf(stderr,"TFmdm_set_cosm(): Redshift < -1 is illegal.\n");
|
|---|
| 147 | exit(1);
|
|---|
| 148 | } else if (redshift>99.0) {
|
|---|
| 149 | fprintf(stderr,
|
|---|
| 150 | "TFmdm_set_cosm(): Large redshift entered. TF may be inaccurate.\n");
|
|---|
| 151 | qwarn = 1;
|
|---|
| 152 | }
|
|---|
| 153 | if (degen_hdm<1) degen_hdm=1;
|
|---|
| 154 | num_degen_hdm = (double) degen_hdm;
|
|---|
| 155 | /* Have to save this for TFmdm_onek_mpc() */
|
|---|
| 156 | /* This routine would crash if baryons or neutrinos were zero,
|
|---|
| 157 | so don't allow that */
|
|---|
| 158 | if (omega_baryon<=0) omega_baryon=1e-5;
|
|---|
| 159 | if (omega_hdm<=0) omega_hdm=1e-5;
|
|---|
| 160 |
|
|---|
| 161 | omega_curv = 1.0-omega_matter-omega_lambda;
|
|---|
| 162 | omhh = omega_matter*SQR(hubble);
|
|---|
| 163 | obhh = omega_baryon*SQR(hubble);
|
|---|
| 164 | onhh = omega_hdm*SQR(hubble);
|
|---|
| 165 | f_baryon = omega_baryon/omega_matter;
|
|---|
| 166 | f_hdm = omega_hdm/omega_matter;
|
|---|
| 167 | f_cdm = 1.0-f_baryon-f_hdm;
|
|---|
| 168 | f_cb = f_cdm+f_baryon;
|
|---|
| 169 | f_bnu = f_baryon+f_hdm;
|
|---|
| 170 |
|
|---|
| 171 | /* Compute the equality scale. */
|
|---|
| 172 | z_equality = 25000.0*omhh/SQR(SQR(theta_cmb)); /* Actually 1+z_eq */
|
|---|
| 173 | k_equality = 0.0746*omhh/SQR(theta_cmb);
|
|---|
| 174 |
|
|---|
| 175 | /* Compute the drag epoch and sound horizon */
|
|---|
| 176 | z_drag_b1 = 0.313*pow(omhh,-0.419)*(1+0.607*pow(omhh,0.674));
|
|---|
| 177 | z_drag_b2 = 0.238*pow(omhh,0.223);
|
|---|
| 178 | z_drag = 1291*pow(omhh,0.251)/(1.0+0.659*pow(omhh,0.828))*
|
|---|
| 179 | (1.0+z_drag_b1*pow(obhh,z_drag_b2));
|
|---|
| 180 | y_drag = z_equality/(1.0+z_drag);
|
|---|
| 181 |
|
|---|
| 182 | sound_horizon_fit = 44.5*log(9.83/omhh)/sqrt(1.0+10.0*pow(obhh,0.75));
|
|---|
| 183 |
|
|---|
| 184 | /* Set up for the free-streaming & infall growth function */
|
|---|
| 185 | p_c = 0.25*(5.0-sqrt(1+24.0*f_cdm));
|
|---|
| 186 | p_cb = 0.25*(5.0-sqrt(1+24.0*f_cb));
|
|---|
| 187 |
|
|---|
| 188 | omega_denom = omega_lambda+SQR(1.0+redshift)*(omega_curv+
|
|---|
| 189 | omega_matter*(1.0+redshift));
|
|---|
| 190 | omega_lambda_z = omega_lambda/omega_denom;
|
|---|
| 191 | omega_matter_z = omega_matter*SQR(1.0+redshift)*(1.0+redshift)/omega_denom;
|
|---|
| 192 | growth_k0 = z_equality/(1.0+redshift)*2.5*omega_matter_z/
|
|---|
| 193 | (pow(omega_matter_z,4.0/7.0)-omega_lambda_z+
|
|---|
| 194 | (1.0+omega_matter_z/2.0)*(1.0+omega_lambda_z/70.0));
|
|---|
| 195 | growth_to_z0 = z_equality*2.5*omega_matter/(pow(omega_matter,4.0/7.0)
|
|---|
| 196 | -omega_lambda + (1.0+omega_matter/2.0)*(1.0+omega_lambda/70.0));
|
|---|
| 197 | growth_to_z0 = growth_k0/growth_to_z0;
|
|---|
| 198 |
|
|---|
| 199 | /* Compute small-scale suppression */
|
|---|
| 200 | alpha_nu = f_cdm/f_cb*(5.0-2.*(p_c+p_cb))/(5.-4.*p_cb)*
|
|---|
| 201 | pow(1+y_drag,p_cb-p_c)*
|
|---|
| 202 | (1+f_bnu*(-0.553+0.126*f_bnu*f_bnu))/
|
|---|
| 203 | (1-0.193*sqrt(f_hdm*num_degen_hdm)+0.169*f_hdm*pow(num_degen_hdm,0.2))*
|
|---|
| 204 | (1+(p_c-p_cb)/2*(1+1/(3.-4.*p_c)/(7.-4.*p_cb))/(1+y_drag));
|
|---|
| 205 | alpha_gamma = sqrt(alpha_nu);
|
|---|
| 206 | beta_c = 1/(1-0.949*f_bnu);
|
|---|
| 207 | /* Done setting scalar variables */
|
|---|
| 208 | hhubble = hubble; /* Need to pass Hubble constant to TFmdm_onek_hmpc() */
|
|---|
| 209 | return qwarn;
|
|---|
| 210 | }
|
|---|
| 211 |
|
|---|
| 212 | /* ---------------------------- TFmdm_onek_mpc() ---------------------- */
|
|---|
| 213 |
|
|---|
| 214 | double TFmdm_onek_mpc(double kk)
|
|---|
| 215 | /* Given a wavenumber in Mpc^-1, return the transfer function for the
|
|---|
| 216 | cosmology held in the global variables. */
|
|---|
| 217 | /* Input: kk -- Wavenumber in Mpc^-1 */
|
|---|
| 218 | /* Output: The following are set as global variables:
|
|---|
| 219 | growth_cb -- the transfer function for density-weighted
|
|---|
| 220 | CDM + Baryon perturbations.
|
|---|
| 221 | growth_cbnu -- the transfer function for density-weighted
|
|---|
| 222 | CDM + Baryon + Massive Neutrino perturbations. */
|
|---|
| 223 | /* The function returns growth_cb */
|
|---|
| 224 | {
|
|---|
| 225 | double tf_sup_L, tf_sup_C;
|
|---|
| 226 | double temp1, temp2;
|
|---|
| 227 |
|
|---|
| 228 | qq = kk/omhh*SQR(theta_cmb);
|
|---|
| 229 |
|
|---|
| 230 | /* Compute the scale-dependent growth functions */
|
|---|
| 231 | y_freestream = 17.2*f_hdm*(1+0.488*pow(f_hdm,-7.0/6.0))*
|
|---|
| 232 | SQR(num_degen_hdm*qq/f_hdm);
|
|---|
| 233 | temp1 = pow(growth_k0, 1.0-p_cb);
|
|---|
| 234 | temp2 = pow(growth_k0/(1+y_freestream),0.7);
|
|---|
| 235 | growth_cb = pow(1.0+temp2, p_cb/0.7)*temp1;
|
|---|
| 236 | growth_cbnu = pow(pow(f_cb,0.7/p_cb)+temp2, p_cb/0.7)*temp1;
|
|---|
| 237 |
|
|---|
| 238 | /* Compute the master function */
|
|---|
| 239 | gamma_eff =omhh*(alpha_gamma+(1-alpha_gamma)/
|
|---|
| 240 | (1+SQR(SQR(kk*sound_horizon_fit*0.43))));
|
|---|
| 241 | qq_eff = qq*omhh/gamma_eff;
|
|---|
| 242 |
|
|---|
| 243 | tf_sup_L = log(2.71828+1.84*beta_c*alpha_gamma*qq_eff);
|
|---|
| 244 | tf_sup_C = 14.4+325/(1+60.5*pow(qq_eff,1.11));
|
|---|
| 245 | tf_sup = tf_sup_L/(tf_sup_L+tf_sup_C*SQR(qq_eff));
|
|---|
| 246 |
|
|---|
| 247 | qq_nu = 3.92*qq*sqrt(num_degen_hdm/f_hdm);
|
|---|
| 248 | max_fs_correction = 1+1.2*pow(f_hdm,0.64)*pow(num_degen_hdm,0.3+0.6*f_hdm)/
|
|---|
| 249 | (pow(qq_nu,-1.6)+pow(qq_nu,0.8));
|
|---|
| 250 | tf_master = tf_sup*max_fs_correction;
|
|---|
| 251 |
|
|---|
| 252 | /* Now compute the CDM+HDM+baryon transfer functions */
|
|---|
| 253 | tf_cb = tf_master*growth_cb/growth_k0;
|
|---|
| 254 | tf_cbnu = tf_master*growth_cbnu/growth_k0;
|
|---|
| 255 | return tf_cb;
|
|---|
| 256 | }
|
|---|
| 257 |
|
|---|
| 258 | /* ---------------------------- TFmdm_onek_hmpc() ---------------------- */
|
|---|
| 259 |
|
|---|
| 260 | double TFmdm_onek_hmpc(double kk)
|
|---|
| 261 | /* Given a wavenumber in h Mpc^-1, return the transfer function for the
|
|---|
| 262 | cosmology held in the global variables. */
|
|---|
| 263 | /* Input: kk -- Wavenumber in h Mpc^-1 */
|
|---|
| 264 | /* Output: The following are set as global variables:
|
|---|
| 265 | growth_cb -- the transfer function for density-weighted
|
|---|
| 266 | CDM + Baryon perturbations.
|
|---|
| 267 | growth_cbnu -- the transfer function for density-weighted
|
|---|
| 268 | CDM + Baryon + Massive Neutrino perturbations. */
|
|---|
| 269 | /* The function returns growth_cb */
|
|---|
| 270 | {
|
|---|
| 271 | return TFmdm_onek_mpc(kk*hhubble);
|
|---|
| 272 | }
|
|---|
| 273 |
|
|---|
| 274 |
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
