There might exist some typos in this guide, please check it carefully with our original source codes. If you encounter any problems during your modifications, please contact us.
::::: subroutine output(EV,y, j,tau,sources) :::::
Add: After integer j
integer nu_j !EFTCAMB MOD
Add: After real(dl) ISW
! EFTCAMB MOD START: definitions of EFTCAMB quantities
! background quantities. Designer approach.
real(dl) EFT_H0,Hdot, Hdotdot, adotdota
! storage for EFT functions.
real(dl) EFTOmegaV, EFTOmegaP, EFTOmegaPP, EFTOmegaPPP
real(dl) EFTAlpha1V, EFTAlpha1P, EFTAlpha2V, EFTAlpha2P, EFTAlpha3V, EFTAlpha3P
real(dl) EFTAlpha4V, EFTAlpha4P, EFTAlpha4PP, EFTAlpha5V, EFTAlpha5P, EFTAlpha6V, EFTAlpha6P
! Background quantities.
real(dl) EFTc, EFTLambda, EFTcdot, EFTLambdadot
real(dl) EFTgrhoq, EFTgpresq, EFTgrhodotq, EFTgpresdotq
real(dl) EFT_grhonu, EFT_gpinu, EFT_grhonudot, EFT_gpinudot, grhormass_t
real(dl) EFT_grhonudot_tot, EFT_grhonu_tot, EFT_gpinu_tot, EFT_gpinudot_tot
real(dl) EFT_gpinudotdot, EFT_gpinudotdot_tot, EFT_gpipinudotdot
! perturbations quantities.
real(dl) EFTpiA, EFTpiB, EFTpiC, EFTpiD, EFTpiE ! pi field equations
real(dl) EFTeomF, EFTeomN, EFTeomX, EFTeomY, EFTeomG, EFTeomU ! equations of motion
real(dl) EFTeomL, EFTeomM, EFTeomV, EFTeomNdot, EFTeomVdot, EFTeomXdot
real(dl) pidotdot
! EFT quantities relevant for sources:
real(dl) EFTISW, EFTLensing, EFTsigmadot, polterdot
! mu and gamma
real(dl) EFT_Psi, EFT_Phi, EFT_function_mu, EFT_function_gamma, EFT_xi
! Debug quantities
real(dl) temp1, temp2, temp3,temp4, temp5, temp6, temp7
! EFTCAMB MOD END
! background quantities. Designer approach.
real(dl) EFT_H0,Hdot, Hdotdot, adotdota
! storage for EFT functions.
real(dl) EFTOmegaV, EFTOmegaP, EFTOmegaPP, EFTOmegaPPP
real(dl) EFTAlpha1V, EFTAlpha1P, EFTAlpha2V, EFTAlpha2P, EFTAlpha3V, EFTAlpha3P
real(dl) EFTAlpha4V, EFTAlpha4P, EFTAlpha4PP, EFTAlpha5V, EFTAlpha5P, EFTAlpha6V, EFTAlpha6P
! Background quantities.
real(dl) EFTc, EFTLambda, EFTcdot, EFTLambdadot
real(dl) EFTgrhoq, EFTgpresq, EFTgrhodotq, EFTgpresdotq
real(dl) EFT_grhonu, EFT_gpinu, EFT_grhonudot, EFT_gpinudot, grhormass_t
real(dl) EFT_grhonudot_tot, EFT_grhonu_tot, EFT_gpinu_tot, EFT_gpinudot_tot
real(dl) EFT_gpinudotdot, EFT_gpinudotdot_tot, EFT_gpipinudotdot
! perturbations quantities.
real(dl) EFTpiA, EFTpiB, EFTpiC, EFTpiD, EFTpiE ! pi field equations
real(dl) EFTeomF, EFTeomN, EFTeomX, EFTeomY, EFTeomG, EFTeomU ! equations of motion
real(dl) EFTeomL, EFTeomM, EFTeomV, EFTeomNdot, EFTeomVdot, EFTeomXdot
real(dl) pidotdot
! EFT quantities relevant for sources:
real(dl) EFTISW, EFTLensing, EFTsigmadot, polterdot
! mu and gamma
real(dl) EFT_Psi, EFT_Phi, EFT_function_mu, EFT_function_gamma, EFT_xi
! Debug quantities
real(dl) temp1, temp2, temp3,temp4, temp5, temp6, temp7
! EFTCAMB MOD END
Replace:
! EFTCAMB MOD START: Computation of DE background density.
if (CP%EFTflag==0) then
grhov_t=grhov*a**(-1-3*w_lam)
else if (CP%EFTflag/=0) then
grhov_t=grhov*EFTw(a,3)
EFT_H0 = (CP%h0/c_EFT)*1000._dl
end if
! Original code:
! grhov_t=grhov*a**(-1-3*w_lam)
! EFTCAMB MOD END
if (CP%EFTflag==0) then
grhov_t=grhov*a**(-1-3*w_lam)
else if (CP%EFTflag/=0) then
grhov_t=grhov*EFTw(a,3)
EFT_H0 = (CP%h0/c_EFT)*1000._dl
end if
! Original code:
! grhov_t=grhov*a**(-1-3*w_lam)
! EFTCAMB MOD END
Replace:
! EFTCAMB MOD START: computation of gpres.
if (CP%EFTflag==0) then ! Normal CAMB code
gpres=(grhog_t+grhor_t)/3+grhov_t*w_lam
else if (CP%EFTflag/=0) then
gpres=(grhog_t+grhor_t)/3._dl +EFTw(a,0)*grhov_t
end if
! Original code:
! gpres=(grhog_t+grhor_t)/3+grhov_t*w_lam
! EFTCAMB MOD END
if (CP%EFTflag==0) then ! Normal CAMB code
gpres=(grhog_t+grhor_t)/3+grhov_t*w_lam
else if (CP%EFTflag/=0) then
gpres=(grhog_t+grhor_t)/3._dl +EFTw(a,0)*grhov_t
end if
! Original code:
! gpres=(grhog_t+grhor_t)/3+grhov_t*w_lam
! EFTCAMB MOD END
Replace:
! EFTCAMB MOD START: initialization of DE perturbations.
if (CP%EFTflag==0) then
if (w_lam /= -1 .and. w_Perturb) then
clxq=y(EV%w_ix)
vq=y(EV%w_ix+1)
dgrho=dgrho + clxq*grhov_t
dgq = dgq + vq*grhov_t*(1+w_lam)
end if
else if (CP%EFTflag/=0.and.EV%EFTCAMBactive) then
clxq=y(EV%w_ix)
vq=y(EV%w_ix+1)
pidotdot = yprime(EV%w_ix+1)
end if
! Original code:
! if (w_lam /= -1 .and. w_Perturb) then
! clxq=y(EV%w_ix)
! vq=y(EV%w_ix+1)
! dgrho=dgrho + clxq*grhov_t
! dgq = dgq + vq*grhov_t*(1+w_lam)
! end if
! EFTCAMB MOD END
if (CP%EFTflag==0) then
if (w_lam /= -1 .and. w_Perturb) then
clxq=y(EV%w_ix)
vq=y(EV%w_ix+1)
dgrho=dgrho + clxq*grhov_t
dgq = dgq + vq*grhov_t*(1+w_lam)
end if
else if (CP%EFTflag/=0.and.EV%EFTCAMBactive) then
clxq=y(EV%w_ix)
vq=y(EV%w_ix+1)
pidotdot = yprime(EV%w_ix+1)
end if
! Original code:
! if (w_lam /= -1 .and. w_Perturb) then
! clxq=y(EV%w_ix)
! vq=y(EV%w_ix+1)
! dgrho=dgrho + clxq*grhov_t
! dgq = dgq + vq*grhov_t*(1+w_lam)
! end if
! EFTCAMB MOD END
Add: After adotoa=sqrt((grho+grhok)/3)
! EFTCAMB MOD START: Computation of background quantities:
! 1) FRW background quantities.
if (CP%EFTflag/=0.and.a>=EFTturnonpi) then
adotdota=(adotoa*adotoa-gpres)/2._dl
Hdot =adotdota-adotoa**2._dl
! Massive neutrinos:
EFT_gpinudot_tot = 0._dl
if (CP%Num_Nu_Massive /= 0) then
do nu_j = 1, CP%Nu_mass_eigenstates
EFT_grhonu = 0._dl
EFT_gpinu = 0._dl
EFT_grhonudot = 0._dl
EFT_gpinudot = 0._dl
EFT_gpinudotdot = 0._dl
grhormass_t=grhormass(nu_j)/a**2
call Nu_background(a*nu_masses(nu_j),EFT_grhonu,EFT_gpinu)
EFT_grhonu_tot = EFT_grhonu_tot + grhormass_t*EFT_grhonu
EFT_gpinu_tot = EFT_gpinu_tot + grhormass_t*EFT_gpinu
EFT_grhonudot_tot = EFT_grhonudot_tot + grhormass_t*(Nu_drho(a*nu_masses(nu_j),adotoa, EFT_grhonu)&
& -4._dl*adotoa*EFT_grhonu)
EFT_gpinudot = Nu_pidot(a*nu_masses(nu_j),adotoa, EFT_gpinu)
EFT_gpipinudotdot = Nu_pidotdot(a*nu_masses(nu_j),adotoa,Hdot,EFT_gpinu,EFT_gpinudot)
EFT_gpinudot_tot = EFT_gpinudot_tot + grhormass_t*(EFT_gpinudot -4._dl*adotoa*EFT_gpinu)
EFT_gpinudotdot_tot = EFT_gpinudotdot_tot + grhormass_t*(EFT_gpipinudotdot &
& -8._dl*adotoa*EFT_gpinudot +4._dl*EFT_gpinu*(+4._dl*adotoa**2-Hdot))
end do
end if
Hdotdot = 2._dl*adotoa*Hdot &
& + 0.5_dl*adotoa*(grhob_t + grhoc_t + 8._dl*(grhog_t+grhor_t)/3._dl)&
& + 0.5_dl*adotoa*grhov_t*((1._dl+EFTw(a,0))*(1._dl+3._dl*EFTw(a,0))-a*EFTw(a,1))&
& + adotoa/6._dl*EFT_grhonu_tot -0.5_dl*adotoa*EFT_gpinu_tot -0.5_dl*EFT_gpinudot_tot
end if
! 2) Compute and write once for all the values of the EFT functions.
if (CP%EFTflag/=0 .and. a>=EFTturnonpi*0.1_dl) then
EFTOmegaV = EFTOmega(a,0)
EFTOmegaP = EFTOmega(a,1)
EFTOmegaPP = EFTOmega(a,2)
EFTOmegaPPP = EFTOmega(a,3)
EFTAlpha1V = EFTAlpha1(a,0)
EFTAlpha1P = EFTAlpha1(a,1)
EFTAlpha2V = EFTAlpha2(a,0)
EFTAlpha2P = EFTAlpha2(a,1)
EFTAlpha3V = EFTAlpha3(a,0)
EFTAlpha3P = EFTAlpha3(a,1)
EFTAlpha4V = EFTAlpha4(a,0)
EFTAlpha4P = EFTAlpha4(a,1)
EFTAlpha4PP = EFTAlpha4(a,2)
EFTAlpha5V = EFTAlpha5(a,0)
EFTAlpha5P = EFTAlpha5(a,1)
EFTAlpha6V = EFTAlpha6(a,0)
EFTAlpha6P = EFTAlpha6(a,1)
end if
! 3) Computation of EFT background quantities.
if (CP%EFTflag/=0.and.a>=EFTturnonpi) then
!8*pi*G*c*a^2
EFTc = (adotoa*adotoa - Hdot)*(EFTOmegaV + a*EFTOmegaP*0.5_dl) &
& - 0.5_dl*a2*adotoa*adotoa*EFTOmegaPP&
& + 0.5_dl*grhov_t*(1._dl+EFTw(a,0))
!8*pi*G*Lambda*a^2
EFTLambda = +EFTw(a,0)*grhov_t &
&-EFTOmegaV*(2._dl*Hdot+adotoa**2._dl) &
&-a*EFTOmegaP*(2._dl*adotoa**2._dl + Hdot) &
&-a2*adotoa**2._dl*EFTOmegaPP
!EFT C DOT: 8*pi*G*cdot*a^2
EFTcdot = -EFTOmegaV*(Hdotdot-4._dl*adotoa*Hdot+2._dl*adotoa*adotoa*adotoa) &
& + 0.5_dl*a*EFTOmegaP*(-Hdotdot+adotoa*Hdot+adotoa*adotoa*adotoa)&
& +0.5_dl*a2*adotoa*EFTOmegaPP*(adotoa*adotoa-3._dl*Hdot)&
& -0.5_dl*a*a2*adotoa*adotoa*adotoa*EFTOmegaPPP&
& +0.5_dl*adotoa*grhov_t*(-3._dl*(1._dl+EFTw(a,0))**2 + a*EFTw(a,1))
!EFT LAMBDA DOT: 8*pi*G*Ldot*a^2
EFTLambdadot = -2._dl*EFTOmegaV*(Hdotdot-adotoa*Hdot-adotoa*adotoa*adotoa)&
& - a*EFTOmegaP*(4._dl*adotoa*Hdot+Hdotdot)&
& -a2*EFTOmegaPP*adotoa*(3._dl*Hdot+2._dl*adotoa*adotoa)&
& -a*a2*EFTOmegaPPP*adotoa*adotoa*adotoa&
& +grhov_t*adotoa*(a*EFTw(a,1)-3._dl*EFTw(a,0)*(1._dl+EFTw(a,0)))
if (CP%EFTflag==2.and.CP%DesignerEFTmodel/=2) then
EFTc = EFTcTemp(a,0)
EFTcdot = EFTcTemp(a,1)
EFTLambda = EFTLambdaTemp(a,0)
EFTLambdadot = EFTLambdaTemp(a,1)
end if
EFTgrhoq = 2._dl*EFTc -EFTLambda -3._dl*a*adotoa*adotoa*EFTOmegaP
EFTgpresq = EFTLambda + a2*adotoa*adotoa*EFTOmegaPP&
& +a*EFTOmegaP*(Hdot + 2._dl*adotoa*adotoa)
EFTgrhodotq = -3._dl*adotoa*(EFTgrhoq+EFTgpresq) + 3._dl*a*adotoa**3._dl*EFTOmegaP
EFTgpresdotq = EFTLambdadot +a2*a*adotoa**3*EFTOmegaPPP + 3._dl*a2*adotoa*Hdot*EFTOmegaPP&
& +a*EFTOmegaP*Hdotdot +3._dl*a*adotoa*Hdot*EFTOmegaP +2._dl*a2*adotoa**3*EFTOmegaPP&
& -2._dl*a*adotoa**3*EFTOmegaP
end if
! Debug code. Print every background EFT quantity.
if (a>= EFTturnonpi*0.1_dl.and.DebugEFTOutput) then
write (1,'(8e15.5)') a, tau, adotoa, Hdot, Hdotdot, grhov_t
write (2,'(25e15.5)') a, tau, EFTOmegaV, EFTOmegaP, EFTOmegaPP, EFTOmegaPPP,&
& EFTAlpha1V, EFTAlpha1P, EFTAlpha2V, EFTAlpha2P, EFTAlpha3V, EFTAlpha3P,&
& EFTAlpha4V, EFTAlpha4P, EFTAlpha4PP, EFTAlpha5V, EFTAlpha5P, EFTAlpha6V, EFTAlpha6P
write (3,'(12e15.6)') a, tau, EFTc, EFTLambda, EFTcdot, EFTLambdadot, EFTgrhoq,&
& EFTgpresq, EFTgrhodotq, EFTgpresdotq
end if
! EFTCAMB MOD END
! EFTCAMB MOD START: computation of Einstein equations factors.
if (CP%EFTflag/=0.and.a>=EFTturnonpi) then
EFTeomF = 1.5_dl/(k*(1._dl+EFTOmegaV))*((EFTgrhoQ+EFTgpresQ)*clxq& ! Background op
& +adotoa*adotoa*a*EFTOmegaP*clxq+a*adotoa*EFTOmegaP*vq)&
& +1.5_dl*a*EFT_H0*EFTAlpha2V**3*(vq+adotoa*clxq)/(k*(1._dl+EFTOmegaV))& ! Alpha2
& -1.5_dl*EFTAlpha3V**2/(1._dl+EFTOmegaV)*(k-3._dl*(Hdot-adotoa**2)/(k))*clxq& ! ALpha3
& +1.5_dl*EFTAlpha4V**2/(1._dl+EFTOmegaV)*k*clxq& ! Alpha4
& -1.5_dl*EFTAlpha4V**2/(1._dl+EFTOmegaV)*(Hdot-adotoa**2)/k*clxq
EFTeomN = -a*adotoa*EFTOmegaP/(1._dl+EFTOmegaV)*k*clxq& ! Background op
& +2._dl*adotoa/(1._dl+EFTOmegaV)*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P)*k*clxq& !Alpha4
& +EFTAlpha4V**2/(1._dl+EFTOmegaV)*k*vq&
& +2._dl*EFTAlpha5V**2*k*(vq+adotoa*clxq)/(1._dl+EFTOmegaV) ! Alpha5
EFTeomNdot = -a*Hdot*EFTOmegaP/(1._dl+EFTOmegaV)*k*clxq& ! Backgrpund op
& -a*adotoa*EFTOmegaP/(1._dl+EFTOmegaV)*k*vq&
& -a*adotoa**2/(1._dl+EFTOmegaV)*(EFTOmegaP+a*EFTOmegaPP&
& -a*EFTOmegaP**2/(1._dl+EFTOmegaV))*k*clxq&
& +EFTAlpha4V**2*k*Pidotdot/(1._dl+EFTOmegaV)+ a*adotoa*k*vq/(1._dl+EFTOmegaV)*&! Alpha4
&(2._dl*EFTAlpha4V*EFTAlpha4P-EFTAlpha4V**2*EFTOmegaP/(1._dl+EFTOmegaV))&
& +2._dl*k/(1._dl+EFTOmegaV)*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P)*(Hdot*clxq+adotoa*vq)&
& +2._dl*a*adotoa**2*k*clxq/(1._dl+EFTOmegaV)*(a*EFTAlpha4V*EFTAlpha4PP+a*EFTAlpha4P**2&
& +3._dl*EFTAlpha4V*EFTAlpha4P-EFTOmegaP/(1._dl+EFTOmegaV)*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P))&
& +2._dl*EFTAlpha5V**2*k/(1._dl+EFTOmegaV)*(Pidotdot + adotoa*vq + Hdot*clxq)& ! Alpha5
& +2._dl*a*k*adotoa/(1._dl+EFTOmegaV)*(vq+adotoa*clxq)*(2._dl*EFTAlpha5V*EFTAlpha5P&
& -EFTAlpha5V**2*EFTOmegaP/(1._dl+EFTOmegaV))
EFTeomX = 1._dl& ! Backgrond op
& -EFTAlpha4V**2/(1._dl+EFTOmegaV) ! Alpha4
EFTeomXdot = -a*adotoa/(1._dl+EFTOmegaV)*(2._dl*EFTAlpha4V*EFTAlpha4P& ! Alpha4
& -EFTAlpha4V**2*EFTOmegaP/(1._dl+EFTOmegaV))
EFTeomY = +0.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)& ! Background op
& +1.5_dl/(1._dl+EFTOmegaV)*(EFTAlpha3V**2+a*EFTAlpha3V*EFTAlpha3P)& ! Alpha3
& +0.5_dl*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P)/(1._dl+EFTOmegaV) ! Alpha4
EFTeomG = +1._dl + 0.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)& ! Background op
& +0.5_dl*a*EFT_H0*EFTAlpha2V**3/adotoa/(1._dl+EFTOmegaV)& ! Alpha2
& +1.5_dl*EFTAlpha3V**2/(1._dl+EFTOmegaV)& ! Alpha3
& +EFTAlpha4V**2/(1._dl+EFTOmegaV) ! Alpha4
EFTeomU = 1._dl& ! Background op
& +1.5_dl*EFTAlpha3V**2/(1._dl+EFTOmegaV)& ! Alpha3
& +0.5_dl*EFTAlpha4V**2/(1._dl+EFTOmegaV) ! Alpha4
EFTeomL = -1.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)*(3._dl*adotoa*adotoa-Hdot)*clxq& ! Background op
& -1.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)*adotoa*vq&
& -0.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)*k2*clxq&
& +0.5_dl*clxq/(adotoa*(1._dl+EFTOmegaV))*EFTgrhodotQ&
& +(vq+adotoa*clxq)/(adotoa*(1+EFTOmegaV))*EFTc&
& +2._dl*a2*EFT_H0**2*EFTAlpha1V**4/adotoa/(1._dl+EFTOmegaV)*(vq+adotoa*clxq)& ! Alpha1
& +1.5_dl*a*EFT_H0*EFTAlpha2V**3/(1._dl+EFTOmegaV)*& ! Alpha2
&((Hdot/adotoa-2._dl*adotoa-k2/(3._dl*adotoa))*clxq-vq)&
& -1.5_dl*EFTAlpha3V**2/(1._dl+EFTOmegaV)*(k2-3._dl*(Hdot-adotoa**2))*clxq& ! Alpha3
& +3._dl*EFTAlpha4V**2/(1._dl+EFTOmegaV)*(Hdot-adotoa**2-k2/3._dl)*clxq& ! Alpha4
& +4._dl*EFTAlpha6V**2*k2/adotoa*(vq+adotoa*clxq)/(1._dl+EFTOmegaV) ! Alpha6
EFTeomM = EFTgpresdotQ*clxq + (EFTgrhoQ+ EFTgpresQ)*(vq+adotoa*clxq)& ! Background op
& +a2*adotoa**2*EFTOmegaPP*vq +a2*adotoa**3*EFTOmegaPP*clxq&
& +a*adotoa*EFTOmegaP*(pidotdot + (Hdot/adotoa+4._dl*adotoa)*vq&
& + (2._dl*Hdot+6._dl*adotoa**2 +2._dl/3._dl*k2)*clxq)&
& +a*EFT_H0*EFTAlpha2V**2*(EFTAlpha2V*pidotdot& ! Alpha2
& +(4._dl*EFTAlpha2V+3._dl*a*EFTAlpha2P)*adotoa*vq +(3._dl*adotoa**2*EFTAlpha2V&
& +Hdot*EFTAlpha2V +3._dl*a*adotoa**2*EFTAlpha2P)*clxq)&
& +EFTAlpha3V**2*(3._dl*adotoa**2-3._dl*Hdot+k2)*vq& ! Alpha3
& +EFTAlpha3V**2*(6._dl*adotoa**3-3._dl*Hdotdot)*clxq&
& +2._dl*adotoa*k2*clxq*(EFTAlpha3V**2+a*EFTAlpha3V*EFTAlpha3P)&
& -6._dl*a*adotoa*(Hdot-adotoa**2)*EFTAlpha3V*EFTAlpha3P*clxq&
& -EFTAlpha4V**2*(Hdot-adotoa**2-k2/3._dl)*vq& ! Alpha4
& -2._dl*adotoa*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P)*(Hdot-adotoa**2-k2/3._dl)*clxq&
& -EFTAlpha4V**2*(Hdotdot-2*adotoa*Hdot)*clxq&
& -4._dl*EFTAlpha5V**2*k2*(vq+adotoa*clxq)/3._dl ! Alpha5
EFTeomV = +0.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)& ! Background op
& -(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P)/(1._dl+EFTOmegaV) ! Alpha4
EFTeomVdot = +0.5_dl*a*adotoa/(1._dl+EFTOmegaV)*(EFTOmegaP+a*EFTOmegaPP& ! Background op
& -a*EFTOmegaP**2/(1._dl+EFTOmegaV))&
& -a*adotoa/(1._dl+EFTOmegaV)*(a*EFTAlpha4V*EFTAlpha4PP+a*EFTAlpha4P**2& ! Alpha4
& +3._dl*EFTAlpha4V*EFTAlpha4P-EFTOmegaP/(1._dl+EFTOmegaV)*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P))
end if
! Debug code. Prints Einstein equations factors.
if (a>= EFTturnonpi .and. DebugEFTOutput) then
write (7,'(15E15.5)') a, tau, EFTeomF, EFTeomN, EFTeomX, EFTeomY, EFTeomG, EFTeomU,&
& EFTeomL, EFTeomM, EFTeomNdot, EFTeomVdot, EFTeomXdot
end if
! EFTCAMB MOD END
! 1) FRW background quantities.
if (CP%EFTflag/=0.and.a>=EFTturnonpi) then
adotdota=(adotoa*adotoa-gpres)/2._dl
Hdot =adotdota-adotoa**2._dl
! Massive neutrinos:
EFT_gpinudot_tot = 0._dl
if (CP%Num_Nu_Massive /= 0) then
do nu_j = 1, CP%Nu_mass_eigenstates
EFT_grhonu = 0._dl
EFT_gpinu = 0._dl
EFT_grhonudot = 0._dl
EFT_gpinudot = 0._dl
EFT_gpinudotdot = 0._dl
grhormass_t=grhormass(nu_j)/a**2
call Nu_background(a*nu_masses(nu_j),EFT_grhonu,EFT_gpinu)
EFT_grhonu_tot = EFT_grhonu_tot + grhormass_t*EFT_grhonu
EFT_gpinu_tot = EFT_gpinu_tot + grhormass_t*EFT_gpinu
EFT_grhonudot_tot = EFT_grhonudot_tot + grhormass_t*(Nu_drho(a*nu_masses(nu_j),adotoa, EFT_grhonu)&
& -4._dl*adotoa*EFT_grhonu)
EFT_gpinudot = Nu_pidot(a*nu_masses(nu_j),adotoa, EFT_gpinu)
EFT_gpipinudotdot = Nu_pidotdot(a*nu_masses(nu_j),adotoa,Hdot,EFT_gpinu,EFT_gpinudot)
EFT_gpinudot_tot = EFT_gpinudot_tot + grhormass_t*(EFT_gpinudot -4._dl*adotoa*EFT_gpinu)
EFT_gpinudotdot_tot = EFT_gpinudotdot_tot + grhormass_t*(EFT_gpipinudotdot &
& -8._dl*adotoa*EFT_gpinudot +4._dl*EFT_gpinu*(+4._dl*adotoa**2-Hdot))
end do
end if
Hdotdot = 2._dl*adotoa*Hdot &
& + 0.5_dl*adotoa*(grhob_t + grhoc_t + 8._dl*(grhog_t+grhor_t)/3._dl)&
& + 0.5_dl*adotoa*grhov_t*((1._dl+EFTw(a,0))*(1._dl+3._dl*EFTw(a,0))-a*EFTw(a,1))&
& + adotoa/6._dl*EFT_grhonu_tot -0.5_dl*adotoa*EFT_gpinu_tot -0.5_dl*EFT_gpinudot_tot
end if
! 2) Compute and write once for all the values of the EFT functions.
if (CP%EFTflag/=0 .and. a>=EFTturnonpi*0.1_dl) then
EFTOmegaV = EFTOmega(a,0)
EFTOmegaP = EFTOmega(a,1)
EFTOmegaPP = EFTOmega(a,2)
EFTOmegaPPP = EFTOmega(a,3)
EFTAlpha1V = EFTAlpha1(a,0)
EFTAlpha1P = EFTAlpha1(a,1)
EFTAlpha2V = EFTAlpha2(a,0)
EFTAlpha2P = EFTAlpha2(a,1)
EFTAlpha3V = EFTAlpha3(a,0)
EFTAlpha3P = EFTAlpha3(a,1)
EFTAlpha4V = EFTAlpha4(a,0)
EFTAlpha4P = EFTAlpha4(a,1)
EFTAlpha4PP = EFTAlpha4(a,2)
EFTAlpha5V = EFTAlpha5(a,0)
EFTAlpha5P = EFTAlpha5(a,1)
EFTAlpha6V = EFTAlpha6(a,0)
EFTAlpha6P = EFTAlpha6(a,1)
end if
! 3) Computation of EFT background quantities.
if (CP%EFTflag/=0.and.a>=EFTturnonpi) then
!8*pi*G*c*a^2
EFTc = (adotoa*adotoa - Hdot)*(EFTOmegaV + a*EFTOmegaP*0.5_dl) &
& - 0.5_dl*a2*adotoa*adotoa*EFTOmegaPP&
& + 0.5_dl*grhov_t*(1._dl+EFTw(a,0))
!8*pi*G*Lambda*a^2
EFTLambda = +EFTw(a,0)*grhov_t &
&-EFTOmegaV*(2._dl*Hdot+adotoa**2._dl) &
&-a*EFTOmegaP*(2._dl*adotoa**2._dl + Hdot) &
&-a2*adotoa**2._dl*EFTOmegaPP
!EFT C DOT: 8*pi*G*cdot*a^2
EFTcdot = -EFTOmegaV*(Hdotdot-4._dl*adotoa*Hdot+2._dl*adotoa*adotoa*adotoa) &
& + 0.5_dl*a*EFTOmegaP*(-Hdotdot+adotoa*Hdot+adotoa*adotoa*adotoa)&
& +0.5_dl*a2*adotoa*EFTOmegaPP*(adotoa*adotoa-3._dl*Hdot)&
& -0.5_dl*a*a2*adotoa*adotoa*adotoa*EFTOmegaPPP&
& +0.5_dl*adotoa*grhov_t*(-3._dl*(1._dl+EFTw(a,0))**2 + a*EFTw(a,1))
!EFT LAMBDA DOT: 8*pi*G*Ldot*a^2
EFTLambdadot = -2._dl*EFTOmegaV*(Hdotdot-adotoa*Hdot-adotoa*adotoa*adotoa)&
& - a*EFTOmegaP*(4._dl*adotoa*Hdot+Hdotdot)&
& -a2*EFTOmegaPP*adotoa*(3._dl*Hdot+2._dl*adotoa*adotoa)&
& -a*a2*EFTOmegaPPP*adotoa*adotoa*adotoa&
& +grhov_t*adotoa*(a*EFTw(a,1)-3._dl*EFTw(a,0)*(1._dl+EFTw(a,0)))
if (CP%EFTflag==2.and.CP%DesignerEFTmodel/=2) then
EFTc = EFTcTemp(a,0)
EFTcdot = EFTcTemp(a,1)
EFTLambda = EFTLambdaTemp(a,0)
EFTLambdadot = EFTLambdaTemp(a,1)
end if
EFTgrhoq = 2._dl*EFTc -EFTLambda -3._dl*a*adotoa*adotoa*EFTOmegaP
EFTgpresq = EFTLambda + a2*adotoa*adotoa*EFTOmegaPP&
& +a*EFTOmegaP*(Hdot + 2._dl*adotoa*adotoa)
EFTgrhodotq = -3._dl*adotoa*(EFTgrhoq+EFTgpresq) + 3._dl*a*adotoa**3._dl*EFTOmegaP
EFTgpresdotq = EFTLambdadot +a2*a*adotoa**3*EFTOmegaPPP + 3._dl*a2*adotoa*Hdot*EFTOmegaPP&
& +a*EFTOmegaP*Hdotdot +3._dl*a*adotoa*Hdot*EFTOmegaP +2._dl*a2*adotoa**3*EFTOmegaPP&
& -2._dl*a*adotoa**3*EFTOmegaP
end if
! Debug code. Print every background EFT quantity.
if (a>= EFTturnonpi*0.1_dl.and.DebugEFTOutput) then
write (1,'(8e15.5)') a, tau, adotoa, Hdot, Hdotdot, grhov_t
write (2,'(25e15.5)') a, tau, EFTOmegaV, EFTOmegaP, EFTOmegaPP, EFTOmegaPPP,&
& EFTAlpha1V, EFTAlpha1P, EFTAlpha2V, EFTAlpha2P, EFTAlpha3V, EFTAlpha3P,&
& EFTAlpha4V, EFTAlpha4P, EFTAlpha4PP, EFTAlpha5V, EFTAlpha5P, EFTAlpha6V, EFTAlpha6P
write (3,'(12e15.6)') a, tau, EFTc, EFTLambda, EFTcdot, EFTLambdadot, EFTgrhoq,&
& EFTgpresq, EFTgrhodotq, EFTgpresdotq
end if
! EFTCAMB MOD END
! EFTCAMB MOD START: computation of Einstein equations factors.
if (CP%EFTflag/=0.and.a>=EFTturnonpi) then
EFTeomF = 1.5_dl/(k*(1._dl+EFTOmegaV))*((EFTgrhoQ+EFTgpresQ)*clxq& ! Background op
& +adotoa*adotoa*a*EFTOmegaP*clxq+a*adotoa*EFTOmegaP*vq)&
& +1.5_dl*a*EFT_H0*EFTAlpha2V**3*(vq+adotoa*clxq)/(k*(1._dl+EFTOmegaV))& ! Alpha2
& -1.5_dl*EFTAlpha3V**2/(1._dl+EFTOmegaV)*(k-3._dl*(Hdot-adotoa**2)/(k))*clxq& ! ALpha3
& +1.5_dl*EFTAlpha4V**2/(1._dl+EFTOmegaV)*k*clxq& ! Alpha4
& -1.5_dl*EFTAlpha4V**2/(1._dl+EFTOmegaV)*(Hdot-adotoa**2)/k*clxq
EFTeomN = -a*adotoa*EFTOmegaP/(1._dl+EFTOmegaV)*k*clxq& ! Background op
& +2._dl*adotoa/(1._dl+EFTOmegaV)*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P)*k*clxq& !Alpha4
& +EFTAlpha4V**2/(1._dl+EFTOmegaV)*k*vq&
& +2._dl*EFTAlpha5V**2*k*(vq+adotoa*clxq)/(1._dl+EFTOmegaV) ! Alpha5
EFTeomNdot = -a*Hdot*EFTOmegaP/(1._dl+EFTOmegaV)*k*clxq& ! Backgrpund op
& -a*adotoa*EFTOmegaP/(1._dl+EFTOmegaV)*k*vq&
& -a*adotoa**2/(1._dl+EFTOmegaV)*(EFTOmegaP+a*EFTOmegaPP&
& -a*EFTOmegaP**2/(1._dl+EFTOmegaV))*k*clxq&
& +EFTAlpha4V**2*k*Pidotdot/(1._dl+EFTOmegaV)+ a*adotoa*k*vq/(1._dl+EFTOmegaV)*&! Alpha4
&(2._dl*EFTAlpha4V*EFTAlpha4P-EFTAlpha4V**2*EFTOmegaP/(1._dl+EFTOmegaV))&
& +2._dl*k/(1._dl+EFTOmegaV)*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P)*(Hdot*clxq+adotoa*vq)&
& +2._dl*a*adotoa**2*k*clxq/(1._dl+EFTOmegaV)*(a*EFTAlpha4V*EFTAlpha4PP+a*EFTAlpha4P**2&
& +3._dl*EFTAlpha4V*EFTAlpha4P-EFTOmegaP/(1._dl+EFTOmegaV)*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P))&
& +2._dl*EFTAlpha5V**2*k/(1._dl+EFTOmegaV)*(Pidotdot + adotoa*vq + Hdot*clxq)& ! Alpha5
& +2._dl*a*k*adotoa/(1._dl+EFTOmegaV)*(vq+adotoa*clxq)*(2._dl*EFTAlpha5V*EFTAlpha5P&
& -EFTAlpha5V**2*EFTOmegaP/(1._dl+EFTOmegaV))
EFTeomX = 1._dl& ! Backgrond op
& -EFTAlpha4V**2/(1._dl+EFTOmegaV) ! Alpha4
EFTeomXdot = -a*adotoa/(1._dl+EFTOmegaV)*(2._dl*EFTAlpha4V*EFTAlpha4P& ! Alpha4
& -EFTAlpha4V**2*EFTOmegaP/(1._dl+EFTOmegaV))
EFTeomY = +0.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)& ! Background op
& +1.5_dl/(1._dl+EFTOmegaV)*(EFTAlpha3V**2+a*EFTAlpha3V*EFTAlpha3P)& ! Alpha3
& +0.5_dl*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P)/(1._dl+EFTOmegaV) ! Alpha4
EFTeomG = +1._dl + 0.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)& ! Background op
& +0.5_dl*a*EFT_H0*EFTAlpha2V**3/adotoa/(1._dl+EFTOmegaV)& ! Alpha2
& +1.5_dl*EFTAlpha3V**2/(1._dl+EFTOmegaV)& ! Alpha3
& +EFTAlpha4V**2/(1._dl+EFTOmegaV) ! Alpha4
EFTeomU = 1._dl& ! Background op
& +1.5_dl*EFTAlpha3V**2/(1._dl+EFTOmegaV)& ! Alpha3
& +0.5_dl*EFTAlpha4V**2/(1._dl+EFTOmegaV) ! Alpha4
EFTeomL = -1.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)*(3._dl*adotoa*adotoa-Hdot)*clxq& ! Background op
& -1.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)*adotoa*vq&
& -0.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)*k2*clxq&
& +0.5_dl*clxq/(adotoa*(1._dl+EFTOmegaV))*EFTgrhodotQ&
& +(vq+adotoa*clxq)/(adotoa*(1+EFTOmegaV))*EFTc&
& +2._dl*a2*EFT_H0**2*EFTAlpha1V**4/adotoa/(1._dl+EFTOmegaV)*(vq+adotoa*clxq)& ! Alpha1
& +1.5_dl*a*EFT_H0*EFTAlpha2V**3/(1._dl+EFTOmegaV)*& ! Alpha2
&((Hdot/adotoa-2._dl*adotoa-k2/(3._dl*adotoa))*clxq-vq)&
& -1.5_dl*EFTAlpha3V**2/(1._dl+EFTOmegaV)*(k2-3._dl*(Hdot-adotoa**2))*clxq& ! Alpha3
& +3._dl*EFTAlpha4V**2/(1._dl+EFTOmegaV)*(Hdot-adotoa**2-k2/3._dl)*clxq& ! Alpha4
& +4._dl*EFTAlpha6V**2*k2/adotoa*(vq+adotoa*clxq)/(1._dl+EFTOmegaV) ! Alpha6
EFTeomM = EFTgpresdotQ*clxq + (EFTgrhoQ+ EFTgpresQ)*(vq+adotoa*clxq)& ! Background op
& +a2*adotoa**2*EFTOmegaPP*vq +a2*adotoa**3*EFTOmegaPP*clxq&
& +a*adotoa*EFTOmegaP*(pidotdot + (Hdot/adotoa+4._dl*adotoa)*vq&
& + (2._dl*Hdot+6._dl*adotoa**2 +2._dl/3._dl*k2)*clxq)&
& +a*EFT_H0*EFTAlpha2V**2*(EFTAlpha2V*pidotdot& ! Alpha2
& +(4._dl*EFTAlpha2V+3._dl*a*EFTAlpha2P)*adotoa*vq +(3._dl*adotoa**2*EFTAlpha2V&
& +Hdot*EFTAlpha2V +3._dl*a*adotoa**2*EFTAlpha2P)*clxq)&
& +EFTAlpha3V**2*(3._dl*adotoa**2-3._dl*Hdot+k2)*vq& ! Alpha3
& +EFTAlpha3V**2*(6._dl*adotoa**3-3._dl*Hdotdot)*clxq&
& +2._dl*adotoa*k2*clxq*(EFTAlpha3V**2+a*EFTAlpha3V*EFTAlpha3P)&
& -6._dl*a*adotoa*(Hdot-adotoa**2)*EFTAlpha3V*EFTAlpha3P*clxq&
& -EFTAlpha4V**2*(Hdot-adotoa**2-k2/3._dl)*vq& ! Alpha4
& -2._dl*adotoa*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P)*(Hdot-adotoa**2-k2/3._dl)*clxq&
& -EFTAlpha4V**2*(Hdotdot-2*adotoa*Hdot)*clxq&
& -4._dl*EFTAlpha5V**2*k2*(vq+adotoa*clxq)/3._dl ! Alpha5
EFTeomV = +0.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)& ! Background op
& -(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P)/(1._dl+EFTOmegaV) ! Alpha4
EFTeomVdot = +0.5_dl*a*adotoa/(1._dl+EFTOmegaV)*(EFTOmegaP+a*EFTOmegaPP& ! Background op
& -a*EFTOmegaP**2/(1._dl+EFTOmegaV))&
& -a*adotoa/(1._dl+EFTOmegaV)*(a*EFTAlpha4V*EFTAlpha4PP+a*EFTAlpha4P**2& ! Alpha4
& +3._dl*EFTAlpha4V*EFTAlpha4P-EFTOmegaP/(1._dl+EFTOmegaV)*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P))
end if
! Debug code. Prints Einstein equations factors.
if (a>= EFTturnonpi .and. DebugEFTOutput) then
write (7,'(15E15.5)') a, tau, EFTeomF, EFTeomN, EFTeomX, EFTeomY, EFTeomG, EFTeomU,&
& EFTeomL, EFTeomM, EFTeomNdot, EFTeomVdot, EFTeomXdot
end if
! EFTCAMB MOD END
Replace:
! EFTCAMB MOD START: compute z,dz before loading radiation and photon
if (CP%EFTflag==0.or.a
z=(0.5_dl*dgrho/k + etak)/adotoa
dz= -adotoa*z - 0.5_dl*dgrho/k
else if (CP%EFTflag/=0.and.a>=EFTturnonpi) then
! RSA approximation.
z= 1._dl/EFTeomG*(etak/adotoa + 0.5_dl*dgrho/(k*adotoa*(1._dl+EFTOmegaV)) + EFTeomL/(k*EFT_H0))
dz= 1._dl/EFTeomU*(-2._dl*adotoa*z*(1._dl+EFTeomY-0.5_dl*EFTeomG) &
& - 0.5_dl*dgrho/k/(1._dl+EFTOmegaV) -adotoa*EFTeomL/(k*EFT_H0) -1.5_dl/k/(1._dl+EFTOmegaV)*EFTeomM/EFT_H0)
end if
!equations of motion 1.
if (EV%no_nu_multpoles) then
clxr=-4*dz/k
qr=-4._dl/3*z
pir=0
pirdot=0
else
clxr=y(EV%r_ix)
qr =y(EV%r_ix+1)
pir =y(EV%r_ix+2)
pirdot=yprime(EV%r_ix+2)
end if
! Original code:
! if (EV%no_nu_multpoles) then
! z=(0.5_dl*dgrho/k + etak)/adotoa
! dz= -adotoa*z - 0.5_dl*dgrho/k
! clxr=-4*dz/k
! qr=-4._dl/3*z
! pir=0
! pirdot=0
! else
! clxr=y(EV%r_ix)
! qr =y(EV%r_ix+1)
! pir =y(EV%r_ix+2)
! pirdot=yprime(EV%r_ix+2)
! end if
! EFTCAMB MOD END
if (CP%EFTflag==0.or.a
dz= -adotoa*z - 0.5_dl*dgrho/k
else if (CP%EFTflag/=0.and.a>=EFTturnonpi) then
! RSA approximation.
z= 1._dl/EFTeomG*(etak/adotoa + 0.5_dl*dgrho/(k*adotoa*(1._dl+EFTOmegaV)) + EFTeomL/(k*EFT_H0))
dz= 1._dl/EFTeomU*(-2._dl*adotoa*z*(1._dl+EFTeomY-0.5_dl*EFTeomG) &
& - 0.5_dl*dgrho/k/(1._dl+EFTOmegaV) -adotoa*EFTeomL/(k*EFT_H0) -1.5_dl/k/(1._dl+EFTOmegaV)*EFTeomM/EFT_H0)
end if
!equations of motion 1.
if (EV%no_nu_multpoles) then
clxr=-4*dz/k
qr=-4._dl/3*z
pir=0
pirdot=0
else
clxr=y(EV%r_ix)
qr =y(EV%r_ix+1)
pir =y(EV%r_ix+2)
pirdot=yprime(EV%r_ix+2)
end if
! if (EV%no_nu_multpoles) then
! z=(0.5_dl*dgrho/k + etak)/adotoa
! dz= -adotoa*z - 0.5_dl*dgrho/k
! clxr=-4*dz/k
! qr=-4._dl/3*z
! pir=0
! pirdot=0
! else
! clxr=y(EV%r_ix)
! qr =y(EV%r_ix+1)
! pir =y(EV%r_ix+2)
! pirdot=yprime(EV%r_ix+2)
! end if
! EFTCAMB MOD END
Comment off: After if (EV%no_phot_multpoles) then
! z=(0.5_dl*dgrho/k + etak)/adotoa !EFTCAMB MOD
! dz= -adotoa*z - 0.5_dl*dgrho/k !EFTCAMB MOD
! dz= -adotoa*z - 0.5_dl*dgrho/k !EFTCAMB MOD
Replace:
! EFTCAMB MOD START: equation of motion 3.
! Get sigma (shear) and z from the constraints
! have to get z from eta for numerical stability
if (CP%EFTflag==0 .or. a
z=(0.5_dl*dgrho/k + etak)/adotoa
sigma=(z+1.5_dl*dgq/k2)/EV%Kf(1)
else if (CP%EFTflag/=0 .and. a>=EFTturnonpi) then
z= 1._dl/EFTeomG*(etak/adotoa + 0.5_dl*dgrho/(k*adotoa*(1._dl+EFTOmegaV)) + EFTeomL/(k*EFT_H0))
dz=1._dl/EFTeomU*(-2*adotoa*(1._dl+EFTeomY)*z +etak -0.5_dl/k/(1._dl+EFTOmegaV)*(grhog_t*clxg+grhor_t*clxr)&
& -1.5_dl/k/(1._dl+EFTOmegaV)*EFTeomM/EFT_H0)
sigma= 1._dl/EFTeomX*(z*EFTeomU +1.5_dl*dgq/(k2*(1._dl+EFTOmegaV)) +EFTeomF/EFT_H0)
end if
! Original code:
! z=(0.5_dl*dgrho/k + etak)/adotoa
! sigma=(z+1.5_dl*dgq/k2)/EV%Kf(1)
! EFTCAMB MOD END
! Get sigma (shear) and z from the constraints
! have to get z from eta for numerical stability
if (CP%EFTflag==0 .or. a
sigma=(z+1.5_dl*dgq/k2)/EV%Kf(1)
else if (CP%EFTflag/=0 .and. a>=EFTturnonpi) then
z= 1._dl/EFTeomG*(etak/adotoa + 0.5_dl*dgrho/(k*adotoa*(1._dl+EFTOmegaV)) + EFTeomL/(k*EFT_H0))
dz=1._dl/EFTeomU*(-2*adotoa*(1._dl+EFTeomY)*z +etak -0.5_dl/k/(1._dl+EFTOmegaV)*(grhog_t*clxg+grhor_t*clxr)&
& -1.5_dl/k/(1._dl+EFTOmegaV)*EFTeomM/EFT_H0)
sigma= 1._dl/EFTeomX*(z*EFTeomU +1.5_dl*dgq/(k2*(1._dl+EFTOmegaV)) +EFTeomF/EFT_H0)
end if
! z=(0.5_dl*dgrho/k + etak)/adotoa
! sigma=(z+1.5_dl*dgq/k2)/EV%Kf(1)
! EFTCAMB MOD END
Add: Before !Maple's fortran output - see scal_eqs.map
! EFTCAMB MOD START: computation of quantities used to construct the sources.
if (CP%EFTflag/=0 .and. a>=EFTturnonpi) then
EFTsigmadot= 1._dl/EFTeomX*(-2._dl*adotoa*(1._dl+EFTeomV)*sigma +etak&
& -1._dl/k*dgpi/(1._dl+EFTOmegaV) +EFTeomN/EFT_H0)
EFTISW = 1._dl/EFTeomX*(-2._dl*(1._dl+EFTeomV)*(Hdot*sigma+ adotoa*EFTsigmadot)&
& -2._dl*adotoa*sigma*EFTeomVdot + 0.5_dl*dgq*(1._dl+EFTeomX)/(1._dl+EFTOmegaV)&
& +a*adotoa*EFTOmegaP/k*dgpi/(1._dl+EFTOmegaV)**2&
& -1._dl/(1._dl+EFTOmegaV)*(2._dl*adotoa*dgpi/k+diff_rhopi/k)&
& +(1._dl+EFTeomX)*k2/3._dl/EFT_H0*EFTeomF-EFTeomXdot*EFTsigmadot+EFTeomNdot/EFT_H0)
EFTLensing = 1._dl/EFTeomX*(-2._dl*adotoa*(1._dl+EFTeomV)*sigma +(1._dl+EFTeomX)*etak&
& -1._dl/k*dgpi/(1._dl+EFTOmegaV) +EFTeomN/EFT_H0)
end if
! EFTCAMB MOD END
if (CP%EFTflag/=0 .and. a>=EFTturnonpi) then
EFTsigmadot= 1._dl/EFTeomX*(-2._dl*adotoa*(1._dl+EFTeomV)*sigma +etak&
& -1._dl/k*dgpi/(1._dl+EFTOmegaV) +EFTeomN/EFT_H0)
EFTISW = 1._dl/EFTeomX*(-2._dl*(1._dl+EFTeomV)*(Hdot*sigma+ adotoa*EFTsigmadot)&
& -2._dl*adotoa*sigma*EFTeomVdot + 0.5_dl*dgq*(1._dl+EFTeomX)/(1._dl+EFTOmegaV)&
& +a*adotoa*EFTOmegaP/k*dgpi/(1._dl+EFTOmegaV)**2&
& -1._dl/(1._dl+EFTOmegaV)*(2._dl*adotoa*dgpi/k+diff_rhopi/k)&
& +(1._dl+EFTeomX)*k2/3._dl/EFT_H0*EFTeomF-EFTeomXdot*EFTsigmadot+EFTeomNdot/EFT_H0)
EFTLensing = 1._dl/EFTeomX*(-2._dl*adotoa*(1._dl+EFTeomV)*sigma +(1._dl+EFTeomX)*etak&
& -1._dl/k*dgpi/(1._dl+EFTOmegaV) +EFTeomN/EFT_H0)
end if
! EFTCAMB MOD END
Replace:
! EFTCAMB MOD START: TT source
if (CP%EFTflag==0.or.a < EFTturnonpi) then ! Normal CAMB code.
!Maple's fortran output - see scal_eqs.map
!2phi' term (\phi' + \psi' in Newtonian gauge)
ISW = (4.D0/3.D0*k*EV%Kf(1)*sigma+(-2.D0/3.D0*sigma-2.D0/3.D0*etak/adotoa)*k &
-diff_rhopi/k**2-1.D0/adotoa*dgrho/3.D0+(3.D0*gpres+5.D0*grho)*sigma/k/3.D0 &
-2.D0/k*adotoa/EV%Kf(1)*etak)*expmmu(j)
!e.g. to get only late-time ISW
! if (1/a-1 < 30) ISW=0
!The rest, note y(9)->octg, yprime(9)->octgprime (octopoles)
sources(1)= ISW + ((-9.D0/160.D0*pig-27.D0/80.D0*ypol(2))/k**2*opac(j)+(11.D0/10.D0*sigma- &
3.D0/8.D0*EV%Kf(2)*ypol(3)+vb-9.D0/80.D0*EV%Kf(2)*octg+3.D0/40.D0*qg)/k-(- &
180.D0*ypolprime(2)-30.D0*pigdot)/k**2/160.D0)*dvis(j)+(-(9.D0*pigdot+ &
54.D0*ypolprime(2))/k**2*opac(j)/160.D0+pig/16.D0+clxg/4.D0+3.D0/8.D0*ypol(2)+(- &
21.D0/5.D0*adotoa*sigma-3.D0/8.D0*EV%Kf(2)*ypolprime(3)+vbdot+3.D0/40.D0*qgdot- &
9.D0/80.D0*EV%Kf(2)*octgprime)/k+(-9.D0/160.D0*dopac(j)*pig-21.D0/10.D0*dgpi-27.D0/ &
80.D0*dopac(j)*ypol(2))/k**2)*vis(j)+(3.D0/16.D0*ddvis(j)*pig+9.D0/ &
8.D0*ddvis(j)*ypol(2))/k**2+21.D0/10.D0/k/EV%Kf(1)*vis(j)*etak
else if (CP%EFTflag/=0 .and. a>=EFTturnonpi) then
polterdot = 9._dl/15._dl*ypolprime(2) + 0.1_dl*pigdot
ISW = expmmu(j)*(EFTISW)/k
sources(1) = ISW&
& +vis(j)* (clxg/4.D0 + polter/1.6d0 + vbdot/k -9.D0*(polterdot)/k2*opac(j)/16.D0 -9.D0/16.D0*dopac(j)*polter/k2&
& + 21.D0/10.D0*EFTsigmadot/k + 3.D0/40.D0*qgdot/k &
& +(-3.D0/8.D0*EV%Kf(2)*ypolprime(3) - 9.D0/80.D0*EV%Kf(2)*octgprime)/k)&
& +((-9.D0/160.D0*pig-27.D0/80.D0*ypol(2))/k**2*opac(j)+(11.D0/10.D0*sigma- &
& 3.D0/8.D0*EV%Kf(2)*ypol(3)+vb-9.D0/80.D0*EV%Kf(2)*octg+3.D0/40.D0*qg)/k-(- &
& 180.D0*ypolprime(2)-30.D0*pigdot)/k**2/160.D0)*dvis(j)&
& +(3.D0/16.D0*ddvis(j)*pig+9.D0/8.D0*ddvis(j)*ypol(2))/k**2
end if
! EFTCAMB MOD END
if (CP%EFTflag==0.or.a < EFTturnonpi) then ! Normal CAMB code.
!Maple's fortran output - see scal_eqs.map
!2phi' term (\phi' + \psi' in Newtonian gauge)
ISW = (4.D0/3.D0*k*EV%Kf(1)*sigma+(-2.D0/3.D0*sigma-2.D0/3.D0*etak/adotoa)*k &
-diff_rhopi/k**2-1.D0/adotoa*dgrho/3.D0+(3.D0*gpres+5.D0*grho)*sigma/k/3.D0 &
-2.D0/k*adotoa/EV%Kf(1)*etak)*expmmu(j)
!e.g. to get only late-time ISW
! if (1/a-1 < 30) ISW=0
!The rest, note y(9)->octg, yprime(9)->octgprime (octopoles)
sources(1)= ISW + ((-9.D0/160.D0*pig-27.D0/80.D0*ypol(2))/k**2*opac(j)+(11.D0/10.D0*sigma- &
3.D0/8.D0*EV%Kf(2)*ypol(3)+vb-9.D0/80.D0*EV%Kf(2)*octg+3.D0/40.D0*qg)/k-(- &
180.D0*ypolprime(2)-30.D0*pigdot)/k**2/160.D0)*dvis(j)+(-(9.D0*pigdot+ &
54.D0*ypolprime(2))/k**2*opac(j)/160.D0+pig/16.D0+clxg/4.D0+3.D0/8.D0*ypol(2)+(- &
21.D0/5.D0*adotoa*sigma-3.D0/8.D0*EV%Kf(2)*ypolprime(3)+vbdot+3.D0/40.D0*qgdot- &
9.D0/80.D0*EV%Kf(2)*octgprime)/k+(-9.D0/160.D0*dopac(j)*pig-21.D0/10.D0*dgpi-27.D0/ &
80.D0*dopac(j)*ypol(2))/k**2)*vis(j)+(3.D0/16.D0*ddvis(j)*pig+9.D0/ &
8.D0*ddvis(j)*ypol(2))/k**2+21.D0/10.D0/k/EV%Kf(1)*vis(j)*etak
else if (CP%EFTflag/=0 .and. a>=EFTturnonpi) then
polterdot = 9._dl/15._dl*ypolprime(2) + 0.1_dl*pigdot
ISW = expmmu(j)*(EFTISW)/k
sources(1) = ISW&
& +vis(j)* (clxg/4.D0 + polter/1.6d0 + vbdot/k -9.D0*(polterdot)/k2*opac(j)/16.D0 -9.D0/16.D0*dopac(j)*polter/k2&
& + 21.D0/10.D0*EFTsigmadot/k + 3.D0/40.D0*qgdot/k &
& +(-3.D0/8.D0*EV%Kf(2)*ypolprime(3) - 9.D0/80.D0*EV%Kf(2)*octgprime)/k)&
& +((-9.D0/160.D0*pig-27.D0/80.D0*ypol(2))/k**2*opac(j)+(11.D0/10.D0*sigma- &
& 3.D0/8.D0*EV%Kf(2)*ypol(3)+vb-9.D0/80.D0*EV%Kf(2)*octg+3.D0/40.D0*qg)/k-(- &
& 180.D0*ypolprime(2)-30.D0*pigdot)/k**2/160.D0)*dvis(j)&
& +(3.D0/16.D0*ddvis(j)*pig+9.D0/8.D0*ddvis(j)*ypol(2))/k**2
end if
! EFTCAMB MOD END
Replace:
! EFTCAMB MOD START: Lensing source
if (CP%EFTflag==0 .or. a
!phi_lens = Phi - 1/2 kappa (a/k)^2 sum_i rho_i pi_i
!Neglect pi contributions because not accurate at late time anyway
phi = -(dgrho +3*dgq*adotoa/k)/(k2*EV%Kf(1)*2)
! - (grhor_t*pir + grhog_t*pig+ pinu*gpnu_t)/k2
sources(3) = -2*phi*f_K(tau-tau_maxvis)/(f_K(CP%tau0-tau_maxvis)*f_K(CP%tau0-tau))
!sources(3) = -2*phi*(tau-tau_maxvis)/((CP%tau0-tau_maxvis)*(CP%tau0-tau))
!We include the lensing factor of two here
else if (CP%EFTflag/=0 .and. a>=EFTturnonpi) then
sources(3) = -(EFTLensing)/k*f_K(tau-tau_maxvis)/(f_K(CP%tau0-tau_maxvis)*f_K(CP%tau0-tau))
end if
! Original code:
! !phi_lens = Phi - 1/2 kappa (a/k)^2 sum_i rho_i pi_i
! !Neglect pi contributions because not accurate at late time anyway
! phi = -(dgrho +3*dgq*adotoa/k)/(k2*EV%Kf(1)*2)
! ! - dgpi/k2/2
! sources(3) = -2*phi*f_K(tau-tau_maxvis)/(f_K(CP%tau0-tau_maxvis)*f_K(CP%tau0-tau))
! ! sources(3) = -2*phi*(tau-tau_maxvis)/((CP%tau0-tau_maxvis)*(CP%tau0-tau))
! EFTCAMB MOD END
if (CP%EFTflag==0 .or. a
!Neglect pi contributions because not accurate at late time anyway
phi = -(dgrho +3*dgq*adotoa/k)/(k2*EV%Kf(1)*2)
! - (grhor_t*pir + grhog_t*pig+ pinu*gpnu_t)/k2
sources(3) = -2*phi*f_K(tau-tau_maxvis)/(f_K(CP%tau0-tau_maxvis)*f_K(CP%tau0-tau))
!sources(3) = -2*phi*(tau-tau_maxvis)/((CP%tau0-tau_maxvis)*(CP%tau0-tau))
!We include the lensing factor of two here
else if (CP%EFTflag/=0 .and. a>=EFTturnonpi) then
sources(3) = -(EFTLensing)/k*f_K(tau-tau_maxvis)/(f_K(CP%tau0-tau_maxvis)*f_K(CP%tau0-tau))
end if
! !phi_lens = Phi - 1/2 kappa (a/k)^2 sum_i rho_i pi_i
! !Neglect pi contributions because not accurate at late time anyway
! phi = -(dgrho +3*dgq*adotoa/k)/(k2*EV%Kf(1)*2)
! ! - dgpi/k2/2
! sources(3) = -2*phi*f_K(tau-tau_maxvis)/(f_K(CP%tau0-tau_maxvis)*f_K(CP%tau0-tau))
! ! sources(3) = -2*phi*(tau-tau_maxvis)/((CP%tau0-tau_maxvis)*(CP%tau0-tau))
! EFTCAMB MOD END
Add: Before the end
! EFTCAMB MOD START: computation of mu and gamma. These are not used but it's nice to have the possibility to plot them...
if (CP%EFTflag/=0.and.a>=EFTturnonpi) then
EFT_Psi = EFTsigmadot/k + adotoa*sigma/k
EFT_Phi = etak/k - adotoa*sigma/k
EFT_function_mu = -2._dl*k*(EFTsigmadot+adotoa*sigma)/(dgrho)
EFT_function_gamma = (etak - adotoa*sigma)/(EFTsigmadot + adotoa*sigma)
EFT_xi = vq/(adotoa*clxq)
! This is used to print mu and gamma for different k values.
! If you don't want to print all the other quantities then
! turn off DenugEFTO and turn the other to true.
if (DebugEFTOutput.or..false.) then
write (10,'(10e15.5)') a, tau, 1._dl/a-1._dl, k, EFT_function_mu, EFT_function_gamma, EFT_xi
end if
end if
! EFTCAMB MOD END
if (CP%EFTflag/=0.and.a>=EFTturnonpi) then
EFT_Psi = EFTsigmadot/k + adotoa*sigma/k
EFT_Phi = etak/k - adotoa*sigma/k
EFT_function_mu = -2._dl*k*(EFTsigmadot+adotoa*sigma)/(dgrho)
EFT_function_gamma = (etak - adotoa*sigma)/(EFTsigmadot + adotoa*sigma)
EFT_xi = vq/(adotoa*clxq)
! This is used to print mu and gamma for different k values.
! If you don't want to print all the other quantities then
! turn off DenugEFTO and turn the other to true.
if (DebugEFTOutput.or..false.) then
write (10,'(10e15.5)') a, tau, 1._dl/a-1._dl, k, EFT_function_mu, EFT_function_gamma, EFT_xi
end if
end if
! EFTCAMB MOD END
::::: subroutine initial(EV,y, tau) :::::
Add: Before y=0
! EFTCAMB MOD START: compute the time at which turn on EFT.
EV%EFTturnOnTime = DeltaTime(0._dl,EFTturnonpi)
! EFTCAMB MOD END
EV%EFTturnOnTime = DeltaTime(0._dl,EFTturnonpi)
! EFTCAMB MOD END
-2-