/* The following routines implement all of the fitting formulae in Eisenstein \& Hu (1997) */ /* There are two sets of routines here. The first set, TFfit_hmpc(), TFset_parameters(), and TFfit_onek(), calculate the transfer function for an arbitrary CDM+baryon universe using the fitting formula in Section 3 of the paper. The second set, TFsound_horizon_fit(), TFk_peak(), TFnowiggles(), and TFzerobaryon(), calculate other quantities given in Section 4 of the paper. */ #include #include #include void TFset_parameters(double omega0hh, double f_baryon, double Tcmb); double TFfit_onek(double k, double *tf_baryon, double *tf_cdm); void TFfit_hmpc(double omega0, double f_baryon, double hubble, double Tcmb, int numk, double *k, double *tf_full, double *tf_baryon, double *tf_cdm); double TFsound_horizon_fit(double omega0, double f_baryon, double hubble); double TFk_peak(double omega0, double f_baryon, double hubble); double TFnowiggles(double omega0, double f_baryon, double hubble, double Tcmb, double k_hmpc); double TFzerobaryon(double omega0, double hubble, double Tcmb, double k_hmpc); /* ------------------------ DRIVER ROUTINE --------------------------- */ /* The following is an example of a driver routine you might use. */ /* Basically, the driver routine needs to call TFset_parameters() to set all the scalar parameters, and then call TFfit_onek() for each wavenumber k you desire. */ /* While the routines use Mpc^-1 units internally, this driver has been written to take an array of wavenumbers in units of h Mpc^-1. On the other hand, if you want to use Mpc^-1 externally, you can do this by altering the variables you pass to the driver: omega0 -> omega0*hubble*hubble, hubble -> 1.0 */ /* INPUT: omega0 -- the matter density (baryons+CDM) in units of critical f_baryon -- the ratio of baryon density to matter density hubble -- the Hubble constant, in units of 100 km/s/Mpc Tcmb -- the CMB temperature in Kelvin. T<=0 uses the COBE value 2.728. numk -- the length of the following zero-offset array k[] -- the array of wavevectors k[0..numk-1] */ /* INPUT/OUTPUT: There are three output arrays of transfer functions. All are zero-offset and, if used, must have storage [0..numk-1] declared in the calling program. However, if you substitute the NULL pointer for one or more of the arrays, then that particular transfer function won't be outputted. The transfer functions are: tf_full[] -- The full fitting formula, eq. (16), for the matter transfer function. tf_baryon[] -- The baryonic piece of the full fitting formula, eq. 21. tf_cdm[] -- The CDM piece of the full fitting formula, eq. 17. */ /* Again, you can set these pointers to NULL in the function call if you don't want a particular output. */ /* Various intermediate scalar quantities are stored in global variables, so that you might more easily access them. However, this also means that you would be better off not simply #include'ing this file in your programs, but rather compiling it separately, calling only the driver, and using extern declarations to access the intermediate quantities. */ /* ----------------------------- DRIVER ------------------------------- */ void TFfit_hmpc(double omega0, double f_baryon, double hubble, double Tcmb, int numk, double *k, double *tf_full, double *tf_baryon, double *tf_cdm) /* Remember: k[0..numk-1] is in units of h Mpc^-1. */ { int j; double tf_thisk, baryon_piece, cdm_piece; TFset_parameters(omega0*hubble*hubble, f_baryon, Tcmb); for (j=0;j 1 and omega0 -> omega0*hubble^2. */ { double omhh, sound_horizon_fit_mpc; omhh = omega0*hubble*hubble; sound_horizon_fit_mpc = 44.5*log(9.83/omhh)/sqrt(1+10.0*pow(omhh*f_baryon,0.75)); return sound_horizon_fit_mpc*hubble; } double TFk_peak(double omega0, double f_baryon, double hubble) /* Input: omega0 -- CDM density, in units of critical density f_baryon -- Baryon fraction, the ratio of baryon to CDM density. hubble -- Hubble constant, in units of 100 km/s/Mpc */ /* Output: The approximate location of the first baryonic peak, in h Mpc^-1 */ /* Note: If you prefer to have the answer in units of Mpc^-1, use hubble -> 1 and omega0 -> omega0*hubble^2. */ { double omhh, k_peak_mpc; omhh = omega0*hubble*hubble; k_peak_mpc = 2.5*3.14159*(1+0.217*omhh)/TFsound_horizon_fit(omhh,f_baryon,1.0); return k_peak_mpc/hubble; } double TFnowiggles(double omega0, double f_baryon, double hubble, double Tcmb, double k_hmpc) /* Input: omega0 -- CDM density, in units of critical density f_baryon -- Baryon fraction, the ratio of baryon to CDM density. hubble -- Hubble constant, in units of 100 km/s/Mpc Tcmb -- Temperature of the CMB in Kelvin; Tcmb<=0 forces use of COBE FIRAS value of 2.728 K k_hmpc -- Wavenumber in units of (h Mpc^-1). */ /* Output: The value of an approximate transfer function that captures the non-oscillatory part of a partial baryon transfer function. In other words, the baryon oscillations are left out, but the suppression of power below the sound horizon is included. See equations (30) and (31). */ /* Note: If you prefer to use wavenumbers in units of Mpc^-1, use hubble -> 1 and omega0 -> omega0*hubble^2. */ { double k, omhh, theta_cmb, k_equality, q, xx, alpha_gamma, gamma_eff; double q_eff, T_nowiggles_L0, T_nowiggles_C0; k = k_hmpc*hubble; /* Convert to Mpc^-1 */ omhh = omega0*hubble*hubble; if (Tcmb<=0.0) Tcmb=2.728; /* COBE FIRAS */ theta_cmb = Tcmb/2.7; k_equality = 0.0746*omhh/SQR(theta_cmb); q = k/13.41/k_equality; xx = k*TFsound_horizon_fit(omhh, f_baryon, 1.0); alpha_gamma = 1-0.328*log(431.0*omhh)*f_baryon + 0.38*log(22.3*omhh)* SQR(f_baryon); gamma_eff = omhh*(alpha_gamma+(1-alpha_gamma)/(1+POW4(0.43*xx))); q_eff = q*omhh/gamma_eff; T_nowiggles_L0 = log(2.0*2.718282+1.8*q_eff); T_nowiggles_C0 = 14.2 + 731.0/(1+62.5*q_eff); return T_nowiggles_L0/(T_nowiggles_L0+T_nowiggles_C0*SQR(q_eff)); } /* ======================= Zero Baryon Formula =========================== */ double TFzerobaryon(double omega0, double hubble, double Tcmb, double k_hmpc) /* Input: omega0 -- CDM density, in units of critical density hubble -- Hubble constant, in units of 100 km/s/Mpc Tcmb -- Temperature of the CMB in Kelvin; Tcmb<=0 forces use of COBE FIRAS value of 2.728 K k_hmpc -- Wavenumber in units of (h Mpc^-1). */ /* Output: The value of the transfer function for a zero-baryon universe. */ /* Note: If you prefer to use wavenumbers in units of Mpc^-1, use hubble -> 1 and omega0 -> omega0*hubble^2. */ { double k, omhh, theta_cmb, k_equality, q, T_0_L0, T_0_C0; k = k_hmpc*hubble; /* Convert to Mpc^-1 */ omhh = omega0*hubble*hubble; if (Tcmb<=0.0) Tcmb=2.728; /* COBE FIRAS */ theta_cmb = Tcmb/2.7; k_equality = 0.0746*omhh/SQR(theta_cmb); q = k/13.41/k_equality; T_0_L0 = log(2.0*2.718282+1.8*q); T_0_C0 = 14.2 + 731.0/(1+62.5*q); return T_0_L0/(T_0_L0+T_0_C0*q*q); } /* Added by cmv 24/08/2007 */ void TFprint_parameters(void) { printf("omhh=%g obhh=%g theta_cmb=%g\n",omhh,obhh,theta_cmb); printf("z_equality=%g k_equality=%g R_equality=%g\n",z_equality,k_equality,R_equality); printf("z_drag=%g R_drag=%g sound_horizon=%g\n",z_drag,R_drag,sound_horizon); printf("k_silk=%g\n",k_silk); printf("alpha_c=%g beta_c=%g\n",alpha_c,beta_c); printf("alpha_b=%g beta_b=%g\n",alpha_b,beta_b); printf("beta_node=%g\n",beta_node); printf("k_peak=%g sound_horizon_fit=%g\n",k_peak,sound_horizon_fit); printf("alpha_gamma=%g\n",alpha_gamma); }