Replace: After "call MassiveNuVars(EV,ay,a,grho,gpres,dgrho,dgq, wnu_arr)"
! 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:
! 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:
! 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
Add: follows above
! EFTCAMB MOD START: Computation of background quantities:
! 1) Computation of FRW background quantities.
if (CP%EFTflag/=0) then
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))
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
if (CP%EFTflag==1) 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)))
else if (CP%EFTflag==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,'(25e15.5)') a, tau, EFTOmegaV, EFTOmegaP, EFTOmegaPP, EFTOmegaPPP,&
& EFTAlpha1V, EFTAlpha1P, EFTAlpha2V, EFTAlpha2P, EFTAlpha3V, EFTAlpha3P,&
& EFTAlpha4V, EFTAlpha4P, EFTAlpha4PP, EFTAlpha5V, EFTAlpha5P, EFTAlpha6V, EFTAlpha6P
write (2,'(12e15.6)') a, tau, EFTc, EFTLambda, EFTcdot, EFTLambdadot, EFTgrhoq, EFTgpresq, EFTgrhodotq, EFTgpresdotq
write (3,'(8e15.5)') a, tau, adotoa, Hdot, Hdotdot, grhov_t
end if
! EFTCAMB MOD END
! 1) Computation of FRW background quantities.
if (CP%EFTflag/=0) then
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))
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
if (CP%EFTflag==1) 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)))
else if (CP%EFTflag==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,'(25e15.5)') a, tau, EFTOmegaV, EFTOmegaP, EFTOmegaPP, EFTOmegaPPP,&
& EFTAlpha1V, EFTAlpha1P, EFTAlpha2V, EFTAlpha2P, EFTAlpha3V, EFTAlpha3P,&
& EFTAlpha4V, EFTAlpha4P, EFTAlpha4PP, EFTAlpha5V, EFTAlpha5P, EFTAlpha6V, EFTAlpha6P
write (2,'(12e15.6)') a, tau, EFTc, EFTLambda, EFTcdot, EFTLambdadot, EFTgrhoq, EFTgpresq, EFTgrhodotq, EFTgpresdotq
write (3,'(8e15.5)') a, tau, adotoa, Hdot, Hdotdot, grhov_t
end if
! EFTCAMB MOD END
Replace: follows above
! 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:
! 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:
! 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: follows 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 (6,'(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 (6,'(15E15.5)') a, tau, EFTeomF, EFTeomN, EFTeomX, EFTeomY, EFTeomG, EFTeomU,&
& EFTeomL, EFTeomM, EFTeomNdot, EFTeomVdot, EFTeomXdot
end if
! EFTCAMB MOD END
Replace: follows 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
if (EV%no_nu_multpoles) then
!RSA approximation of arXiv:1104.2933, dropping opactity terms in the velocity
!Approximate total density variables with just matter terms
clxr=-4*dz/k
qr=-4._dl/3*z
pir=0
else
! Massless neutrinos
clxr=ay(EV%r_ix)
qr =ay(EV%r_ix+1)
pir =ay(EV%r_ix+2)
end if
if (EV%no_phot_multpoles) then
if (.not. EV%no_nu_multpoles) then
clxg=-4*dz/k-4/k*opacity*(vb+z)
qg=-4._dl/3*z
else
clxg=clxr-4/k*opacity*(vb+z)
qg=qr
end if
pig=0
else
! Photons
clxg=ay(EV%g_ix)
qg=ay(EV%g_ix+1)
if (.not. EV%TightCoupling) pig=ay(EV%g_ix+2)
end if
!... Original:
! if (EV%no_nu_multpoles) then
! !RSA approximation of arXiv:1104.2933, dropping opactity terms in the velocity
! !Approximate total density variables with just matter terms
! 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
! else
! ! Massless neutrinos
! clxr=ay(EV%r_ix)
! qr =ay(EV%r_ix+1)
! pir =ay(EV%r_ix+2)
! endif
! if (EV%no_phot_multpoles) then
! if (.not. EV%no_nu_multpoles) then
! z=(0.5_dl*dgrho/k + etak)/adotoa
! dz= -adotoa*z - 0.5_dl*dgrho/k
! clxg=-4*dz/k-4/k*opacity*(vb+z)
! qg=-4._dl/3*z
! else
! clxg=clxr-4/k*opacity*(vb+z)
! qg=qr
! end if
! pig=0
! else
! ! Photons
! clxg=ay(EV%g_ix)
! qg=ay(EV%g_ix+1)
! if (.not. EV%TightCoupling) pig=ay(EV%g_ix+2)
! 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
if (EV%no_nu_multpoles) then
!RSA approximation of arXiv:1104.2933, dropping opactity terms in the velocity
!Approximate total density variables with just matter terms
clxr=-4*dz/k
qr=-4._dl/3*z
pir=0
else
! Massless neutrinos
clxr=ay(EV%r_ix)
qr =ay(EV%r_ix+1)
pir =ay(EV%r_ix+2)
end if
if (EV%no_phot_multpoles) then
if (.not. EV%no_nu_multpoles) then
clxg=-4*dz/k-4/k*opacity*(vb+z)
qg=-4._dl/3*z
else
clxg=clxr-4/k*opacity*(vb+z)
qg=qr
end if
pig=0
else
! Photons
clxg=ay(EV%g_ix)
qg=ay(EV%g_ix+1)
if (.not. EV%TightCoupling) pig=ay(EV%g_ix+2)
end if
!... Original:
! if (EV%no_nu_multpoles) then
! !RSA approximation of arXiv:1104.2933, dropping opactity terms in the velocity
! !Approximate total density variables with just matter terms
! 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
! else
! ! Massless neutrinos
! clxr=ay(EV%r_ix)
! qr =ay(EV%r_ix+1)
! pir =ay(EV%r_ix+2)
! endif
! if (EV%no_phot_multpoles) then
! if (.not. EV%no_nu_multpoles) then
! z=(0.5_dl*dgrho/k + etak)/adotoa
! dz= -adotoa*z - 0.5_dl*dgrho/k
! clxg=-4*dz/k-4/k*opacity*(vb+z)
! qg=-4._dl/3*z
! else
! clxg=clxr-4/k*opacity*(vb+z)
! qg=qr
! end if
! pig=0
! else
! ! Photons
! clxg=ay(EV%g_ix)
! qg=ay(EV%g_ix+1)
! if (.not. EV%TightCoupling) pig=ay(EV%g_ix+2)
! end if
! EFTCAMB MOD END
Replace: After "ayprime(1)=adotoa*a"
! EFTCAMB MOD START: equation of motion. Now we use the non-RSA expression for dz.
if (CP%EFTflag==0 .or. .not. EV%EFTCAMBactive) then
!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
if (CP%flat) then
sigma=(z+1.5_dl*dgq/k2)
ayprime(2)=0.5_dl*dgq
else
sigma=(z+1.5_dl*dgq/k2)/EV%Kf(1)
ayprime(2)=0.5_dl*dgq + CP%curv*z
end if
else if (CP%EFTflag/=0.and.CP%flat.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*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)
ayprime(2)=0.5_dl*dgq/(1._dl+EFTOmegaV) + k2/3._dl*EFTeomF/EFT_H0
end if
!... Original:
! ! 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
! if (CP%flat) then
! !eta*k equation
! sigma=(z+1.5_dl*dgq/k2)
! ayprime(2)=0.5_dl*dgq
! else
! sigma=(z+1.5_dl*dgq/k2)/EV%Kf(1)
! ayprime(2)=0.5_dl*dgq + CP%curv*z
! end if
! EFTCAMB MOD END
if (CP%EFTflag==0 .or. .not. EV%EFTCAMBactive) then
!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
if (CP%flat) then
sigma=(z+1.5_dl*dgq/k2)
ayprime(2)=0.5_dl*dgq
else
sigma=(z+1.5_dl*dgq/k2)/EV%Kf(1)
ayprime(2)=0.5_dl*dgq + CP%curv*z
end if
else if (CP%EFTflag/=0.and.CP%flat.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*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)
ayprime(2)=0.5_dl*dgq/(1._dl+EFTOmegaV) + k2/3._dl*EFTeomF/EFT_H0
end if
!... Original:
! ! 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
! if (CP%flat) then
! !eta*k equation
! sigma=(z+1.5_dl*dgq/k2)
! ayprime(2)=0.5_dl*dgq
! else
! sigma=(z+1.5_dl*dgq/k2)/EV%Kf(1)
! ayprime(2)=0.5_dl*dgq + CP%curv*z
! end if
! EFTCAMB MOD END
Replace: follows above
! EFTCAMB MOD START: Pi field equation.
! 1) Computation of pi field equation prefactors.
if (CP%EFTflag/=0 .and. a>=EFTturnonpi*0.1_dl) then
EFTpiA = EFTc +2*a2*EFT_H0**2*EFTAlpha1V**4 +3._dl/2._dl*a2*(adotoa*EFTOmegaP+EFT_H0*EFTAlpha2V**3)**2&
&/(2._dl*(1+EFTOmegaV)+EFTAlpha3V**2+EFTAlpha4V**2) +4._dl*a2*k2*EFTAlpha6V**2
!
EFTpiB = EFTcdot +4._dl*adotoa*EFTc +8._dl*a2*adotoa*EFT_H0*(EFTAlpha1V**4+a*EFTAlpha1V**3*EFTAlpha1P)&
& +4._dl*a2*k2*adotoa*(3._dl*EFTAlpha6V**2+a*EFTAlpha6V*EFTAlpha6P) +a*k2*(EFTAlpha4V**2&
& +2._dl*EFTAlpha5V**2)/(2._dl*(1._dl+EFTOmegaV)-2._dl*EFTAlpha4V**2)*(adotoa*EFTOmegaP+EFT_H0*EFTAlpha2V**3)&
& -a*(adotoa*EFTOmegaP+EFT_H0*EFTAlpha2V**3)/(4._dl*(1._dl+EFTOmegaV)+6._dl*EFTAlpha3V**2 +2._dl*EFTAlpha4V**2)*&
&(-3._dl*(EFTgrhoQ + EFTgpresQ) -3._dl*a*adotoa**2*EFTOmegaP*(4._dl +Hdot/(adotoa**2)) -3._dl*a2*adotoa**2*EFTOmegaPP&
& -3*a*adotoa*EFT_H0*(4._dl*EFTAlpha2V**3 +3._dl*a*EFTAlpha2V**2*EFTAlpha2P) -(9._dl*EFTAlpha3V**2-3._dl*EFTAlpha4V**2)*&
&(Hdot-adotoa**2) +k2*(3._dl*EFTAlpha3V**2-EFTAlpha4V**2+4._dl*EFTAlpha5V**2))&
& +1._dl/(1._dl+EFTOmegaV+2._dl*EFTAlpha5V**2)*(a*adotoa*EFTOmegaP+2._dl*adotoa*(EFTAlpha5V**2+2._dl*EFTAlpha5V*EFTAlpha5P)&
& -(1._dl+EFTOmegaV)*(a*adotoa*EFTOmegaP+a*EFT_H0*EFTAlpha2V**3)/(2._dl*(1._dl+EFTOmegaV)+3._dl*EFTAlpha3V**2+EFTAlpha4V**2))*&
&(-EFTc +1.5_dl*a*adotoa**2*EFTOmegaP-2._dl*a2*EFT_H0*EFTAlpha1V**4-4._dl*EFTAlpha6V**2*k2+1.5_dl*a*adotoa*EFT_H0*EFTAlpha2V**3)
!
EFTpiC = +adotoa*EFTcdot +(6._dl*adotoa**2-2._dl*Hdot)*EFTc +1.5_dl*a*adotoa*EFTOmegaP*(Hdotdot-2._dl*adotoa**3) &
& +6._dl*a2*adotoa**2*EFT_H0**2*EFTAlpha1V**4 +2._dl*a2*Hdot*EFT_H0**2*EFTAlpha1V**4 &
& +8._dl*a2*a*adotoa**2*EFT_H0**2*EFTAlpha1V**3*EFTAlpha1P +1.5_dl*(Hdot-adotoa**2)**2*(EFTAlpha4V**2+3._dl*EFTAlpha3V**2)&
& +4.5_dl*adotoa*EFT_H0*a*(Hdot-adotoa**2)*(EFTAlpha2V**3+a*EFTAlpha2V**2*EFTAlpha2P)&
& +0.5_dl*a*EFT_H0*EFTAlpha2V**3*(3._dl*Hdotdot -12._dl*Hdot*adotoa +6._dl*adotoa**3) &
& -a*(adotoa*EFTOmegaP+EFT_H0*EFTAlpha2V**3)/(4._dl*(1._dl+EFTOmegaV)+6._dl*EFTAlpha3V**2 +2._dl*EFTAlpha4V**2)*&
&(-3._dl*EFTgpresdotQ -3._dl*adotoa*(EFTgrhoQ +EFTgpresQ)-3._dl*a*adotoa**3*(a*EFTOmegaPP +6._dl*EFTOmegaP) &
& -6._dl*a*adotoa*Hdot*EFTOmegaP +3._dl*(Hdotdot-2._dl*adotoa*Hdot)*(EFTAlpha4V**2+3._dl*EFTAlpha3V**2)&
& +6._dl*adotoa*(Hdot-adotoa**2)*(3._dl*EFTAlpha3V**2+3._dl*a*EFTAlpha3V*EFTAlpha3P +EFTAlpha4V**2 +a*EFTAlpha4V*EFTAlpha4P)&
& -3._dl*a*EFT_H0*(3._dl*adotoa**2*EFTAlpha2V**3 +Hdot*EFTAlpha2V**3 +3._dl*a*adotoa**2*EFTAlpha2V**2*EFTAlpha2P))&
& +1._dl/(1._dl+EFTOmegaV+2._dl*EFTAlpha5V**2)*(a*adotoa*EFTOmegaP+2._dl*adotoa*(EFTAlpha5V**2+2._dl*a*EFTAlpha5V*EFTAlpha5P)&
& -(1._dl+EFTOmegaV)*(a*adotoa*EFTOmegaP+a*EFT_H0*EFTAlpha2V**3)/(2._dl*(1._dl+EFTOmegaV)+3._dl*EFTAlpha3V**2+EFTAlpha4V**2))*&
&(-0.5*EFTgrhodotQ -adotoa*EFTc +1.5_dl*a*adotoa*EFTOmegaP*(3._dl*adotoa**2-Hdot) -2._dl*a2*adotoa*EFT_H0**2*EFTAlpha1V**4&
& -1.5_dl*a*EFT_H0*EFTAlpha2V**3*(Hdot-2._dl*adotoa**2) -3._dl*adotoa*(Hdot-adotoa**2)*(1.5_dl*EFTAlpha3V**2+EFTAlpha4V**2))
!
EFTpiD = EFTc -0.5_dl*a*adotoa*EFT_H0*(EFTAlpha2V**3+3._dl*a*EFTAlpha2V**2*EFTAlpha2P) +(adotoa**2-Hdot)*(3._dl*EFTAlpha3V**2+EFTAlpha4V**2)&
& +4._dl*a2*(Hdot*EFTAlpha6V**2 +2._dl*adotoa**2*EFTAlpha6V**2 +a*adotoa**2*EFTAlpha6V*EFTAlpha6P)&
& +2._dl*(Hdot*EFTAlpha5V**2 +2._dl*a*adotoa**2*EFTAlpha5V*EFTAlpha5P)&
& -a*(adotoa*EFTOmegaP+EFT_H0*EFTAlpha2V**3)/(4._dl*(1._dl+EFTOmegaV)+6._dl*EFTAlpha3V**2 +2._dl*EFTAlpha4V**2)*&
&(-2._dl*a*adotoa*EFTOmegaP +4._dl*adotoa*EFTAlpha5V**2 -2._dl*adotoa*(3._dl*EFTAlpha3V**2 +3._dl*a*EFTAlpha3V*EFTAlpha3P &
& +EFTAlpha4V**2 +a*EFTAlpha4V*EFTAlpha4P))&
& +1._dl/(1._dl+EFTOmegaV+2._dl*EFTAlpha5V**2)*(a*adotoa*EFTOmegaP+2._dl*adotoa*(EFTAlpha5V**2+2._dl*a*EFTAlpha5V*EFTAlpha5P)&
& -(1._dl+EFTOmegaV)*(a*adotoa*EFTOmegaP+a*EFT_H0*EFTAlpha2V**3)/(2._dl*(1._dl+EFTOmegaV)+3._dl*EFTAlpha3V**2+EFTAlpha4V**2))*&
&(+0.5_dl*a*adotoa*EFTOmegaP -2._dl*adotoa*EFTAlpha5V**2 +0.5_dl*a*EFT_H0*EFTAlpha2V**3 +1.5_dl*adotoa*EFTAlpha3V**2&
& -adotoa*EFTAlpha4V**2 -4._dl*adotoa*EFTAlpha6V**2)&
& +(EFTAlpha4V**2 +2._dl*EFTAlpha5V**2)/(2._dl*(1._dl+EFTOmegaV) -2._dl*EFTAlpha4V**2)*(EFTgrhoQ +EFTgpresQ +a*adotoa**2*EFTOmegaP&
& -EFTAlpha4V**2*(Hdot-adotoa**2) +a*adotoa*EFT_H0*EFTAlpha2V**3 +3._dl*EFTAlpha3V**2*(adotoa**2-Hdot))&
& +k2*(+0.5_dl*EFTAlpha3V**2 +0.5_dl*EFTAlpha4V**2 &
& +(EFTAlpha4V**2 +2._dl*EFTAlpha5V**2)/(2._dl*(1._dl+EFTOmegaV) -2._dl*EFTAlpha4V**2)*(EFTAlpha3V**2+EFTAlpha4V**2))
!
EFTpiE = (EFTc -1.5_dl*a*adotoa**2*EFTOmegaP -0.5_dl*a*adotoa*EFT_H0*(2._dl*EFTAlpha2V**3 +3._dl*a*EFTAlpha2V**2*EFTAlpha2P)&
& +0.5_dl*EFTAlpha3V**2*(k2-3._dl*Hdot+adotoa**2) +0.5_dl*EFTAlpha4V**2*(k2-Hdot+adotoa**2)&
& -a*(adotoa*EFTOmegaP+EFT_H0*EFTAlpha2V**3)/(4._dl*(1._dl+EFTOmegaV)+6._dl*EFTAlpha3V**2 +2._dl*EFTAlpha4V**2)*&
&(-2._dl*adotoa*(a*EFTOmegaP +2._dl*(1._dl+EFTOmegaV)) -2._dl*adotoa*(3._dl*EFTAlpha3V**2+3._dl*a*EFTAlpha3V*EFTAlpha3P&
& +EFTAlpha4V**2 +a*EFTAlpha4V*EFTAlpha4P))&
& +1._dl/(1._dl+EFTOmegaV+2._dl*EFTAlpha5V**2)*(a*adotoa*EFTOmegaP+2._dl*adotoa*(EFTAlpha5V**2+2._dl*a*EFTAlpha5V*EFTAlpha5P)&
& -(1._dl+EFTOmegaV)*(a*adotoa*EFTOmegaP+a*EFT_H0*EFTAlpha2V**3)/(2._dl*(1._dl+EFTOmegaV)+3._dl*EFTAlpha3V**2+EFTAlpha4V**2))*&
&( +adotoa*(1._dl+EFTOmegaV+0.5_dl*a*EFTOmegaP) +0.5_dl*a*EFT_H0*EFTAlpha2V**3 +1.5_dl*adotoa*EFTAlpha3V**2 +adotoa*EFTAlpha4V**2)&
& +(EFTAlpha4V**2 +2._dl*EFTAlpha5V**2)/(2._dl*(1._dl+EFTOmegaV) -2._dl*EFTAlpha4V**2)*k2*(EFTAlpha4V**2+EFTAlpha3V**2))*k*z&
& +1._dl*a*(adotoa*EFTOmegaP+EFT_H0*EFTAlpha2V**3)/(4._dl*(1._dl+EFTOmegaV)+6._dl*EFTAlpha3V**2 +2._dl*EFTAlpha4V**2)*&
&(grhog_t*clxg+grhor_t*clxr) ++(EFTAlpha4V**2 +2._dl*EFTAlpha5V**2)/(2._dl*(1._dl+EFTOmegaV) -2._dl*EFTAlpha4V**2)*k*dgq&
& -0.5_dl/(1._dl+EFTOmegaV+2._dl*EFTAlpha5V**2)*(a*adotoa*EFTOmegaP+2._dl*adotoa*(EFTAlpha5V**2+2._dl*a*EFTAlpha5V*EFTAlpha5P)&
& -(1._dl+EFTOmegaV)*(a*adotoa*EFTOmegaP+a*EFT_H0*EFTAlpha2V**3)/(2._dl*(1._dl+EFTOmegaV)+3._dl*EFTAlpha3V**2+EFTAlpha4V**2))*dgrho
!
end if
! 2) Stability check of the theory: some unstable models may slip the stability check performed at the beginning.
if (CP%EFTflag /=0 .and. a>=EFTturnonpi) then
!1- Positive gravitational constant:
if (1._dl +EFTOmegaV <= 0) stop "Negative gravitational constant."
!2- Ghost instability:
if (EFTpiA <= 0) stop "Ghost instability."
!3- Tachionic perturbations:
!if (EFTpiD < 0) stop "Tachionic perturbations."
!4- Positive mass:
!if (EFTpiC < 0) stop "Negative scalaron mass."
end if
! 5) Solution of the pi field equation of motion.
if (CP%EFTflag==0) then
if (w_lam /= -1 .and. w_Perturb) then
ayprime(EV%w_ix)= -3*adotoa*(cs2_lam-w_lam)*(clxq+3*adotoa*(1+w_lam)*vq/k) &
-(1+w_lam)*k*vq -(1+w_lam)*k*z
ayprime(EV%w_ix+1) = -adotoa*(1-3*cs2_lam)*vq + k*cs2_lam*clxq/(1+w_lam)
end if
else if (CP%EFTflag/=0.and.EV%EFTCAMBactive) then
ayprime(EV%w_ix)= vq
ayprime(EV%w_ix+1)= -EFTpiB*vq/EFTpiA - (EFTpiC+ k*k*EFTpiD)*clxq/EFTpiA - EFT_H0*EFTpiE/EFTpiA
end if
! 5) Debug code: Print of resuls.
if (a>=EFTturnonpi.and. DebugEFT) then
write (4,'(12e15.5)') a, tau, EFTpiA, EFTpiB, EFTpiC, EFTpiD, EFTpiE
write (5,'(6d15.5)') a, tau, clxq, vq, ayprime(EV%w_ix+1)
write (7,'(6e15.5)') a, tau, sigma, z, dz, etak
end if
!... Original:
! if (w_lam /= -1 .and. w_Perturb) then
! ayprime(EV%w_ix)= -3*adotoa*(cs2_lam-w_lam)*(clxq+3*adotoa*(1+w_lam)*vq/k) &
! -(1+w_lam)*k*vq -(1+w_lam)*k*z
! ayprime(EV%w_ix+1) = -adotoa*(1-3*cs2_lam)*vq + k*cs2_lam*clxq/(1+w_lam)
! end if
! EFTCAMB MOD END
! 1) Computation of pi field equation prefactors.
if (CP%EFTflag/=0 .and. a>=EFTturnonpi*0.1_dl) then
EFTpiA = EFTc +2*a2*EFT_H0**2*EFTAlpha1V**4 +3._dl/2._dl*a2*(adotoa*EFTOmegaP+EFT_H0*EFTAlpha2V**3)**2&
&/(2._dl*(1+EFTOmegaV)+EFTAlpha3V**2+EFTAlpha4V**2) +4._dl*a2*k2*EFTAlpha6V**2
!
EFTpiB = EFTcdot +4._dl*adotoa*EFTc +8._dl*a2*adotoa*EFT_H0*(EFTAlpha1V**4+a*EFTAlpha1V**3*EFTAlpha1P)&
& +4._dl*a2*k2*adotoa*(3._dl*EFTAlpha6V**2+a*EFTAlpha6V*EFTAlpha6P) +a*k2*(EFTAlpha4V**2&
& +2._dl*EFTAlpha5V**2)/(2._dl*(1._dl+EFTOmegaV)-2._dl*EFTAlpha4V**2)*(adotoa*EFTOmegaP+EFT_H0*EFTAlpha2V**3)&
& -a*(adotoa*EFTOmegaP+EFT_H0*EFTAlpha2V**3)/(4._dl*(1._dl+EFTOmegaV)+6._dl*EFTAlpha3V**2 +2._dl*EFTAlpha4V**2)*&
&(-3._dl*(EFTgrhoQ + EFTgpresQ) -3._dl*a*adotoa**2*EFTOmegaP*(4._dl +Hdot/(adotoa**2)) -3._dl*a2*adotoa**2*EFTOmegaPP&
& -3*a*adotoa*EFT_H0*(4._dl*EFTAlpha2V**3 +3._dl*a*EFTAlpha2V**2*EFTAlpha2P) -(9._dl*EFTAlpha3V**2-3._dl*EFTAlpha4V**2)*&
&(Hdot-adotoa**2) +k2*(3._dl*EFTAlpha3V**2-EFTAlpha4V**2+4._dl*EFTAlpha5V**2))&
& +1._dl/(1._dl+EFTOmegaV+2._dl*EFTAlpha5V**2)*(a*adotoa*EFTOmegaP+2._dl*adotoa*(EFTAlpha5V**2+2._dl*EFTAlpha5V*EFTAlpha5P)&
& -(1._dl+EFTOmegaV)*(a*adotoa*EFTOmegaP+a*EFT_H0*EFTAlpha2V**3)/(2._dl*(1._dl+EFTOmegaV)+3._dl*EFTAlpha3V**2+EFTAlpha4V**2))*&
&(-EFTc +1.5_dl*a*adotoa**2*EFTOmegaP-2._dl*a2*EFT_H0*EFTAlpha1V**4-4._dl*EFTAlpha6V**2*k2+1.5_dl*a*adotoa*EFT_H0*EFTAlpha2V**3)
!
EFTpiC = +adotoa*EFTcdot +(6._dl*adotoa**2-2._dl*Hdot)*EFTc +1.5_dl*a*adotoa*EFTOmegaP*(Hdotdot-2._dl*adotoa**3) &
& +6._dl*a2*adotoa**2*EFT_H0**2*EFTAlpha1V**4 +2._dl*a2*Hdot*EFT_H0**2*EFTAlpha1V**4 &
& +8._dl*a2*a*adotoa**2*EFT_H0**2*EFTAlpha1V**3*EFTAlpha1P +1.5_dl*(Hdot-adotoa**2)**2*(EFTAlpha4V**2+3._dl*EFTAlpha3V**2)&
& +4.5_dl*adotoa*EFT_H0*a*(Hdot-adotoa**2)*(EFTAlpha2V**3+a*EFTAlpha2V**2*EFTAlpha2P)&
& +0.5_dl*a*EFT_H0*EFTAlpha2V**3*(3._dl*Hdotdot -12._dl*Hdot*adotoa +6._dl*adotoa**3) &
& -a*(adotoa*EFTOmegaP+EFT_H0*EFTAlpha2V**3)/(4._dl*(1._dl+EFTOmegaV)+6._dl*EFTAlpha3V**2 +2._dl*EFTAlpha4V**2)*&
&(-3._dl*EFTgpresdotQ -3._dl*adotoa*(EFTgrhoQ +EFTgpresQ)-3._dl*a*adotoa**3*(a*EFTOmegaPP +6._dl*EFTOmegaP) &
& -6._dl*a*adotoa*Hdot*EFTOmegaP +3._dl*(Hdotdot-2._dl*adotoa*Hdot)*(EFTAlpha4V**2+3._dl*EFTAlpha3V**2)&
& +6._dl*adotoa*(Hdot-adotoa**2)*(3._dl*EFTAlpha3V**2+3._dl*a*EFTAlpha3V*EFTAlpha3P +EFTAlpha4V**2 +a*EFTAlpha4V*EFTAlpha4P)&
& -3._dl*a*EFT_H0*(3._dl*adotoa**2*EFTAlpha2V**3 +Hdot*EFTAlpha2V**3 +3._dl*a*adotoa**2*EFTAlpha2V**2*EFTAlpha2P))&
& +1._dl/(1._dl+EFTOmegaV+2._dl*EFTAlpha5V**2)*(a*adotoa*EFTOmegaP+2._dl*adotoa*(EFTAlpha5V**2+2._dl*a*EFTAlpha5V*EFTAlpha5P)&
& -(1._dl+EFTOmegaV)*(a*adotoa*EFTOmegaP+a*EFT_H0*EFTAlpha2V**3)/(2._dl*(1._dl+EFTOmegaV)+3._dl*EFTAlpha3V**2+EFTAlpha4V**2))*&
&(-0.5*EFTgrhodotQ -adotoa*EFTc +1.5_dl*a*adotoa*EFTOmegaP*(3._dl*adotoa**2-Hdot) -2._dl*a2*adotoa*EFT_H0**2*EFTAlpha1V**4&
& -1.5_dl*a*EFT_H0*EFTAlpha2V**3*(Hdot-2._dl*adotoa**2) -3._dl*adotoa*(Hdot-adotoa**2)*(1.5_dl*EFTAlpha3V**2+EFTAlpha4V**2))
!
EFTpiD = EFTc -0.5_dl*a*adotoa*EFT_H0*(EFTAlpha2V**3+3._dl*a*EFTAlpha2V**2*EFTAlpha2P) +(adotoa**2-Hdot)*(3._dl*EFTAlpha3V**2+EFTAlpha4V**2)&
& +4._dl*a2*(Hdot*EFTAlpha6V**2 +2._dl*adotoa**2*EFTAlpha6V**2 +a*adotoa**2*EFTAlpha6V*EFTAlpha6P)&
& +2._dl*(Hdot*EFTAlpha5V**2 +2._dl*a*adotoa**2*EFTAlpha5V*EFTAlpha5P)&
& -a*(adotoa*EFTOmegaP+EFT_H0*EFTAlpha2V**3)/(4._dl*(1._dl+EFTOmegaV)+6._dl*EFTAlpha3V**2 +2._dl*EFTAlpha4V**2)*&
&(-2._dl*a*adotoa*EFTOmegaP +4._dl*adotoa*EFTAlpha5V**2 -2._dl*adotoa*(3._dl*EFTAlpha3V**2 +3._dl*a*EFTAlpha3V*EFTAlpha3P &
& +EFTAlpha4V**2 +a*EFTAlpha4V*EFTAlpha4P))&
& +1._dl/(1._dl+EFTOmegaV+2._dl*EFTAlpha5V**2)*(a*adotoa*EFTOmegaP+2._dl*adotoa*(EFTAlpha5V**2+2._dl*a*EFTAlpha5V*EFTAlpha5P)&
& -(1._dl+EFTOmegaV)*(a*adotoa*EFTOmegaP+a*EFT_H0*EFTAlpha2V**3)/(2._dl*(1._dl+EFTOmegaV)+3._dl*EFTAlpha3V**2+EFTAlpha4V**2))*&
&(+0.5_dl*a*adotoa*EFTOmegaP -2._dl*adotoa*EFTAlpha5V**2 +0.5_dl*a*EFT_H0*EFTAlpha2V**3 +1.5_dl*adotoa*EFTAlpha3V**2&
& -adotoa*EFTAlpha4V**2 -4._dl*adotoa*EFTAlpha6V**2)&
& +(EFTAlpha4V**2 +2._dl*EFTAlpha5V**2)/(2._dl*(1._dl+EFTOmegaV) -2._dl*EFTAlpha4V**2)*(EFTgrhoQ +EFTgpresQ +a*adotoa**2*EFTOmegaP&
& -EFTAlpha4V**2*(Hdot-adotoa**2) +a*adotoa*EFT_H0*EFTAlpha2V**3 +3._dl*EFTAlpha3V**2*(adotoa**2-Hdot))&
& +k2*(+0.5_dl*EFTAlpha3V**2 +0.5_dl*EFTAlpha4V**2 &
& +(EFTAlpha4V**2 +2._dl*EFTAlpha5V**2)/(2._dl*(1._dl+EFTOmegaV) -2._dl*EFTAlpha4V**2)*(EFTAlpha3V**2+EFTAlpha4V**2))
!
EFTpiE = (EFTc -1.5_dl*a*adotoa**2*EFTOmegaP -0.5_dl*a*adotoa*EFT_H0*(2._dl*EFTAlpha2V**3 +3._dl*a*EFTAlpha2V**2*EFTAlpha2P)&
& +0.5_dl*EFTAlpha3V**2*(k2-3._dl*Hdot+adotoa**2) +0.5_dl*EFTAlpha4V**2*(k2-Hdot+adotoa**2)&
& -a*(adotoa*EFTOmegaP+EFT_H0*EFTAlpha2V**3)/(4._dl*(1._dl+EFTOmegaV)+6._dl*EFTAlpha3V**2 +2._dl*EFTAlpha4V**2)*&
&(-2._dl*adotoa*(a*EFTOmegaP +2._dl*(1._dl+EFTOmegaV)) -2._dl*adotoa*(3._dl*EFTAlpha3V**2+3._dl*a*EFTAlpha3V*EFTAlpha3P&
& +EFTAlpha4V**2 +a*EFTAlpha4V*EFTAlpha4P))&
& +1._dl/(1._dl+EFTOmegaV+2._dl*EFTAlpha5V**2)*(a*adotoa*EFTOmegaP+2._dl*adotoa*(EFTAlpha5V**2+2._dl*a*EFTAlpha5V*EFTAlpha5P)&
& -(1._dl+EFTOmegaV)*(a*adotoa*EFTOmegaP+a*EFT_H0*EFTAlpha2V**3)/(2._dl*(1._dl+EFTOmegaV)+3._dl*EFTAlpha3V**2+EFTAlpha4V**2))*&
&( +adotoa*(1._dl+EFTOmegaV+0.5_dl*a*EFTOmegaP) +0.5_dl*a*EFT_H0*EFTAlpha2V**3 +1.5_dl*adotoa*EFTAlpha3V**2 +adotoa*EFTAlpha4V**2)&
& +(EFTAlpha4V**2 +2._dl*EFTAlpha5V**2)/(2._dl*(1._dl+EFTOmegaV) -2._dl*EFTAlpha4V**2)*k2*(EFTAlpha4V**2+EFTAlpha3V**2))*k*z&
& +1._dl*a*(adotoa*EFTOmegaP+EFT_H0*EFTAlpha2V**3)/(4._dl*(1._dl+EFTOmegaV)+6._dl*EFTAlpha3V**2 +2._dl*EFTAlpha4V**2)*&
&(grhog_t*clxg+grhor_t*clxr) ++(EFTAlpha4V**2 +2._dl*EFTAlpha5V**2)/(2._dl*(1._dl+EFTOmegaV) -2._dl*EFTAlpha4V**2)*k*dgq&
& -0.5_dl/(1._dl+EFTOmegaV+2._dl*EFTAlpha5V**2)*(a*adotoa*EFTOmegaP+2._dl*adotoa*(EFTAlpha5V**2+2._dl*a*EFTAlpha5V*EFTAlpha5P)&
& -(1._dl+EFTOmegaV)*(a*adotoa*EFTOmegaP+a*EFT_H0*EFTAlpha2V**3)/(2._dl*(1._dl+EFTOmegaV)+3._dl*EFTAlpha3V**2+EFTAlpha4V**2))*dgrho
!
end if
! 2) Stability check of the theory: some unstable models may slip the stability check performed at the beginning.
if (CP%EFTflag /=0 .and. a>=EFTturnonpi) then
!1- Positive gravitational constant:
if (1._dl +EFTOmegaV <= 0) stop "Negative gravitational constant."
!2- Ghost instability:
if (EFTpiA <= 0) stop "Ghost instability."
!3- Tachionic perturbations:
!if (EFTpiD < 0) stop "Tachionic perturbations."
!4- Positive mass:
!if (EFTpiC < 0) stop "Negative scalaron mass."
end if
! 5) Solution of the pi field equation of motion.
if (CP%EFTflag==0) then
if (w_lam /= -1 .and. w_Perturb) then
ayprime(EV%w_ix)= -3*adotoa*(cs2_lam-w_lam)*(clxq+3*adotoa*(1+w_lam)*vq/k) &
-(1+w_lam)*k*vq -(1+w_lam)*k*z
ayprime(EV%w_ix+1) = -adotoa*(1-3*cs2_lam)*vq + k*cs2_lam*clxq/(1+w_lam)
end if
else if (CP%EFTflag/=0.and.EV%EFTCAMBactive) then
ayprime(EV%w_ix)= vq
ayprime(EV%w_ix+1)= -EFTpiB*vq/EFTpiA - (EFTpiC+ k*k*EFTpiD)*clxq/EFTpiA - EFT_H0*EFTpiE/EFTpiA
end if
! 5) Debug code: Print of resuls.
if (a>=EFTturnonpi.and. DebugEFT) then
write (4,'(12e15.5)') a, tau, EFTpiA, EFTpiB, EFTpiC, EFTpiD, EFTpiE
write (5,'(6d15.5)') a, tau, clxq, vq, ayprime(EV%w_ix+1)
write (7,'(6e15.5)') a, tau, sigma, z, dz, etak
end if
!... Original:
! if (w_lam /= -1 .and. w_Perturb) then
! ayprime(EV%w_ix)= -3*adotoa*(cs2_lam-w_lam)*(clxq+3*adotoa*(1+w_lam)*vq/k) &
! -(1+w_lam)*k*vq -(1+w_lam)*k*z
! ayprime(EV%w_ix+1) = -adotoa*(1-3*cs2_lam)*vq + k*cs2_lam*clxq/(1+w_lam)
! end if
! EFTCAMB MOD END
-2-