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.
Comment off: After !Approximate total density variables with just matter terms
!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
Comment off: After if (.not. EV%no_nu_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. 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 code:
! 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 code:
! 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:
! 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 +3._dl*EFT_dgpnu) +(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 (8,'(10e15.5)') a, tau, EFTpiA, EFTpiB, EFTpiC, EFTpiD, EFTpiE
write (4,'(7d15.5)') a, tau, clxq, vq, ayprime(EV%w_ix+1)
write (5,'(7e15.5)') a, tau, sigma, z, dz, etak
write (6,'(7e15.5)') a, tau, clxc, clxb, clxg, clxr
end if
! Original code:
! 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 +3._dl*EFT_dgpnu) +(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 (8,'(10e15.5)') a, tau, EFTpiA, EFTpiB, EFTpiC, EFTpiD, EFTpiE
write (4,'(7d15.5)') a, tau, clxq, vq, ayprime(EV%w_ix+1)
write (5,'(7e15.5)') a, tau, sigma, z, dz, etak
write (6,'(7e15.5)') a, tau, clxc, clxb, clxg, clxr
end if
! Original code:
! 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
Replace:
! EFTCAMB MOD START:
if (CP%EFTflag==0) then
gpres=gpres+ (grhog_t+grhor_t)/3 +grhov_t*w_lam
adotdota=(adotoa*adotoa-gpres)/2
! Original code:
! ! ddota/a
! gpres=gpres+ (grhog_t+grhor_t)/3 +grhov_t*w_lam
! adotdota=(adotoa*adotoa-gpres)/2
end if
! EFTCAMB MOD END
if (CP%EFTflag==0) then
gpres=gpres+ (grhog_t+grhor_t)/3 +grhov_t*w_lam
adotdota=(adotoa*adotoa-gpres)/2
! Original code:
! ! ddota/a
! gpres=gpres+ (grhog_t+grhor_t)/3 +grhov_t*w_lam
! adotdota=(adotoa*adotoa-gpres)/2
end if
! EFTCAMB MOD END
Replace:
! EFTCAMB MOD START: shear
! Define shear derivative to first order
if (CP%EFTflag==0 .or. .not. EV%EFTCAMBactive) then
sigmadot = -2*adotoa*sigma-dgs/k+etak
else if (CP%EFTflag/=0 .and. EV%EFTCAMBactive) then
sigmadot = 1._dl/EFTeomX*(-2._dl*adotoa*(1._dl+EFTeomV)*sigma +etak&
& -dgs/k/(1._dl+EFTOmegaV) +EFTeomN/EFT_H0)
end if
! Original code:
! sigmadot = -2*adotoa*sigma-dgs/k+etak
! EFTCAMB MOD END
! Define shear derivative to first order
if (CP%EFTflag==0 .or. .not. EV%EFTCAMBactive) then
sigmadot = -2*adotoa*sigma-dgs/k+etak
else if (CP%EFTflag/=0 .and. EV%EFTCAMBactive) then
sigmadot = 1._dl/EFTeomX*(-2._dl*adotoa*(1._dl+EFTeomV)*sigma +etak&
& -dgs/k/(1._dl+EFTOmegaV) +EFTeomN/EFT_H0)
end if
! Original code:
! sigmadot = -2*adotoa*sigma-dgs/k+etak
! EFTCAMB MOD END
-4-