*** A Routine for Computing Various Cosmological Quantities *** March 17, 2012: E. Komatsu This routine, "secondary.f90", returns various "secondary parameters/quantities/functions" such as the redshift of equality, the comoving angular diameter distance, etc. This routine was used for deriving the secondary parameters from the WMAP 5-year and 7-year data. *************** IMPORTANT NOTES *************** 1. The equality epoch is computed from a_EQ=Omega_R/Omega_M, where Omega_R contains photons and relativistic neutrinos. This routine computes a_EQ assuming that all of Neff neutrino species were relativistic at the equality epoch. This should be an excellent approximation given the current upper bound on the neutrino masses. (As the equality redshift is zeq~3150, neutrinos were relativistic at zeq if the INDIVIDUAL neutrino mass is much less than 1.8eV.) 2. The photon decoupling redshift and the BAO drag redshift are calculated using the fitting formulae developed by Hu&Sugiyama and Eisenstein&Hu. These formulae were derived for Neff=3.04, and thus they may not be very accurate when Neff is not equal to 3.04. 3. Massless-massive neutrino transition is included in this routine, following Section 3.3 of the WMAP 7-year paper [Komatsu et al., ApJS, 192, 18 (2011)]. ALL of Neff neutrinos (not just three neutrinos) are assumed to have the same mass, and thus ALL of Neff neutrinos become non-relativistic at the same redshift. Be careful when varying both Neff and omega_nu*h^2 simultaneously. For example, you may wish to assume that only 3 neutrino species are massive and Neff-3.04 species are completely massless at all redshifts. For such a case you need to modify the code ("DOUBLE PRECISION FUNCTION H2(x)", which follows Eq.27 of Komatsu et al.) accordingly. --------------------- << Core Parameters >> --------------------- In order to use this routine, you must specify the following "core parameters": - om0 : Omega_m at present [This does NOT include massive neutrinos!] - ob0 : Omega_b at present - ode0 : Omega_de at present - onuh2: Omega_massiveneutrino * h^2 at present - h : H0/100(km/s/Mpc) - w : equation of state parameter - wa : derivative of equation of state parameter w.r.t 1-a - Neff : effective number of massless neutrinos These are double-precision variables declared in "MODULE cosmo_parameters". Therefore, you need to write "USE cosmo_parameters" at the beginning of your code. Once you specify the core parameters, you need to call the routine by "CALL get_rest_of_parameters". Then the code will calculate the secondary parameters as well as various cosmological functions. ------------------------------------- << Secondary Parameters/Quantities >> ------------------------------------- The secondary parameters/quantities are double-precision variables declared in "MODULE cosmo_parameters". They include: - aeq : scale factor at equality, assuming that all Neff neutrinos were still relativistic at the equality epoch - eta10 : baryon-to-photon ratio times 10^{10} - dt1608: D_t for B1608+656 in Mpc [Suyu et al., ApJ, 711,201 (2010)] - zdEH : Eisenstein&Hu's fitting formula for the drag redshift - zsHS : Hu&Sugiyama's fitting formula for the decoupling redshift ** CAUTION: these formulae for zdEH and zsHS were derived assuming Neff=3.04 ** - H0 : Hubble's constant in units of 1/Mpc - ok0 : 1-Omega_m-Omega_de at present [radiation contribution is ignored] - omh2 : Omega_m*h^2 at present [This does NOT include massive neutrinos!] - obh2 : Omega_b*h^2 at present - okh2 : (1-Omega_m-Omega_de)*h^2 at present - odeh2 : Omega_de*h^2 at present - og0 : Omega_photon at present - ogh2 : Omega_photon*h^2 at present - or0 : Omega_radiation = Omega_photon + Omega_relativisticneutrino at present, assuming that all Neff neutrinos are still relativistic - orh2 : Omega_radiation*h^2, assuming that all Neff neutrinos are still relativistic ------------------------- << Secondary Functions >> ------------------------- The secondary functions are double-precision functions, which are NOT declared in "MODULE cosmo_parameters"; thus, you need to declare these functions to be double-precision before you use them. - da(a) : comoving angular diameter distance in Mpc - da2(a1,a2): comoving angular diameter distance between a1 & a2 in Mpc - rs(a) : comoving sound horizon in Mpc - theta(a): rs(a)/da(a) in radians - la(a) : acoustic scale, pi/theta - age(a) : age in years; age(1d0) is the present age - dv(a) : D_v in Mpc [Eisenstein et al., ApJ, 633, 650 (2005)] - chi(a) : comoving distance in Mpc; chi(0d0) is the particle horizon - H2(a) : H^2(a) in 1/Mpc^2 - weff(a) : effective equation of state parameter, given by weff = (-1/3)*dln(rho_de)/dln(a) - 1 * a is the scale factor. ----------- << Usage >> ----------- (1) Put "USE cosmo_parameters" at the beginning of the program (2) Declare the secondary functions that you would like to use to be double-precision (3) Specify the core cosmological parameters (4) CALL get_rest_of_parameters Or, you may simply cut-and-paste the following: !=============================================================== ! In order to use "secondary.f90", cut-and-paster from here... USE cosmo_parameters IMPLICIT none ! Available functions: DOUBLE PRECISION :: da ! d_A(a), comoving angular diameter distance in Mpc DOUBLE PRECISION :: da2 ! d_A between a1 & a2 in Mpc DOUBLE PRECISION :: rs ! r_s(a), sound horizon in Mpc DOUBLE PRECISION :: theta ! r_s(a)/d_A(a), acoustic scale in radians DOUBLE PRECISION :: la ! l_A(a), acoustic wavenumber, pi/theta DOUBLE PRECISION :: age ! t(a), cosmic age in years DOUBLE PRECISION :: dv ! D_V(a), distance measured by BAO DOUBLE PRECISION :: chi ! chi(a), comoving distance in Mpc DOUBLE PRECISION :: H2 ! H^2(a), hubble function squared in 1/Mpc^2 DOUBLE PRECISION :: weff ! weff(a) = (-1/3)*dln(rhoDE)/dln(a) - 1 ! Other parameters (not functions, declared double-precision already): ! - aeq : scale factor at equality ! - eta10 : baryon-to-photon ratio times 10^{10} ! - dt1608 : D_t for B1608+656 in Mpc [Suyu et al., ApJ, 711,201 (2010)] ! - zdEH : Eisenstein&Hu's fitting formula for the drag redshift ! - zsHS : Hu&Sugiyama's fitting formula for the decoupling redshift ! * CAUTION: zdEH and zsHS were derived assuming Neff=3.04. ! ! "Core parameters" that you need to specify. ! The following examples are the WMAP7+BAO+H0 Maximum Likelihood parameters. om0 = 0.272d0 ! Omega_m at present [This does NOT include massive neutrinos!] ob0 = 0.0455d0! Omega_b at present ode0 = 0.728d0 ! Omega_de at present onuh2= 0d0 ! Omega_massiveneutrino * h^2 at present h = 0.704d0 ! H0/100(km/s/Mpc) w = -1d0 ! equation of state parameter wa = 0d0 ! derivative of equation of state parameter w.r.t 1-a Neff = 3.04d0 ! effective number of massless neutrinos ! *** IMPORTANT: CALL get_rest_of_parameters TO GET THE DERIVED PARAMETERS *** call get_rest_of_parameters ! ...to here. !=============================================================== --------------- << Test Code >> --------------- For your convenience, a test routine, "test.f90", is included in this package. To run this code, edit Makefile and make. It will generate an executable called "test", which you can run by "./test".