! Exponential potential ! ! Equations module allowing for fairly general quintessence models ! ! by Antony Lewis (http://cosmologist.info/) ! This version Oct 2004. !!FIX March 2005: corrected update to next treatment of tight coupling ! Notes at http://antonylewis.com/notes/CAMB.ps !This module is not well tested, use at your own risk! !Need to specify Vofphi function, and also initial_phi in Quint_init_background !You may also need to change other things to get it to work with different types of quintessence model !It works backwards, in that it assumes Omega_lambda is Omega_Q today, then does a binary search on the !initial conditions to find what is required to give that Omega_Q today after evolution. !!!This is a dummy for the default constant parameterization module LambdaGeneral use precision real(dl) :: w_lam = -1, cs2_lam=1 !p/rho for the dark energy (assumed constant) end module LambdaGeneral module Quint use precision use ModelParams use LambdaGeneral implicit none integer, parameter :: NumPoints = 2000, NumPointsEx=NumPoints+2 !We actually calulate NumPointEx-NumPoints points after a=1 to get good interpolation at a=1 real(dl) phi_a(NumPointsEx),phidot_a(NumPointsEx) real(dl) ddphi_a(NumPointsEx),ddphidot_a(NumPointsEx) real(dl), parameter :: Mpc = 3.085678e22_dl, G=6.6726e-11_dl, kappa=2*fourpi*G, & c = 2.99792458e8_dl real(dl) adot logical, parameter :: has_extra = .true. !If true Om_Lambda is interpreted as quintessence, otherwise just Lambda real(dl) :: initial_phi, da real(dl), parameter :: m =100 !5d-7 real(dl) :: sigma = 1.2247 ! real(dl) :: sigma = 3 contains function Vofphi(phi,deriv) !Returns (8*Pi*G)^(1-deriv/2)*d^{deriv}V(psi)/d^{deriv}psi evaluated at psi !times (Mpc/c)^2 to get units in 1/Mpc^2 !The input variable phi is sqrt(8*Pi*G)*psi implicit none real(dl) phi,Vofphi integer deriv real(dl), parameter :: root8pi = 5.0132565492620010048_dl real(dl), parameter :: Tpl=5.391d-44*root8pi !Tpl = sqrt(8 pi G hbar/c^5) real(dl) lambda, norm lambda = sigma !dimensionless normalization kappa*V*T_pl^2 ! norm= 1 norm= 1d-122 norm = norm * (Mpc/c)**2 /Tpl**2 !convert to units of 1/Mpc^2 ! if (deriv==0) then ! !norm multiplies kappa*V Vofphi= norm*exp(-lambda*phi) ! else if (deriv ==1) then ! Vofphi=-norm*lambda*exp(-lambda*phi) ! else if (deriv ==2) then ! Vofphi=norm*lambda**2*exp(-lambda*phi) ! else ! stop 'Invalid deriv in Vofphi' ! end if ! return !Example linear potential ! sigma= sqrt(3./2.) if (deriv==0) then Vofphi= norm*m*exp(-sigma*phi) else if (deriv ==1) then Vofphi=-norm*m*sigma*exp(-sigma*phi) else if (deriv ==2) then Vofphi=norm*m*sigma**2*exp(-sigma*phi) else stop 'Invalid deriv in Vofphi' end if end function Vofphi subroutine EvolveBackground(dum,num,a,y,yprime) implicit none real dum integer num real(dl) tau, y(num),yprime(num) real(dl) a, a2, tot real(dl) phi,dphida, tmp, phidot a2=a**2 phi = y(1) phidot = y(2)/a2 tmp=a2*(0.5d0*phidot**2 + a2*Vofphi(phi,0)) tot=(grhoc+grhob)*a+grhog+grhornomass +grhok*a2 + tmp adot=sqrt(tot/3.0d0) dphida=phidot/adot yprime(1)=dphida yprime(2)= -a2**2*Vofphi(phi,1)/adot end subroutine EvolveBackground function Quint_GetOmegaFromInitial(astart,phi,phidot,atol) !Get CP%omegav today given particular conditions phi and phidot at a=astart real(dl), intent(IN) :: astart, phi,phidot, atol real(dl) Quint_GetOmegaFromInitial integer, parameter :: NumEqs=2 real(dl) c(24),w(NumEqs,9), y(NumEqs), ast integer ind, i real dum ast=astart ind=1 y(1)=phi y(2)=phidot*astart**2 !Fixed Dec 02 call dverk(dum,NumEqs,EvolveBackground,ast,y,1._dl,atol,ind,c,NumEqs,w) call EvolveBackground(dum,NumEqs,1._dl,y,w(:,1)) Quint_GetOmegaFromInitial=(0.5d0*y(2)**2 + Vofphi(y(1),0))/(3*adot**2) end function Quint_GetOmegaFromInitial function Quint_phidot_start(phi) real(dl) :: phi, Quint_phidot_start !!! Quint_phidot_start = 0 !sqrt(4*VofPhi(phi,0)) end function Quint_phidot_start subroutine Quint_init_background real(dl) astart,aend, afrom integer, parameter :: NumEqs=2 real(dl) c(24),w(NumEqs,9), y(NumEqs),atol integer ind, i, iter real(dl) aVals(NumPointsEx), splZero real(dl) om1,om2, initial_phi2, phi, deltaphi, om logical OK real dum real(dl) effint, omint, p, rhofrac real(dl) initial_phidot !Make interpolation table, etc, !At this point massive neutrinos have been initialized !so nu1 can be used to get their density and pressure at scale factor a !Other built-in components have density and pressure scaling trivially with a !These two must bracket the correct value to give CP%omegav today !Assume that higher initial phi gives higher CP%omegav today !Can fix initial_phi to correct value initial_phi = 7 ! 0.3*grhom/m**3 initial_phi2 =-4 ! 6*grhom/m**3 ! initial_phi = 65 ! 0.3*grhom/m**3 ! initial_phi2 = 65 ! 6*grhom/m**3 astart=1d-9 !See if initial conditions are giving correct CP%omegav now atol=1d-8 initial_phidot = astart*Quint_phidot_start(initial_phi) om1= Quint_GetOmegaFromInitial(astart,initial_phi,initial_phidot, atol) !!! if (.true.) then if (abs(om1-CP%omegav) > 1d-3) then !if not, do binary search in the interval OK=.false. initial_phidot = astart*Quint_phidot_start(initial_phi2) om2= Quint_GetOmegaFromInitial(astart,initial_phi2,initial_phidot, atol) if (om1 > CP%omegav .or. om2 < CP%omegav) then write (*,*) 'initial phi values must bracket required value. ' write (*,*) 'om1, om2 = ', real(om1), real(om2) stop end if do iter=1,100 deltaphi=initial_phi2-initial_phi phi =initial_phi + deltaphi/2 initial_phidot = astart*Quint_phidot_start(phi) om = Quint_GetOmegaFromInitial(astart,phi,initial_phidot,atol) if (om < CP%omegav) then om1=om initial_phi=phi else om2=om initial_phi2=phi end if if (om2-om1 < 1d-3) then OK=.true. initial_phi = (initial_phi2+initial_phi)/2 if (FeedbackLevel > 0) write(*,*) 'phi_initial = ',initial_phi exit end if end do !iterations if (.not. OK) stop 'Search for good intial conditions did not converge' !this shouldn't happen end if !Find initial end if atol=1d-8 ind=1 !y(1)=astart/adotrad y(1)=initial_phi initial_phidot = astart*Quint_phidot_start(initial_phi) y(2)= initial_phidot*astart**2 phi_a(1)=y(1) phidot_a(1)=y(2)/astart**2 afrom=astart da=1._dl/(NumPoints-1) !phi_a(NumPoints) is value now aVals(1)=astart !!Might need to be more careful about first point effint = 0 ; omint = 0 do i=1, NumPointsEx-1 aend = da*i aVals(i+1)=aend call dverk(dum,NumEqs,EvolveBackground,afrom,y,aend,atol,ind,c,NumEqs,w) call EvolveBackground(dum,NumEqs,aend,y,w(:,1)) phi_a(i+1)=y(1) phidot_a(i+1)=y(2)/aend**2 p = phidot_a(i+1)**2/2 - aend**2*Vofphi(phi_a(i+1),0) rhofrac = 3*(adot/aend)**2 effint = effint + p/rhofrac omint = omint + (phidot_a(i+1)**2 -p)/rhofrac if (i==NumPoints-1) then if (FeedbackLevel > 0) then write(*,*) 'Omega_Q_0=',real((0.5d0*phidot_a(i+1)**2 + Vofphi(phi_a(i+1),0))/(3*adot**2)), & ' w_0=',real((0.5d0*phidot_a(i+1)**2 - Vofphi(phi_a(i+1),0))/(0.5d0*phidot_a(i+1)**2 + & Vofphi(phi_a(i+1),0))) write (*,*) 'w_effective = ', real(effint/omint) end if end if end do splZero=0 call spline(aVals,phi_a,NumPointsEx,splZero,splZero,ddphi_a) call spline(aVals,phidot_a,NumPointsEx,splZero,splZero,ddphidot_a) end subroutine Quint_init_background subroutine Quint_ValsAta(a,aphi,aphidot) implicit none !Do interpolation for phi and phidot at a real(dl) a, aphi, aphidot real(dl) a0,b0,ho2o6,a03,b03 integer ix ix = int(a/da)+1 a0 = (ix*da -a)/da b0 = (a-(ix-1)*da)/da ho2o6 = da**2/6._dl a03=(a0**3-a0) b03=(b0**3-b0) aphi=a0*phi_a(ix)+b0*phi_a(ix+1)+(a03*ddphi_a(ix)+b03*ddphi_a(ix+1))*ho2o6 aphidot=a0*phidot_a(ix)+b0*phidot_a(ix+1)+(a03*ddphidot_a(ix)+b03*ddphidot_a(ix+1))*ho2o6 end subroutine Quint_ValsAta end module Quint !Return OmegaK - modify this if you add extra fluid components function GetOmegak() use precision use Quint use ModelParams real(dl) GetOmegak GetOmegak = 1 - (CP%omegab+CP%omegac+CP%omegav+CP%omegan) end function GetOmegak subroutine init_background use Quint !This is only called once per model, and is a good point to do any extra initialization. !It is called before first call to dtauda, but after !massive neutrinos are initialized and after GetOmegak if (has_extra) call Quint_init_background end subroutine init_background !Background evolution function dtauda(a) !get d tau / d a use precision use ModelParams use MassiveNu use Quint implicit none real(dl) dtauda integer nu_i real(dl), intent(IN) :: a real(dl) rhonu,pnu,grhoa2, a2, phi, phidot a2=a**2 ! 8*pi*G*rho*a**4. grhoa2=grhok*a2 +(grhoc+grhob)*a+grhog+grhornomass if (CP%Num_Nu_massive /= 0) then !Get massive neutrino density relative to massless do nu_i = 1, CP%nu_mass_eigenstates call Nu_rho(a*nu_masses(nu_i),rhonu) grhoa2=grhoa2+rhonu*grhormass(nu_i) end do end if !Get any extra contributions if (has_extra) then call Quint_ValsAta(a,phi,phidot) grhoa2=grhoa2 + a2*(phidot**2/2 + a2*Vofphi(phi,0)) else grhoa2=grhoa2 +grhov*a2**2 end if dtauda=sqrt(3/grhoa2) end function dtauda !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !Gauge-dependent perturbation equations module GaugeInterface use precision use ModelParams use MassiveNu use Quint use LambdaGeneral implicit none public !Description of this file. Change if you make modifications. character(LEN=*), parameter :: Eqns_name = 'gauge_quint' logical :: DoTensorNeutrinos = .false. logical :: DoLateRadTruncation = .false. real(dl) vec_sig0, Magnetic !(dummies) integer, parameter :: max_l_evolve = 50 !Maximum l we are ever likely to propagate !Supported scalar initial condition flags integer, parameter :: initial_adiabatic=1, initial_iso_CDM=2, & initial_iso_baryon=3, initial_iso_neutrino=4, initial_iso_neutrino_vel=5, initial_iso_quint=6,initial_vector = 0 integer, parameter :: initial_nummodes = initial_iso_quint type EvolutionVars real(dl) q, q2 real(dl) k_buf,k2_buf ! set in initial integer w_ix !Index of two quintessence equations integer q_ix !index into q_evolve array that gives the value q logical TransferOnly ! nvar (t) - number of scalar (tensor) equations for this k integer nvar,nvarv, nvart !Max_l for the various hierarchies integer lmaxg,lmaxnr,lmaxnu,lmaxgpol,MaxlNeeded integer lmaxnrt, lmaxnut, lmaxt, lmaxpolt integer lmaxnrv, lmaxv, lmaxpolv integer polind !index into scalar array of polarization hierarchy !array indices for massive neutrino equations integer iq0,iq1,iq2 !array index for tensor massive neutrino equations integer iqt !Buffers for non-flat vars real(dl) Kf(max_l_evolve),Kft(max_l_evolve) !Initial values for massive neutrino v*3 variables calculated when switching !to non-relativistic approx real(dl) G11(max_nu),G30(max_nu) real(dl) w_nu !equation of state parameter for massive neutrinos !True when using non-relativistic approximation logical MassiveNuApprox !Massive neutrino scheme being used at the moment integer NuMethod !Tru when using tight-coupling approximation (required for stability at early times) logical TightCoupling !Numer of scalar equations we are propagating integer ScalEqsToPropagate !beta > l for closed models integer FirstZerolForBeta !Tensor vars real(dl) tenspigdot, aux_buf real(dl) pig !For tight coupling real(dl) poltruncfac end type EvolutionVars !precalculated arrays real(dl) polfac(max_l_evolve),tensfac(max_l_evolve),tensfacpol(max_l_evolve), & denl(max_l_evolve) real(dl), parameter :: ep0=1.0d-2 real(dl) epsw integer debug public GetNumEqns,output,outputt,initial,initialt,SwitchToMassiveNuApprox public outtransf,GaugeInterface_EvolveScal contains subroutine GaugeInterface_EvolveScal(EV,tau,y,tauend,tol1,ind,c,w) use ThermoData type(EvolutionVars) EV real(dl) c(24),w(EV%nvar,9), y(EV%nvar), tol1, tau, tauend integer ind real(dl) ep, tau_switch, tau_check real(dl) cs2, opacity !Evolve equations from tau to tauend, performing switch-off of !tight coupling if neccessary. !In principle the tight coupling evolution routine could be different !which could probably be used to improve the speed a bit !It's possible due to instabilities that changing later is actually more !accurate, so cannot expect accuracy to vary monotonically with the switch if (EV%TightCoupling) then tau_switch = tight_tau !when 1/(opacity*tau) = 0.01ish !The numbers here are a bit of guesswork !The high k increase saves time for very small loss of accuracy !The lower k ones are more delicate. Nead avoid instabilities at same time !as ensuring tight coupling is accurate enough if (EV%k_buf > epsw) then if (EV%k_buf > epsw*5) then ep=ep0*5/AccuracyBoost else ep=ep0 end if else ep=ep0 end if !Check the k/opacity criterion tau_check = min(tauend, tau_switch) call thermo(tau_check,cs2,opacity) if (EV%k_buf/opacity > ep) then !so need to switch in this time interval tau_switch = Thermo_OpacityToTime(EV%k_buf/ep) end if if (tauend > tau_switch) then if (tau_switch > tau) then if (CP%flat) then call dverk(EV,EV%ScalEqsToPropagate,fderivs,tau,y,tau_switch,tol1,ind,c,EV%nvar,w) else call dverk(EV,EV%ScalEqsToPropagate,derivs,tau,y,tau_switch,tol1,ind,c,EV%nvar,w) end if end if !Set up variables with their tight coupling values y(8) = EV%pig y(EV%polind+2) = EV%pig/4 EV%TightCoupling = .false. end if end if if (CP%flat) then call dverk(EV,EV%ScalEqsToPropagate,fderivs,tau,y,tauend,tol1,ind,c,EV%nvar,w) else call dverk(EV,EV%ScalEqsToPropagate, derivs,tau,y,tauend,tol1,ind,c,EV%nvar,w) end if end subroutine GaugeInterface_EvolveScal subroutine GaugeInterface_Init !Precompute various arrays integer j epsw = 100/CP%tau0 if (CP%WantScalars) then do j=2,max_l_evolve polfac(j)=real((j+3)*(j-1),dl)/(j+1) end do end if if (CP%WantTensors) then do j=2,max_l_evolve tensfac(j)=real((j+3)*(j-1),dl)/(j+1) tensfacpol(j)=tensfac(j)**2/(j+1) end do end if do j=1,max_l_evolve denl(j)=1._dl/(2*j+1) end do end subroutine GaugeInterface_Init subroutine GetNumEqns(EV) use MassiveNu !Set the numer of equations in each hierarchy, and get total number of equations for this k type(EvolutionVars) EV real(dl) scal if (CP%Num_Nu_massive == 0) then EV%lmaxnu=0 else EV%NuMethod = CP%MassiveNuMethod if (EV%NuMethod == Nu_Best) then if (all(nu_masses(1:CP%Nu_mass_eigenstates) < 1800._dl/AccuracyBoost) .and. & .not. CP%WantTransfer .and. .not. CP%DoLensing) then !If light then approx is very good for CMB EV%NuMethod = Nu_approx else EV%NuMethod = Nu_Trunc end if end if !l_max for massive neutrinos !if relativistic at recombination set to massless lmaxnr below EV%lmaxnu=max(3,nint(6*lAccuracyBoost)) if (CP%Transfer%high_precision) EV%lmaxnu=nint(25*lAccuracyBoost) end if if (CP%WantScalars) then EV%lmaxg = max(nint(8*lAccuracyBoost),3) EV%lmaxnr = max(nint(7*lAccuracyBoost),3) EV%lmaxgpol = EV%lmaxg if (.not.CP%AccuratePolarization) EV%lmaxgpol=max(nint(4*lAccuracyBoost),3) if (EV%q < 0.05) then !Large scales need fewer equations scal = 1 if (CP%AccuratePolarization) scal = 4 !But need more to get polarization right EV%lmaxgpol=max(3,nint(min(8,nint(scal* 150* EV%q))*lAccuracyBoost)) EV%lmaxnr=max(3,nint(min(7,nint(sqrt(scal)* 150 * EV%q))*lAccuracyBoost)) EV%lmaxg=max(3,nint(min(8,nint(sqrt(scal) *300 * EV%q))*lAccuracyBoost)) end if if (CP%Transfer%high_precision) EV%lmaxnr=max(nint(25*lAccuracyBoost),3) EV%nvar=5+ (EV%lmaxg+1) + EV%lmaxgpol-1 +(EV%lmaxnr+1) if (has_extra) then EV%w_ix = EV%nvar+1 EV%nvar=EV%nvar+2 end if if (CP%Num_Nu_massive /= 0) then if (CP%MassiveNuMethod == Nu_approx) then EV%iq0=EV%nvar + 1 EV%nvar= EV%nvar+(EV%lmaxnu+1) else EV%iq0=EV%nvar + 1 EV%iq1=EV%iq0+nqmax EV%iq2=EV%iq1+nqmax EV%nvar=EV%nvar+ nqmax*(EV%lmaxnu+1) endif end if EV%MaxlNeeded=max(EV%lmaxg,EV%lmaxnr,EV%lmaxgpol,EV%lmaxnu) if (EV%MaxlNeeded > max_l_evolve) stop 'Need to increase max_l_evolve' EV%poltruncfac=real(EV%lmaxgpol,dl)/(EV%lmaxgpol-2) EV%polind = 6+EV%lmaxnr+EV%lmaxg EV%lmaxt=0 else EV%nvar=0 end if if (CP%WantTensors) then EV%lmaxt=max(3,nint(8*lAccuracyBoost)) EV%lmaxpolt = max(3,nint(5*lAccuracyBoost)) if (EV%q < 0.05) then !Large scales need fewer equations scal = 1 if (CP%AccuratePolarization) scal = 4 !But need more to get polarization right EV%lmaxt=max(3,nint(min(8,nint(scal *150 * EV%q))*lAccuracyBoost)) EV%lmaxpolt=max(3,nint(min(5,nint( scal*150 * EV%q))*lAccuracyBoost)) end if EV%nvart=(EV%lmaxt-1)+(EV%lmaxpolt-1)*2+3 if (DoTensorNeutrinos) then EV%lmaxnrt=nint(6*lAccuracyBoost) if (EV%q < 0.05) EV%lmaxnrt = max(3,nint(min(6,nint(150 * EV%q))*lAccuracyBoost)) EV%lmaxnut=nint(4*lAccuracyBoost) EV%nvart=EV%nvart+EV%lmaxnrt-1 if (CP%Num_Nu_massive /= 0 ) then EV%iqt=EV%nvart + 1 EV%nvart=EV%nvart+ nqmax*(EV%lmaxnut-1) end if end if else EV%nvart=0 end if end subroutine GetNumEqns !cccccccccccccccccccccccccccccccccc subroutine SwitchToMassiveNuApprox(EV,y) !When the neutrinos are no longer highly relativistic we use a truncated !energy-integrated hierarchy going up to third order in velocity dispersion type(EvolutionVars) EV real(dl) a,a2,pnu,clxnu,dpnu,pinu,rhonu real(dl) shearnu, qnu real(dl) y(EV%nvar) integer nu_i, off_ix, qoff_ix a=y(1) a2=a*a do nu_i = 1, CP%Nu_mass_eigenstates off_ix = (nu_i-1)*4 qoff_ix = (nu_i-1)*nqmax*(EV%lmaxnu+1) !Get density and pressure as ratio to massles by interpolation from table call Nu_background(a*nu_masses(nu_i),rhonu,pnu) !Integrate over q call Nu_Integrate(a*nu_masses(nu_i),clxnu,qnu,dpnu,shearnu, & y(EV%iq0+qoff_ix),y(EV%iq1+qoff_ix),y(EV%iq2+qoff_ix)) !clxnu_here = rhonu*clxnu, qnu_here = qnu*rhonu !Could save time by only calculating shearnu in output dpnu=dpnu/rhonu qnu=qnu/rhonu clxnu = clxnu/rhonu pinu=1.5_dl*shearnu/rhonu EV%MassiveNuApprox=.true. y(EV%iq0+off_ix)=clxnu y(EV%iq0+off_ix+1)=dpnu y(EV%iq0+off_ix+2)=qnu y(EV%iq0+off_ix+3)=pinu call Nu_Intvsq(a*nu_masses(nu_i),EV%G11(nu_i),EV%G30(nu_i),y(EV%iq1+qoff_ix),y(EV%iq0+qoff_ix+3*nqmax)) !Analytic solution for higher moments, proportional to a^{-3} EV%G11(nu_i)=EV%G11(nu_i)*a2*a/rhonu EV%G30(nu_i)=EV%G30(nu_i)*a2*a/rhonu end do EV%ScalEqsToPropagate=(EV%iq0-1)+4*CP%Nu_mass_eigenstates end subroutine SwitchToMassiveNuApprox subroutine MassiveNuVarsOut(EV,y,yprime,a,grho,gpres,dgrho,dgq,dgpi, gdpi_diff,pidot_sum) implicit none type(EvolutionVars) EV real(dl) :: y(EV%nvar), yprime(EV%nvar),a, grho,gpres,dgrho,dgq,dgpi, gdpi_diff,pidot_sum !grho = a^2 kappa rho !gpres = a^2 kappa p !dgrho = a^2 kappa \delta\rho !dgp = a^2 kappa \delta p !dgq = a^2 kappa q (heat flux) !dgpi = a^2 kappa pi (anisotropic stress) !dgpi_diff = a^2 kappa (3*p -rho)*pi integer nu_i, off_ix real(dl) pinudot,grhormass_t, rhonu, pnu, shearnu, rhonudot real(dl) shearnudot, adotoa, grhonu_t,gpnu_t real(dl) clxnu, qnu, pinu, dpnu real(dl) dtauda EV%w_nu = 0 do nu_i = 1, CP%Nu_mass_eigenstates grhormass_t=grhormass(nu_i)/a**2 !Get density and pressure as ratio to massless by interpolation from table call Nu_background(a*nu_masses(nu_i),rhonu,pnu) if (EV%NuMethod == Nu_approx) then off_ix = EV%iq0+(nu_i-1)*(EV%lmaxnu+1) clxnu=y(off_ix) qnu=y(off_ix+1) pinu=y(off_ix+2) pinudot=yprime(off_ix+2) else if (EV%MassiveNuApprox) then off_ix = (nu_i-1)*4 clxnu=y(EV%iq0+off_ix) !dpnu = y(EV%iq0+1+off_ix) qnu=y(EV%iq0+off_ix+2) pinu=y(EV%iq0+off_ix+3) pinudot=yprime(EV%iq0+off_ix+3) else off_ix = (nu_i-1)*nqmax*(EV%lmaxnu+1) !Integrate over q call Nu_Integrate(a*nu_masses(nu_i),clxnu,qnu,dpnu,shearnu, & y(EV%iq0+off_ix),y(EV%iq1+off_ix),y(EV%iq2+off_ix)) !clxnu_here = rhonu*clxnu, qnu_here = qnu*rhonu !dpnu=dpnu/rhonu qnu=qnu/rhonu clxnu = clxnu/rhonu pinu=1.5_dl*shearnu/rhonu adotoa = 1/(a*dtauda(a)) call Nu_derivs(a*nu_masses(nu_i),adotoa,rhonu,rhonudot,shearnudot, & y(EV%iq2+off_ix),yprime(EV%iq2+off_ix)) pinudot=1.5_dl*shearnudot/rhonu - rhonudot/rhonu*pinu endif end if grhonu_t=grhormass_t*rhonu gpnu_t=grhormass_t*pnu grho = grho + grhonu_t gpres= gpres + gpnu_t EV%w_nu = max(EV%w_nu,pnu/rhonu) dgrho= dgrho + grhonu_t*clxnu dgq = dgq + grhonu_t*qnu dgpi = dgpi + grhonu_t*pinu gdpi_diff = gdpi_diff + pinu*(3*gpnu_t-grhonu_t) pidot_sum = pidot_sum + grhonu_t*pinudot end do end subroutine MassiveNuVarsOut subroutine MassiveNuVars(EV,y,a,grho,gpres,dgrho,dgq, wnu_arr) implicit none type(EvolutionVars) EV real(dl) :: y(EV%nvar), a, grho,gpres,dgrho,dgq real(dl), intent(out), optional :: wnu_arr(max_nu) !grho = a^2 kappa rho !gpres = a^2 kappa p !dgrho = a^2 kappa \delta\rho !dgp = a^2 kappa \delta p !dgq = a^2 kappa q (heat flux) integer nu_i, off_ix real(dl) grhormass_t, rhonu, qnu, clxnu, grhonu_t, gpnu_t, pnu do nu_i = 1, CP%Nu_mass_eigenstates grhormass_t=grhormass(nu_i)/a**2 !Get density and pressure as ratio to massless by interpolation from table call Nu_background(a*nu_masses(nu_i),rhonu,pnu) if (EV%NuMethod == Nu_approx) then off_ix = EV%iq0+(nu_i-1)*(EV%lmaxnu+1) clxnu=y(off_ix) qnu=y(off_ix+1) else if (EV%MassiveNuApprox) then off_ix = (nu_i-1)*4 clxnu=y(EV%iq0+off_ix) qnu=y(EV%iq0+off_ix+2) else off_ix = (nu_i-1)*nqmax*(EV%lmaxnu+1) !Integrate over q call Nu_Integrate01(a*nu_masses(nu_i),clxnu,qnu,y(EV%iq0+off_ix),y(EV%iq1+off_ix)) !clxnu_here = rhonu*clxnu, qnu_here = qnu*rhonu qnu=qnu/rhonu clxnu = clxnu/rhonu endif end if grhonu_t=grhormass_t*rhonu gpnu_t=grhormass_t*pnu grho = grho + grhonu_t gpres= gpres + gpnu_t dgrho= dgrho + grhonu_t*clxnu dgq = dgq + grhonu_t*qnu if (present(wnu_arr)) then wnu_arr(nu_i) =pnu/rhonu end if end do end subroutine MassiveNuVars !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine output(EV,y, n,j,tau,sources) use ThermoData use lvalues use ModelData implicit real(dl) (s-t) implicit logical (a-r,u-z) integer n,j type(EvolutionVars) EV real(dl), target :: y(n),yprime(n) real(dl), dimension(:),pointer :: ypol,ypolprime real(dl) dgq,grhob_t,grhor_t,grhoc_t,grhog_t,grhov_t,grhonu_t,gpnu_t,sigma,polter real(dl) qgdot,pigdot,pirdot,vbdot,dgrho real(dl) a,a2,z,clxc,clxb,vb,clxg,qg,pig,clxr,qr,pir,rhonu,pnu,pinu,clxnu,dpnu,qnu real(dl) tau,x,divfac,pinudot,rhonudot,shearnudot, grhormass_t real(dl) shearnu real(dl) k,k2 ,adotoa, grho, gpres,etak,phi, dgpi real(dl) clv, clvdot, phidot, grhoex_t real(dl) sources(CTransScal%NumSources) if (CP%flat) then call fderivs(EV,EV%ScalEqsToPropagate,tau,y,yprime) else call derivs(EV,EV%ScalEqsToPropagate,tau,y,yprime) end if ypolprime => yprime(EV%polind+1:) ypol => y(EV%polind+1:) k=EV%k_buf k2=EV%k2_buf if (EV%TightCoupling) then y(8) = EV%pig ypol(2) = EV%pig/4 end if a =y(1) a2 =a*a etak=y(2) clxc=y(3) clxb=y(4) vb =y(5) clxg=y(6) qg =y(7) pig =y(8) clxr=y(7+EV%lmaxg) qr =y(8+EV%lmaxg) pir =y(9+EV%lmaxg) pigdot=yprime(8) pirdot=yprime(EV%lmaxg+9) vbdot =yprime(5) qgdot =yprime(7) ! Compute expansion rate from: grho 8*pi*rho*a**2 grhob_t=grhob/a grhoc_t=grhoc/a grhor_t=grhornomass/a2 grhog_t=grhog/a2 grhov_t=grhov*a2 grhonu_t=0 gpnu_t=0 grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhonu_t gpres=(grhog_t+grhor_t)/3 +gpnu_t if (has_extra) then call Quint_ValsAta(a,phi,phidot) grhoex_t=phidot**2/2 + a2*Vofphi(phi,0) grho=grho+grhoex_t gpres=gpres - grhoex_t+phidot**2 else grho=grho +grhov_t gpres=gpres - grhov_t end if if (CP%Num_Nu_Massive /= 0) then stop 'doesn''t support massive neutrinos' ! call MassiveNuVarsOut(EV,y,yprime,a,grho,gpres,dgrho,dgq,dgpi, dgpi_diff,pidot_sum) end if adotoa=sqrt((grho+grhok)/3) ! 8*pi*a*a*SUM[rho_i*clx_i] dgrho=grhob_t*clxb+grhoc_t*clxc + grhog_t*clxg+grhor_t*clxr ! 8*pi*a*a*SUM[(rho_i+p_i)*v_i] dgq=grhob_t*vb+grhog_t*qg+grhor_t*qr if (has_extra) then clv=y(EV%w_ix) clvdot=y(EV%w_ix+1) dgrho=dgrho + phidot*clvdot +clv*a2*Vofphi(phi,1) dgq=dgq + k*phidot*clv end if ! Get sigma (shear) and z from the constraints ! have to get z from eta for numerical stability z=(0.5_dl*dgrho/k + etak)/adotoa sigma=(z+1.5_dl*dgq/k2)/EV%Kf(1) polter = 0.1_dl*pig+9._dl/15._dl*ypol(2) if (CP%flat) then x=k*(CP%tau0-tau) divfac=x*x else x=(CP%tau0-tau)/CP%r divfac=(CP%r*rofChi(x))**2*k2 end if if (CP%Num_Nu_Massive == 0) then pinudot=0._dl else end if pinu=0 if (EV%TightCoupling) then dgpi = grhor_t*pir + grhonu_t*pinu + grhog_t*pig pigdot = -dopac(j)/opac(j)*pig + 32._dl/45*k/opac(j)*(-2*adotoa*sigma & +etak/EV%Kf(1)- dgpi/k +vbdot ) ypolprime(2)= pigdot/4 end if t5 = 1/adotoa t14 = 1/k t20 = grhor_t*pir t22 = 0 t30 = k**2 t31 = 1/t30 s1 = 4.D0/3.D0*k*EV%Kf(1)*expmmu(j)*sigma+(-t5*dgrho/3+2.D0/3.D0*(- & sigma-t5*etak)*k+(gpres+5.D0/3.D0*grho)*sigma*t14+((3*pinu*gpnu_t+4* & grhog_t*pig+4*t20+3*t22)*adotoa-grhor_t*pirdot-grhonu_t*pinudot- & grhog_t*pigdot)*t31)*expmmu(j)+(clxg/4+3.D0/8.D0*ypol(2))*vis(j) t101 = s1+(-21.D0/5.D0*adotoa*vis(j)*sigma+(3.D0/40.D0*qgdot-9.D0/ & 80.D0*EV%Kf(2)*yprime(9)+vbdot-3.D0/8.D0*EV%Kf(2)*ypolprime(3))*vis(j)+ & 11.D0/10.D0*sigma*dvis(j)+(vb+3.D0/40.D0*qg-3.D0/8.D0*EV%Kf(2)*ypol(3) & -9.D0/80.D0*EV%Kf(2)*y(9))*dvis(j))*t14+vis(j)*pig/16+((-9.D0/160.D0* & dvis(j)*opac(j)+(-21.D0/10.D0*grhog_t-9.D0/160.D0*dopac(j))*vis(j)& +3.D0/16.D0*ddvis(j))*pig+(3.D0/16.D0*pigdot-27.D0/80.D0*opac(j)* & ypol(2)+9.D0/8.D0*ypolprime(2))*dvis(j)+((-9.D0/160.D0*pigdot-27.D0 & /80.D0*ypolprime(2))*opac(j)-21.D0/10.D0*t22-21.D0/10.D0*t20-27.D0& /80.D0*dopac(j)*ypol(2))*vis(j)+9.D0/8.D0*ddvis(j)*ypol(2))*t31+(-& 2*etak*adotoa*t14*expmmu(j)+21.D0/10.D0*etak*vis(j)*t14)/EV%Kf(1) sources(1)=t101 if (x > 0._dl) then !E polarization source sources(2)=vis(j)*polter*(15._dl/8._dl)/divfac !factor of four because no 1/16 later else sources(2)=0 end if if (CTransScal%NumSources > 2) then !Get lensing sources !Can modify this here if you want to get power spectra for other tracer if (tau>taurend .and. CP%tau0-tau > 0.1_dl) then !phi_lens = phi - 1/2 kappa (a/k)^2 sum_i rho_i pi_i phi = -(dgrho +3*dgq*adotoa/k)/(k2*EV%Kf(1)*2) - (grhor_t*pirdot + grhog_t*pig+ pinu*gpnu_t)/k2 sources(3) = -2*expmmu(j)*phi*(tau-tau_maxvis)/((CP%tau0-tau_maxvis)*(CP%tau0-tau)) !We include the lensing factor of two here else sources(3) = 0 end if end if end subroutine output !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine outputt(EV,yt,n,j,tau,dt,dte,dtb) !calculate the tensor sources for open and closed case use ThermoData implicit none integer j,n type(EvolutionVars) :: EV real(dl), target :: yt(n), ytprime(n) real(dl) tau,dt,dte,dtb,x,polterdot,polterddot,prefac real(dl) pig, pigdot, aux, polter, shear real(dl) sinhxr,cothxor real(dl) k,k2 real(dl), dimension(:),pointer :: E,Bprime,Eprime if (CP%flat) then call fderivst(EV,EV%nvart,tau,yt,ytprime) else call derivst(EV,EV%nvart,tau,yt,ytprime) end if k2=EV%k2_buf k=EV%k_buf pigdot=EV%tenspigdot pig=yt(4) aux=EV%aux_buf shear = yt(3) x=(CP%tau0-tau)/CP%r sinhxr=rofChi(x)*CP%r if (EV%q*sinhxr > 1.e-8_dl) then prefac=sqrt(EV%q2*CP%r*CP%r-CP%Ksign) cothxor=cosfunc(x)/sinhxr E => yt(EV%lmaxt+3-1:) Eprime=> ytprime(EV%lmaxt+3-1:) Bprime => Eprime(EV%lmaxpolt:) polter = 0.1_dl*pig + 9._dl/15._dl*E(2) polterdot=9._dl/15._dl*Eprime(2) + 0.1_dl*pigdot polterddot = 9._dl/15._dl*(-dopac(j)*(E(2)-polter)-opac(j)*( & Eprime(2)-polterdot) + k*(2._dl/3._dl*Bprime(2)*aux - & 5._dl/27._dl*Eprime(3)*EV%Kft(2))) & +0.1_dl*(k*(-ytprime(5)*EV%Kft(2)/3._dl + 8._dl/15._dl*ytprime(3)) - & dopac(j)*(pig - polter) - opac(j)*(pigdot-polterdot)) dt=(shear*expmmu(j) + (15._dl/8._dl)*polter*vis(j)/k)*CP%r/sinhxr**2/prefac dte=CP%r*15._dl/8._dl/k/prefac* & ((ddvis(j)*polter + 2._dl*dvis(j)*polterdot + vis(j)*polterddot) & + 4._dl*cothxor*(dvis(j)*polter + vis(j)*polterdot) - & vis(j)*polter*(k2 -6*cothxor**2)) dtb=15._dl/4._dl*EV%q*CP%r/k/prefac*(vis(j)*(2._dl*cothxor*polter + polterdot) + dvis(j)*polter) else dt=0._dl dte=0._dl dtb=0._dl end if end subroutine outputt !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine initial(EV,y, tau) ! Initial conditions. implicit none type(EvolutionVars) EV real(dl) y(EV%nvar) real(dl) Rv,Rg,Rp15,tau,x,x2,x3,aj3r,om,omtau,eta,Rc,Rb,grhonu,chi,phi,iqg real(dl) k,k2 real(dl) a,a2,clxc,clxg,clxr,clxb,qg,qr,vb,pir integer l,i integer, parameter :: i_clxg=1,i_clxr=2,i_clxc=3, i_clxb=4, & i_qg=5,i_qr=6,i_vb=7,i_pir=8, i_eta=9, i_aj3r=10,i_clv=11,i_clvdot=12 integer, parameter :: i_max = i_clvdot real(dl) initv(6,1:i_max), initvec(1:i_max) EV%ScalEqsToPropagate=EV%nvar if (CP%flat) then EV%k_buf=EV%q EV%k2_buf=EV%q2 EV%Kf(1:EV%MaxlNeeded)=1._dl !initialize for CP%flat case else EV%k2_buf=EV%q2-CP%curv EV%k_buf=sqrt(EV%k2_buf) do l=1,EV%MaxlNeeded EV%Kf(l)=1._dl-CP%curv*(l*(l+2))/EV%k2_buf end do end if k=EV%k_buf k2=EV%k2_buf if (CP%closed) then EV%FirstZerolForBeta = nint(EV%q*CP%r) else EV%FirstZerolForBeta=l0max !a large number end if ! k*tau, (k*tau)**2, (k*tau)**3 x=k*tau x2=x*x x3=x2*x grhonu=grhornomass om = (grhob+grhoc)/sqrt(3*(grhog+grhonu)) omtau=om*tau Rv=grhonu/(grhonu+grhog) Rg = 1-Rv Rc=CP%omegac/(CP%omegac+CP%omegab) Rb=1-Rc Rp15=4*Rv+15 if (CP%Scalar_initial_condition > initial_nummodes) & stop 'Invalid initial condition for scalar modes' a=tau*adotrad*(1+omtau/4) a2=a*a initv=0 ! Set adiabatic initial conditions chi=1 !Get transfer function for chi initv(1,i_clxg)=-chi*EV%Kf(1)/3*x2*(1-omtau/5) initv(1,i_clxr)= initv(1,i_clxg) initv(1,i_clxb)=0.75_dl*initv(1,i_clxg) initv(1,i_clxc)=initv(1,i_clxb) initv(1,i_qg)=initv(1,i_clxg)*x/9._dl initv(1,i_qr)=-chi*EV%Kf(1)*(4*Rv+23)/Rp15*x3/27 initv(1,i_vb)=0.75_dl*initv(1,i_qg) initv(1,i_pir)=chi*4._dl/3*x2/Rp15*(1+omtau/4*(4*Rv-5)/(2*Rv+15)) initv(1,i_aj3r)=chi*4/21._dl/Rp15*x3 initv(1,i_eta)=-chi*2*EV%Kf(1)*(1 - x2/12*(-10._dl/Rp15 + EV%Kf(1))) if (CP%Scalar_initial_condition/= initial_adiabatic) then !CDM isocurvature initv(2,i_clxg)= Rc*omtau*(-2._dl/3 + omtau/4) initv(2,i_clxr)=initv(2,i_clxg) initv(2,i_clxb)=initv(2,i_clxg)*0.75_dl initv(2,i_clxc)=1+initv(2,i_clxb) initv(2,i_qg)=-Rc/9*omtau*x initv(2,i_qr)=initv(2,i_qg) initv(2,i_vb)=0.75_dl*initv(2,i_qg) initv(2,i_pir)=-Rc*omtau*x2/3/(2*Rv+15._dl) initv(2,i_eta)= Rc*omtau*(1._dl/3 - omtau/8)*EV%Kf(1) initv(2,i_aj3r)=0 !Baryon isocurvature if (Rc==0) stop 'Isocurvature initial conditions assume non-zero dark matter' initv(3,:) = initv(2,:)*(Rb/Rc) initv(3,i_clxc) = initv(3,i_clxb) initv(3,i_clxb)= initv(3,i_clxb)+1 !neutrino isocurvature density mode initv(4,i_clxg)=Rv/Rg*(-1 + x2/6) initv(4,i_clxr)=1-x2/6 initv(4,i_clxc)=-omtau*x2/80*Rv*Rb/Rg initv(4,i_clxb)= Rv/Rg/8*x2 iqg = - Rv/Rg*(x/3 - Rb/4/Rg*omtau*x) initv(4,i_qg) =iqg initv(4,i_qr) = x/3 initv(4,i_vb)=0.75_dl*iqg initv(4,i_pir)=x2/Rp15 initv(4,i_eta)=EV%Kf(1)*Rv/Rp15/3*x2 initv(6,i_aj3r)=0 !neutrino isocurvature velocity mode initv(5,i_clxg)=Rv/Rg*x - 2*x*omtau/16*Rb*(2+Rg)/Rg**2 initv(5,i_clxr)=-x -3*x*omtau*Rb/16/Rg initv(5,i_clxc)=-9*omtau*x/64*Rv*Rb/Rg initv(5,i_clxb)= 3*Rv/4/Rg*x - 9*omtau*x/64*Rb*(2+Rg)/Rg**2 iqg = Rv/Rg*(-1 + 3*Rb/4/Rg*omtau+x2/6 +3*omtau**2/16*Rb/Rg**2*(Rg-3*Rb)) initv(5,i_qg) =iqg initv(5,i_qr) = 1 - x2/6*(1+4*EV%Kf(1)/(4*Rv+5)) initv(5,i_vb)=0.75_dl*iqg initv(5,i_pir)=2*x/(4*Rv+5)+omtau*x*6/Rp15/(4*Rv+5) initv(5,i_eta)=2*EV%Kf(1)*x*Rv/(4*Rv+5) + omtau*x*3*EV%Kf(1)*Rv/32*(Rb/Rg - 80/Rp15/(4*Rv+5)) initv(5,i_aj3r) = 3._dl/7*x2/(4*Rv+5) !quintessence isocurvature mode initv(6,i_clv) = 1 - x2/6+omtau*x2/72 initv(6,i_clvdot) = k*(-x/3 + omtau*x/24) end if if (CP%Scalar_initial_condition==initial_vector) then InitVec = 0 do i=1,initial_nummodes InitVec = InitVec+ initv(i,:)*CP%InitialConditionVector(i) end do else InitVec = initv(CP%Scalar_initial_condition,:) if (CP%Scalar_initial_condition==initial_adiabatic) InitVec = -InitVec !So we start with chi=-1 as before end if y(1)=a y(2)= -InitVec(i_eta)*k/2 !get eta_s*k, where eta_s is synchronous gauge variable ! CDM y(3)=InitVec(i_clxc) ! Baryons y(4)=InitVec(i_clxb) y(5)=InitVec(i_vb) ! Photons y(6)=InitVec(i_clxg) y(7)=InitVec(i_qg) y(8:EV%nvar)=0._dl if (has_extra) then y(EV%w_ix) = InitVec(i_clv) y(EV%w_ix+1) = InitVec(i_clvdot) end if ! Neutrinos y(7+EV%lmaxg)=InitVec(i_clxr) y(8+EV%lmaxg)=InitVec(i_qr) y(9+EV%lmaxg)=InitVec(i_pir) if (EV%FirstZerolForBeta==3) then y(10+EV%lmaxg)=0._dl else y(10+EV%lmaxg)=InitVec(i_aj3r) endif EV%TightCoupling=.true. EV%MassiveNuApprox=.false. if (CP%Num_Nu_massive == 0) return if (CP%MassiveNuMethod == Nu_approx) then !Values are just same as massless neutrino ones since highly relativistic y(EV%iq0)=InitVec(i_clxr) y(EV%iq0+1)=InitVec(i_qr) y(EV%iq0+2)=InitVec(i_pir) if (EV%FirstZerolForBeta/=3) then y(EV%iq0+3)=InitVec(i_aj3r) endif else do i=1,nqmax y(EV%iq0+i-1)=-0.25_dl*dlfdlq(i)*InitVec(i_clxr) y(EV%iq1+i-1)=-0.25_dl*dlfdlq(i)*InitVec(i_qr) y(EV%iq2+i-1)=-0.25_dl*dlfdlq(i)*InitVec(i_pir) end do end if end subroutine initial !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine initialt(EV,yt,tau) ! Initial conditions for tensors implicit none real(dl) bigR,tau,x,aj3r,elec, pir integer l type(EvolutionVars) EV real(dl) k,k2 ,a real(dl) yt(EV%nvart) if (CP%flat) then EV%aux_buf=1._dl EV%k2_buf=EV%q2 EV%k_buf=EV%q EV%Kft(1:EV%lmaxt)=1._dl!initialize for CP%flat case else EV%k2_buf=EV%q2-3*CP%curv EV%k_buf=sqrt(EV%k2_buf) EV%aux_buf=sqrt(1._dl+3*CP%curv/EV%k2_buf) do l=1,EV%lmaxt EV%Kft(l)=1._dl-CP%curv*((l+1)**2-3)/EV%k2_buf end do endif k=EV%k_buf k2=EV%k2_buf if (CP%closed) then EV%FirstZerolForBeta = nint(EV%q*CP%r) else EV%FirstZerolForBeta=l0max !a large number end if a=tau*adotrad if (DoTensorNeutrinos) then bigR = (grhornomass)/(grhornomass+grhog) else bigR = 0._dl end if yt(1)=a yt(2)= 1._dl elec=-(1+2*CP%curv/k2)*(2*bigR+10)/(4*bigR+15) !elec, with H=1 x=k*tau !shear yt(3)=-5._dl/2/(bigR+5)*x*elec yt(4:EV%nvart)=0._dl ! Neutrinos if (DoTensorNeutrinos) then pir=-2._dl/3._dl/(bigR+5)*x**2*elec aj3r= -2._dl/21._dl/(bigR+5)*x**3*elec yt((EV%lmaxt-1)+(EV%lmaxt-1)*2+3+1)=pir yt((EV%lmaxt-1)+(EV%lmaxt-1)*2+3+2)=aj3r end if end subroutine initialt subroutine outtransf(EV, y, Arr) !write out clxc, clxb, clxg, clxn use Transfer implicit none type(EvolutionVars) EV real(dl) clxc, clxb, clxg, clxr, clxnu, clxtot,k,k2 real Arr(Transfer_max) real(dl) y(EV%nvar) real(dl) qnu, a, rhonu, pnu a = y(1) clxc = y(3) clxb = y(4) clxg = y(6) clxr = y(7+EV%lmaxg) k = EV%k_buf k2 = EV%k2_buf clxnu = 0 rhonu = 0 clxtot = (clxc*grhoc/a + clxb*grhob/a ) clxtot=clxtot/(grhoc/a+grhob/a) Arr(Transfer_kh) = k/(CP%h0/100._dl) Arr(Transfer_cdm) = clxc/k2 Arr(Transfer_b) = clxb/k2 Arr(Transfer_g) = clxg/k2 Arr(Transfer_r) = clxr/k2 Arr(Transfer_nu) = clxnu/k2 Arr(Transfer_tot) = clxtot/k2 end subroutine outtransf !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine fderivs(EV,n,tau,ay,ayprime) ! Evaluate the time derivatives of the perturbations, flat case ! ayprime is not neccessarily GaugeInterface.yprime, so keep them distinct use ThermoData use MassiveNu implicit none type(EvolutionVars) EV integer n real(dl) ay(n),ayprime(n) real(dl) tau,ep,w real(dl) k,k2 ! Internal variables. real(dl) opacity !!FIX March 2005 real(dl) E2 real(dl) photbar,cs2,pb43,grho,slip,clxgdot,grhormass_t, & clxcdot,clxbdot,drag,adotdota,gpres,clxrdot,etak real(dl) q,aq,v,akv(nqmax0) real(dl) G11_t,G30_t real(dl) dgq,grhob_t,grhor_t,grhoc_t,grhog_t,grhov_t,grhonu_t,gpnu_t,sigma,polter real(dl) qgdot,qrdot,pigdot,pirdot,vbdot,dgrho,adotoa real(dl) a,a2,z,clxc,clxb,vb,clxg,qg,pig,clxr,qr,pir real(dl) rhonu,pnu,clxnu,qnu,f, dopacity real(dl) phi,phidot, clv, clvdot, grhoex_t integer l,i,ind logical tcp k=EV%k_buf k2=EV%k2_buf a=ay(1) etak=ay(2) ! CDM variables clxc=ay(3) ! Baryon variables clxb=ay(4) vb=ay(5) ! Photons clxg=ay(6) qg=ay(7) pig=ay(8) ! Massless neutrinos clxr=ay(7+EV%lmaxg) qr=ay(8+EV%lmaxg) pir=ay(9+EV%lmaxg) !!FIX March 2005 E2=ay(EV%polind+2) a2=a*a ! Get sound speed and ionisation fraction. if (EV%TightCoupling) then call thermo(tau,cs2,opacity,dopacity) else call thermo(tau,cs2,opacity) end if ! Compute expansion rate from: grho 8*pi*rho*a**2 grhob_t=grhob/a grhoc_t=grhoc/a grhor_t=grhornomass/a2 grhog_t=grhog/a2 grhov_t=grhov*a2 grhormass_t=0 if (CP%Num_Nu_Massive == 0) then rhonu=1 pnu=1._dl/3._dl clxnu=0 qnu = 0 else end if grhonu_t=grhormass_t*rhonu gpnu_t=grhormass_t*pnu grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhonu_t gpres=(grhog_t+grhor_t)/3._dl +gpnu_t if (has_extra) then call Quint_ValsAta(a,phi,phidot) grhoex_t=phidot**2/2 + a2*Vofphi(phi,0) grho=grho+grhoex_t gpres=gpres-grhoex_t+phidot**2 ! write(*,'(4E15.5)') a, tau, grhoex_t/grho, (-grhoex_t+phidot**2)/grhoex_t !if (a>0.99) stop else grho=grho +grhov_t gpres=gpres - grhov_t end if adotoa=sqrt(grho/3) ayprime(1)=adotoa*a ! Photon mass density over baryon mass density photbar=grhog_t/grhob_t pb43=4._dl/3._dl*photbar ! Evaluate the time derivative of sz (gradient of expansion) ! 8*pi*a*a*SUM[rho_i*clx_i] dgrho=grhob_t*clxb+grhoc_t*clxc + grhog_t*clxg+grhor_t*clxr +grhonu_t*clxnu ! 8*pi*a*a*SUM[(rho_i+p_i)*v_i] dgq=grhob_t*vb+grhog_t*qg+grhor_t*qr + grhonu_t*qnu if (has_extra) then clv=ay(EV%w_ix) clvdot=ay(EV%w_ix+1) dgrho=dgrho + phidot*clvdot +clv*a2*Vofphi(phi,1) dgq=dgq + k*phidot*clv end if !!! if (k>1e-5) then !tau a w Om_Q delta_Q ! write(*,'(6E15.5)') tau,a,(-grhoex_t+phidot**2)/grhoex_t,grhoex_t/grho, & ! (phidot*clvdot +clv*a2*Vofphi(phi,1))/grhoex_t,clxc ! end if ! Get sigma (shear) and z from the constraints ! have to get z from eta for numerical stability z=(0.5_dl*dgrho/k + etak)/adotoa sigma=z+1.5_dl*dgq/k2 ! ddota/a adotdota=0.5_dl*(adotoa*adotoa-gpres) if (has_extra) then ayprime(EV%w_ix)= clvdot ayprime(EV%w_ix+1) = - 2*adotoa*clvdot - k*z*phidot - k2*clv - clv*Vofphi(phi,2) end if !eta*k equation ayprime(2)=0.5_dl*dgq ! CDM equation of motion clxcdot=-k*z ayprime(3)=clxcdot ! Baryon equation of motion. clxbdot=-k*(z+vb) ayprime(4)=clxbdot ! Photon equation of motion clxgdot=-k*(4._dl/3._dl*z+qg) ! Drag: a*n_e*sigma_T*(4/3*vb-qg) drag=opacity*(4._dl/3._dl*vb-qg) ! use tight coupling approx? if (k > epsw) then ep=ep0 else ep=0.5_dl*ep0 end if ! There was an instability for low k when we ! switched off tight coupling. Minor change in ! criteria ends the problem. (Thanks to David Spergel) tcp = ((k/opacity > ep).or.(1._dl/(opacity*tau) > ep .and. k/opacity > 1d-4)) ! Use explicit equation for vb if appropriate if (EV%TightCoupling) then !!FIX March 2005 pig = 32._dl/45/opacity*k*(sigma+vb) E2 = pig/4 EV%pig = pig ! Use tight-coupling approximation for vb ! zeroth order (in t_c) approximation to vbdot vbdot=(-adotoa*vb+cs2*k*clxb & +0.25_dl*k*pb43*(clxg-2._dl*pig))/(1._dl+pb43) slip = - (2*adotoa/(1+pb43) + dopacity/opacity)*(vb-3._dl/4*qg) & +(-adotdota*vb-k/2*adotoa*clxg +k*(cs2*clxbdot-clxgdot/4))/(opacity*(1+pb43)) ! First-order approximation to vbdot vbdot=vbdot+pb43/(1._dl+pb43)*slip else vbdot=-adotoa*vb+cs2*k*clxb-photbar*opacity*(4._dl/3*vb-qg) end if ayprime(5)=vbdot ! Photon equations of motion ayprime(6)=clxgdot qgdot=4._dl/3._dl*(-vbdot-adotoa*vb+cs2*k*clxb)/pb43 & +k*(clxg-2._dl*pig)/3._dl ayprime(7)=qgdot !!FIX March 2005 polter = pig/10+9._dl/15*E2 ! Use explicit equations for photon moments if appropriate if (EV%tightcoupling) then ! Use tight-coupling approximation where moments are zero for l>1 pigdot=0._dl ayprime(8:EV%lmaxg+6)=0._dl else pigdot=0.4_dl*k*qg-0.6_dl*k*ay(9)-opacity*(pig - polter) & +8._dl/15._dl*k*sigma ayprime(8)=pigdot do l=3,EV%lmaxg-1 ayprime(l+6)=k*denl(l)*(l*ay(l+5)-(l+1)*ay(l+7))-opacity*ay(l+6) end do ! Truncate the photon moment expansion ayprime(EV%lmaxg+6)=k*ay(EV%lmaxg+5)-(EV%lmaxg+1)/tau*ay(EV%lmaxg+6) & -opacity*ay(EV%lmaxg+6) end if ! Massless neutrino equations of motion. clxrdot=-k*(4._dl/3._dl*z+qr) ayprime(EV%lmaxg+7)=clxrdot qrdot=k*(clxr-2._dl*pir)/3._dl ayprime(EV%lmaxg+8)=qrdot pirdot=k*(0.4_dl*qr-0.6_dl*ay(EV%lmaxg+10)+8._dl/15._dl*sigma) ayprime(EV%lmaxg+9)=pirdot do l=3,EV%lmaxnr-1 ayprime(l+EV%lmaxg+7)=k*denl(l)*(l*ay(l+EV%lmaxg+6) -(l+1)*ay(l+EV%lmaxg+8)) end do ! Truncate the neutrino expansion ayprime(EV%lmaxnr+EV%lmaxg+7)=k*ay(EV%lmaxnr+EV%lmaxg+6)- & (EV%lmaxnr+1)/tau*ay(EV%lmaxnr+EV%lmaxg+7) ! Polarization if (EV%TightCoupling) then ayprime(EV%polind+2:EV%polind+EV%lmaxgpol)=0._dl else !l=2 (defined to be zero for l<2) ayprime(EV%polind+2) = -opacity*(ay(EV%polind+2) - polter) - k/3._dl*ay(EV%polind+3) !and the rest do l=3,EV%lmaxgpol-1 ayprime(EV%polind+l)=-opacity*ay(EV%polind+l) + k*denl(l)*(l*ay(EV%polind+l-1) -& polfac(l)*ay(EV%polind+l+1)) end do !truncate ayprime(EV%polind+EV%lmaxgpol)=-opacity*ay(EV%polind+EV%lmaxgpol) + & k*EV%poltruncfac*ay(EV%polind+EV%lmaxgpol-1)-(EV%lmaxgpol+3)*ay(EV%polind+EV%lmaxgpol)/tau end if ! Massive neutrino equations of motion. if (CP%Num_Nu_massive == 0) return end subroutine fderivs subroutine fderivsv(EV,n,tau,yv,yvprime) ! Evaluate the time derivatives of the vector perturbations, flat case implicit none type(EvolutionVars) EV integer n real(dl), target :: yv(n),yvprime(n) real(dl) tau stop 'vectors not implemented in this module' end subroutine fderivsv subroutine initialv(EV,yv,tau) ! Initial conditions for vectors implicit none type(EvolutionVars) EV real(dl) yv(EV%nvarv) real(dl) tau stop 'vectors not implemented in this module' end subroutine initialv subroutine outputv(EV,yv,n,j,tau,dt,dte,dtb) !calculate the vector sources use ThermoData implicit none integer n,j type(EvolutionVars) :: EV real(dl), target :: yv(n), yvprime(n) real(dl) tau,dt,dte,dtb stop 'vectors not implemented in this module' end subroutine outputv !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine fderivst(EV,n,tau,ayt,aytprime) ! Evaluate the time derivatives of the tensor perturbations, flat case use ThermoData use MassiveNu implicit none type(EvolutionVars) EV integer n,l,i,ind real(dl), target :: ayt(n),aytprime(n) real(dl) ep,tau,grho,rhopi,cs2,opacity,gpres,pirdt,grhormass_t logical tcp real(dl), dimension(:),pointer :: neut,neutprime,E,B,Eprime,Bprime real(dl) q,aq,v,akv(nqmax0),shearnu real(dl) grhob_t,grhor_t,grhoc_t,grhog_t,grhov_t,grhonu_t,gpnu_t,polter real(dl) Hchi,shear, pig real(dl) k,k2,a,a2 real(dl) pinu,pir,pnu,adotoa, rhonu real(dl) phi,phidot, grhoex_t k2=EV%k2_buf k=EV%k_buf !E and B start at l=2. Set up pointers accordingly to fill in y arrays E => ayt(EV%lmaxt+3-1:) Eprime=> aytprime(EV%lmaxt+3-1:) B => E(EV%lmaxpolt:) Bprime => Eprime(EV%lmaxpolt:) a=ayt(1) !Hchi = metric perturbation variable, conserved on large scales Hchi=ayt(2) !shear (Hchi' = - k*shear) shear=ayt(3) ! Photon anisotropic stress pig=ayt(4) a2=a*a ! Get sound speed and opacity, and see if should use tight-coupling call thermo(tau,cs2,opacity) if (k > 0.06_dl*epsw) then ep=ep0 else ep=0.2_dl*ep0 end if tcp = ((k/opacity > ep).or.(1._dl/(opacity*tau) > ep)) !Do massive neutrinos if (CP%Num_Nu_Massive ==0) then rhonu=1._dl pnu=1._dl/3._dl else stop end if ! Compute expansion rate from: grho=8*pi*rho*a**2 ! Also calculate gpres: 8*pi*p*a**2 grhob_t=grhob/a grhoc_t=grhoc/a grhor_t=grhornomass/a2 grhog_t=grhog/a2 grhov_t=grhov*a2 grhormass_t=0 grhonu_t=grhormass_t*rhonu gpnu_t=grhormass_t*pnu grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhonu_t gpres=(grhog_t+grhor_t)/3._dl+gpnu_t if (has_extra) then call Quint_ValsAta(a,phi,phidot) grhoex_t=phidot**2/2 + a2*Vofphi(phi,0) grho=grho+grhoex_t gpres=gpres-grhoex_t+phidot**2 else grho=grho +grhov_t gpres=gpres - grhov_t end if adotoa=sqrt(grho/3._dl) aytprime(1)=adotoa*a polter = 0.1_dl*pig + 9._dl/15._dl*E(2) if (tcp) then ! Use explicit equations: ! Equation for the photon anisotropic stress aytprime(4)=k*(-1._dl/3._dl*ayt(5)+8._dl/15._dl*shear) & -opacity*(pig - polter) ! And for the moments do l=3,EV%lmaxt-1 aytprime(l+2)=k*denl(l)*(l*ayt(l+1)- & tensfac(l)*ayt(l+3))-opacity*ayt(l+2) end do ! Truncate the hierarchy aytprime(EV%lmaxt+2)=k*EV%lmaxt/(EV%lmaxt-2._dl)*ayt(EV%lmaxt+1)- & (EV%lmaxt+3._dl)*ayt(EV%lmaxt+2)/tau-opacity*ayt(EV%lmaxt+2) !E equations Eprime(2) = - opacity*(E(2) - polter) + k*(4._dl/6._dl*B(2) - & 5._dl/27._dl*E(3)) do l=3,EV%lmaxpolt-1 Eprime(l) =-opacity*E(l) + k*(denl(l)*(l*E(l-1) - & tensfacpol(l)*E(l+1)) + 4._dl/(l*(l+1))*B(l)) end do !truncate Eprime(EV%lmaxpolt)=0._dl !B-bar equations do l=2,EV%lmaxpolt-1 Bprime(l) =-opacity*B(l) + k*(denl(l)*(l*B(l-1) - & tensfacpol(l)*B(l+1)) - 4._dl/(l*(l+1))*E(l)) end do !truncate Bprime(EV%lmaxpolt)=0._dl else ayt(4)=32._dl/45._dl*k/opacity*shear E(2)=ayt(4)/4._dl ! Set the derivatives to zero aytprime(4:n)=0._dl endif ! Neutrino equations: ! Anisotropic stress if (DoTensorNeutrinos) then neutprime => Bprime(EV%lmaxpolt:) neut => B(EV%lmaxpolt:) ! Massless neutrino anisotropic stress pir=neut(2) pirdt=-1._dl/3._dl*k*neut(3)+ 8._dl/15._dl*k*shear neutprime(2)=pirdt ! And for the moments do l=3,EV%lmaxnrt-1 neutprime(l)=k*denl(l)*(l*neut(l-1)- tensfac(l)*neut(l+1)) end do ! Truncate the hierarchy neutprime(EV%lmaxnrt)=k*EV%lmaxnrt/(EV%lmaxnrt-2._dl)*neut(EV%lmaxnrt-1)- & (EV%lmaxnrt+3._dl)*neut(EV%lmaxnrt)/tau ! Massive neutrino equations of motion. if (CP%Num_Nu_massive /= 0) then else pinu=0._dl end if else pinu=0 pir=0 pirdt=0 end if ! Get the propagation equation for the shear rhopi=grhog_t*pig+grhor_t*pir+grhonu_t*pinu aytprime(3)=-2*adotoa*shear+k*Hchi-rhopi/k if (tcp) then ! Use the full expression for pigdt EV%tenspigdot=aytprime(4) else ! Use the tight-coupling approximation EV%tenspigdot= 32._dl/45._dl*k/opacity*(2._dl*adotoa*shear+aytprime(3)) endif aytprime(2)=-k*shear end subroutine fderivst !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine derivs(EV,n,tau,ay,ayprime) ! Evaluate the time derivatives of the perturbations. !ayprime is not neccessarily GaugeInterface.yprime, so keep them distinct use ThermoData use MassiveNu implicit none type(EvolutionVars) EV integer n real(dl) ay(n),ayprime(n) real(dl) tau real(dl) dgq,grhob_t,grhor_t,grhoc_t,grhog_t,grhov_t,grhonu_t,gpnu_t,sigma,polter real(dl) qgdot,qrdot,pigdot,pirdot,vbdot,dgrho,adotoa real(dl) a,a2,z,clxc,clxb,vb,clxg,qg,pig,clxr,qr,pir,rhonu,pnu,clxnu,qnu real(dl) k,k2 ! Internal variables. real(dl) opacity,ep,cothxor,w real(dl) photbar,cs2,pb43,grho,slip,clxgdot,etak, & clxcdot,clxbdot,drag,adotdota,gpres,clxrdot real(dl) q,aq,v,akv(nqmax0),grhormass_t integer l,i,ind logical tcp !false if can use tight coupling approx real(dl) G11_t,G30_t, f, E2 real(dl) clv,clvdot, phidot, phi,grhoex_t stop 'Non-flat routines out of date' k2=EV%k2_buf k=EV%k_buf cothxor=1._dl/tanfunc(tau/CP%r)/CP%r a=ay(1) !z=ay(2)/a etak=ay(2) ! CDM variables clxc=ay(3) ! Baryon variables clxb=ay(4) vb=ay(5) ! Photons clxg=ay(6) qg=ay(7) pig=ay(8) E2=ay(EV%polind+2) ! Massless neutrinos clxr=ay(7+EV%lmaxg) qr=ay(8+EV%lmaxg) pir=ay(9+EV%lmaxg) a2=a*a ! Get baryon temperature, sound speed and ionisation fraction. call thermo(tau,cs2,opacity) ! Compute expansion rate from: grho 8*pi*rho*a**2 grhob_t=grhob/a grhoc_t=grhoc/a grhor_t=grhornomass/a2 grhog_t=grhog/a2 grhov_t=grhov*a2 grhormass_t=0 if (CP%Num_Nu_Massive == 0) then rhonu=1._dl pnu=1._dl/3._dl clxnu=0._dl qnu = 0 else end if grhonu_t=grhormass_t*rhonu gpnu_t=grhormass_t*pnu grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhonu_t gpres=(grhog_t+grhor_t)/3._dl +gpnu_t if (has_extra) then call Quint_ValsAta(a,phi,phidot) grhoex_t=phidot**2/2 + a2*Vofphi(phi,0) grho=grho+grhoex_t gpres=gpres-grhoex_t+phidot**2 else grho=grho +grhov_t gpres=gpres - grhov_t end if adotoa=sqrt((grho+grhok)/3._dl) ayprime(1)=adotoa*a ! Photon mass density over baryon mass density photbar=grhog_t/grhob_t pb43=4._dl/3._dl*photbar ! 8*pi*a*a*SUM[rho_i*clx_i] dgrho=grhob_t*clxb+grhoc_t*clxc + grhog_t*clxg+grhor_t*clxr +grhonu_t*clxnu ! 8*pi*a*a*SUM[(rho_i+p_i)*v_i] dgq=grhob_t*vb+grhog_t*qg+grhor_t*qr+ grhonu_t*qnu if (has_extra) then clv=ay(EV%w_ix) clvdot=ay(EV%w_ix+1) dgrho=dgrho + phidot*clvdot +clv*a2*Vofphi(phi,1) dgq=dgq + k*phidot*clv end if ! Get sigma (shear) and z from the constraints ! have to get z from eta for numerical stability z=(0.5_dl*dgrho/k + etak)/adotoa sigma=(z+1.5_dl*dgq/k2)/EV%Kf(1) ! ddota/a adotdota=0.5_dl*(adotoa*adotoa-gpres) if (has_extra) then ayprime(EV%w_ix)= clvdot ayprime(EV%w_ix+1) = - 2*adotoa*clvdot - k*z*phidot - k2*clv - clv*Vofphi(phi,2) end if ! eta*k equation ayprime(2)=0.5_dl*dgq + CP%curv*z ! CDM equation of motion clxcdot=-k*z ayprime(3)=clxcdot ! Baryon equation of motion. clxbdot=-k*(z+vb) ayprime(4)=clxbdot ! Photon equation of motion clxgdot=-k*(4._dl/3._dl*z+qg) ! Drag: a*n_e*sigma_T*(4/3*vb-qg) drag=opacity*(4._dl/3._dl*vb-qg) ! use tight coupling approx? if (k > 0.06_dl*epsw) then ep=ep0 else ep=1.1_dl*ep0 end if ! There was an instability for low k when we ! switched off tight coupling. Minor change in ! criteria ends the problem. (Thanks to David Spergel) tcp = ((EV%q/opacity > ep).or.(1._dl/(opacity*tau) > ep .and. EV%q/opacity > 1d-4)) ! Use explicit equation for vb if appropriate if (tcp) then vbdot=-adotoa*vb+cs2*k*clxb-photbar*drag else ! Use tight-coupling approximation for vb ! zeroth order (in t_c) approximation to vbdot vbdot=(-adotoa*vb+cs2*k*clxb & +0.25_dl*k*pb43*(clxg-2._dl*EV%Kf(1)*pig))/(1._dl+pb43) ! First-order (in t_c) approximation to baryon-photon splip slip=2._dl*pb43/(1._dl+pb43)*adotoa*(vb-3._dl/4._dl*qg) & +1._dl/opacity*(-adotdota*vb-0.5_dl*k*adotoa*clxg & +k*(cs2*clxbdot-0.25_dl*clxgdot))/(1._dl+pb43) ! First-order approximation to vbdot vbdot=vbdot+pb43/(1._dl+pb43)*slip end if ayprime(5)=vbdot ! Photon equations of motion ayprime(6)=clxgdot qgdot=4._dl/3._dl*(-vbdot-adotoa*vb+cs2*k*clxb)/pb43 & +k*(clxg-2._dl*pig*EV%Kf(1))/3._dl ayprime(7)=qgdot if (EV%FirstZerolForBeta <= EV%MaxlNeeded) ayprime(8:EV%ScalEqsToPropagate)=0 !bit lazy this ! Use explicit equations for photon moments if appropriate polter = 0.1_dl*pig+9._dl/15._dl*ay(EV%polind+2) !2/15*(3/4 pig + 9/2 E2) if (tcp) then pigdot=0.4_dl*k*qg-0.6_dl*k*EV%Kf(2)*ay(9)-opacity*(pig - polter) & +8._dl/15._dl*k*sigma ayprime(8)=pigdot do l=3,min(EV%FirstZerolForBeta,EV%lmaxg)-1 ayprime(l+6)=k*denl(l)*(l*ay(l+5)-(l+1)*EV%Kf(l)*ay(l+7))-opacity*ay(l+6) end do ! Truncate the photon moment expansion if (EV%lmaxg/=EV%FirstZerolForBeta) then ayprime(EV%lmaxg+6)=k*ay(EV%lmaxg+5)-(EV%lmaxg+1)*cothxor*ay(EV%lmaxg+6) & -opacity*ay(EV%lmaxg+6) end if else ! Use tight-coupling approximation where moments are zero for l>1 pigdot=0._dl ayprime(8:EV%lmaxg+6)=0._dl end if ! Massless neutrino equations of motion. clxrdot=-k*(4._dl/3._dl*z+qr) ayprime(EV%lmaxg+7)=clxrdot qrdot=k*(clxr-2._dl*pir*EV%Kf(1))/3._dl ayprime(EV%lmaxg+8)=qrdot pirdot=k*(0.4_dl*qr-0.6_dl*EV%Kf(2)*ay(EV%lmaxg+10)+8._dl/15._dl*sigma) ayprime(EV%lmaxg+9)=pirdot do l=3,min(EV%FirstZerolForBeta,EV%lmaxnr)-1 ayprime(l+EV%lmaxg+7)=k*denl(l)*(l*ay(l+EV%lmaxg+6) -(l+1)*EV%Kf(l)*ay(l+EV%lmaxg+8)) end do ! Truncate the neutrino expansion if (EV%lmaxnr/=EV%FirstZerolForBeta) then ayprime(EV%lmaxnr+EV%lmaxg+7)=k*ay(EV%lmaxnr+EV%lmaxg+6)- & (EV%lmaxnr+1)*cothxor*ay(EV%lmaxnr+EV%lmaxg+7) end if ! Polarization if (tcp) then !l=2 (defined to be zero for l<2) ayprime(EV%polind+2) = -opacity*(ay(EV%polind+2) - polter) - k/3._dl*EV%Kf(2)*ay(EV%polind+3) !and the rest do l=3,min(EV%FirstZerolForBeta,EV%lmaxgpol)-1 ayprime(EV%polind+l)=-opacity*ay(EV%polind+l) + k*denl(l)*(l*ay(EV%polind+l-1) -& polfac(l)*EV%Kf(l)*ay(EV%polind+l+1)) end do !truncate if (EV%lmaxgpol/=EV%FirstZerolForBeta) then ayprime(EV%polind+EV%lmaxgpol)=-opacity*ay(EV%polind+EV%lmaxgpol) + & k*EV%poltruncfac*ay(EV%polind+EV%lmaxgpol-1)-(EV%lmaxgpol+3)*cothxor*ay(EV%polind+EV%lmaxgpol) end if else ayprime(EV%polind+2:EV%polind+EV%lmaxgpol)=0._dl end if ! Massive neutrino equations of motion. if (CP%Num_Nu_massive == 0) return end subroutine derivs !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine derivst(EV,n,tau,ayt,aytprime) ! Evaluate the time derivatives of the tensor perturbations. use ThermoData use MassiveNu implicit none type(EvolutionVars) EV integer n,l,i,ind real(dl), target :: ayt(n),aytprime(n) real(dl) ep,tau,grho,rhopi,cs2,opacity,gpres,pirdt,cothxor,grhormass_t logical tcp real(dl), dimension(:),pointer :: neut,neutprime,E,B,Eprime,Bprime real(dl) q,aq,v,akv(nqmax0),shearnu real(dl) grhob_t,grhor_t,grhoc_t,grhog_t,grhov_t,grhonu_t,gpnu_t,polter real(dl) Hchi,shear,aux, pig real(dl) k,k2,a,a2 real(dl) pinu,pir,pnu,adotoa, rhonu real(dl) grhoex_t,phi,phidot k2=EV%k2_buf k= EV%k_buf aux=EV%aux_buf !E and B start at l=2. Set up pointers accordingly to fill in ayt arrays E => ayt(EV%lmaxt+3-1:) Eprime=> aytprime(EV%lmaxt+3-1:) B => E(EV%lmaxpolt:) Bprime => Eprime(EV%lmaxpolt:) cothxor=1._dl/tanfunc(tau/CP%r)/CP%r a=ayt(1) ! Electric part of Weyl tensor and shear tensor ! elec=ayt(2) Hchi=ayt(2) shear=ayt(3) ! Photon anisotropic stress pig=ayt(4) a2=a*a ! Get sound speed and opacity, and see if should use tight-coupling call thermo(tau,cs2,opacity) if (k > 0.06_dl*epsw) then ep=ep0 else ep=0.2_dl*ep0 end if tcp = ((k/opacity > ep).or.(1._dl/(opacity*tau) > ep)) if (CP%Num_Nu_Massive == 0) then rhonu=1._dl pnu=1._dl/3._dl else ! call Nu_background(a*amnu,rhonu,pnu) end if ! Compute expansion rate from: grho=8*pi*rho*a**2 ! Also calculate gpres: 8*pi*p*a**2 grhob_t=grhob/a grhoc_t=grhoc/a grhor_t=grhornomass/a2 grhog_t=grhog/a2 grhov_t=grhov*a2 grhormass_t=0 grhonu_t=grhormass_t*rhonu gpnu_t=grhormass_t*pnu grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhonu_t gpres=(grhog_t+grhor_t)/3._dl +gpnu_t if (has_extra) then call Quint_ValsAta(a,phi,phidot) grhoex_t=phidot**2/2 + a2*Vofphi(phi,0) grho=grho+grhoex_t gpres=gpres-grhoex_t+phidot**2 else grho=grho +grhov_t gpres=gpres - grhov_t end if adotoa=sqrt((grho+grhok)/3._dl) aytprime(1)=adotoa*a polter = 0.1_dl*pig + 9._dl/15._dl*E(2) if (tcp) then ! Don't use tight coupling approx - use explicit equations: ! Equation for the photon anisotropic stress aytprime(4)=k*(-1._dl/3._dl*EV%Kft(2)*ayt(5)+8._dl/15._dl*shear) & -opacity*(pig - polter) do l=3,min(EV%FirstZerolForBeta,EV%lmaxt)-1 aytprime(l+2)=k*denl(l)*(l*ayt(l+1)- & tensfac(l)*EV%Kft(l)*ayt(l+3))-opacity*ayt(l+2) end do !Truncate the hierarchy if (EV%lmaxt/=EV%FirstZerolForBeta) then aytprime(EV%lmaxt+2)=k*EV%lmaxt/(EV%lmaxt-2._dl)*ayt(EV%lmaxt+1)- & (EV%lmaxt+3._dl)*cothxor*ayt(EV%lmaxt+2)-opacity*ayt(EV%lmaxt+2) end if !E and B-bar equations Eprime(2) = - opacity*(E(2) - polter) + k*(4._dl/6._dl*aux*B(2) - & 5._dl/27._dl*EV%Kft(2)*E(3)) do l=3,min(EV%FirstZerolForBeta,EV%lmaxpolt)-1 Eprime(l) =-opacity*E(l) + k*(denl(l)*(l*E(l-1) - & tensfacpol(l)*EV%Kft(l)*E(l+1)) + 4._dl/(l*(l+1))*aux*B(l)) end do !truncate: difficult, but zetting to zero seems to work OK Eprime(EV%lmaxpolt)=0._dl do l=2,min(EV%FirstZerolForBeta,EV%lmaxpolt)-1 Bprime(l) =-opacity*B(l) + k*(denl(l)*(l*B(l-1) - & tensfacpol(l)*EV%Kft(l)*B(l+1)) - 4._dl/(l*(l+1))*aux*E(l)) end do !truncate Bprime(EV%lmaxpolt)=0._dl else !Tight coupling ayt(4)=32._dl/45._dl*k/opacity*shear E(2)=ayt(4)/4._dl aytprime(4:n)=0._dl endif ! Neutrino equations: ! Anisotropic stress if (DoTensorNeutrinos) then neutprime => Bprime(EV%lmaxpolt:) neut => B(EV%lmaxpolt:) ! Massless neutrino anisotropic stress pir=neut(2) pirdt=k*(-1._dl/3._dl*EV%Kft(2)*neut(3)+ 8._dl/15._dl*shear) neutprime(2)=pirdt ! And for the moments do l=3,min(EV%FirstZerolForBeta,EV%lmaxnrt)-1 neutprime(l)=k*denl(l)*(l*neut(l-1)- tensfac(l)*EV%Kft(l)*neut(l+1)) end do ! Truncate the hierarchy if (EV%lmaxnrt/=EV%FirstZerolForBeta) then neutprime(EV%lmaxnrt)=k*EV%lmaxnrt/(EV%lmaxnrt-2._dl)*neut(EV%lmaxnrt-1)- & (EV%lmaxnrt+3._dl)*cothxor*neut(EV%lmaxnrt) endif end if ! Massive neutrino equations of motion. end subroutine derivst !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc end module GaugeInterface