recfast.f90 File Source

Recombination module for CAMB, using RECFAST


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C Integrator for Cosmic Recombination of Hydrogen and Helium,
C developed by Douglas Scott (dscott@astro.ubc.ca)
C based on calculations in the paper Seager, Sasselov & Scott
C (ApJ, 523, L1, 1999).
and "fudge" updates in Wong, Moss & Scott (2008).
C
C Permission to use, copy, modify and distribute without fee or royalty at
C any tier, this software and its documentation, for any purpose and without
C fee or royalty is hereby granted, provided that you agree to comply with
C the following copyright notice and statements, including the disclaimer,
C and that the same appear on ALL copies of the software and documentation,
C including modifications that you make for internal use or for distribution:
C
C Copyright 1999-2010 by University of British Columbia.  All rights reserved.
C
C THIS SOFTWARE IS PROVIDED "AS IS", AND U.B.C. MAKES NO
C REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.
C BY WAY OF EXAMPLE, BUT NOT LIMITATION,
c U.B.C. MAKES NO REPRESENTATIONS OR WARRANTIES OF
C MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT
C THE USE OF THE LICENSED SOFTWARE OR DOCUMENTATION WILL NOT INFRINGE
C ANY THIRD PARTY PATENTS, COPYRIGHTS, TRADEMARKS OR OTHER RIGHTS.
C
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

CN     Name:        RECFAST
CV     Version: 1.5.2
C
CP     Purpose:  Calculate ionised fraction as a function of redshift.
CP            Solves for H and He simultaneously, and includes
CP           H "fudge factor" for low z effect, as well as
CP           HeI fudge factor.
C
CD     Description: Solves for ionisation history since recombination
CD     using the equations in Seager, Sasselov & Scott (ApJ, 1999).
CD     The Cosmological model can be flat or open.
CD  The matter temperature is also followed, with an update from
CD  Scott & Scott (2009).
CD  The values for \alpha_B for H are from Hummer (1994).
CD  The singlet HeI coefficient is a fit from the full code.
CD  Additional He "fudge factors" are as described in Wong, Moss
CD  and Scott (2008).
CD  Extra fitting function included (in optical depth) to account
CD  for extra H physics described in Rubino-Martin et al. (2010).
CD  Care is taken to use the most accurate constants.
C
CA     Arguments:
CA     Name, Description
CA     real(dl) throughout
CA
CA     z is redshift - W is sqrt(1+z), like conformal time
CA     x is total ionised fraction, relative to H
CA     x_H is ionized fraction of H - y(1) in R-K routine
CA     x_He is ionized fraction of He - y(2) in R-K routine
CA       (note that x_He=n_He+/n_He here and not n_He+/n_H)
CA     Tmat is matter temperature - y(3) in R-K routine
CA     f's are the derivatives of the Y's
CA     alphaB is case B recombination rate
CA     alpHe is the singlet only HeII recombination rate
CA     a_PPB is Pequignot, Petitjean & Boisson fitting parameter for Hydrogen
CA     b_PPB is Pequignot, Petitjean & Boisson fitting parameter for Hydrogen
CA     c_PPB is Pequignot, Petitjean & Boisson fitting parameter for Hydrogen
CA     d_PPB is Pequignot, Petitjean & Boisson fitting parameter for Hydrogen
CA     a_VF is Verner and Ferland type fitting parameter for Helium
CA     b_VF is Verner and Ferland type fitting parameter for Helium
CA     T_0 is Verner and Ferland type fitting parameter for Helium
CA     T_1 is Verner and Ferland type fitting parameter for Helium
CA     Tnow is the observed CMB temperature today
CA     Yp is the primordial helium abundace
CA     fHe is He/H number ratio = Yp/4(1-Yp)
CA     Trad and Tmat are radiation and matter temperatures
CA     epsilon is the approximate difference (=Trad-Tmat) at high z
CA     OmegaB is Omega in baryons today
CA     H is Hubble constant in units of 100 km/s/Mpc
CA     HO is Hubble constant in SI units
CA     bigH is 100 km/s/Mpc in SI units
CA     Hz is the value of H at the specific z (in ION)
CA     G is grvitational constant
CA     n is number density of hydrogen
CA     Nnow is number density today
CA     x0 is initial ionized fraction
CA     x_H0 is initial ionized fraction of Hydrogen
CA     x_He0 is initial ionized fraction of Helium
CA     rhs is dummy for calculating x0
CA     zinitial and zfinal are starting and ending redshifts
CA     zeq is the redshift of matter-radiation equality
CA     zstart and zend are for each pass to the integrator
CA     C,k_B,h_P: speed of light, Boltzmann's and Planck's constants
CA     m_e,m_H: electron mass and mass of H atom in SI
CA     not4: ratio of 4He atomic mass to 1H atomic mass
CA     sigma: Thomson cross-section
CA     a_rad: radiation constant for u=aT^4
CA     Lambda: 2s-1s two photon rate for Hydrogen
CA     Lambda_He: 2s-1s two photon rate for Helium
CA     DeltaB: energy of first excited state from continuum = 3.4eV
CA     DeltaB_He: energy of first excited state from cont. for He = 3.4eV
CA     L_H_ion: level for H ionization in m^-1
CA     L_H_alpha: level for H Ly alpha in m^-1
CA     L_He1_ion: level for HeI ionization
CA     L_He2_ion: level for HeII ionization
CA     L_He_2s: level for HeI 2s
CA     L_He_2p: level for HeI 2p (21P1-11S0) in m^-1
CA     Lalpha: Ly alpha wavelength in SI
CA     Lalpha_He: Helium I 2p-1s wavelength in SI
CA     mu_H,mu_T: mass per H atom and mass per particle
CA     H_frac: follow Tmat when t_Compton / t_Hubble > H_frac
CA     CDB=DeltaB/k_B                     Constants derived from B1,B2,R
CA     CDB_He=DeltaB_He/k_B  n=2-infinity for He in Kelvin
CA     CB1=CDB*4.         Lalpha and sigma_Th, calculated
CA     CB1_He1: CB1 for HeI ionization potential
CA     CB1_He2: CB1 for HeII ionization potential
CA     CR=2*Pi*(m_e/h_P)*(k_B/h_P)  once and passed in a common block
CA     CK=Lalpha**3/(8.*Pi)
CA     CK_He=Lalpha_He**3/(8.*Pi)
CA     CL=C*h_P/(k_B*Lalpha)
CA     CL_He=C*h_P/(k_B*Lalpha_He)
CA     CT=(8./3.)*(sigma/(m_e*C))*a
CA     Bfact=exp((E_2p-E_2s)/kT)    Extra Boltzmann factor
CA b_He= "fudge factor" for HeI, to approximate higher z behaviour
CA Heswitch=integer for modifying HeI recombination
CA Parameters and quantities to describe the extra triplet states
CA  and also the continuum opacity of H, with a fitting function
CA  suggested by KIV, astro-ph/0703438
CA a_trip: used to fit HeI triplet recombination rate
CA b_trip: used to fit HeI triplet recombination rate
CA L_He_2Pt: level for 23P012-11S0 in m^-1
CA L_He_2St: level for 23S1-11S0 in m^-1
CA L_He2St_ion: level for 23S1-continuum in m^-1
CA A2P_s: Einstein A coefficient for He 21P1-11S0
CA A2P_t: Einstein A coefficient for He 23P1-11S0
CA sigma_He_2Ps: H ionization x-section at HeI 21P1-11S0 freq. in m^2
CA sigma_He_2Pt: H ionization x-section at HeI 23P1-11S0 freq. in m^2
CA CL_PSt = h_P*C*(L_He_2Pt - L_He_2st)/k_B
CA CfHe_t: triplet statistical correction
CA Hswitch is an boolean for modifying the H recombination
CA AGauss1 is the amplitude of the 1st Gaussian for the H fudging
CA AGauss2 is the amplitude of the 2nd Gaussian for the H fudging
CA zGauss1 is the ln(1+z) central value of the 1st Gaussian
CA zGauss2 is the ln(1+z) central value of the 2nd Gaussian
CA wGauss1 is the width of the 1st Gaussian
CA wGauss2 is the width of the 2nd Gaussian




CA     tol: tolerance for the integrator
CA     cw(24),w(3,9): work space for DVERK
CA     Ndim: number of d.e.'s to solve (integer)
CA     Nz: number of output redshitf (integer)
CA     I: loop index (integer)
CA     ind,nw: work-space for DVERK (integer)
C
CF     File & device access:
CF     Unit /I,IO,O  /Name (if known)
C
CM     Modules called:
CM     DVERK (numerical integrator)
CM     GET_INIT (initial values for ionization fractions)
CM     ION (ionization and Temp derivatives)
C
CC     Comments:
CC     none
C
CH     History:
CH     CREATED            (simplest version) 19th March 1989
CH     RECREATED    11th January 1995
CH               includes variable Cosmology
CH               uses DVERK integrator
CH               initial conditions are Saha
CH     TESTED              a bunch, well, OK, not really
CH     MODIFIED     January 1995 (include Hummer's 1994 alpha table)
CH               January 1995 (include new value for 2s-1s rate)
CH               January 1995 (expand comments)
CH               March 1995 (add Saha for Helium)
CH               August 1997 (add HeII alpha table)
CH               July 1998 (include OmegaT correction and H fudge factor)
CH               Nov 1998 (change Trad to Tmat in Rup)
CH               Jan 1999 (tidied up for public consumption)
CH               Sept 1999 (switch to formula for alpha's, fix glitch)
CH                  Sept 1999 modified to CMBFAST by US & MZ
CH                     Nov 1999 modified for F90 and CAMB (AML)
CH                     Aug 2000 modified to prevent overflow erorr in He_Boltz (AML)
CH                     Feb 2001 corrected fix of Aug 2000 (AML)
CH                     Oct 2001 fixed error in hubble parameter, now uses global function (AML)
March 2003 fixed bugs reported by savita gahlaut
March 2005 added option for corrections from astro-ph/0501672.
thanks to V.K.Dubrovich, S.I.Grachev
June 2006 defined RECFAST_fudge as free parameter (AML)
October 2006 (included new value for G)
October 2006 (improved m_He/m_H to be "not4")
October 2006 (fixed error, x for x_H in part of f(1))
CH              January 2008 (improved HeI recombination effects,
CH                       including HeI rec. fudge factor)
Feb 2008   Recfast 1.4 changes above added (AML)
removed Dubrovich option (wrong anyway)
CH              Sept 2008 (added extra term to make transition, smoother for Tmat evolution)
Sept 2008 Recfast 1.4.2 changes above added (AML)
General recombination module structure, fix to make He x_e smooth also in recfast (AML)
CH      Jan 2010 (added fitting function to modify K
CH              to match x_e(z) for new H physics)
AL             June 2012 updated fudge parameters to match HyRec and CosmoRec (AML)
AL             Sept 2012 changes now in public recfast, version number changed to match Recfast 1.5.2.

Modules

Dependencies