Replace: Before "pig = 32._dl/45/opacity*k*(sigma+vb)"
! EFTCAMB MOD START:
if (CP%EFTflag==0) then
gpres=gpres+ (grhog_t+grhor_t)/3 +grhov_t*w_lam
adotdota=(adotoa*adotoa-gpres)/2
end if
!... Original:
! gpres=gpres+ (grhog_t+grhor_t)/3 +grhov_t*w_lam
! adotdota=(adotoa*adotoa-gpres)/2
! EFTCAMB MOD END
if (CP%EFTflag==0) then
gpres=gpres+ (grhog_t+grhor_t)/3 +grhov_t*w_lam
adotdota=(adotoa*adotoa-gpres)/2
end if
!... Original:
! gpres=gpres+ (grhog_t+grhor_t)/3 +grhov_t*w_lam
! adotdota=(adotoa*adotoa-gpres)/2
! EFTCAMB MOD END
Replace: Before "!Once know slip, recompute qgdot, pig, pigdot"
! 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:
! ! Define shear derivative to first order
! 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:
! ! Define shear derivative to first order
! sigmadot = -2*adotoa*sigma-dgs/k+etak
! EFTCAMB MOD END
Add: after "end subroutine derivst"
! EFTCAMB MOD START: subroutine which sets initial conditions for pi.
subroutine EFTpiInitialConditions(y,EV,tau)
use ThermoData
implicit none
type(EvolutionVars) EV
real(dl) :: y(EV%nvar)
real(dl) k, k2,a, a2, etak, clxc, clxb, vb, clxr, clxg, &
grhob_t,grhor_t,grhoc_t,grhog_t,grhov_t, gpres, grho, adotoa, &
dgrho, z, dz, adotdota, opacity, cs2, dopacity, tau, dgq, qr, qg
! background quantities.
real(dl) EFT_H0,Hdot, Hdotdot, PiFieldScale
! 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
! perturbations quantities.
real(dl) EFTpiA, EFTpiB, EFTpiC, EFTpiD, EFTpiE, EFTpiAdot, EFTpiEdot ! pi field equation.
!
real(dl) temp1, temp2
k=EV%k_buf
k2=EV%k2_buf
a=y(1)
a2=a*a
etak=y(2)
clxc=y(3)
clxb=y(4)
vb=y(5)
! Compute expansion rate from: grho 8*pi*rho*a**2
grhob_t=grhob/a
grhoc_t=grhoc/a
grhor_t=grhornomass/a2
grhog_t=grhog/a2
grhov_t=grhov*EFTw(a,2)
! Get sound speed and ionisation fraction.
if (EV%TightCoupling) then
call thermo(tau,cs2,opacity,dopacity)
else
call thermo(tau,cs2,opacity)
end if
gpres=0
grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhov_t
adotoa=sqrt(grho/3)
dgrho=grhob_t*clxb+grhoc_t*clxc
z=(0.5_dl*dgrho/k + etak)/adotoa
dz= -2._dl*adotoa*z + etak
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
else
! Massless neutrinos
clxr=y(EV%r_ix)
qr =y(EV%r_ix+1)
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
else
! Photons
clxg=y(EV%g_ix)
qg=y(EV%g_ix+1)
end if
dgrho=grhob_t*clxb+grhoc_t*clxc + grhog_t*clxg+grhor_t*clxr
dgq=grhob_t*vb + grhog_t*qg+grhor_t*qr
z=(0.5_dl*dgrho/k + etak)/adotoa
dz= -2._dl*adotoa*z + etak
! EFT initial stuff
EFT_H0 = (CP%h0/c_EFT)*1000._dl
! Compute FRW background quantities
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))
! Compute EFT background quantities
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)
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
! Compute E which is the only one that we cannot compute outside.
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
! Compute Edot. This is just for background operators.
if ((CP%EFTflag==1.and.CP%PureEFTmodelAlpha1==0.and.CP%PureEFTmodelAlpha2==0.and.&
&CP%PureEFTmodelAlpha3==0.and.CP%PureEFTmodelAlpha4==0.and.CP%PureEFTmodelAlpha5==0.and.&
&CP%PureEFTmodelAlpha6==0).or.(CP%EFTflag==2)) then
EFTpiEdot = (EFTcdot +0.75_dl*(2._dl*adotoa*Hdot*a2*EFTOmegaP**2/(1._dl+EFTOmegaV)&
& +2._dl*adotoa**3*a**3*EFTOmegaP*EFTOmegaPP/(1._dl + EFTOmegaV)-a**3*adotoa**3*EFTOmegaP**3/(1._dl+EFTOmegaV)**2))*k*z &
& +(EFTc+0.75_dl*adotoa**2*(a2*EFTOmegaP**2)/(1._dl+EFTOmegaV))*k*dz +2._dl*adotoa*EFTpiE&
& +(-dgrho+(grhog_t*clxg+grhor_t*clxr))/4._dl*(+a*adotoa**2*EFTOmegaP +a*Hdot*EFTOmegaP+a2*adotoa**2*EFTOmegaPP &
& -a2*adotoa**2*EFTOmegaP**2/(1._dl+EFTOmegaV))/(1._dl+EFTOmegaV)&
& -adotoa/4._dl*a*EFTOmegaP/(1._dl+EFTOmegaV)*(grhob_t*(-k*(z+vb)-3._dl*adotoa*clxb) + grhoc_t*(-k*z -3._dl*adotoa*clxc))
else
EFTpiEdot = 0._dl
end if
! Write initial conditions for pi and pidot
! 1- Test initial conditions
y(EV%w_ix)= 0._dl
y(EV%w_ix+1)= 0._dl
! 2- Initial conditions on the source
if (EFTpiCfunction(a,k)+k*k*EFTpiDfunction(a,k)/=0) then
y(EV%w_ix)= -EFT_H0*EFTpiE/(EFTpiCfunction(a,k)+k*k*EFTpiDfunction(a,k))
y(EV%w_ix+1)= EFT_H0*(EFTpiE/(EFTpiCfunction(a,k)+k*k*EFTpiDfunction(a,k))**2*a*adotoa*(EFTpiCdotFunction(a,k)+k2*EFTpiDdotFunction(a,k))&
& - EFTpiEdot/(EFTpiCfunction(a,k)+k2*EFTpiDfunction(a,k)))
else
y(EV%w_ix)= 0._dl
y(EV%w_ix+1)= 0._dl
end if
end subroutine EFTpiInitialConditions
function EFTpiAfunction(a,k)
use ModelParams
implicit none
real(dl) EFTpiAfunction
real(dl), intent(IN) :: a, k
real(dl) a2,k2, grhob_t,grhor_t,grhoc_t,grhog_t,grhov_t, gpres, grho, adotoa, adotdota
! background quantities.
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, EFTpiA
a2 = a*a
k2 = k*k
EFT_H0 = (CP%h0/c_EFT)*1000._dl
! Compute expansion rate from: grho 8*pi*rho*a**2
grhob_t=grhob/a
grhoc_t=grhoc/a
grhor_t=grhornomass/a2
grhog_t=grhog/a2
grhov_t=grhov*EFTw(a,2)
grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhov_t
! Compute FRW background quantities
adotoa=sqrt(grho/3)
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))
! Compute EFT background quantities
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)
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))
else if (CP%EFTflag==2) then
EFTc = EFTcTemp(a,0)
end if
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
EFTpiAfunction = EFTpiA
return
end function EFTpiAfunction
function EFTpiBfunction(a,k)
use ModelParams
implicit none
real(dl) EFTpiBfunction
real(dl), intent(IN) :: a,k
real(dl) a2,k2,grhob_t,grhor_t,grhoc_t,grhog_t,grhov_t, gpres, grho, adotoa,adotdota
! background quantities.
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
! perturbations quantities.
real(dl) EFTpiB
a2 = a*a
k2 = k*k
EFT_H0 = (CP%h0/c_EFT)*1000._dl
! Compute expansion rate from: grho 8*pi*rho*a**2
grhob_t=grhob/a
grhoc_t=grhoc/a
grhor_t=grhornomass/a2
grhog_t=grhog/a2
grhov_t=grhov*EFTw(a,2)
grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhov_t
! Compute FRW background quantities
adotoa=sqrt(grho/3)
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))
! Compute EFT background quantities
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)
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
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)
EFTpiBfunction = EFTpiB
return
end function EFTpiBfunction
function EFTpiCfunction(a,k)
use ModelParams
implicit none
real(dl) EFTpiCfunction
real(dl), intent(IN) :: a,k
real(dl) a2,k2,grhob_t,grhor_t,grhoc_t,grhog_t,grhov_t, gpres, grho, adotoa,adotdota
! background quantities.
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
! perturbations quantities.
real(dl) EFTpiC
a2 = a*a
k2 = k*k
EFT_H0 = (CP%h0/c_EFT)*1000._dl
! Compute expansion rate from: grho 8*pi*rho*a**2
grhob_t=grhob/a
grhoc_t=grhoc/a
grhor_t=grhornomass/a2
grhog_t=grhog/a2
grhov_t=grhov*EFTw(a,2)
grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhov_t
! Compute FRW background quantities
adotoa=sqrt(grho/3)
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))
! Compute EFT background quantities
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)
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
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))
EFTpiCfunction = EFTpiC
return
end function EFTpiCfunction
function EFTpiDfunction(a,k)
use ModelParams
implicit none
real(dl) EFTpiDfunction
real(dl), intent(IN) :: a,k
real(dl) a2,k2,grhob_t,grhor_t,grhoc_t,grhog_t,grhov_t, gpres, grho, adotoa,adotdota
! background quantities.
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
! perturbations quantities.
real(dl) EFTpiD
a2 = a*a
k2 = k*k
EFT_H0 = (CP%h0/c_EFT)*1000._dl
! Compute expansion rate from: grho 8*pi*rho*a**2
grhob_t=grhob/a
grhoc_t=grhoc/a
grhor_t=grhornomass/a2
grhog_t=grhog/a2
grhov_t=grhov*EFTw(a,2)
grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhov_t
! Compute FRW background quantities
adotoa=sqrt(grho/3)
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))
! Compute EFT background quantities
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)
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
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))
EFTpiDfunction = EFTpiD
return
end function EFTpiDfunction
! EFTCAMB: now we take the numerical derivative of the functions appearing in the pi field equation
! to set up initial conditions. This part of code is called just once also because the accuracy
! of initial conditions is not an issue.
! EFTCAMB: numerical derivative of the function A
function EFTpiAdotFunction(a,k)
implicit none
real(dl) EFTpiAdotFunction, temp, err
real(dl), intent(IN) :: a,k
EFTpiAdotfunction = dfridr(EFTpiAfunction,a,k,0.1_dl*EFTturnonpi,err)
return
end function EFTpiAdotFunction
! EFTCAMB: numerical derivative of the function B
function EFTpiBdotFunction(a,k)
implicit none
real(dl) EFTpiBdotFunction, temp, err
real(dl), intent(IN) :: a,k
EFTpiBdotfunction = dfridr(EFTpiBfunction,a,k,0.1_dl*EFTturnonpi,err)
return
end function EFTpiBdotFunction
! EFTCAMB: numerical derivative of the function C
function EFTpiCdotFunction(a,k)
implicit none
real(dl) EFTpiCdotFunction, temp, err
real(dl), intent(IN) :: a,k
EFTpiCdotfunction = dfridr(EFTpiCfunction,a,k,0.1_dl*EFTturnonpi,err)
return
end function EFTpiCdotFunction
! EFTCAMB: numerical derivative of the function D
function EFTpiDdotFunction(a,k)
implicit none
real(dl) EFTpiDdotFunction, temp, err
real(dl), intent(IN) :: a,k
EFTpiDdotfunction = dfridr(EFTpiDfunction,a,k,0.1_dl*EFTturnonpi,err)
return
end function EFTpiDdotFunction
! EFTCAMB: algorithm to compute the numerical derivative. Modified to accept two function arguments.
function dfridr(func,x,k,h,err)
implicit none
integer,parameter :: ntab = 100
real(dl) dfridr,err,h,x,k,func
external func
real(dl), parameter :: CON=1.4_dl ! decrease of the stepsize.
real(dl), parameter :: CON2=CON*CON
real(dl), parameter :: BIG=1.d+30
real(dl), parameter :: SAFE=2._dl
integer i,j
real(dl) errt, fac, hh
real(dl), dimension(ntab,ntab) :: a
if (h.eq.0._dl) h = 1.d-8
hh=h
a(1,1)=(func(x+hh,k)-func(x-hh,k))/(2.0_dl*hh)
err=BIG
do 12 i=2,NTAB
hh=hh/CON
a(1,i)=(func(x+hh,k)-func(x-hh,k))/(2.0_dl*hh)
fac=CON2
do 11 j=2,i
a(j,i)=(a(j-1,i)*fac-a(j-1,i-1))/(fac-1._dl)
fac=CON2*fac
errt=max(abs(a(j,i)-a(j-1,i)),abs(a(j,i)-a(j-1,i-1)))
if (errt.le.err) then
err=errt
dfridr=a(j,i)
endif
11 continue
if(abs(a(i,i)-a(i-1,i-1)).ge.SAFE*err)return
12 continue
return
end function dfridr
! EFTCAMB MOD END
subroutine EFTpiInitialConditions(y,EV,tau)
use ThermoData
implicit none
type(EvolutionVars) EV
real(dl) :: y(EV%nvar)
real(dl) k, k2,a, a2, etak, clxc, clxb, vb, clxr, clxg, &
grhob_t,grhor_t,grhoc_t,grhog_t,grhov_t, gpres, grho, adotoa, &
dgrho, z, dz, adotdota, opacity, cs2, dopacity, tau, dgq, qr, qg
! background quantities.
real(dl) EFT_H0,Hdot, Hdotdot, PiFieldScale
! 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
! perturbations quantities.
real(dl) EFTpiA, EFTpiB, EFTpiC, EFTpiD, EFTpiE, EFTpiAdot, EFTpiEdot ! pi field equation.
!
real(dl) temp1, temp2
k=EV%k_buf
k2=EV%k2_buf
a=y(1)
a2=a*a
etak=y(2)
clxc=y(3)
clxb=y(4)
vb=y(5)
! Compute expansion rate from: grho 8*pi*rho*a**2
grhob_t=grhob/a
grhoc_t=grhoc/a
grhor_t=grhornomass/a2
grhog_t=grhog/a2
grhov_t=grhov*EFTw(a,2)
! Get sound speed and ionisation fraction.
if (EV%TightCoupling) then
call thermo(tau,cs2,opacity,dopacity)
else
call thermo(tau,cs2,opacity)
end if
gpres=0
grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhov_t
adotoa=sqrt(grho/3)
dgrho=grhob_t*clxb+grhoc_t*clxc
z=(0.5_dl*dgrho/k + etak)/adotoa
dz= -2._dl*adotoa*z + etak
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
else
! Massless neutrinos
clxr=y(EV%r_ix)
qr =y(EV%r_ix+1)
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
else
! Photons
clxg=y(EV%g_ix)
qg=y(EV%g_ix+1)
end if
dgrho=grhob_t*clxb+grhoc_t*clxc + grhog_t*clxg+grhor_t*clxr
dgq=grhob_t*vb + grhog_t*qg+grhor_t*qr
z=(0.5_dl*dgrho/k + etak)/adotoa
dz= -2._dl*adotoa*z + etak
! EFT initial stuff
EFT_H0 = (CP%h0/c_EFT)*1000._dl
! Compute FRW background quantities
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))
! Compute EFT background quantities
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)
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
! Compute E which is the only one that we cannot compute outside.
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
! Compute Edot. This is just for background operators.
if ((CP%EFTflag==1.and.CP%PureEFTmodelAlpha1==0.and.CP%PureEFTmodelAlpha2==0.and.&
&CP%PureEFTmodelAlpha3==0.and.CP%PureEFTmodelAlpha4==0.and.CP%PureEFTmodelAlpha5==0.and.&
&CP%PureEFTmodelAlpha6==0).or.(CP%EFTflag==2)) then
EFTpiEdot = (EFTcdot +0.75_dl*(2._dl*adotoa*Hdot*a2*EFTOmegaP**2/(1._dl+EFTOmegaV)&
& +2._dl*adotoa**3*a**3*EFTOmegaP*EFTOmegaPP/(1._dl + EFTOmegaV)-a**3*adotoa**3*EFTOmegaP**3/(1._dl+EFTOmegaV)**2))*k*z &
& +(EFTc+0.75_dl*adotoa**2*(a2*EFTOmegaP**2)/(1._dl+EFTOmegaV))*k*dz +2._dl*adotoa*EFTpiE&
& +(-dgrho+(grhog_t*clxg+grhor_t*clxr))/4._dl*(+a*adotoa**2*EFTOmegaP +a*Hdot*EFTOmegaP+a2*adotoa**2*EFTOmegaPP &
& -a2*adotoa**2*EFTOmegaP**2/(1._dl+EFTOmegaV))/(1._dl+EFTOmegaV)&
& -adotoa/4._dl*a*EFTOmegaP/(1._dl+EFTOmegaV)*(grhob_t*(-k*(z+vb)-3._dl*adotoa*clxb) + grhoc_t*(-k*z -3._dl*adotoa*clxc))
else
EFTpiEdot = 0._dl
end if
! Write initial conditions for pi and pidot
! 1- Test initial conditions
y(EV%w_ix)= 0._dl
y(EV%w_ix+1)= 0._dl
! 2- Initial conditions on the source
if (EFTpiCfunction(a,k)+k*k*EFTpiDfunction(a,k)/=0) then
y(EV%w_ix)= -EFT_H0*EFTpiE/(EFTpiCfunction(a,k)+k*k*EFTpiDfunction(a,k))
y(EV%w_ix+1)= EFT_H0*(EFTpiE/(EFTpiCfunction(a,k)+k*k*EFTpiDfunction(a,k))**2*a*adotoa*(EFTpiCdotFunction(a,k)+k2*EFTpiDdotFunction(a,k))&
& - EFTpiEdot/(EFTpiCfunction(a,k)+k2*EFTpiDfunction(a,k)))
else
y(EV%w_ix)= 0._dl
y(EV%w_ix+1)= 0._dl
end if
end subroutine EFTpiInitialConditions
function EFTpiAfunction(a,k)
use ModelParams
implicit none
real(dl) EFTpiAfunction
real(dl), intent(IN) :: a, k
real(dl) a2,k2, grhob_t,grhor_t,grhoc_t,grhog_t,grhov_t, gpres, grho, adotoa, adotdota
! background quantities.
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, EFTpiA
a2 = a*a
k2 = k*k
EFT_H0 = (CP%h0/c_EFT)*1000._dl
! Compute expansion rate from: grho 8*pi*rho*a**2
grhob_t=grhob/a
grhoc_t=grhoc/a
grhor_t=grhornomass/a2
grhog_t=grhog/a2
grhov_t=grhov*EFTw(a,2)
grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhov_t
! Compute FRW background quantities
adotoa=sqrt(grho/3)
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))
! Compute EFT background quantities
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)
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))
else if (CP%EFTflag==2) then
EFTc = EFTcTemp(a,0)
end if
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
EFTpiAfunction = EFTpiA
return
end function EFTpiAfunction
function EFTpiBfunction(a,k)
use ModelParams
implicit none
real(dl) EFTpiBfunction
real(dl), intent(IN) :: a,k
real(dl) a2,k2,grhob_t,grhor_t,grhoc_t,grhog_t,grhov_t, gpres, grho, adotoa,adotdota
! background quantities.
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
! perturbations quantities.
real(dl) EFTpiB
a2 = a*a
k2 = k*k
EFT_H0 = (CP%h0/c_EFT)*1000._dl
! Compute expansion rate from: grho 8*pi*rho*a**2
grhob_t=grhob/a
grhoc_t=grhoc/a
grhor_t=grhornomass/a2
grhog_t=grhog/a2
grhov_t=grhov*EFTw(a,2)
grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhov_t
! Compute FRW background quantities
adotoa=sqrt(grho/3)
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))
! Compute EFT background quantities
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)
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
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)
EFTpiBfunction = EFTpiB
return
end function EFTpiBfunction
function EFTpiCfunction(a,k)
use ModelParams
implicit none
real(dl) EFTpiCfunction
real(dl), intent(IN) :: a,k
real(dl) a2,k2,grhob_t,grhor_t,grhoc_t,grhog_t,grhov_t, gpres, grho, adotoa,adotdota
! background quantities.
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
! perturbations quantities.
real(dl) EFTpiC
a2 = a*a
k2 = k*k
EFT_H0 = (CP%h0/c_EFT)*1000._dl
! Compute expansion rate from: grho 8*pi*rho*a**2
grhob_t=grhob/a
grhoc_t=grhoc/a
grhor_t=grhornomass/a2
grhog_t=grhog/a2
grhov_t=grhov*EFTw(a,2)
grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhov_t
! Compute FRW background quantities
adotoa=sqrt(grho/3)
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))
! Compute EFT background quantities
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)
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
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))
EFTpiCfunction = EFTpiC
return
end function EFTpiCfunction
function EFTpiDfunction(a,k)
use ModelParams
implicit none
real(dl) EFTpiDfunction
real(dl), intent(IN) :: a,k
real(dl) a2,k2,grhob_t,grhor_t,grhoc_t,grhog_t,grhov_t, gpres, grho, adotoa,adotdota
! background quantities.
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
! perturbations quantities.
real(dl) EFTpiD
a2 = a*a
k2 = k*k
EFT_H0 = (CP%h0/c_EFT)*1000._dl
! Compute expansion rate from: grho 8*pi*rho*a**2
grhob_t=grhob/a
grhoc_t=grhoc/a
grhor_t=grhornomass/a2
grhog_t=grhog/a2
grhov_t=grhov*EFTw(a,2)
grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhov_t
! Compute FRW background quantities
adotoa=sqrt(grho/3)
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))
! Compute EFT background quantities
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)
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
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))
EFTpiDfunction = EFTpiD
return
end function EFTpiDfunction
! EFTCAMB: now we take the numerical derivative of the functions appearing in the pi field equation
! to set up initial conditions. This part of code is called just once also because the accuracy
! of initial conditions is not an issue.
! EFTCAMB: numerical derivative of the function A
function EFTpiAdotFunction(a,k)
implicit none
real(dl) EFTpiAdotFunction, temp, err
real(dl), intent(IN) :: a,k
EFTpiAdotfunction = dfridr(EFTpiAfunction,a,k,0.1_dl*EFTturnonpi,err)
return
end function EFTpiAdotFunction
! EFTCAMB: numerical derivative of the function B
function EFTpiBdotFunction(a,k)
implicit none
real(dl) EFTpiBdotFunction, temp, err
real(dl), intent(IN) :: a,k
EFTpiBdotfunction = dfridr(EFTpiBfunction,a,k,0.1_dl*EFTturnonpi,err)
return
end function EFTpiBdotFunction
! EFTCAMB: numerical derivative of the function C
function EFTpiCdotFunction(a,k)
implicit none
real(dl) EFTpiCdotFunction, temp, err
real(dl), intent(IN) :: a,k
EFTpiCdotfunction = dfridr(EFTpiCfunction,a,k,0.1_dl*EFTturnonpi,err)
return
end function EFTpiCdotFunction
! EFTCAMB: numerical derivative of the function D
function EFTpiDdotFunction(a,k)
implicit none
real(dl) EFTpiDdotFunction, temp, err
real(dl), intent(IN) :: a,k
EFTpiDdotfunction = dfridr(EFTpiDfunction,a,k,0.1_dl*EFTturnonpi,err)
return
end function EFTpiDdotFunction
! EFTCAMB: algorithm to compute the numerical derivative. Modified to accept two function arguments.
function dfridr(func,x,k,h,err)
implicit none
integer,parameter :: ntab = 100
real(dl) dfridr,err,h,x,k,func
external func
real(dl), parameter :: CON=1.4_dl ! decrease of the stepsize.
real(dl), parameter :: CON2=CON*CON
real(dl), parameter :: BIG=1.d+30
real(dl), parameter :: SAFE=2._dl
integer i,j
real(dl) errt, fac, hh
real(dl), dimension(ntab,ntab) :: a
if (h.eq.0._dl) h = 1.d-8
hh=h
a(1,1)=(func(x+hh,k)-func(x-hh,k))/(2.0_dl*hh)
err=BIG
do 12 i=2,NTAB
hh=hh/CON
a(1,i)=(func(x+hh,k)-func(x-hh,k))/(2.0_dl*hh)
fac=CON2
do 11 j=2,i
a(j,i)=(a(j-1,i)*fac-a(j-1,i-1))/(fac-1._dl)
fac=CON2*fac
errt=max(abs(a(j,i)-a(j-1,i)),abs(a(j,i)-a(j-1,i-1)))
if (errt.le.err) then
err=errt
dfridr=a(j,i)
endif
11 continue
if(abs(a(i,i)-a(i-1,i-1)).ge.SAFE*err)return
12 continue
return
end function dfridr
! EFTCAMB MOD END
-3-