There might exist some typos in this guide, please check it carefully with our original source codes. If you encounter any problems during your modifications, please contact us.
::::: subroutine derivst(EV,n,tau,ayt,aytprime) :::::
Add: After real(dl) pir, adotoa, rhonu, shear
! EFTCAMB MOD START: define quantities:
real(dl) EFT_H0, EFTOmegaV, EFTOmegaP, EFTAlpha4V, EFTAlpha4P, EFTAT, EFTBT, EFTDT
! EFTCAMB MOD END.
real(dl) EFT_H0, EFTOmegaV, EFTOmegaP, EFTAlpha4V, EFTAlpha4P, EFTAT, EFTBT, EFTDT
! EFTCAMB MOD END.
Replace:
! EFTCAMB MOD START: Compute the DE background density.
if (CP%EFTflag==0) then
if (w_lam==-1._dl) then
grhov_t=grhov*a2
else
grhov_t=grhov*a**(-1._dl-3._dl*w_lam)
end if
else if (CP%EFTflag/=0) then
grhov_t=grhov*EFTw(a,3)
EFT_H0 = (CP%h0/c_EFT)*1000._dl
end if
! Original code:
! if (w_lam==-1._dl) then
! grhov_t=grhov*a2
! else
! grhov_t=grhov*a**(-1-3*w_lam)
! end if
! EFTCAMB MOD END.
if (CP%EFTflag==0) then
if (w_lam==-1._dl) then
grhov_t=grhov*a2
else
grhov_t=grhov*a**(-1._dl-3._dl*w_lam)
end if
else if (CP%EFTflag/=0) then
grhov_t=grhov*EFTw(a,3)
EFT_H0 = (CP%h0/c_EFT)*1000._dl
end if
! Original code:
! if (w_lam==-1._dl) then
! grhov_t=grhov*a2
! else
! grhov_t=grhov*a**(-1-3*w_lam)
! end if
! EFTCAMB MOD END.
Replace:
! EFTCAMB MOD START: modified tensor equation of motion:
if (CP%EFTflag==0.or.a < EFTturnonpi) then
if (CP%flat) then
aytprime(3)=-2*adotoa*shear+k*Hchi-rhopi/k
else
aytprime(3)=-2*adotoa*shear+k*Hchi*(1+2*CP%curv/k2)-rhopi/k
endif
else if (CP%EFTflag/=0.and.a>=EFTturnonpi) then
EFTOmegaV = EFTOmega(a,0)
EFTOmegaP = EFTOmega(a,1)
EFTAlpha4V = EFTAlpha4(a,0)
EFTAlpha4P = EFTAlpha4(a,1)
EFTAT = 1._dl + EFTOmegaV - EFTAlpha4V**2
EFTBT = 2._dl*adotoa*(1._dl + EFTOmegaV - EFTAlpha4V**2 +0.5_dl*a*EFTOmegaP -a*EFTAlpha4V*EFTAlpha4P)
EFTDT = 1._dl + EFTOmegaV
aytprime(3)= -EFTBT/EFTAT*shear +EFTDT/EFTAT*k*Hchi -rhopi/k/EFTAT
end if
! Original code:
! if (CP%flat) then
! aytprime(3)=-2*adotoa*shear+k*Hchi-rhopi/k
! else
! aytprime(3)=-2*adotoa*shear+k*Hchi*(1+2*CP%curv/k2)-rhopi/k
! endif
! EFTCAMB MOD END.
if (CP%EFTflag==0.or.a < EFTturnonpi) then
if (CP%flat) then
aytprime(3)=-2*adotoa*shear+k*Hchi-rhopi/k
else
aytprime(3)=-2*adotoa*shear+k*Hchi*(1+2*CP%curv/k2)-rhopi/k
endif
else if (CP%EFTflag/=0.and.a>=EFTturnonpi) then
EFTOmegaV = EFTOmega(a,0)
EFTOmegaP = EFTOmega(a,1)
EFTAlpha4V = EFTAlpha4(a,0)
EFTAlpha4P = EFTAlpha4(a,1)
EFTAT = 1._dl + EFTOmegaV - EFTAlpha4V**2
EFTBT = 2._dl*adotoa*(1._dl + EFTOmegaV - EFTAlpha4V**2 +0.5_dl*a*EFTOmegaP -a*EFTAlpha4V*EFTAlpha4P)
EFTDT = 1._dl + EFTOmegaV
aytprime(3)= -EFTBT/EFTAT*shear +EFTDT/EFTAT*k*Hchi -rhopi/k/EFTAT
end if
! Original code:
! if (CP%flat) then
! aytprime(3)=-2*adotoa*shear+k*Hchi-rhopi/k
! else
! aytprime(3)=-2*adotoa*shear+k*Hchi*(1+2*CP%curv/k2)-rhopi/k
! endif
! 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
real(dl) wnu_arr(max_nu)
integer nu_i
! 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
real(dl) EFT_grhonu, EFT_gpinu, EFT_grhonudot, EFT_gpinudot, grhormass_t
real(dl) EFT_grhonudot_tot, EFT_grhonu_tot, EFT_gpinu_tot, EFT_gpinudot_tot
real(dl) EFT_dgpnu
! perturbations quantities.
real(dl) EFTpiA, EFTpiB, EFTpiC, EFTpiD, EFTpiE, 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,3)
! 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
dgrho=grhob_t*clxb+grhoc_t*clxc
dgq=grhob_t*vb
EFT_dgpnu = 0._dl
if (CP%Num_Nu_Massive > 0) then
call MassiveNuVars(EV,y,a,grho,gpres,dgrho,dgq, wnu_arr=wnu_arr, dgp=EFT_dgpnu)
end if
adotoa=sqrt(grho/3)
gpres=gpres + (grhog_t+grhor_t)/3._dl +EFTw(a,0)*grhov_t
adotdota=(adotoa*adotoa-gpres)/2._dl
Hdot =adotdota-adotoa**2._dl
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= dgrho + grhog_t*clxg+grhor_t*clxr
dgq= dgq + 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
! Massive neutrinos:
EFT_gpinudot_tot = 0._dl
if (CP%Num_Nu_Massive /= 0) then
do nu_i = 1, CP%Nu_mass_eigenstates
EFT_grhonu = 0._dl
EFT_gpinu = 0._dl
EFT_grhonudot = 0._dl
EFT_gpinudot = 0._dl
grhormass_t=grhormass(nu_i)/a**2
call Nu_background(a*nu_masses(nu_i),EFT_grhonu,EFT_gpinu)
EFT_grhonudot_tot = EFT_grhonudot_tot + grhormass_t*(Nu_drho(a*nu_masses(nu_i),adotoa, EFT_grhonu)&
& -4._dl*adotoa*EFT_grhonu)
EFT_gpinudot_tot = EFT_gpinudot_tot + grhormass_t*(Nu_pidot(a*nu_masses(nu_i),adotoa, EFT_gpinu)&
& -4._dl*adotoa*EFT_gpinu)
EFT_grhonu_tot = EFT_grhonu_tot + grhormass_t*EFT_grhonu
EFT_gpinu_tot = EFT_gpinu_tot + grhormass_t*EFT_gpinu
end do
end if
Hdotdot = 2._dl*adotoa*Hdot &
& + 0.5_dl*adotoa*(grhob_t + grhoc_t + 8._dl*(grhog_t+grhor_t)/3._dl)&
& + 0.5_dl*adotoa*grhov_t*((1._dl+EFTw(a,0))*(1._dl+3._dl*EFTw(a,0))-a*EFTw(a,1))&
& + adotoa/6._dl*EFT_grhonu_tot -0.5_dl*adotoa*EFT_gpinu_tot -0.5_dl*EFT_gpinudot_tot
! 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)
!8*pi*G*c*a^2
EFTc = (adotoa*adotoa - Hdot)*(EFTOmegaV + a*EFTOmegaP*0.5_dl) &
& - 0.5_dl*a2*adotoa*adotoa*EFTOmegaPP&
& + 0.5_dl*grhov_t*(1._dl+EFTw(a,0))
!8*pi*G*Lambda*a^2
EFTLambda = +EFTw(a,0)*grhov_t &
&-EFTOmegaV*(2._dl*Hdot+adotoa**2._dl) &
&-a*EFTOmegaP*(2._dl*adotoa**2._dl + Hdot) &
&-a2*adotoa**2._dl*EFTOmegaPP
!EFT C DOT: 8*pi*G*cdot*a^2
EFTcdot = -EFTOmegaV*(Hdotdot-4._dl*adotoa*Hdot+2._dl*adotoa*adotoa*adotoa) &
& + 0.5_dl*a*EFTOmegaP*(-Hdotdot+adotoa*Hdot+adotoa*adotoa*adotoa)&
& +0.5_dl*a2*adotoa*EFTOmegaPP*(adotoa*adotoa-3._dl*Hdot)&
& -0.5_dl*a*a2*adotoa*adotoa*adotoa*EFTOmegaPPP&
& +0.5_dl*adotoa*grhov_t*(-3._dl*(1._dl+EFTw(a,0))**2 + a*EFTw(a,1))
!EFT LAMBDA DOT: 8*pi*G*Ldot*a^2
EFTLambdadot = -2._dl*EFTOmegaV*(Hdotdot-adotoa*Hdot-adotoa*adotoa*adotoa)&
& - a*EFTOmegaP*(4._dl*adotoa*Hdot+Hdotdot)&
& -a2*EFTOmegaPP*adotoa*(3._dl*Hdot+2._dl*adotoa*adotoa)&
& -a*a2*EFTOmegaPPP*adotoa*adotoa*adotoa&
& +grhov_t*adotoa*(a*EFTw(a,1)-3._dl*EFTw(a,0)*(1._dl+EFTw(a,0)))
if (CP%EFTflag==2.and.CP%DesignerEFTmodel/=2) then
EFTc = EFTcTemp(a,0)
EFTcdot = EFTcTemp(a,1)
EFTLambda = EFTLambdaTemp(a,0)
EFTLambdadot = EFTLambdaTemp(a,1)
end if
EFTgrhoq = 2._dl*EFTc -EFTLambda -3._dl*a*adotoa*adotoa*EFTOmegaP
EFTgpresq = EFTLambda + a2*adotoa*adotoa*EFTOmegaPP&
& +a*EFTOmegaP*(Hdot + 2._dl*adotoa*adotoa)
EFTgrhodotq = -3._dl*adotoa*(EFTgrhoq+EFTgpresq) + 3._dl*a*adotoa**3._dl*EFTOmegaP
EFTgpresdotq = EFTLambdadot +a2*a*adotoa**3*EFTOmegaPPP + 3._dl*a2*adotoa*Hdot*EFTOmegaPP&
& +a*EFTOmegaP*Hdotdot +3._dl*a*adotoa*Hdot*EFTOmegaP +2._dl*a2*adotoa**3*EFTOmegaPP&
& -2._dl*a*adotoa**3*EFTOmegaP
! 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 +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
! 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
real(dl) EFT_grhonu, EFT_gpinu, EFT_grhonudot, EFT_gpinudot, grhormass_t
real(dl) EFT_grhonudot_tot, EFT_grhonu_tot, EFT_gpinu_tot, EFT_gpinudot_tot
! 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
integer nu_i
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,3)
grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhov_t
! Massive neutrinos
if (CP%Num_Nu_Massive /= 0) then
do nu_i = 1, CP%Nu_mass_eigenstates
EFT_grhonu = 0._dl
EFT_gpinu = 0._dl
EFT_grhonudot = 0._dl
EFT_gpinudot = 0._dl
grhormass_t=grhormass(nu_i)/a**2
call Nu_background(a*nu_masses(nu_i),EFT_grhonu,EFT_gpinu)
EFT_grhonu_tot = EFT_grhonu_tot + grhormass_t*EFT_grhonu
EFT_gpinu_tot = EFT_gpinu_tot + grhormass_t*EFT_gpinu
end do
end if
grho = grho + EFT_grhonu_tot
gpres = +EFT_gpinu_tot +(grhog_t+grhor_t)/3._dl +EFTw(a,0)*grhov_t
! 3) Hubble and derivatives:
adotoa=sqrt(grho/3)
adotdota=(adotoa*adotoa-gpres)/2._dl
Hdot =adotdota-adotoa**2._dl
EFT_gpinu_tot = 0._dl
EFT_grhonu_tot = 0._dl
EFT_gpinudot_tot = 0._dl
if (CP%Num_Nu_Massive /= 0) then
do nu_i = 1, CP%Nu_mass_eigenstates
EFT_grhonu = 0._dl
EFT_gpinu = 0._dl
EFT_grhonudot = 0._dl
EFT_gpinudot = 0._dl
grhormass_t=grhormass(nu_i)/a**2
call Nu_background(a*nu_masses(nu_i),EFT_grhonu,EFT_gpinu)
EFT_grhonu_tot = EFT_grhonu_tot + grhormass_t*EFT_grhonu
EFT_gpinu_tot = EFT_gpinu_tot + grhormass_t*EFT_gpinu
EFT_gpinudot_tot = EFT_gpinudot_tot + grhormass_t*(Nu_pidot(a*nu_masses(nu_i),adotoa, EFT_gpinu)&
& -4._dl*adotoa*EFT_gpinu)
end do
end if
Hdotdot = 2._dl*adotoa*Hdot &
& + 0.5_dl*adotoa*(grhob_t + grhoc_t + 8._dl*(grhog_t+grhor_t)/3._dl)&
& + 0.5_dl*adotoa*grhov_t*((1._dl+EFTw(a,0))*(1._dl+3._dl*EFTw(a,0))-a*EFTw(a,1))&
& + adotoa/6._dl*EFT_grhonu_tot -0.5_dl*adotoa*EFT_gpinu_tot -0.5_dl*EFT_gpinudot_tot
! 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)
!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))
if (CP%EFTflag==2.and.CP%DesignerEFTmodel/=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
real(dl) EFT_grhonu, EFT_gpinu, EFT_grhonudot, EFT_gpinudot, grhormass_t
real(dl) EFT_grhonudot_tot, EFT_grhonu_tot, EFT_gpinu_tot, EFT_gpinudot_tot
! 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
integer nu_i
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,3)
grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhov_t
! Massive neutrinos
if (CP%Num_Nu_Massive /= 0) then
do nu_i = 1, CP%Nu_mass_eigenstates
EFT_grhonu = 0._dl
EFT_gpinu = 0._dl
EFT_grhonudot = 0._dl
EFT_gpinudot = 0._dl
grhormass_t=grhormass(nu_i)/a**2
call Nu_background(a*nu_masses(nu_i),EFT_grhonu,EFT_gpinu)
EFT_grhonu_tot = EFT_grhonu_tot + grhormass_t*EFT_grhonu
EFT_gpinu_tot = EFT_gpinu_tot + grhormass_t*EFT_gpinu
end do
end if
grho = grho + EFT_grhonu_tot
gpres = +EFT_gpinu_tot +(grhog_t+grhor_t)/3._dl +EFTw(a,0)*grhov_t
! 3) Hubble and derivatives:
adotoa=sqrt(grho/3)
adotdota=(adotoa*adotoa-gpres)/2._dl
Hdot =adotdota-adotoa**2._dl
EFT_gpinu_tot = 0._dl
EFT_grhonu_tot = 0._dl
EFT_gpinudot_tot = 0._dl
if (CP%Num_Nu_Massive /= 0) then
do nu_i = 1, CP%Nu_mass_eigenstates
EFT_grhonu = 0._dl
EFT_gpinu = 0._dl
EFT_grhonudot = 0._dl
EFT_gpinudot = 0._dl
grhormass_t=grhormass(nu_i)/a**2
call Nu_background(a*nu_masses(nu_i),EFT_grhonu,EFT_gpinu)
EFT_grhonu_tot = EFT_grhonu_tot + grhormass_t*EFT_grhonu
EFT_gpinu_tot = EFT_gpinu_tot + grhormass_t*EFT_gpinu
EFT_gpinudot_tot = EFT_gpinudot_tot + grhormass_t*(Nu_pidot(a*nu_masses(nu_i),adotoa, EFT_gpinu)&
& -4._dl*adotoa*EFT_gpinu)
end do
end if
Hdotdot = 2._dl*adotoa*Hdot &
& + 0.5_dl*adotoa*(grhob_t + grhoc_t + 8._dl*(grhog_t+grhor_t)/3._dl)&
& + 0.5_dl*adotoa*grhov_t*((1._dl+EFTw(a,0))*(1._dl+3._dl*EFTw(a,0))-a*EFTw(a,1))&
& + adotoa/6._dl*EFT_grhonu_tot -0.5_dl*adotoa*EFT_gpinu_tot -0.5_dl*EFT_gpinudot_tot
! 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)
!8*pi*G*c*a^2
EFTc = (adotoa*adotoa - Hdot)*(EFTOmegaV + a*EFTOmegaP*0.5_dl) &
& - 0.5_dl*a2*adotoa*adotoa*EFTOmegaPP&
& + 0.5_dl*grhov_t*(1._dl+EFTw(a,0))
!8*pi*G*Lambda*a^2
EFTLambda = +EFTw(a,0)*grhov_t &
&-EFTOmegaV*(2._dl*Hdot+adotoa**2._dl) &
&-a*EFTOmegaP*(2._dl*adotoa**2._dl + Hdot) &
&-a2*adotoa**2._dl*EFTOmegaPP
!EFT C DOT: 8*pi*G*cdot*a^2
EFTcdot = -EFTOmegaV*(Hdotdot-4._dl*adotoa*Hdot+2._dl*adotoa*adotoa*adotoa) &
& + 0.5_dl*a*EFTOmegaP*(-Hdotdot+adotoa*Hdot+adotoa*adotoa*adotoa)&
& +0.5_dl*a2*adotoa*EFTOmegaPP*(adotoa*adotoa-3._dl*Hdot)&
& -0.5_dl*a*a2*adotoa*adotoa*adotoa*EFTOmegaPPP&
& +0.5_dl*adotoa*grhov_t*(-3._dl*(1._dl+EFTw(a,0))**2 + a*EFTw(a,1))
!EFT LAMBDA DOT: 8*pi*G*Ldot*a^2
EFTLambdadot = -2._dl*EFTOmegaV*(Hdotdot-adotoa*Hdot-adotoa*adotoa*adotoa)&
& - a*EFTOmegaP*(4._dl*adotoa*Hdot+Hdotdot)&
& -a2*EFTOmegaPP*adotoa*(3._dl*Hdot+2._dl*adotoa*adotoa)&
& -a*a2*EFTOmegaPPP*adotoa*adotoa*adotoa&
& +grhov_t*adotoa*(a*EFTw(a,1)-3._dl*EFTw(a,0)*(1._dl+EFTw(a,0)))
if (CP%EFTflag==2.and.CP%DesignerEFTmodel/=2) then
EFTc = EFTcTemp(a,0)
EFTcdot = EFTcTemp(a,1)
EFTLambda = EFTLambdaTemp(a,0)
EFTLambdadot = EFTLambdaTemp(a,1)
end if
EFTgrhoq = 2._dl*EFTc -EFTLambda -3._dl*a*adotoa*adotoa*EFTOmegaP
EFTgpresq = EFTLambda + a2*adotoa*adotoa*EFTOmegaPP&
& +a*EFTOmegaP*(Hdot + 2._dl*adotoa*adotoa)
EFTgrhodotq = -3._dl*adotoa*(EFTgrhoq+EFTgpresq) + 3._dl*a*adotoa**3._dl*EFTOmegaP
EFTgpresdotq = EFTLambdadot +a2*a*adotoa**3*EFTOmegaPPP + 3._dl*a2*adotoa*Hdot*EFTOmegaPP&
& +a*EFTOmegaP*Hdotdot +3._dl*a*adotoa*Hdot*EFTOmegaP +2._dl*a2*adotoa**3*EFTOmegaPP&
& -2._dl*a*adotoa**3*EFTOmegaP
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
real(dl) EFT_grhonu, EFT_gpinu, EFT_grhonudot, EFT_gpinudot, grhormass_t
real(dl) EFT_grhonudot_tot, EFT_grhonu_tot, EFT_gpinu_tot, EFT_gpinudot_tot
! 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
integer nu_i
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,3)
grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhov_t
! Massive neutrinos
if (CP%Num_Nu_Massive /= 0) then
do nu_i = 1, CP%Nu_mass_eigenstates
EFT_grhonu = 0._dl
EFT_gpinu = 0._dl
EFT_grhonudot = 0._dl
EFT_gpinudot = 0._dl
grhormass_t=grhormass(nu_i)/a**2
call Nu_background(a*nu_masses(nu_i),EFT_grhonu,EFT_gpinu)
EFT_grhonu_tot = EFT_grhonu_tot + grhormass_t*EFT_grhonu
EFT_gpinu_tot = EFT_gpinu_tot + grhormass_t*EFT_gpinu
end do
end if
grho = grho + EFT_grhonu_tot
gpres = +EFT_gpinu_tot +(grhog_t+grhor_t)/3._dl +EFTw(a,0)*grhov_t
! 3) Hubble and derivatives:
adotoa=sqrt(grho/3)
adotdota=(adotoa*adotoa-gpres)/2._dl
Hdot =adotdota-adotoa**2._dl
EFT_gpinu_tot = 0._dl
EFT_grhonu_tot = 0._dl
EFT_gpinudot_tot = 0._dl
if (CP%Num_Nu_Massive /= 0) then
do nu_i = 1, CP%Nu_mass_eigenstates
EFT_grhonu = 0._dl
EFT_gpinu = 0._dl
EFT_grhonudot = 0._dl
EFT_gpinudot = 0._dl
grhormass_t=grhormass(nu_i)/a**2
call Nu_background(a*nu_masses(nu_i),EFT_grhonu,EFT_gpinu)
EFT_grhonu_tot = EFT_grhonu_tot + grhormass_t*EFT_grhonu
EFT_gpinu_tot = EFT_gpinu_tot + grhormass_t*EFT_gpinu
EFT_gpinudot_tot = EFT_gpinudot_tot + grhormass_t*(Nu_pidot(a*nu_masses(nu_i),adotoa, EFT_gpinu)&
& -4._dl*adotoa*EFT_gpinu)
end do
end if
Hdotdot = 2._dl*adotoa*Hdot &
& + 0.5_dl*adotoa*(grhob_t + grhoc_t + 8._dl*(grhog_t+grhor_t)/3._dl)&
& + 0.5_dl*adotoa*grhov_t*((1._dl+EFTw(a,0))*(1._dl+3._dl*EFTw(a,0))-a*EFTw(a,1))&
& + adotoa/6._dl*EFT_grhonu_tot -0.5_dl*adotoa*EFT_gpinu_tot -0.5_dl*EFT_gpinudot_tot
! 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)
!8*pi*G*c*a^2
EFTc = (adotoa*adotoa - Hdot)*(EFTOmegaV + a*EFTOmegaP*0.5_dl) &
& - 0.5_dl*a2*adotoa*adotoa*EFTOmegaPP&
& + 0.5_dl*grhov_t*(1._dl+EFTw(a,0))
!8*pi*G*Lambda*a^2
EFTLambda = +EFTw(a,0)*grhov_t &
&-EFTOmegaV*(2._dl*Hdot+adotoa**2._dl) &
&-a*EFTOmegaP*(2._dl*adotoa**2._dl + Hdot) &
&-a2*adotoa**2._dl*EFTOmegaPP
!EFT C DOT: 8*pi*G*cdot*a^2
EFTcdot = -EFTOmegaV*(Hdotdot-4._dl*adotoa*Hdot+2._dl*adotoa*adotoa*adotoa) &
& + 0.5_dl*a*EFTOmegaP*(-Hdotdot+adotoa*Hdot+adotoa*adotoa*adotoa)&
& +0.5_dl*a2*adotoa*EFTOmegaPP*(adotoa*adotoa-3._dl*Hdot)&
& -0.5_dl*a*a2*adotoa*adotoa*adotoa*EFTOmegaPPP&
& +0.5_dl*adotoa*grhov_t*(-3._dl*(1._dl+EFTw(a,0))**2 + a*EFTw(a,1))
!EFT LAMBDA DOT: 8*pi*G*Ldot*a^2
EFTLambdadot = -2._dl*EFTOmegaV*(Hdotdot-adotoa*Hdot-adotoa*adotoa*adotoa)&
& - a*EFTOmegaP*(4._dl*adotoa*Hdot+Hdotdot)&
& -a2*EFTOmegaPP*adotoa*(3._dl*Hdot+2._dl*adotoa*adotoa)&
& -a*a2*EFTOmegaPPP*adotoa*adotoa*adotoa&
& +grhov_t*adotoa*(a*EFTw(a,1)-3._dl*EFTw(a,0)*(1._dl+EFTw(a,0)))
if (CP%EFTflag==2.and.CP%DesignerEFTmodel/=2) then
EFTc = EFTcTemp(a,0)
EFTcdot = EFTcTemp(a,1)
EFTLambda = EFTLambdaTemp(a,0)
EFTLambdadot = EFTLambdaTemp(a,1)
end if
EFTgrhoq = 2._dl*EFTc -EFTLambda -3._dl*a*adotoa*adotoa*EFTOmegaP
EFTgpresq = EFTLambda + a2*adotoa*adotoa*EFTOmegaPP&
& +a*EFTOmegaP*(Hdot + 2._dl*adotoa*adotoa)
EFTgrhodotq = -3._dl*adotoa*(EFTgrhoq+EFTgpresq) + 3._dl*a*adotoa**3._dl*EFTOmegaP
EFTgpresdotq = EFTLambdadot +a2*a*adotoa**3*EFTOmegaPPP + 3._dl*a2*adotoa*Hdot*EFTOmegaPP&
& +a*EFTOmegaP*Hdotdot +3._dl*a*adotoa*Hdot*EFTOmegaP +2._dl*a2*adotoa**3*EFTOmegaPP&
& -2._dl*a*adotoa**3*EFTOmegaP
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
real(dl) EFT_grhonu, EFT_gpinu, EFT_grhonudot, EFT_gpinudot, grhormass_t
real(dl) EFT_grhonudot_tot, EFT_grhonu_tot, EFT_gpinu_tot, EFT_gpinudot_tot
! 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
integer nu_i
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,3)
grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhov_t
! Massive neutrinos
if (CP%Num_Nu_Massive /= 0) then
do nu_i = 1, CP%Nu_mass_eigenstates
EFT_grhonu = 0._dl
EFT_gpinu = 0._dl
EFT_grhonudot = 0._dl
EFT_gpinudot = 0._dl
grhormass_t=grhormass(nu_i)/a**2
call Nu_background(a*nu_masses(nu_i),EFT_grhonu,EFT_gpinu)
EFT_grhonu_tot = EFT_grhonu_tot + grhormass_t*EFT_grhonu
EFT_gpinu_tot = EFT_gpinu_tot + grhormass_t*EFT_gpinu
end do
end if
grho = grho + EFT_grhonu_tot
gpres = +EFT_gpinu_tot +(grhog_t+grhor_t)/3._dl +EFTw(a,0)*grhov_t
! 3) Hubble and derivatives:
adotoa=sqrt(grho/3)
adotdota=(adotoa*adotoa-gpres)/2._dl
Hdot =adotdota-adotoa**2._dl
EFT_gpinu_tot = 0._dl
EFT_grhonu_tot = 0._dl
EFT_gpinudot_tot = 0._dl
if (CP%Num_Nu_Massive /= 0) then
do nu_i = 1, CP%Nu_mass_eigenstates
EFT_grhonu = 0._dl
EFT_gpinu = 0._dl
EFT_grhonudot = 0._dl
EFT_gpinudot = 0._dl
grhormass_t=grhormass(nu_i)/a**2
call Nu_background(a*nu_masses(nu_i),EFT_grhonu,EFT_gpinu)
EFT_grhonu_tot = EFT_grhonu_tot + grhormass_t*EFT_grhonu
EFT_gpinu_tot = EFT_gpinu_tot + grhormass_t*EFT_gpinu
EFT_gpinudot_tot = EFT_gpinudot_tot + grhormass_t*(Nu_pidot(a*nu_masses(nu_i),adotoa, EFT_gpinu)&
& -4._dl*adotoa*EFT_gpinu)
end do
end if
Hdotdot = 2._dl*adotoa*Hdot &
& + 0.5_dl*adotoa*(grhob_t + grhoc_t + 8._dl*(grhog_t+grhor_t)/3._dl)&
& + 0.5_dl*adotoa*grhov_t*((1._dl+EFTw(a,0))*(1._dl+3._dl*EFTw(a,0))-a*EFTw(a,1))&
& + adotoa/6._dl*EFT_grhonu_tot -0.5_dl*adotoa*EFT_gpinu_tot -0.5_dl*EFT_gpinudot_tot
! 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)
!8*pi*G*c*a^2
EFTc = (adotoa*adotoa - Hdot)*(EFTOmegaV + a*EFTOmegaP*0.5_dl) &
& - 0.5_dl*a2*adotoa*adotoa*EFTOmegaPP&
& + 0.5_dl*grhov_t*(1._dl+EFTw(a,0))
!8*pi*G*Lambda*a^2
EFTLambda = +EFTw(a,0)*grhov_t &
&-EFTOmegaV*(2._dl*Hdot+adotoa**2._dl) &
&-a*EFTOmegaP*(2._dl*adotoa**2._dl + Hdot) &
&-a2*adotoa**2._dl*EFTOmegaPP
!EFT C DOT: 8*pi*G*cdot*a^2
EFTcdot = -EFTOmegaV*(Hdotdot-4._dl*adotoa*Hdot+2._dl*adotoa*adotoa*adotoa) &
& + 0.5_dl*a*EFTOmegaP*(-Hdotdot+adotoa*Hdot+adotoa*adotoa*adotoa)&
& +0.5_dl*a2*adotoa*EFTOmegaPP*(adotoa*adotoa-3._dl*Hdot)&
& -0.5_dl*a*a2*adotoa*adotoa*adotoa*EFTOmegaPPP&
& +0.5_dl*adotoa*grhov_t*(-3._dl*(1._dl+EFTw(a,0))**2 + a*EFTw(a,1))
!EFT LAMBDA DOT: 8*pi*G*Ldot*a^2
EFTLambdadot = -2._dl*EFTOmegaV*(Hdotdot-adotoa*Hdot-adotoa*adotoa*adotoa)&
& - a*EFTOmegaP*(4._dl*adotoa*Hdot+Hdotdot)&
& -a2*EFTOmegaPP*adotoa*(3._dl*Hdot+2._dl*adotoa*adotoa)&
& -a*a2*EFTOmegaPPP*adotoa*adotoa*adotoa&
& +grhov_t*adotoa*(a*EFTw(a,1)-3._dl*EFTw(a,0)*(1._dl+EFTw(a,0)))
if (CP%EFTflag==2.and.CP%DesignerEFTmodel/=2) then
EFTc = EFTcTemp(a,0)
EFTcdot = EFTcTemp(a,1)
EFTLambda = EFTLambdaTemp(a,0)
EFTLambdadot = EFTLambdaTemp(a,1)
end if
EFTgrhoq = 2._dl*EFTc -EFTLambda -3._dl*a*adotoa*adotoa*EFTOmegaP
EFTgpresq = EFTLambda + a2*adotoa*adotoa*EFTOmegaPP&
& +a*EFTOmegaP*(Hdot + 2._dl*adotoa*adotoa)
EFTgrhodotq = -3._dl*adotoa*(EFTgrhoq+EFTgpresq) + 3._dl*a*adotoa**3._dl*EFTOmegaP
EFTgpresdotq = EFTLambdadot +a2*a*adotoa**3*EFTOmegaPPP + 3._dl*a2*adotoa*Hdot*EFTOmegaPP&
& +a*EFTOmegaP*Hdotdot +3._dl*a*adotoa*Hdot*EFTOmegaP +2._dl*a2*adotoa**3*EFTOmegaPP&
& -2._dl*a*adotoa**3*EFTOmegaP
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
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
real(dl) wnu_arr(max_nu)
integer nu_i
! 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
real(dl) EFT_grhonu, EFT_gpinu, EFT_grhonudot, EFT_gpinudot, grhormass_t
real(dl) EFT_grhonudot_tot, EFT_grhonu_tot, EFT_gpinu_tot, EFT_gpinudot_tot
real(dl) EFT_dgpnu
! perturbations quantities.
real(dl) EFTpiA, EFTpiB, EFTpiC, EFTpiD, EFTpiE, 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,3)
! 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
dgrho=grhob_t*clxb+grhoc_t*clxc
dgq=grhob_t*vb
EFT_dgpnu = 0._dl
if (CP%Num_Nu_Massive > 0) then
call MassiveNuVars(EV,y,a,grho,gpres,dgrho,dgq, wnu_arr=wnu_arr, dgp=EFT_dgpnu)
end if
adotoa=sqrt(grho/3)
gpres=gpres + (grhog_t+grhor_t)/3._dl +EFTw(a,0)*grhov_t
adotdota=(adotoa*adotoa-gpres)/2._dl
Hdot =adotdota-adotoa**2._dl
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= dgrho + grhog_t*clxg+grhor_t*clxr
dgq= dgq + 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
! Massive neutrinos:
EFT_gpinudot_tot = 0._dl
if (CP%Num_Nu_Massive /= 0) then
do nu_i = 1, CP%Nu_mass_eigenstates
EFT_grhonu = 0._dl
EFT_gpinu = 0._dl
EFT_grhonudot = 0._dl
EFT_gpinudot = 0._dl
grhormass_t=grhormass(nu_i)/a**2
call Nu_background(a*nu_masses(nu_i),EFT_grhonu,EFT_gpinu)
EFT_grhonudot_tot = EFT_grhonudot_tot + grhormass_t*(Nu_drho(a*nu_masses(nu_i),adotoa, EFT_grhonu)&
& -4._dl*adotoa*EFT_grhonu)
EFT_gpinudot_tot = EFT_gpinudot_tot + grhormass_t*(Nu_pidot(a*nu_masses(nu_i),adotoa, EFT_gpinu)&
& -4._dl*adotoa*EFT_gpinu)
EFT_grhonu_tot = EFT_grhonu_tot + grhormass_t*EFT_grhonu
EFT_gpinu_tot = EFT_gpinu_tot + grhormass_t*EFT_gpinu
end do
end if
Hdotdot = 2._dl*adotoa*Hdot &
& + 0.5_dl*adotoa*(grhob_t + grhoc_t + 8._dl*(grhog_t+grhor_t)/3._dl)&
& + 0.5_dl*adotoa*grhov_t*((1._dl+EFTw(a,0))*(1._dl+3._dl*EFTw(a,0))-a*EFTw(a,1))&
& + adotoa/6._dl*EFT_grhonu_tot -0.5_dl*adotoa*EFT_gpinu_tot -0.5_dl*EFT_gpinudot_tot
! 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)
!8*pi*G*c*a^2
EFTc = (adotoa*adotoa - Hdot)*(EFTOmegaV + a*EFTOmegaP*0.5_dl) &
& - 0.5_dl*a2*adotoa*adotoa*EFTOmegaPP&
& + 0.5_dl*grhov_t*(1._dl+EFTw(a,0))
!8*pi*G*Lambda*a^2
EFTLambda = +EFTw(a,0)*grhov_t &
&-EFTOmegaV*(2._dl*Hdot+adotoa**2._dl) &
&-a*EFTOmegaP*(2._dl*adotoa**2._dl + Hdot) &
&-a2*adotoa**2._dl*EFTOmegaPP
!EFT C DOT: 8*pi*G*cdot*a^2
EFTcdot = -EFTOmegaV*(Hdotdot-4._dl*adotoa*Hdot+2._dl*adotoa*adotoa*adotoa) &
& + 0.5_dl*a*EFTOmegaP*(-Hdotdot+adotoa*Hdot+adotoa*adotoa*adotoa)&
& +0.5_dl*a2*adotoa*EFTOmegaPP*(adotoa*adotoa-3._dl*Hdot)&
& -0.5_dl*a*a2*adotoa*adotoa*adotoa*EFTOmegaPPP&
& +0.5_dl*adotoa*grhov_t*(-3._dl*(1._dl+EFTw(a,0))**2 + a*EFTw(a,1))
!EFT LAMBDA DOT: 8*pi*G*Ldot*a^2
EFTLambdadot = -2._dl*EFTOmegaV*(Hdotdot-adotoa*Hdot-adotoa*adotoa*adotoa)&
& - a*EFTOmegaP*(4._dl*adotoa*Hdot+Hdotdot)&
& -a2*EFTOmegaPP*adotoa*(3._dl*Hdot+2._dl*adotoa*adotoa)&
& -a*a2*EFTOmegaPPP*adotoa*adotoa*adotoa&
& +grhov_t*adotoa*(a*EFTw(a,1)-3._dl*EFTw(a,0)*(1._dl+EFTw(a,0)))
if (CP%EFTflag==2.and.CP%DesignerEFTmodel/=2) then
EFTc = EFTcTemp(a,0)
EFTcdot = EFTcTemp(a,1)
EFTLambda = EFTLambdaTemp(a,0)
EFTLambdadot = EFTLambdaTemp(a,1)
end if
EFTgrhoq = 2._dl*EFTc -EFTLambda -3._dl*a*adotoa*adotoa*EFTOmegaP
EFTgpresq = EFTLambda + a2*adotoa*adotoa*EFTOmegaPP&
& +a*EFTOmegaP*(Hdot + 2._dl*adotoa*adotoa)
EFTgrhodotq = -3._dl*adotoa*(EFTgrhoq+EFTgpresq) + 3._dl*a*adotoa**3._dl*EFTOmegaP
EFTgpresdotq = EFTLambdadot +a2*a*adotoa**3*EFTOmegaPPP + 3._dl*a2*adotoa*Hdot*EFTOmegaPP&
& +a*EFTOmegaP*Hdotdot +3._dl*a*adotoa*Hdot*EFTOmegaP +2._dl*a2*adotoa**3*EFTOmegaPP&
& -2._dl*a*adotoa**3*EFTOmegaP
! 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 +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
! 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
real(dl) EFT_grhonu, EFT_gpinu, EFT_grhonudot, EFT_gpinudot, grhormass_t
real(dl) EFT_grhonudot_tot, EFT_grhonu_tot, EFT_gpinu_tot, EFT_gpinudot_tot
! 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
integer nu_i
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,3)
grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhov_t
! Massive neutrinos
if (CP%Num_Nu_Massive /= 0) then
do nu_i = 1, CP%Nu_mass_eigenstates
EFT_grhonu = 0._dl
EFT_gpinu = 0._dl
EFT_grhonudot = 0._dl
EFT_gpinudot = 0._dl
grhormass_t=grhormass(nu_i)/a**2
call Nu_background(a*nu_masses(nu_i),EFT_grhonu,EFT_gpinu)
EFT_grhonu_tot = EFT_grhonu_tot + grhormass_t*EFT_grhonu
EFT_gpinu_tot = EFT_gpinu_tot + grhormass_t*EFT_gpinu
end do
end if
grho = grho + EFT_grhonu_tot
gpres = +EFT_gpinu_tot +(grhog_t+grhor_t)/3._dl +EFTw(a,0)*grhov_t
! 3) Hubble and derivatives:
adotoa=sqrt(grho/3)
adotdota=(adotoa*adotoa-gpres)/2._dl
Hdot =adotdota-adotoa**2._dl
EFT_gpinu_tot = 0._dl
EFT_grhonu_tot = 0._dl
EFT_gpinudot_tot = 0._dl
if (CP%Num_Nu_Massive /= 0) then
do nu_i = 1, CP%Nu_mass_eigenstates
EFT_grhonu = 0._dl
EFT_gpinu = 0._dl
EFT_grhonudot = 0._dl
EFT_gpinudot = 0._dl
grhormass_t=grhormass(nu_i)/a**2
call Nu_background(a*nu_masses(nu_i),EFT_grhonu,EFT_gpinu)
EFT_grhonu_tot = EFT_grhonu_tot + grhormass_t*EFT_grhonu
EFT_gpinu_tot = EFT_gpinu_tot + grhormass_t*EFT_gpinu
EFT_gpinudot_tot = EFT_gpinudot_tot + grhormass_t*(Nu_pidot(a*nu_masses(nu_i),adotoa, EFT_gpinu)&
& -4._dl*adotoa*EFT_gpinu)
end do
end if
Hdotdot = 2._dl*adotoa*Hdot &
& + 0.5_dl*adotoa*(grhob_t + grhoc_t + 8._dl*(grhog_t+grhor_t)/3._dl)&
& + 0.5_dl*adotoa*grhov_t*((1._dl+EFTw(a,0))*(1._dl+3._dl*EFTw(a,0))-a*EFTw(a,1))&
& + adotoa/6._dl*EFT_grhonu_tot -0.5_dl*adotoa*EFT_gpinu_tot -0.5_dl*EFT_gpinudot_tot
! 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)
!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))
if (CP%EFTflag==2.and.CP%DesignerEFTmodel/=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
real(dl) EFT_grhonu, EFT_gpinu, EFT_grhonudot, EFT_gpinudot, grhormass_t
real(dl) EFT_grhonudot_tot, EFT_grhonu_tot, EFT_gpinu_tot, EFT_gpinudot_tot
! 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
integer nu_i
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,3)
grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhov_t
! Massive neutrinos
if (CP%Num_Nu_Massive /= 0) then
do nu_i = 1, CP%Nu_mass_eigenstates
EFT_grhonu = 0._dl
EFT_gpinu = 0._dl
EFT_grhonudot = 0._dl
EFT_gpinudot = 0._dl
grhormass_t=grhormass(nu_i)/a**2
call Nu_background(a*nu_masses(nu_i),EFT_grhonu,EFT_gpinu)
EFT_grhonu_tot = EFT_grhonu_tot + grhormass_t*EFT_grhonu
EFT_gpinu_tot = EFT_gpinu_tot + grhormass_t*EFT_gpinu
end do
end if
grho = grho + EFT_grhonu_tot
gpres = +EFT_gpinu_tot +(grhog_t+grhor_t)/3._dl +EFTw(a,0)*grhov_t
! 3) Hubble and derivatives:
adotoa=sqrt(grho/3)
adotdota=(adotoa*adotoa-gpres)/2._dl
Hdot =adotdota-adotoa**2._dl
EFT_gpinu_tot = 0._dl
EFT_grhonu_tot = 0._dl
EFT_gpinudot_tot = 0._dl
if (CP%Num_Nu_Massive /= 0) then
do nu_i = 1, CP%Nu_mass_eigenstates
EFT_grhonu = 0._dl
EFT_gpinu = 0._dl
EFT_grhonudot = 0._dl
EFT_gpinudot = 0._dl
grhormass_t=grhormass(nu_i)/a**2
call Nu_background(a*nu_masses(nu_i),EFT_grhonu,EFT_gpinu)
EFT_grhonu_tot = EFT_grhonu_tot + grhormass_t*EFT_grhonu
EFT_gpinu_tot = EFT_gpinu_tot + grhormass_t*EFT_gpinu
EFT_gpinudot_tot = EFT_gpinudot_tot + grhormass_t*(Nu_pidot(a*nu_masses(nu_i),adotoa, EFT_gpinu)&
& -4._dl*adotoa*EFT_gpinu)
end do
end if
Hdotdot = 2._dl*adotoa*Hdot &
& + 0.5_dl*adotoa*(grhob_t + grhoc_t + 8._dl*(grhog_t+grhor_t)/3._dl)&
& + 0.5_dl*adotoa*grhov_t*((1._dl+EFTw(a,0))*(1._dl+3._dl*EFTw(a,0))-a*EFTw(a,1))&
& + adotoa/6._dl*EFT_grhonu_tot -0.5_dl*adotoa*EFT_gpinu_tot -0.5_dl*EFT_gpinudot_tot
! 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)
!8*pi*G*c*a^2
EFTc = (adotoa*adotoa - Hdot)*(EFTOmegaV + a*EFTOmegaP*0.5_dl) &
& - 0.5_dl*a2*adotoa*adotoa*EFTOmegaPP&
& + 0.5_dl*grhov_t*(1._dl+EFTw(a,0))
!8*pi*G*Lambda*a^2
EFTLambda = +EFTw(a,0)*grhov_t &
&-EFTOmegaV*(2._dl*Hdot+adotoa**2._dl) &
&-a*EFTOmegaP*(2._dl*adotoa**2._dl + Hdot) &
&-a2*adotoa**2._dl*EFTOmegaPP
!EFT C DOT: 8*pi*G*cdot*a^2
EFTcdot = -EFTOmegaV*(Hdotdot-4._dl*adotoa*Hdot+2._dl*adotoa*adotoa*adotoa) &
& + 0.5_dl*a*EFTOmegaP*(-Hdotdot+adotoa*Hdot+adotoa*adotoa*adotoa)&
& +0.5_dl*a2*adotoa*EFTOmegaPP*(adotoa*adotoa-3._dl*Hdot)&
& -0.5_dl*a*a2*adotoa*adotoa*adotoa*EFTOmegaPPP&
& +0.5_dl*adotoa*grhov_t*(-3._dl*(1._dl+EFTw(a,0))**2 + a*EFTw(a,1))
!EFT LAMBDA DOT: 8*pi*G*Ldot*a^2
EFTLambdadot = -2._dl*EFTOmegaV*(Hdotdot-adotoa*Hdot-adotoa*adotoa*adotoa)&
& - a*EFTOmegaP*(4._dl*adotoa*Hdot+Hdotdot)&
& -a2*EFTOmegaPP*adotoa*(3._dl*Hdot+2._dl*adotoa*adotoa)&
& -a*a2*EFTOmegaPPP*adotoa*adotoa*adotoa&
& +grhov_t*adotoa*(a*EFTw(a,1)-3._dl*EFTw(a,0)*(1._dl+EFTw(a,0)))
if (CP%EFTflag==2.and.CP%DesignerEFTmodel/=2) then
EFTc = EFTcTemp(a,0)
EFTcdot = EFTcTemp(a,1)
EFTLambda = EFTLambdaTemp(a,0)
EFTLambdadot = EFTLambdaTemp(a,1)
end if
EFTgrhoq = 2._dl*EFTc -EFTLambda -3._dl*a*adotoa*adotoa*EFTOmegaP
EFTgpresq = EFTLambda + a2*adotoa*adotoa*EFTOmegaPP&
& +a*EFTOmegaP*(Hdot + 2._dl*adotoa*adotoa)
EFTgrhodotq = -3._dl*adotoa*(EFTgrhoq+EFTgpresq) + 3._dl*a*adotoa**3._dl*EFTOmegaP
EFTgpresdotq = EFTLambdadot +a2*a*adotoa**3*EFTOmegaPPP + 3._dl*a2*adotoa*Hdot*EFTOmegaPP&
& +a*EFTOmegaP*Hdotdot +3._dl*a*adotoa*Hdot*EFTOmegaP +2._dl*a2*adotoa**3*EFTOmegaPP&
& -2._dl*a*adotoa**3*EFTOmegaP
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
real(dl) EFT_grhonu, EFT_gpinu, EFT_grhonudot, EFT_gpinudot, grhormass_t
real(dl) EFT_grhonudot_tot, EFT_grhonu_tot, EFT_gpinu_tot, EFT_gpinudot_tot
! 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
integer nu_i
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,3)
grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhov_t
! Massive neutrinos
if (CP%Num_Nu_Massive /= 0) then
do nu_i = 1, CP%Nu_mass_eigenstates
EFT_grhonu = 0._dl
EFT_gpinu = 0._dl
EFT_grhonudot = 0._dl
EFT_gpinudot = 0._dl
grhormass_t=grhormass(nu_i)/a**2
call Nu_background(a*nu_masses(nu_i),EFT_grhonu,EFT_gpinu)
EFT_grhonu_tot = EFT_grhonu_tot + grhormass_t*EFT_grhonu
EFT_gpinu_tot = EFT_gpinu_tot + grhormass_t*EFT_gpinu
end do
end if
grho = grho + EFT_grhonu_tot
gpres = +EFT_gpinu_tot +(grhog_t+grhor_t)/3._dl +EFTw(a,0)*grhov_t
! 3) Hubble and derivatives:
adotoa=sqrt(grho/3)
adotdota=(adotoa*adotoa-gpres)/2._dl
Hdot =adotdota-adotoa**2._dl
EFT_gpinu_tot = 0._dl
EFT_grhonu_tot = 0._dl
EFT_gpinudot_tot = 0._dl
if (CP%Num_Nu_Massive /= 0) then
do nu_i = 1, CP%Nu_mass_eigenstates
EFT_grhonu = 0._dl
EFT_gpinu = 0._dl
EFT_grhonudot = 0._dl
EFT_gpinudot = 0._dl
grhormass_t=grhormass(nu_i)/a**2
call Nu_background(a*nu_masses(nu_i),EFT_grhonu,EFT_gpinu)
EFT_grhonu_tot = EFT_grhonu_tot + grhormass_t*EFT_grhonu
EFT_gpinu_tot = EFT_gpinu_tot + grhormass_t*EFT_gpinu
EFT_gpinudot_tot = EFT_gpinudot_tot + grhormass_t*(Nu_pidot(a*nu_masses(nu_i),adotoa, EFT_gpinu)&
& -4._dl*adotoa*EFT_gpinu)
end do
end if
Hdotdot = 2._dl*adotoa*Hdot &
& + 0.5_dl*adotoa*(grhob_t + grhoc_t + 8._dl*(grhog_t+grhor_t)/3._dl)&
& + 0.5_dl*adotoa*grhov_t*((1._dl+EFTw(a,0))*(1._dl+3._dl*EFTw(a,0))-a*EFTw(a,1))&
& + adotoa/6._dl*EFT_grhonu_tot -0.5_dl*adotoa*EFT_gpinu_tot -0.5_dl*EFT_gpinudot_tot
! 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)
!8*pi*G*c*a^2
EFTc = (adotoa*adotoa - Hdot)*(EFTOmegaV + a*EFTOmegaP*0.5_dl) &
& - 0.5_dl*a2*adotoa*adotoa*EFTOmegaPP&
& + 0.5_dl*grhov_t*(1._dl+EFTw(a,0))
!8*pi*G*Lambda*a^2
EFTLambda = +EFTw(a,0)*grhov_t &
&-EFTOmegaV*(2._dl*Hdot+adotoa**2._dl) &
&-a*EFTOmegaP*(2._dl*adotoa**2._dl + Hdot) &
&-a2*adotoa**2._dl*EFTOmegaPP
!EFT C DOT: 8*pi*G*cdot*a^2
EFTcdot = -EFTOmegaV*(Hdotdot-4._dl*adotoa*Hdot+2._dl*adotoa*adotoa*adotoa) &
& + 0.5_dl*a*EFTOmegaP*(-Hdotdot+adotoa*Hdot+adotoa*adotoa*adotoa)&
& +0.5_dl*a2*adotoa*EFTOmegaPP*(adotoa*adotoa-3._dl*Hdot)&
& -0.5_dl*a*a2*adotoa*adotoa*adotoa*EFTOmegaPPP&
& +0.5_dl*adotoa*grhov_t*(-3._dl*(1._dl+EFTw(a,0))**2 + a*EFTw(a,1))
!EFT LAMBDA DOT: 8*pi*G*Ldot*a^2
EFTLambdadot = -2._dl*EFTOmegaV*(Hdotdot-adotoa*Hdot-adotoa*adotoa*adotoa)&
& - a*EFTOmegaP*(4._dl*adotoa*Hdot+Hdotdot)&
& -a2*EFTOmegaPP*adotoa*(3._dl*Hdot+2._dl*adotoa*adotoa)&
& -a*a2*EFTOmegaPPP*adotoa*adotoa*adotoa&
& +grhov_t*adotoa*(a*EFTw(a,1)-3._dl*EFTw(a,0)*(1._dl+EFTw(a,0)))
if (CP%EFTflag==2.and.CP%DesignerEFTmodel/=2) then
EFTc = EFTcTemp(a,0)
EFTcdot = EFTcTemp(a,1)
EFTLambda = EFTLambdaTemp(a,0)
EFTLambdadot = EFTLambdaTemp(a,1)
end if
EFTgrhoq = 2._dl*EFTc -EFTLambda -3._dl*a*adotoa*adotoa*EFTOmegaP
EFTgpresq = EFTLambda + a2*adotoa*adotoa*EFTOmegaPP&
& +a*EFTOmegaP*(Hdot + 2._dl*adotoa*adotoa)
EFTgrhodotq = -3._dl*adotoa*(EFTgrhoq+EFTgpresq) + 3._dl*a*adotoa**3._dl*EFTOmegaP
EFTgpresdotq = EFTLambdadot +a2*a*adotoa**3*EFTOmegaPPP + 3._dl*a2*adotoa*Hdot*EFTOmegaPP&
& +a*EFTOmegaP*Hdotdot +3._dl*a*adotoa*Hdot*EFTOmegaP +2._dl*a2*adotoa**3*EFTOmegaPP&
& -2._dl*a*adotoa**3*EFTOmegaP
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
real(dl) EFT_grhonu, EFT_gpinu, EFT_grhonudot, EFT_gpinudot, grhormass_t
real(dl) EFT_grhonudot_tot, EFT_grhonu_tot, EFT_gpinu_tot, EFT_gpinudot_tot
! 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
integer nu_i
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,3)
grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhov_t
! Massive neutrinos
if (CP%Num_Nu_Massive /= 0) then
do nu_i = 1, CP%Nu_mass_eigenstates
EFT_grhonu = 0._dl
EFT_gpinu = 0._dl
EFT_grhonudot = 0._dl
EFT_gpinudot = 0._dl
grhormass_t=grhormass(nu_i)/a**2
call Nu_background(a*nu_masses(nu_i),EFT_grhonu,EFT_gpinu)
EFT_grhonu_tot = EFT_grhonu_tot + grhormass_t*EFT_grhonu
EFT_gpinu_tot = EFT_gpinu_tot + grhormass_t*EFT_gpinu
end do
end if
grho = grho + EFT_grhonu_tot
gpres = +EFT_gpinu_tot +(grhog_t+grhor_t)/3._dl +EFTw(a,0)*grhov_t
! 3) Hubble and derivatives:
adotoa=sqrt(grho/3)
adotdota=(adotoa*adotoa-gpres)/2._dl
Hdot =adotdota-adotoa**2._dl
EFT_gpinu_tot = 0._dl
EFT_grhonu_tot = 0._dl
EFT_gpinudot_tot = 0._dl
if (CP%Num_Nu_Massive /= 0) then
do nu_i = 1, CP%Nu_mass_eigenstates
EFT_grhonu = 0._dl
EFT_gpinu = 0._dl
EFT_grhonudot = 0._dl
EFT_gpinudot = 0._dl
grhormass_t=grhormass(nu_i)/a**2
call Nu_background(a*nu_masses(nu_i),EFT_grhonu,EFT_gpinu)
EFT_grhonu_tot = EFT_grhonu_tot + grhormass_t*EFT_grhonu
EFT_gpinu_tot = EFT_gpinu_tot + grhormass_t*EFT_gpinu
EFT_gpinudot_tot = EFT_gpinudot_tot + grhormass_t*(Nu_pidot(a*nu_masses(nu_i),adotoa, EFT_gpinu)&
& -4._dl*adotoa*EFT_gpinu)
end do
end if
Hdotdot = 2._dl*adotoa*Hdot &
& + 0.5_dl*adotoa*(grhob_t + grhoc_t + 8._dl*(grhog_t+grhor_t)/3._dl)&
& + 0.5_dl*adotoa*grhov_t*((1._dl+EFTw(a,0))*(1._dl+3._dl*EFTw(a,0))-a*EFTw(a,1))&
& + adotoa/6._dl*EFT_grhonu_tot -0.5_dl*adotoa*EFT_gpinu_tot -0.5_dl*EFT_gpinudot_tot
! 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)
!8*pi*G*c*a^2
EFTc = (adotoa*adotoa - Hdot)*(EFTOmegaV + a*EFTOmegaP*0.5_dl) &
& - 0.5_dl*a2*adotoa*adotoa*EFTOmegaPP&
& + 0.5_dl*grhov_t*(1._dl+EFTw(a,0))
!8*pi*G*Lambda*a^2
EFTLambda = +EFTw(a,0)*grhov_t &
&-EFTOmegaV*(2._dl*Hdot+adotoa**2._dl) &
&-a*EFTOmegaP*(2._dl*adotoa**2._dl + Hdot) &
&-a2*adotoa**2._dl*EFTOmegaPP
!EFT C DOT: 8*pi*G*cdot*a^2
EFTcdot = -EFTOmegaV*(Hdotdot-4._dl*adotoa*Hdot+2._dl*adotoa*adotoa*adotoa) &
& + 0.5_dl*a*EFTOmegaP*(-Hdotdot+adotoa*Hdot+adotoa*adotoa*adotoa)&
& +0.5_dl*a2*adotoa*EFTOmegaPP*(adotoa*adotoa-3._dl*Hdot)&
& -0.5_dl*a*a2*adotoa*adotoa*adotoa*EFTOmegaPPP&
& +0.5_dl*adotoa*grhov_t*(-3._dl*(1._dl+EFTw(a,0))**2 + a*EFTw(a,1))
!EFT LAMBDA DOT: 8*pi*G*Ldot*a^2
EFTLambdadot = -2._dl*EFTOmegaV*(Hdotdot-adotoa*Hdot-adotoa*adotoa*adotoa)&
& - a*EFTOmegaP*(4._dl*adotoa*Hdot+Hdotdot)&
& -a2*EFTOmegaPP*adotoa*(3._dl*Hdot+2._dl*adotoa*adotoa)&
& -a*a2*EFTOmegaPPP*adotoa*adotoa*adotoa&
& +grhov_t*adotoa*(a*EFTw(a,1)-3._dl*EFTw(a,0)*(1._dl+EFTw(a,0)))
if (CP%EFTflag==2.and.CP%DesignerEFTmodel/=2) then
EFTc = EFTcTemp(a,0)
EFTcdot = EFTcTemp(a,1)
EFTLambda = EFTLambdaTemp(a,0)
EFTLambdadot = EFTLambdaTemp(a,1)
end if
EFTgrhoq = 2._dl*EFTc -EFTLambda -3._dl*a*adotoa*adotoa*EFTOmegaP
EFTgpresq = EFTLambda + a2*adotoa*adotoa*EFTOmegaPP&
& +a*EFTOmegaP*(Hdot + 2._dl*adotoa*adotoa)
EFTgrhodotq = -3._dl*adotoa*(EFTgrhoq+EFTgpresq) + 3._dl*a*adotoa**3._dl*EFTOmegaP
EFTgpresdotq = EFTLambdadot +a2*a*adotoa**3*EFTOmegaPPP + 3._dl*a2*adotoa*Hdot*EFTOmegaPP&
& +a*EFTOmegaP*Hdotdot +3._dl*a*adotoa*Hdot*EFTOmegaP +2._dl*a2*adotoa**3*EFTOmegaPP&
& -2._dl*a*adotoa**3*EFTOmegaP
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
----------------------------------------------------------------------------------------------------------------------------
Modify: subroutines.f90Add: After end subroutine dverk
! EFTCAMB MOD START
! -------------------------------------------------------------------------------------------------
! EFTCAMB numerical subroutines.
! Includes several algorithms that are used by the EFTCAMB code.
! They are written in a general form and can be used at different purposes.
! -------------------------------------------------------------------------------------------------
! 1) Brent root finding algorithm.
! This is used to solve numerically the equation func=funcZero by means of the Brent method.
function zbrent(func,x1,x2,tol,funcZero,succes)
use precision
implicit none
real(dl) :: tol
real(dl) :: x1,x2
real(dl) :: funcZero
logical :: succes
real(dl) :: func
real(dl) zbrent
external func
real(dl),parameter :: EPS = 3.d-15
integer ,parameter :: ITMAX = 2000 !Max Number of iterations.
integer iter
real(dl) a,b,c,d,e,fa,fb,fc,p,q,r,s,tol1,xm
succes = .true.
a=x1
b=x2
fa=func(a)-funcZero
fb=func(b)-funcZero
if((fa.gt.0..and.fb.gt.0.).or.(fa.lt.0..and.fb.lt.0.)) then
succes=.false.
return
end if
c=b
fc=fb
do 11 iter=1,ITMAX
if((fb.gt.0..and.fc.gt.0.).or.(fb.lt.0..and.fc.lt.0.))then
c=a
fc=fa
d=b-a
e=d
endif
if(abs(fc).lt.abs(fb)) then
a=b
b=c
c=a
fa=fb
fb=fc
fc=fa
endif
tol1=2.*EPS*abs(b)+0.5*tol
xm=.5*(c-b)
if(abs(xm).le.tol1 .or. fb.eq.0.)then
zbrent=b
return
endif
if(abs(e).ge.tol1 .and. abs(fa).gt.abs(fb)) then
s=fb/fa
if(a.eq.c) then
p=2.*xm*s
q=1.-s
else
q=fa/fc
r=fb/fc
p=s*(2.*xm*q*(q-r)-(b-a)*(r-1.))
q=(q-1.)*(r-1.)*(s-1.)
endif
if(p.gt.0.) q=-q
p=abs(p)
if(2.*p .lt. min(3.*xm*q-abs(tol1*q),abs(e*q))) then
e=d
d=p/q
else
d=xm
e=d
endif
else
d=xm
e=d
endif
a=b
fa=fb
if(abs(d) .gt. tol1) then
b=b+d
else
b=b+sign(tol1,xm)
endif
fb=func(b)-funcZero
11 continue
succes = .false.
zbrent=b
return
end function zbrent
! -------------------------------------------------------------------------------------------------
! 2) Hunting algorithm.
! This is used to efficiently search an ordered table by means of a hunting and
! a bisection algorithm.
subroutine hunt(xx,n,x,jlo)
! xx = the table to be searched (in);
! n = the length of the table (in);
! x = the requested value (in);
! jlo = the closest (from below) entry of the table (out).
use precision
implicit none
integer jlo,n
real(dl) x,xx(n)
integer inc,jhi,jm
logical ascnd
ascnd=xx(n).ge.xx(1)
if(jlo.le.0.or.jlo.gt.n)then
jlo=0
jhi=n+1
goto 3
endif
inc=1
if(x.ge.xx(jlo).eqv.ascnd)then
1 jhi=jlo+inc
if(jhi.gt.n)then
jhi=n+1
else if(x.ge.xx(jhi).eqv.ascnd)then
jlo=jhi
inc=inc+inc
goto 1
endif
else
jhi=jlo
2 jlo=jhi-inc
if(jlo.lt.1)then
jlo=0
else if(x.lt.xx(jlo).eqv.ascnd)then
jhi=jlo
inc=inc+inc
goto 2
endif
endif
3 if(jhi-jlo.eq.1)then
if(x.eq.xx(n))jlo=n-1
if(x.eq.xx(1))jlo=1
return
endif
jm=(jhi+jlo)/2
if(x.ge.xx(jm).eqv.ascnd)then
jlo=jm
else
jhi=jm
endif
goto 3
end subroutine hunt
! -------------------------------------------------------------------------------------------------
! 3) Bracketing subroutine.
! This subroutine does a outward search for the smallest intervall containing a root
! of the equation func=funcZero
subroutine zbrac(func,x1,x2,succes,funcZero)
! func = function to find the root of (in);
! x1 = lower bound of the interval in which to find the root (in);
! x2 = upper bound of the interval in which to find the root.
! The intervall must bracket the root (in);
! success = true if the algorithm succeeds false if not (inout);
! funcZero = the value desired func = funcZero (in).
use precision
implicit none
real(dl) func, funcZero
real(dl) x1,x2
logical succes
integer, parameter :: NTRY = 1000
real(dl), parameter :: FACTOR = 5._dl
real(dl) delta, temp
real(dl) f1,f2
integer j
external func
if (x1.eq.x2) stop 'you have to guess an initial range in zbrac'
if (x2
temp=x2
x2=x1
x1=temp
end if
f1=func(x1)-funcZero
f2=func(x2)-funcZero
succes=.true.
do 11 j=1,NTRY
if(f1*f2.lt.0.)return
delta=ABS(x2-x1)
x1=x1-FACTOR*delta
x2=x2+FACTOR*delta
f1=func(x1)-funcZero
f2=func(x2)-funcZero
11 continue
succes=.false.
return
end subroutine zbrac
! -------------------------------------------------------------------------------------------------
! 4) Neville interpolator.
! This is used to interpolate the EFT functions once the designer/mapping
! code has sampled them
subroutine Polint(n,xa,ya,xpl,ypl,dypl)
! n = number of points in the table (in)
! xa = first coordinate of the points to be interpolated (in)
! ya = second coordinate of the points to be interpolated (in)
! xpl = requested value of x (in)
! ypl = value of the interpolated function at xpl (out)
! dypl =
use precision
implicit none
integer , intent(in) :: n ! Length of the table of points
real(dl), intent(in) :: xa(n),ya(n) ! The tables of points to be interpolated as y=f(x)
real(dl) :: dypl,xpl,ypl
integer :: i,m,ns
real(dl) :: den,dif,dift,ho,hp,wpl,cc(n),d(n)
ns=1
dif=abs(xpl-xa(1))
do i=1,n
dift=abs(xpl-xa(i))
if (dift
ns=i
dif=dift
endif
cc(i)=ya(i)
d(i)=ya(i)
end do
ypl=ya(ns)
ns=ns-1
do m=1,n-1
do i=1,n-m
ho=xa(i)-xpl
hp=xa(i+m)-xpl
wpl=cc(i+1)-d(i)
den=ho-hp
if (den==0.) then
write(*,*) 'failure in polint'
stop
end if
den=wpl/den
d(i)=hp*den
cc(i)=ho*den
end do
if (2*ns
dypl=cc(ns+1)
else
dypl=d(ns)
ns=ns-1
endif
ypl=ypl+dypl
end do
return
end subroutine Polint
! -------------------------------------------------------------------------------------------------
! 5) Fourth order Runge-Kutta.
! This is a very simple algorithm that is used to solve the designer equation.
!
subroutine EFT_rk4(n, y, dydx, x, h, yout, deriv)
! n = dimensionality of the problem;
! y = 'position' at t=x;
! dydx = 'velocity' at t=x;
! x = initial time;
! h = time step;
! yout = 'position' at t=x+h computed using fourth order Runge-Kutta;
! deriv = name of the subroutine that computes dydx.
use precision
implicit none
real(dl) :: x, h
integer :: n
real(dl), dimension(n) :: y, dydx, yout
interface
subroutine deriv(n, x, y, dydx)
use precision
implicit none
real(dl) :: x
integer :: n
real(dl), dimension(n) :: y
real(dl), dimension(n) :: dydx
end subroutine deriv
end interface
real(dl), dimension(n) :: yt, dyt,dym
real(dl) :: hh,h6,xh
integer :: i
hh=h*0.5_dl
h6=h/6._dl
xh=x+hh
yt=y+hh*dydx
call deriv(n, xh, yt, dyt)
yt=y+hh*dyt
call deriv(n, xh, yt, dym)
yt =y+h*dym
dym =dyt+dym
call deriv(n, x+h, yt, dyt)
yout=y+h6*(dydx+dyt+2.0*dym)
return
end subroutine EFT_rk4
! EFTCAMB MOD END
! -------------------------------------------------------------------------------------------------
! EFTCAMB numerical subroutines.
! Includes several algorithms that are used by the EFTCAMB code.
! They are written in a general form and can be used at different purposes.
! -------------------------------------------------------------------------------------------------
! 1) Brent root finding algorithm.
! This is used to solve numerically the equation func=funcZero by means of the Brent method.
function zbrent(func,x1,x2,tol,funcZero,succes)
use precision
implicit none
real(dl) :: tol
real(dl) :: x1,x2
real(dl) :: funcZero
logical :: succes
real(dl) :: func
real(dl) zbrent
external func
real(dl),parameter :: EPS = 3.d-15
integer ,parameter :: ITMAX = 2000 !Max Number of iterations.
integer iter
real(dl) a,b,c,d,e,fa,fb,fc,p,q,r,s,tol1,xm
succes = .true.
a=x1
b=x2
fa=func(a)-funcZero
fb=func(b)-funcZero
if((fa.gt.0..and.fb.gt.0.).or.(fa.lt.0..and.fb.lt.0.)) then
succes=.false.
return
end if
c=b
fc=fb
do 11 iter=1,ITMAX
if((fb.gt.0..and.fc.gt.0.).or.(fb.lt.0..and.fc.lt.0.))then
c=a
fc=fa
d=b-a
e=d
endif
if(abs(fc).lt.abs(fb)) then
a=b
b=c
c=a
fa=fb
fb=fc
fc=fa
endif
tol1=2.*EPS*abs(b)+0.5*tol
xm=.5*(c-b)
if(abs(xm).le.tol1 .or. fb.eq.0.)then
zbrent=b
return
endif
if(abs(e).ge.tol1 .and. abs(fa).gt.abs(fb)) then
s=fb/fa
if(a.eq.c) then
p=2.*xm*s
q=1.-s
else
q=fa/fc
r=fb/fc
p=s*(2.*xm*q*(q-r)-(b-a)*(r-1.))
q=(q-1.)*(r-1.)*(s-1.)
endif
if(p.gt.0.) q=-q
p=abs(p)
if(2.*p .lt. min(3.*xm*q-abs(tol1*q),abs(e*q))) then
e=d
d=p/q
else
d=xm
e=d
endif
else
d=xm
e=d
endif
a=b
fa=fb
if(abs(d) .gt. tol1) then
b=b+d
else
b=b+sign(tol1,xm)
endif
fb=func(b)-funcZero
11 continue
succes = .false.
zbrent=b
return
end function zbrent
! -------------------------------------------------------------------------------------------------
! 2) Hunting algorithm.
! This is used to efficiently search an ordered table by means of a hunting and
! a bisection algorithm.
subroutine hunt(xx,n,x,jlo)
! xx = the table to be searched (in);
! n = the length of the table (in);
! x = the requested value (in);
! jlo = the closest (from below) entry of the table (out).
use precision
implicit none
integer jlo,n
real(dl) x,xx(n)
integer inc,jhi,jm
logical ascnd
ascnd=xx(n).ge.xx(1)
if(jlo.le.0.or.jlo.gt.n)then
jlo=0
jhi=n+1
goto 3
endif
inc=1
if(x.ge.xx(jlo).eqv.ascnd)then
1 jhi=jlo+inc
if(jhi.gt.n)then
jhi=n+1
else if(x.ge.xx(jhi).eqv.ascnd)then
jlo=jhi
inc=inc+inc
goto 1
endif
else
jhi=jlo
2 jlo=jhi-inc
if(jlo.lt.1)then
jlo=0
else if(x.lt.xx(jlo).eqv.ascnd)then
jhi=jlo
inc=inc+inc
goto 2
endif
endif
3 if(jhi-jlo.eq.1)then
if(x.eq.xx(n))jlo=n-1
if(x.eq.xx(1))jlo=1
return
endif
jm=(jhi+jlo)/2
if(x.ge.xx(jm).eqv.ascnd)then
jlo=jm
else
jhi=jm
endif
goto 3
end subroutine hunt
! -------------------------------------------------------------------------------------------------
! 3) Bracketing subroutine.
! This subroutine does a outward search for the smallest intervall containing a root
! of the equation func=funcZero
subroutine zbrac(func,x1,x2,succes,funcZero)
! func = function to find the root of (in);
! x1 = lower bound of the interval in which to find the root (in);
! x2 = upper bound of the interval in which to find the root.
! The intervall must bracket the root (in);
! success = true if the algorithm succeeds false if not (inout);
! funcZero = the value desired func = funcZero (in).
use precision
implicit none
real(dl) func, funcZero
real(dl) x1,x2
logical succes
integer, parameter :: NTRY = 1000
real(dl), parameter :: FACTOR = 5._dl
real(dl) delta, temp
real(dl) f1,f2
integer j
external func
if (x1.eq.x2) stop 'you have to guess an initial range in zbrac'
if (x2
x2=x1
x1=temp
end if
f1=func(x1)-funcZero
f2=func(x2)-funcZero
succes=.true.
do 11 j=1,NTRY
if(f1*f2.lt.0.)return
delta=ABS(x2-x1)
x1=x1-FACTOR*delta
x2=x2+FACTOR*delta
f1=func(x1)-funcZero
f2=func(x2)-funcZero
11 continue
succes=.false.
return
end subroutine zbrac
! -------------------------------------------------------------------------------------------------
! 4) Neville interpolator.
! This is used to interpolate the EFT functions once the designer/mapping
! code has sampled them
subroutine Polint(n,xa,ya,xpl,ypl,dypl)
! n = number of points in the table (in)
! xa = first coordinate of the points to be interpolated (in)
! ya = second coordinate of the points to be interpolated (in)
! xpl = requested value of x (in)
! ypl = value of the interpolated function at xpl (out)
! dypl =
use precision
implicit none
integer , intent(in) :: n ! Length of the table of points
real(dl), intent(in) :: xa(n),ya(n) ! The tables of points to be interpolated as y=f(x)
real(dl) :: dypl,xpl,ypl
integer :: i,m,ns
real(dl) :: den,dif,dift,ho,hp,wpl,cc(n),d(n)
ns=1
dif=abs(xpl-xa(1))
do i=1,n
dift=abs(xpl-xa(i))
if (dift
dif=dift
endif
cc(i)=ya(i)
d(i)=ya(i)
end do
ypl=ya(ns)
ns=ns-1
do m=1,n-1
do i=1,n-m
ho=xa(i)-xpl
hp=xa(i+m)-xpl
wpl=cc(i+1)-d(i)
den=ho-hp
if (den==0.) then
write(*,*) 'failure in polint'
stop
end if
den=wpl/den
d(i)=hp*den
cc(i)=ho*den
end do
if (2*ns
else
dypl=d(ns)
ns=ns-1
endif
ypl=ypl+dypl
end do
return
end subroutine Polint
! -------------------------------------------------------------------------------------------------
! 5) Fourth order Runge-Kutta.
! This is a very simple algorithm that is used to solve the designer equation.
!
subroutine EFT_rk4(n, y, dydx, x, h, yout, deriv)
! n = dimensionality of the problem;
! y = 'position' at t=x;
! dydx = 'velocity' at t=x;
! x = initial time;
! h = time step;
! yout = 'position' at t=x+h computed using fourth order Runge-Kutta;
! deriv = name of the subroutine that computes dydx.
use precision
implicit none
real(dl) :: x, h
integer :: n
real(dl), dimension(n) :: y, dydx, yout
interface
subroutine deriv(n, x, y, dydx)
use precision
implicit none
real(dl) :: x
integer :: n
real(dl), dimension(n) :: y
real(dl), dimension(n) :: dydx
end subroutine deriv
end interface
real(dl), dimension(n) :: yt, dyt,dym
real(dl) :: hh,h6,xh
integer :: i
hh=h*0.5_dl
h6=h/6._dl
xh=x+hh
yt=y+hh*dydx
call deriv(n, xh, yt, dyt)
yt=y+hh*dyt
call deriv(n, xh, yt, dym)
yt =y+h*dym
dym =dyt+dym
call deriv(n, x+h, yt, dyt)
yout=y+h6*(dydx+dyt+2.0*dym)
return
end subroutine EFT_rk4
! EFTCAMB MOD END
Re-compile and Make it.
-5-