/* ------ File: -------------------- evalH3.c ------------ --- */ /* The equations for recombination are computed here. * Specifically, the ode for H recombination, He recombination, * and matter temperature. * * For use in a stiff integrator, the derivatives * with respect to redshift, and the Jacobian are computed. * * This code is in MKS, but the output is dimensionless. * * Sara Seager * Last modified: May 7 1999. */ #include #include "atom.h" #include "constant.h" #include "integrate.h" #include "nrutil.h" #include "input.h" #define MAX 1.0e+300 #define MIN 1.0e-300 /* --- global variables --------------------------------------------------- */ extern struct Input input; /* ------ begin -------------------- evaluate.c -------------------------- */ void evaluate(int choice, double z, double y[], double fvec[], double dfdz[], double **fjac, int n) { int i, j; double xe, xp, xHep, nH1, nHe1; double TR, TM, Comp; double Const; double Hz, Zeq, DD, dDD_dz, C, e, Fnu; double dH_dz; double PHI2, B_12, B_21; double nHTot, DnHTot; double KH, DKH_dz, KHe, dKHe_dz,fHe; double alphaH, dalphaH_dTM, betaH, dbetaH_dTM; double alphaHe, betaHe, dalphaHe_dTM, dbetaHe_dTM; double BH, BHe, dBH_dTM, dBHe_dTM, BHe2p2s, dBHe2p2s_dTM; double SH, SHe, dSH_dTM, dSHe_dTM; double CH, CHe, dCH_dz, dCHe_dz, dCH_dnH1, dCHe_dnHe1; double dCH_dTM, dCHe_dTM; double a1, a2, a3, a4; double T4, T0, T1; /* --- Definition of variables: * Cosmological parameters: * z: redshift * TM: matter temperature * TR: radiation temperature * TNOW: present day CMB temperature (input) * dt_dz: derivative of time wrt dz * Hz: Hubble parameter * dH_dHz: the derivative wrt z * nHTot: total hydrogen nuclei * dnHTot_dz: derivative wrt z * fHe: ratio of He/H by number * * Atomic parameters (H or He appended), defined in atom.h: * Lam2s1s: two photon process * Eion2s: ionization energy in J of n=2s * E2s1s: energy of n=2s from the ground state * E2p2sHe: energy difference between 2p and 2s * Parameters in the matter temperature equation: * sigmaT: Thomson scattering cross section * Comp: part of Compton cooling equation * * Parameters in the integrands: * xe = ne/nHTot = fraction of electrons * xp = np/nHTot = fraction of protons * nH1 = ground state of H atom * xHep = He+ = nHe+/nHTot * We use fractional number densities nx/nHTot, to avoid * evolving volme. Choice of normalization is arbitrary * (nx/nHTot) is arbitrary, could use nx/nHetot, etc. * alphaH: hydrogen case B recombination coefficient * alphaHe: helium case B recombination coefficient * CH: H inhibition factor (See Peebles, 1968, 1993) * CHe: He inhibition factor * KH: H cosmological redshifting factor (See Peebles, 1968, 1993) * KHe: He cosmological redshifting factor * BH the boltzmann equilibrium relation between n=2s and n=1s it comes out of the ode derivation from getting rid of n=2s * SH the Saha equilibrium relation between n=2s and the continuum, used for deriving beta from alpha * a1, a2, a3, a4, T0, T1: fitting coefficients for rec coefficients * * Constants: * see include file constant.h */ /* -- Calculate H(z) ----- ---------------------------------------- */ H_z(z, &Hz, &dH_dz); nHTot = NHnow(z); DnHTot = nHTot*3.0/(1.0+z); fHe = YP/(4.0*(1.0-YP)); /* --- Set variables --------------------------------------------- */ xHep = y[1]; xp = y[2]; xe = xHep + xp; TM = y[3]; nH1 = (1.0-xp)*nHTot; nHe1 = nHTot*(fHe-xHep); /* --- Radiation temperature -------------------------------- */ TR = To*(1.0+z); /* --- Calculate matter temperature equation parameters ---------------- */ /* Comp = 8.0*sigmaT*aRad*pow(TR,4.0)/(3.0*Hz*(1.0+z)*mElect*cLight); */ Const = 1.473883e-21; Comp = Const*pow(TR,4.0)/(3.0*Hz*(1.0+z)); /* --- Calculate the Saha and Boltzmann equilibrium relations--------- --- */ SHe = SahaBoltz(1.0,2.0,1.0, EionHe2s, TM); dSHe_dTM = SHe*(-1.0/TM*(1.5 + (EionHe2s)/(kBoltz*TM))); BHe = Boltzmann(1.0, 1.0, E2s1sHe, TM); dBHe_dTM = BHe * (E2s1sHe)/(kBoltz*TM*TM); SH = SahaBoltz(2.0,1.0,1.0, EionH2s, TM); dSH_dTM = SH*(-1.0/TM*(1.5 + (EionH2s)/(kBoltz*TM))); BH = Boltzmann(1.0, 1.0, E2s1sH, TM); dBH_dTM = BH * E2s1sH/(kBoltz*TM*TM); /* --- For HeI, the energy levels of 2s and 2p are quite different. Therefore the ratio b should be a boltzmann factor, and this will also have to be incorporated into the derivatives wrt TM ---- */ BHe2p2s = Boltzmann(3.0, 1.0, E2p2sHe, TM); dBHe2p2s_dTM = BHe2p2s * E2p2sHe/(kBoltz*TM*TM); /* --- The Saha and Boltzmann relations must be set to MIN and MAX accordingly, due to numerical problems at low z. --- */ if (BHe2p2s < 1.0e-300) BHe2p2s = MIN; if (dBHe2p2s_dTM < 1.0e-300) dBHe2p2s_dTM = MIN; if (SHe > 1.0e+300) SHe = MAX; if (fabs(dSHe_dTM) > 1.0e+300) dSHe_dTM = MAX; if (SH > 1.0e+300) SH = MAX; if (fabs(dSH_dTM) > 1.0e+300) dSH_dTM = MAX; if (BH < 1.0e-300) BH = MIN; if (dBH_dTM < 1.0e-300) dBH_dTM = MIN; if (BHe < 1.0e-300) BHe = MIN; if (dBHe_dTM < 1.0e-300) dBHe_dTM = MIN; /* --- H recombination coefficient alphaBH and photoionization coefficient betaH. --- */ a1 =4.309; a2 = -0.6166; a3 = 0.6703; a4 = 0.5300; T4 = TM/1.0e+04; alphaH = 1.0e-13*a1*pow(T4,a2)/(1.0 + a3*pow(T4,a4))*1.0e-06; /* --- Multiply alpha by the fudge factor input.F. This step is what makes this simple ode method approximate the full multi-level calculation for H. Note that because beta is derived from alpha, it is also affected by F. --- */ alphaH *= input.F; dalphaH_dTM = alphaH*(a2/TM - a3*a4*pow(T4,a4)/(TM*(1.0+ a3*pow(T4,a4)))); betaH = alphaH/SH; dbetaH_dTM = dalphaH_dTM/SH - alphaH/SQ(SH)*dSH_dTM; /* --- He case B recombination and photoionization coefficient No need to multiply a fudge factor. --- */ a1 = 1.691e-12; a2 = 1.519; T0=3.0; T1=3.2026e+04; alphaHe = (sqrt(TM/T0)*pow((1.0+sqrt(TM/T0)),1.0-a2) *pow(1.0+sqrt(TM/T1),1.0+a2)); alphaHe = a1/alphaHe/1.0e+06; dalphaHe_dTM = (0.5*sqrt(T0/TM)*pow((1.0+sqrt(TM/T0)),1.0-a2) + sqrt(TM/T0)*(1.0-a2)*pow((1.0+sqrt(TM/T0)),-a2)*0.5*sqrt(T0/TM)) *pow(1.0+sqrt(TM/T1),1.0+a2) + sqrt(TM/T0)*pow((1.0+sqrt(TM/T0)),1.0-a2)*(1.0+a2)* pow(1.0+sqrt(TM/T1),a2)*0.5*sqrt(T1/TM); dalphaHe_dTM = -1.0e+06/a1*SQ(alphaHe)*dalphaHe_dTM; betaHe = alphaHe/SHe; dbetaHe_dTM = (dalphaHe_dTM/SHe - alphaHe/SQ(SHe)*dSHe_dTM); /* --- Cosmological redshifting --- */ KHe = CUBE(LyalphaHe)/(Hz*8.0*PI); dKHe_dz = -KHe*dH_dz/Hz; KH = CUBE(LyalphaH)/(Hz*8.0*PI); DKH_dz = -KH*dH_dz/Hz; /* --- Inhibition factors --- */ CHe = (1.0+KHe*Lam2s1sHe*nHe1*BHe2p2s)/ (1.0+KHe*(Lam2s1sHe+betaHe)*nHe1*BHe2p2s); dCHe_dz = CHe *((Lam2s1sHe*nHe1*(dKHe_dz*BHe2p2s))/ (1.0+KHe*Lam2s1sHe*nHe1*BHe2p2s) - (nHe1*((Lam2s1sHe+betaHe)*(dKHe_dz*BHe2p2s)))/ (1.0+KHe*(Lam2s1sHe+betaHe)*nHe1*BHe2p2s)); dCHe_dTM = CHe *((Lam2s1sHe*nHe1*(KHe*dBHe2p2s_dTM))/(1.0+KHe*Lam2s1sHe*nHe1*BHe2p2s) - (nHe1*((Lam2s1sHe+betaHe)*(KHe*dBHe2p2s_dTM) +KHe*BHe2p2s*betaHe)) /(1.0+KHe*(Lam2s1sHe+betaHe)*nHe1*BHe2p2s)); dCHe_dnHe1 = CHe*((-KHe*Lam2s1sHe*nHTot*BHe2p2s)/ (1.0+KHe*Lam2s1sHe*nHe1*BHe2p2s) + KHe*(Lam2s1sHe+betaHe)*nHTot*BHe2p2s /(1.0+KHe*(Lam2s1sHe+betaHe)*nHe1*BHe2p2s)); CH = (1.0+KH*Lam2s1sH*nH1)/(1.0+KH*(Lam2s1sH+betaH)*nH1); dCH_dz = CH *((DKH_dz*Lam2s1sH*nH1)/(1.0+KH*Lam2s1sH*nH1) - (DKH_dz*(Lam2s1sH+betaH)*nH1)/(1.0+KH*(Lam2s1sH+betaH)*nH1)); dCH_dTM = -CH/(1.0+KH*(Lam2s1sH+betaH)*nH1)*KH*nH1*dbetaH_dTM; dCH_dnH1 = CH*((-KH*Lam2s1sH*nHTot)/(1.0+KH*Lam2s1sH*nH1) + KH*(Lam2s1sH+betaH)*nHTot/(1.0+KH*(Lam2s1sH+betaH)*nH1)); /* --- Finally calculate the rate equations. 1 = He, 2 = H, 3 = TM -------------- */ fvec[1] = (-alphaHe*xe*xHep*nHTot + betaHe*(fHe-xHep)*BHe)*CHe; fvec[2] = (-alphaH*xe*xp*nHTot + betaH*(1.0-xp)*BH)*CH; fvec[3] = Comp*xe/(1.0+xe+fHe)*(TM - TR) + 2.0*TM/(1.0+z); /* --- Calculate derivatives wrt redshift ---------------------- */ if (choice ==1) { dfdz[1] = -xe*xHep*(alphaHe*DnHTot)*CHe + fvec[1]*dCHe_dz/CHe; dfdz[2] = (-xe*xp*alphaH*DnHTot)*CH + fvec[2]*dCH_dz/CH; dfdz[3] = Comp*xe/(1.0+xe+fHe)*((3.0/(1.0+z) - dH_dz/Hz)* (TM-TR) -To) - 2.0*TM/SQ(1.0+z); /* --- Calculate the jacobian ------------------------------------------ */ fjac[1][1] = (-(alphaHe*(xp + xe)*nHTot) + betaHe*BHe*(-1.0))*CHe + fvec[1]*dCHe_dnHe1/CHe; fjac[1][2] = -alphaHe*xHep*nHTot*CHe; fjac[1][3] = ((-xe*xHep*(dalphaHe_dTM*nHTot) + (fHe-xHep)*(dbetaHe_dTM*BHe + betaHe*dBHe_dTM)))*CHe + fvec[1]*dCHe_dTM/CHe; fjac[2][1] = -alphaH*xp*nHTot*CH; fjac[2][2] = (-(alphaH*(xe+xp)*nHTot) + betaH*BH*(-1.0))*CH + fvec[2]*dCH_dnH1/CH; fjac[2][3] = (-xe*xp*dalphaH_dTM*nHTot + (1.0-xp)*(dbetaH_dTM*BH + betaH*dBH_dTM))*CH + fvec[2]*dCH_dTM/CH; fjac[3][1] = Comp*xe/(1.0+xe+fHe)*(-1.0/(1.0+xe+fHe) + 1.0/xe)*(TM - TR); fjac[3][2] = Comp*xe/(1.0+xe+fHe)*(-1.0/(1.0+xe+fHe) + 1.0/xe)*(TM - TR); fjac[3][3] = Comp*xe/(1.0+xe) + 2.0/(1.0+z); } /* --- Convert from dn/dt to dn/dz by multiplying all population equations by DD=dt/dz --- */ DD = -1.0/((1.0+z)*Hz); C = 0.5/(1.0+z+1.0+Zeq) + 1.5/(1.0+z); dDD_dz = -DD*(1.0/(1.0+z) + C); for (i=1; i<=2; i++) { if (choice == 1) { dfdz[i] = dfdz[i]*DD + fvec[i]*dDD_dz; for (j=1; j<=3; j++) fjac[i][j] *= DD; } fvec[i] *= DD; } } #undef MAX #undef MIN