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 derivs(EV,n,tau,ay,ayprime) :::::
Add: After real(dl) cothxor !1/tau in flat case
! EFTCAMB MOD START: definition of variables
! background quantities. Designer approach.
real(dl) EFT_H0,Hdot, Hdotdot
! 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_dgpnu
! 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
! Debug quantities
real(dl) temp1, temp2, temp3,temp4, temp5, temp6, temp7
! EFTCAMB MOD END
! background quantities. Designer approach.
real(dl) EFT_H0,Hdot, Hdotdot
! 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_dgpnu
! 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
! Debug quantities
real(dl) temp1, temp2, temp3,temp4, temp5, temp6, temp7
! EFTCAMB MOD END
Replace:
! EFTCAMB MOD START: Compute the DE background density.
if (CP%EFTflag==0) then
if (w_lam==-1._dl) then
grhov_t=grhov*a2
else
grhov_t=grhov*a**(-1._dl-3._dl*w_lam)
end if
else if (CP%EFTflag/=0) then
grhov_t=grhov*EFTw(a,3)
EFT_H0 = (CP%h0/c_EFT)*1000._dl
end if
! Original code:
! if (w_lam==-1._dl) then
! grhov_t=grhov*a2
! else
! grhov_t=grhov*a**(-1-3*w_lam)
! end if
! EFTCAMB MOD END
if (CP%EFTflag==0) then
if (w_lam==-1._dl) then
grhov_t=grhov*a2
else
grhov_t=grhov*a**(-1._dl-3._dl*w_lam)
end if
else if (CP%EFTflag/=0) then
grhov_t=grhov*EFTw(a,3)
EFT_H0 = (CP%h0/c_EFT)*1000._dl
end if
! Original code:
! if (w_lam==-1._dl) then
! grhov_t=grhov*a2
! else
! grhov_t=grhov*a**(-1-3*w_lam)
! end if
! EFTCAMB MOD END
Replace:
! EFTCAMB MOD START: compatibility with massive neutrinos.
EFT_dgpnu = 0._dl
if (CP%Num_Nu_Massive > 0) then
call MassiveNuVars(EV,ay,a,grho,gpres,dgrho,dgq, wnu_arr=wnu_arr, dgp=EFT_dgpnu)
end if
! Original code:
! if (CP%Num_Nu_Massive > 0) then
! call MassiveNuVars(EV,ay,a,grho,gpres,dgrho,dgq, wnu_arr)
! end if
! EFTCAMB MOD END.
EFT_dgpnu = 0._dl
if (CP%Num_Nu_Massive > 0) then
call MassiveNuVars(EV,ay,a,grho,gpres,dgrho,dgq, wnu_arr=wnu_arr, dgp=EFT_dgpnu)
end if
! Original code:
! if (CP%Num_Nu_Massive > 0) then
! call MassiveNuVars(EV,ay,a,grho,gpres,dgrho,dgq, wnu_arr)
! end if
! EFTCAMB MOD END.
Replace:
! EFTCAMB MOD START: protection against non flatness.
if (CP%flat) then
adotoa=sqrt(grho/3)
cothxor=1._dl/tau
else if (CP%EFTflag==0) then
adotoa=sqrt((grho+grhok)/3._dl)
cothxor=1._dl/tanfunc(tau/CP%r)/CP%r
else
stop "EFTcamb is not working for non-flat universes."
end if
! Original code:
! if (CP%flat) then
! adotoa=sqrt(grho/3)
! cothxor=1._dl/tau
! else
! adotoa=sqrt((grho+grhok)/3._dl)
! cothxor=1._dl/tanfunc(tau/CP%r)/CP%r
! end if
! EFTCAMB MOD END
if (CP%flat) then
adotoa=sqrt(grho/3)
cothxor=1._dl/tau
else if (CP%EFTflag==0) then
adotoa=sqrt((grho+grhok)/3._dl)
cothxor=1._dl/tanfunc(tau/CP%r)/CP%r
else
stop "EFTcamb is not working for non-flat universes."
end if
! Original code:
! if (CP%flat) then
! adotoa=sqrt(grho/3)
! cothxor=1._dl/tau
! else
! adotoa=sqrt((grho+grhok)/3._dl)
! cothxor=1._dl/tanfunc(tau/CP%r)/CP%r
! end if
! EFTCAMB MOD END
! EFTCAMB MOD START: Computation of background quantities:
! 1) Computation of FRW background quantities.
if (CP%EFTflag/=0) then
EFT_grhonudot_tot = 0._dl
EFT_gpinudot_tot = 0._dl
EFT_grhonu_tot = 0._dl
EFT_gpinu_tot = 0._dl
! Massive neutrinos mod:
if (CP%Num_Nu_Massive /= 0) then
do nu_i = 1, CP%Nu_mass_eigenstates
EFT_grhonu = 0._dl
EFT_gpinu = 0._dl
EFT_grhonudot = 0._dl
EFT_gpinudot = 0._dl
grhormass_t=grhormass(nu_i)/a**2
call Nu_background(a*nu_masses(nu_i),EFT_grhonu,EFT_gpinu)
EFT_grhonudot_tot = EFT_grhonudot_tot + grhormass_t*(Nu_drho(a*nu_masses(nu_i),adotoa, EFT_grhonu)&
& -4._dl*adotoa*EFT_grhonu)
EFT_gpinudot_tot = EFT_gpinudot_tot + grhormass_t*(Nu_pidot(a*nu_masses(nu_i),adotoa, EFT_gpinu)&
& -4._dl*adotoa*EFT_gpinu)
EFT_grhonu_tot = EFT_grhonu_tot + grhormass_t*EFT_grhonu
EFT_gpinu_tot = EFT_gpinu_tot + grhormass_t*EFT_gpinu
end do
end if
gpres=gpres + (grhog_t+grhor_t)/3._dl +EFTw(a,0)*grhov_t
adotdota=(adotoa*adotoa-gpres)/2._dl
Hdot =adotdota-adotoa**2._dl
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*0.1_dl) 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.DebugEFT) 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
! 1) Computation of FRW background quantities.
if (CP%EFTflag/=0) then
EFT_grhonudot_tot = 0._dl
EFT_gpinudot_tot = 0._dl
EFT_grhonu_tot = 0._dl
EFT_gpinu_tot = 0._dl
! Massive neutrinos mod:
if (CP%Num_Nu_Massive /= 0) then
do nu_i = 1, CP%Nu_mass_eigenstates
EFT_grhonu = 0._dl
EFT_gpinu = 0._dl
EFT_grhonudot = 0._dl
EFT_gpinudot = 0._dl
grhormass_t=grhormass(nu_i)/a**2
call Nu_background(a*nu_masses(nu_i),EFT_grhonu,EFT_gpinu)
EFT_grhonudot_tot = EFT_grhonudot_tot + grhormass_t*(Nu_drho(a*nu_masses(nu_i),adotoa, EFT_grhonu)&
& -4._dl*adotoa*EFT_grhonu)
EFT_gpinudot_tot = EFT_gpinudot_tot + grhormass_t*(Nu_pidot(a*nu_masses(nu_i),adotoa, EFT_gpinu)&
& -4._dl*adotoa*EFT_gpinu)
EFT_grhonu_tot = EFT_grhonu_tot + grhormass_t*EFT_grhonu
EFT_gpinu_tot = EFT_gpinu_tot + grhormass_t*EFT_gpinu
end do
end if
gpres=gpres + (grhog_t+grhor_t)/3._dl +EFTw(a,0)*grhov_t
adotdota=(adotoa*adotoa-gpres)/2._dl
Hdot =adotdota-adotoa**2._dl
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*0.1_dl) 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.DebugEFT) 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
Replace:
! EFTCAMB MOD START: initialization of DE perturbations.
if (CP%EFTflag==0) then
if (w_lam /= -1 .and. w_Perturb) then
clxq=ay(EV%w_ix)
vq=ay(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=ay(EV%w_ix)
vq=ay(EV%w_ix+1)
pidotdot = ayprime(EV%w_ix+1)
end if
! Original code:
! if (w_lam /= -1 .and. w_Perturb) then
! clxq=ay(EV%w_ix)
! vq=ay(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=ay(EV%w_ix)
vq=ay(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=ay(EV%w_ix)
vq=ay(EV%w_ix+1)
pidotdot = ayprime(EV%w_ix+1)
end if
! Original code:
! if (w_lam /= -1 .and. w_Perturb) then
! clxq=ay(EV%w_ix)
! vq=ay(EV%w_ix+1)
! dgrho=dgrho + clxq*grhov_t
! dgq = dgq + vq*grhov_t*(1+w_lam)
! end if
! EFTCAMB MOD END
Add: Following the above
! EFTCAMB MOD START: computation of Einstein equations factors.
if (CP%EFTflag/=0.and.a>=EFTturnonpi*0.1_dl) 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 EFT Einstein equations factors.
if (a>=EFTturnonpi*0.1_dl.and.DebugEFT) then
write (7,'(15E15.5)') a, tau, EFTeomF, EFTeomN, EFTeomX, EFTeomY, EFTeomG, EFTeomU,&
& EFTeomL, EFTeomM, EFTeomNdot, EFTeomVdot, EFTeomXdot
end if
! EFTCAMB MOD END
if (CP%EFTflag/=0.and.a>=EFTturnonpi*0.1_dl) 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 EFT Einstein equations factors.
if (a>=EFTturnonpi*0.1_dl.and.DebugEFT) then
write (7,'(15E15.5)') a, tau, EFTeomF, EFTeomN, EFTeomX, EFTeomY, EFTeomG, EFTeomU,&
& EFTeomL, EFTeomM, EFTeomNdot, EFTeomVdot, EFTeomXdot
end if
! EFTCAMB MOD END
Add: Following the above
! EFTCAMB MOD START: compute z,dz before loading radiation and photon
if (CP%EFTflag==0.or..not.EV%EFTCAMBactive) then
z=(0.5_dl*dgrho/k + etak)/adotoa
dz= -adotoa*z - 0.5_dl*dgrho/k
else if (CP%EFTflag/=0.and.EV%EFTCAMBactive) 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._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
! EFTCAMB MOD END
if (CP%EFTflag==0.or..not.EV%EFTCAMBactive) then
z=(0.5_dl*dgrho/k + etak)/adotoa
dz= -adotoa*z - 0.5_dl*dgrho/k
else if (CP%EFTflag/=0.and.EV%EFTCAMBactive) 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._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
! EFTCAMB MOD END
-3-