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.
----------------------------------------------------------------------------------------------------------------------------
Add: EFT_main.f90----------------------------------------------------------------------------------------------------------------------------
Add: EFTstabilitySpace.f90----------------------------------------------------------------------------------------------------------------------------
Modify: Makefile_main
Replace:
EQUATIONS ?= equations_EFT
#... Original:
#EQUATIONS ?= equations
#... Original:
#EQUATIONS ?= equations
Add: After
"#DRIVER ?= tester.f90"
# EFTCAMB: optional program to test stability.
#DRIVER ?= EFTstabilitySpace.f90
#DRIVER ?= EFTstabilitySpace.f90
Replace:
# EFTCAMB MOD START
CAMBOBJ = constants.o utils.o $(EXTCAMBFILES) subroutines.o inifile.o $(POWERSPECTRUM).o $(RECOMBINATION).o \
$(REIONIZATION).o modules.o bessels.o EFT_main.o $(EQUATIONS).o $(NONLINEAR).o lensing.o $(BISPECTRUM).o cmbmain.o camb.o
#... Original:
#CAMBOBJ = constants.o utils.o $(EXTCAMBFILES) subroutines.o inifile.o $(POWERSPECTRUM).o $(RECOMBINATION).o \
# $(REIONIZATION).o modules.o bessels.o $(EQUATIONS).o $(NONLINEAR).o lensing.o $(BISPECTRUM).o cmbmain.o camb.o
# EFTCAMB MOD END
CAMBOBJ = constants.o utils.o $(EXTCAMBFILES) subroutines.o inifile.o $(POWERSPECTRUM).o $(RECOMBINATION).o \
$(REIONIZATION).o modules.o bessels.o EFT_main.o $(EQUATIONS).o $(NONLINEAR).o lensing.o $(BISPECTRUM).o cmbmain.o camb.o
#... Original:
#CAMBOBJ = constants.o utils.o $(EXTCAMBFILES) subroutines.o inifile.o $(POWERSPECTRUM).o $(RECOMBINATION).o \
# $(REIONIZATION).o modules.o bessels.o $(EQUATIONS).o $(NONLINEAR).o lensing.o $(BISPECTRUM).o cmbmain.o camb.o
# EFTCAMB MOD END
----------------------------------------------------------------------------------------------------------------------------
Modify: params.ini
Add: Before
"#Parameters for CAMB"
####### Model selection flags for EFTCAMB #######
# EFT flags: set up in which mode EFTCAMB is used.
# We refer to the documentation (EFTCAMB:numerical notes) for a thorough
# explanation of the effect of these flags.
#
# 1) Main EFT flag:
# EFTflag = 0 : GR code. Every EFT modification is ignored.
# EFTflag = 1 : Pure EFT code.
# EFTflag = 2 : Designer matching EFT.
# EFTflag = 3 : Matching EFT.
EFTflag = 0
# 2) Background Dark Energy equation of state flag:
# EFTwDE = 0 : Cosmological constant
# EFTwDE = 1 : DE with constant Eos determined by EFTw0
# EFTwDE = 2 : CPL parametrization
EFTwDE = 1
# 3) Pure EFT model selection flag:
# The following structure applies to every operator and can be specified for every operatori
# separately.
#
# PureEFTmodel___ = 0 : Zero (operator ignored)
# PureEFTmodel___ = 1 : Constant model
# PureEFTmodel___ = 2 : Linear model
# PureEFTmodel___ = 3 : Power law model
# PureEFTmodel___ = 4 : Exponential model
# PureEFTmodel___ = 5 : User defined
PureEFTmodelOmega = 2
PureEFTmodelAlpha1 = 0
PureEFTmodelAlpha2 = 0
PureEFTmodelAlpha3 = 0
PureEFTmodelAlpha4 = 0
PureEFTmodelAlpha5 = 0
PureEFTmodelAlpha6 = 0
# 4) Designer matching EFT model selection flag:
# DesignerEFTmodel = 1 : designer f(R)
DesignerEFTmodel = 1
# 5) Matching EFT model selection flag:
# No model implemented so far.
MatchingEFTmodel = 1
####### Model parameters for EFTCAMB #######
# Notice that if the model is not selected via the model selection flags then
# the values of the parameters are automatically ignored.
# 1) Background Dark Energy equation of state parameters:
EFTw0 = -1.0001
EFTwa = -0.3
# 2) Pure EFT parameters:
EFTOmega0 = 0.001
EFTOmegaExp = 0
EFTAlpha10 = 0.1
EFTAlpha1Exp = 0
EFTAlpha20 = 0.1
EFTAlpha2Exp = 0
EFTAlpha30 = 0.001
EFTAlpha3Exp = 0
EFTAlpha40 = 0.001
EFTAlpha4Exp = 0
EFTAlpha50 = 0.1
EFTAlpha5Exp = 0
EFTAlpha60 = 0.1
EFTAlpha6Exp = 0
# 3) Designer matching parameters:
# Model 1: designer f(R) theories
EFTB0 = 1
####### Parameters for CAMB #######
# EFT flags: set up in which mode EFTCAMB is used.
# We refer to the documentation (EFTCAMB:numerical notes) for a thorough
# explanation of the effect of these flags.
#
# 1) Main EFT flag:
# EFTflag = 0 : GR code. Every EFT modification is ignored.
# EFTflag = 1 : Pure EFT code.
# EFTflag = 2 : Designer matching EFT.
# EFTflag = 3 : Matching EFT.
EFTflag = 0
# 2) Background Dark Energy equation of state flag:
# EFTwDE = 0 : Cosmological constant
# EFTwDE = 1 : DE with constant Eos determined by EFTw0
# EFTwDE = 2 : CPL parametrization
EFTwDE = 1
# 3) Pure EFT model selection flag:
# The following structure applies to every operator and can be specified for every operatori
# separately.
#
# PureEFTmodel___ = 0 : Zero (operator ignored)
# PureEFTmodel___ = 1 : Constant model
# PureEFTmodel___ = 2 : Linear model
# PureEFTmodel___ = 3 : Power law model
# PureEFTmodel___ = 4 : Exponential model
# PureEFTmodel___ = 5 : User defined
PureEFTmodelOmega = 2
PureEFTmodelAlpha1 = 0
PureEFTmodelAlpha2 = 0
PureEFTmodelAlpha3 = 0
PureEFTmodelAlpha4 = 0
PureEFTmodelAlpha5 = 0
PureEFTmodelAlpha6 = 0
# 4) Designer matching EFT model selection flag:
# DesignerEFTmodel = 1 : designer f(R)
DesignerEFTmodel = 1
# 5) Matching EFT model selection flag:
# No model implemented so far.
MatchingEFTmodel = 1
####### Model parameters for EFTCAMB #######
# Notice that if the model is not selected via the model selection flags then
# the values of the parameters are automatically ignored.
# 1) Background Dark Energy equation of state parameters:
EFTw0 = -1.0001
EFTwa = -0.3
# 2) Pure EFT parameters:
EFTOmega0 = 0.001
EFTOmegaExp = 0
EFTAlpha10 = 0.1
EFTAlpha1Exp = 0
EFTAlpha20 = 0.1
EFTAlpha2Exp = 0
EFTAlpha30 = 0.001
EFTAlpha3Exp = 0
EFTAlpha40 = 0.001
EFTAlpha4Exp = 0
EFTAlpha50 = 0.1
EFTAlpha5Exp = 0
EFTAlpha60 = 0.1
EFTAlpha6Exp = 0
# 3) Designer matching parameters:
# Model 1: designer f(R) theories
EFTB0 = 1
####### Parameters for CAMB #######
----------------------------------------------------------------------------------------------------------------------------
Modify: cmbmain.f90
::::: module CAMBmain :::::
Add: after "use Errors"
Add: after "use Errors"
! EFTCAMB MOD START
use EFTdef
use EFTSTability
use EFTreturntoGR
! EFTCAMB MOD END
use EFTdef
use EFTSTability
use EFTreturntoGR
! EFTCAMB MOD END
::::: subroutine cmbmain :::::
Add: after "integer q_ix"
Add: after "integer q_ix"
! EFTCAMB MOD START
integer EFTstabilityResult
logical DesignerSucces
! EFTCAMB MOD END
integer EFTstabilityResult
logical DesignerSucces
! EFTCAMB MOD END
Add: after "WantLateTime = CP%DoLensing .or. num_redshiftwindows > 0"
! EFTCAMB MOD START
if (EFTCAMBuseinCOSMOMC==1) then
! 1) Call designer codes.
if (CP%EFTflag==2.and.CP%DesignerEFTmodel==1) then
write(*,*) 'EFTCAMB: Calling f(R) background designer code'
call Designer_fR_Background(DesignerSucces)
if (.not.DesignerSucces) then
write(*,*) 'EFTCAMB: wrong background DE parameters.'
stop
end if
end if
! 2) Check for the stability of the theory.
if (CP%EFTflag/=0) then
write(*,*) 'EFTCAMB: Checking the stability of the theory.'
call EFTCheck_Stability(EFTStabilityResult)
if (EFTStabilityResult==0) stop 'EFTCAMB: Unstable theory.'
end if
! 3) Checks the return to GR of the theory.
if (CP%EFTflag/=0) then
call EFTCheckReturnToGR
if (DebugMsgs .and. Feedbacklevel > 0) write(*,*) 'EFTCAMB: a_pi=', EFTturnonpi
end if
end if
! EFTCAMB MOD END
if (EFTCAMBuseinCOSMOMC==1) then
! 1) Call designer codes.
if (CP%EFTflag==2.and.CP%DesignerEFTmodel==1) then
write(*,*) 'EFTCAMB: Calling f(R) background designer code'
call Designer_fR_Background(DesignerSucces)
if (.not.DesignerSucces) then
write(*,*) 'EFTCAMB: wrong background DE parameters.'
stop
end if
end if
! 2) Check for the stability of the theory.
if (CP%EFTflag/=0) then
write(*,*) 'EFTCAMB: Checking the stability of the theory.'
call EFTCheck_Stability(EFTStabilityResult)
if (EFTStabilityResult==0) stop 'EFTCAMB: Unstable theory.'
end if
! 3) Checks the return to GR of the theory.
if (CP%EFTflag/=0) then
call EFTCheckReturnToGR
if (DebugMsgs .and. Feedbacklevel > 0) write(*,*) 'EFTCAMB: a_pi=', EFTturnonpi
end if
end if
! EFTCAMB MOD END
::::: subroutine CalcScalarSources(EV,taustart) :::::
Replace: after "ind=1"
Replace: after "ind=1"
! EFTCAMB MOD START: It's always nice to plot the evolution of perturbations in DE/MG models!
! In order to print the user needs to choose the scale by means of fixq and then turno on the
! flag that prints quantities from the output routine.
! It is suggested to run the code on just one core.
if (fixq/=0._dl) then
tol1=tol/exp(AccuracyBoost-1)
write(*,*) 'EFTCAMB: start printing.'
call CreateTxtFile('Results/Debug_Output/Files/1_EFTfunctions.dat',1)
call CreateTxtFile('Results/Debug_Output/Files/2_EFTBackground.dat',2)
call CreateTxtFile('Results/Debug_Output/Files/3_FRW.dat',3)
call CreateTxtFile('Results/Debug_Output/Files/4_PiFieldCoeff.dat',4)
call CreateTxtFile('Results/Debug_Output/Files/5_PiFieldSol.dat',5)
call CreateTxtFile('Results/Debug_Output/Files/6_EinsteinEq.dat',6)
call CreateTxtFile('Results/Debug_Output/Files/7_MetricSol.dat',7)
call CreateTxtFile('Results/Debug_Output/Files/8_Sources.dat',8)
call CreateTxtFile('Results/Debug_Output/Files/9_MuGamma.dat',9)
call CreateTxtFile('Results/Debug_File/Debug.dat',10)
do j=1,10000
tauend = taustart+(j-1)*(CP%tau0-taustart)/10000
call GaugeInterface_EvolveScal(EV,tau,y,tauend,tol1,ind,c,w)
yprime = 0
call derivs(EV,EV%ScalEqsToPropagate,tau,y,yprime)
call output(EV,y,j,tau,sources)
end do
close(1)
close(2)
close(3)
close(4)
close(5)
close(6)
close(7)
close(8)
close(9)
close(10)
write(*,*) 'Stop printing'
stop
end if
! Original code:
!!Example code for plotting out variable evolution
!!if (fixq/=0._dl) then
!! tol1=tol/exp(AccuracyBoost-1)
!! call CreateTxtFile('evolve_q005.txt',1)
!! do j=1,1000
!! tauend = taustart+(j-1)*(CP%tau0-taustart)/1000
!! call GaugeInterface_EvolveScal(EV,tau,y,tauend,tol1,ind,c,w)
!! yprime = 0
!! call derivs(EV,EV%ScalEqsToPropagate,tau,y,yprime)
!! adotoa = 1/(y(1)*dtauda(y(1)))
!! ddelta= (yprime(3)*grhoc+yprime(4)*grhob)/(grhob+grhoc)
!! delta=(grhoc*y(3)+grhob*y(4))/(grhob+grhoc)
!! growth= ddelta/delta/adotoa
!! write (1,'(7E15.5)') tau, delta, growth, y(3), y(4), y(EV%g_ix), y(1)
!! end do
!! close(1)
!! stop
!!end if
! EFTCAMB MOD END
! In order to print the user needs to choose the scale by means of fixq and then turno on the
! flag that prints quantities from the output routine.
! It is suggested to run the code on just one core.
if (fixq/=0._dl) then
tol1=tol/exp(AccuracyBoost-1)
write(*,*) 'EFTCAMB: start printing.'
call CreateTxtFile('Results/Debug_Output/Files/1_EFTfunctions.dat',1)
call CreateTxtFile('Results/Debug_Output/Files/2_EFTBackground.dat',2)
call CreateTxtFile('Results/Debug_Output/Files/3_FRW.dat',3)
call CreateTxtFile('Results/Debug_Output/Files/4_PiFieldCoeff.dat',4)
call CreateTxtFile('Results/Debug_Output/Files/5_PiFieldSol.dat',5)
call CreateTxtFile('Results/Debug_Output/Files/6_EinsteinEq.dat',6)
call CreateTxtFile('Results/Debug_Output/Files/7_MetricSol.dat',7)
call CreateTxtFile('Results/Debug_Output/Files/8_Sources.dat',8)
call CreateTxtFile('Results/Debug_Output/Files/9_MuGamma.dat',9)
call CreateTxtFile('Results/Debug_File/Debug.dat',10)
do j=1,10000
tauend = taustart+(j-1)*(CP%tau0-taustart)/10000
call GaugeInterface_EvolveScal(EV,tau,y,tauend,tol1,ind,c,w)
yprime = 0
call derivs(EV,EV%ScalEqsToPropagate,tau,y,yprime)
call output(EV,y,j,tau,sources)
end do
close(1)
close(2)
close(3)
close(4)
close(5)
close(6)
close(7)
close(8)
close(9)
close(10)
write(*,*) 'Stop printing'
stop
end if
! Original code:
!!Example code for plotting out variable evolution
!!if (fixq/=0._dl) then
!! tol1=tol/exp(AccuracyBoost-1)
!! call CreateTxtFile('evolve_q005.txt',1)
!! do j=1,1000
!! tauend = taustart+(j-1)*(CP%tau0-taustart)/1000
!! call GaugeInterface_EvolveScal(EV,tau,y,tauend,tol1,ind,c,w)
!! yprime = 0
!! call derivs(EV,EV%ScalEqsToPropagate,tau,y,yprime)
!! adotoa = 1/(y(1)*dtauda(y(1)))
!! ddelta= (yprime(3)*grhoc+yprime(4)*grhob)/(grhob+grhoc)
!! delta=(grhoc*y(3)+grhob*y(4))/(grhob+grhoc)
!! growth= ddelta/delta/adotoa
!! write (1,'(7E15.5)') tau, delta, growth, y(3), y(4), y(EV%g_ix), y(1)
!! end do
!! close(1)
!! stop
!!end if
! EFTCAMB MOD END
----------------------------------------------------------------------------------------------------------------------------
Modify: Modify: modules.f90
::::: type CAMBparams :::::
Add: after "integer :: Nu_mass_numbers(max_nu) !physical number per eigenstate"
Add: after "integer :: Nu_mass_numbers(max_nu) !physical number per eigenstate"
! EFTCAMB MOD START
! 1) Definition of flags:
integer :: EFTflag
integer :: EFTwDE
integer :: PureEFTmodelOmega
integer :: PureEFTmodelAlpha1, PureEFTmodelAlpha2, PureEFTmodelAlpha3
integer :: PureEFTmodelAlpha4, PureEFTmodelAlpha5, PureEFTmodelAlpha6
integer :: DesignerEFTmodel
integer :: MatchingEFTmodel
! 2) Definition of model parameters:
real(dl) :: EFTw0, EFTwa
real(dl) :: EFTOmega0, EFTOmegaExp
real(dl) :: EFTAlpha10, EFTAlpha1Exp, EFTAlpha20, EFTAlpha2Exp
real(dl) :: EFTAlpha30, EFTAlpha3Exp, EFTAlpha40, EFTAlpha4Exp
real(dl) :: EFTAlpha50, EFTAlpha5Exp, EFTAlpha60, EFTAlpha6Exp
real(dl) :: EFTB0
! EFTCAMB MOD END
! 1) Definition of flags:
integer :: EFTflag
integer :: EFTwDE
integer :: PureEFTmodelOmega
integer :: PureEFTmodelAlpha1, PureEFTmodelAlpha2, PureEFTmodelAlpha3
integer :: PureEFTmodelAlpha4, PureEFTmodelAlpha5, PureEFTmodelAlpha6
integer :: DesignerEFTmodel
integer :: MatchingEFTmodel
! 2) Definition of model parameters:
real(dl) :: EFTw0, EFTwa
real(dl) :: EFTOmega0, EFTOmegaExp
real(dl) :: EFTAlpha10, EFTAlpha1Exp, EFTAlpha20, EFTAlpha2Exp
real(dl) :: EFTAlpha30, EFTAlpha3Exp, EFTAlpha40, EFTAlpha4Exp
real(dl) :: EFTAlpha50, EFTAlpha5Exp, EFTAlpha60, EFTAlpha6Exp
real(dl) :: EFTB0
! EFTCAMB MOD END
----------------------------------------------------------------------------------------------------------------------------
Modify: inidriver.F90
::::: program driver :::::
Add: After "use F90_UNIX"
Add: After "use F90_UNIX"
! EFTCAMB MOD START
use EFTdef
! EFTCAMB MOD END
use EFTdef
! EFTCAMB MOD END
Add: After "P%tcmb = Ini_Read_Double('temp_cmb',COBE_CMBTemp)"
! EFTCAMB MOD START
! 1) Initialization of EFTCAMB flags.
P%EFTflag = Ini_Read_Int('EFTflag',0)
P%EFTwDE = Ini_Read_Int('EFTwDE',0)
P%PureEFTmodelOmega = Ini_Read_Int('PureEFTmodelOmega',0)
P%PureEFTmodelAlpha1 = Ini_Read_Int('PureEFTmodelAlpha1',0)
P%PureEFTmodelAlpha2 = Ini_Read_Int('PureEFTmodelAlpha2',0)
P%PureEFTmodelAlpha3 = Ini_Read_Int('PureEFTmodelAlpha3',0)
P%PureEFTmodelAlpha4 = Ini_Read_Int('PureEFTmodelAlpha4',0)
P%PureEFTmodelAlpha5 = Ini_Read_Int('PureEFTmodelAlpha5',0)
P%PureEFTmodelAlpha6 = Ini_Read_Int('PureEFTmodelAlpha6',0)
P%DesignerEFTmodel = Ini_Read_Int('DesignerEFTmodel',1)
P%MatchingEFTmodel = Ini_Read_Int('MatchingEFTmodel',1)
! 2) Initialization of EFTCAMB model parameters.
P%EFTw0 = Ini_Read_Double('EFTw0')
P%EFTwa = Ini_Read_Double('EFTwa')
P%EFTOmega0 = Ini_Read_Double('EFTOmega0')
P%EFTOmegaExp = Ini_Read_Double('EFTOmegaExp')
P%EFTAlpha10 = Ini_Read_Double('EFTAlpha10')
P%EFTAlpha1Exp = Ini_Read_Double('EFTAlpha1Exp')
P%EFTAlpha20 = Ini_Read_Double('EFTAlpha20')
P%EFTAlpha2Exp = Ini_Read_Double('EFTAlpha2Exp')
P%EFTAlpha30 = Ini_Read_Double('EFTAlpha30')
P%EFTAlpha3Exp = Ini_Read_Double('EFTAlpha3Exp')
P%EFTAlpha40 = Ini_Read_Double('EFTAlpha40')
P%EFTAlpha4Exp = Ini_Read_Double('EFTAlpha4Exp')
P%EFTAlpha50 = Ini_Read_Double('EFTAlpha50')
P%EFTAlpha5Exp = Ini_Read_Double('EFTAlpha5Exp')
P%EFTAlpha60 = Ini_Read_Double('EFTAlpha60')
P%EFTAlpha6Exp = Ini_Read_Double('EFTAlpha6Exp')
P%EFTB0 = Ini_Read_Double('EFTB0')
! EFTCAMB MOD END
! 1) Initialization of EFTCAMB flags.
P%EFTflag = Ini_Read_Int('EFTflag',0)
P%EFTwDE = Ini_Read_Int('EFTwDE',0)
P%PureEFTmodelOmega = Ini_Read_Int('PureEFTmodelOmega',0)
P%PureEFTmodelAlpha1 = Ini_Read_Int('PureEFTmodelAlpha1',0)
P%PureEFTmodelAlpha2 = Ini_Read_Int('PureEFTmodelAlpha2',0)
P%PureEFTmodelAlpha3 = Ini_Read_Int('PureEFTmodelAlpha3',0)
P%PureEFTmodelAlpha4 = Ini_Read_Int('PureEFTmodelAlpha4',0)
P%PureEFTmodelAlpha5 = Ini_Read_Int('PureEFTmodelAlpha5',0)
P%PureEFTmodelAlpha6 = Ini_Read_Int('PureEFTmodelAlpha6',0)
P%DesignerEFTmodel = Ini_Read_Int('DesignerEFTmodel',1)
P%MatchingEFTmodel = Ini_Read_Int('MatchingEFTmodel',1)
! 2) Initialization of EFTCAMB model parameters.
P%EFTw0 = Ini_Read_Double('EFTw0')
P%EFTwa = Ini_Read_Double('EFTwa')
P%EFTOmega0 = Ini_Read_Double('EFTOmega0')
P%EFTOmegaExp = Ini_Read_Double('EFTOmegaExp')
P%EFTAlpha10 = Ini_Read_Double('EFTAlpha10')
P%EFTAlpha1Exp = Ini_Read_Double('EFTAlpha1Exp')
P%EFTAlpha20 = Ini_Read_Double('EFTAlpha20')
P%EFTAlpha2Exp = Ini_Read_Double('EFTAlpha2Exp')
P%EFTAlpha30 = Ini_Read_Double('EFTAlpha30')
P%EFTAlpha3Exp = Ini_Read_Double('EFTAlpha3Exp')
P%EFTAlpha40 = Ini_Read_Double('EFTAlpha40')
P%EFTAlpha4Exp = Ini_Read_Double('EFTAlpha4Exp')
P%EFTAlpha50 = Ini_Read_Double('EFTAlpha50')
P%EFTAlpha5Exp = Ini_Read_Double('EFTAlpha5Exp')
P%EFTAlpha60 = Ini_Read_Double('EFTAlpha60')
P%EFTAlpha6Exp = Ini_Read_Double('EFTAlpha6Exp')
P%EFTB0 = Ini_Read_Double('EFTB0')
! EFTCAMB MOD END
----------------------------------------------------------------------------------------------------------------------------
Add: equations_EFT.f90or Modify: equations.f90 to get equations_EFT.f90
::::: function dtauda(a) :::::
Add: after " use LambdaGeneral"
Add: after " use LambdaGeneral"
! EFTCAMB MOD START
use EFTdef
use EFTfunctions
! EFTCAMB MOD END
use EFTdef
use EFTfunctions
! EFTCAMB MOD END
Replace: after "grhoa2=grhok*a2+(grhoc+grhob)*a+grhog+grhornomass"
! EFTCAMB MOD START
if (CP%EFTflag==0) then ! Original CAMB code.
if (w_lam == -1._dl) then
grhoa2=grhoa2+grhov*a2**2
else
grhoa2=grhoa2+grhov*a**(1-3*w_lam)
end if
else if (CP%EFTflag/=0) then
grhoa2=grhoa2+grhov*EFTw(a,2)*a2
end if
!... original:
! if (w_lam == -1._dl) then
! grhoa2=grhoa2+grhov*a2**2
! else
! grhoa2=grhoa2+grhov*a**(1-3*w_lam)
! end if
! EFTCAMB MOD END
if (CP%EFTflag==0) then ! Original CAMB code.
if (w_lam == -1._dl) then
grhoa2=grhoa2+grhov*a2**2
else
grhoa2=grhoa2+grhov*a**(1-3*w_lam)
end if
else if (CP%EFTflag/=0) then
grhoa2=grhoa2+grhov*EFTw(a,2)*a2
end if
!... original:
! if (w_lam == -1._dl) then
! grhoa2=grhoa2+grhov*a2**2
! else
! grhoa2=grhoa2+grhov*a**(1-3*w_lam)
! end if
! EFTCAMB MOD END
::::: module GaugeInterface :::::
Add: After "use Errors"
Add: After "use Errors"
! EFTCAMB MOD START
use EFTdef
use EFTfunctions
! EFTCAMB MOD END
use EFTdef
use EFTfunctions
! EFTCAMB MOD END
Add: After "real(dl) denlkt(4,max_l_evolve),Kft(max_l_evolve)"
! EFTCAMB MOD START
logical EFTCAMBactive
real(dl) EFTturnOnTime
! EFTCAMB MOD END
logical EFTCAMBactive
real(dl) EFTturnOnTime
! EFTCAMB MOD END
::::: recursive subroutine GaugeInterface_EvolveScal(EV,tau,y,tauend,tol1,ind,c,w) :::::
Add: After "tau_switch_nu_massive= noSwitch"
Add: After "tau_switch_nu_massive= noSwitch"
! EFTCAMB MOD START: add the dark energy switch.
if (CP%EFTflag==0) EV%EFTturnOnTime = noSwitch
! EFTCAMB MOD END
if (CP%EFTflag==0) EV%EFTturnOnTime = noSwitch
! EFTCAMB MOD END
Replace: After "tau_switch_no_phot_multpoles =max(15/EV%k_buf,taurend)*AccuracyBoost"
! EFTCAMB MOD START: add the dark energy switch.
next_switch = min(tau_switch_ktau, tau_switch_nu_massless,EV%TightSwitchoffTime, tau_switch_nu_massive, &
tau_switch_no_nu_multpoles, tau_switch_no_phot_multpoles, tau_switch_nu_nonrel,noSwitch, EV%EFTturnOnTime)
!... Original:
! next_switch = min(tau_switch_ktau, tau_switch_nu_massless,EV%TightSwitchoffTime, tau_switch_nu_massive, &
! tau_switch_no_nu_multpoles, tau_switch_no_phot_multpoles, tau_switch_nu_nonrel,noSwitch)
! EFTCAMB MOD END
next_switch = min(tau_switch_ktau, tau_switch_nu_massless,EV%TightSwitchoffTime, tau_switch_nu_massive, &
tau_switch_no_nu_multpoles, tau_switch_no_phot_multpoles, tau_switch_nu_nonrel,noSwitch, EV%EFTturnOnTime)
!... Original:
! next_switch = min(tau_switch_ktau, tau_switch_nu_massless,EV%TightSwitchoffTime, tau_switch_nu_massive, &
! tau_switch_no_nu_multpoles, tau_switch_no_phot_multpoles, tau_switch_nu_nonrel,noSwitch)
! EFTCAMB MOD END
Add: After "else if (next_switch==tau_switch_no_phot_multpoles) then"
! EFTCAMB MOD START: effect of the DE switch.
else if (next_switch==EV%EFTturnOnTime) then
ind=1
EVout%EFTCAMBactive = .true.
EVout%EFTturnOnTime = noSwitch
call SetupScalarArrayIndices(EVout)
call CopyScalarVariableArray(y,yout, EV, EVout)
y=yout
EV=EVout
call EFTpiInitialConditions(y,EV,tau)
! EFTCAMB MOD STOP
else if (next_switch==EV%EFTturnOnTime) then
ind=1
EVout%EFTCAMBactive = .true.
EVout%EFTturnOnTime = noSwitch
call SetupScalarArrayIndices(EVout)
call CopyScalarVariableArray(y,yout, EV, EVout)
y=yout
EV=EVout
call EFTpiInitialConditions(y,EV,tau)
! EFTCAMB MOD STOP
::::: subroutine SetupScalarArrayIndices(EV, max_num_eqns) :::::
Replace: After "maxeq = maxeq + (EV%lmaxg+1)+(EV%lmaxnr+1)+EV%lmaxgpol-1"
Replace: After "maxeq = maxeq + (EV%lmaxg+1)+(EV%lmaxnr+1)+EV%lmaxgpol-1"
! EFTCAMB MOD START
if ((w_lam /= -1 .and. w_Perturb).or.CP%EFTflag/=0) then
EV%w_ix = neq+1
neq=neq+2
maxeq=maxeq+2
else
EV%w_ix=0
end if
!... Original:
! !Dark energy
! if (w_lam /= -1 .and. w_Perturb) then
! EV%w_ix = neq+1
! neq=neq+2
! maxeq=maxeq+2
! else
! EV%w_ix=0
! end if
! EFTCAMB MOD END
if ((w_lam /= -1 .and. w_Perturb).or.CP%EFTflag/=0) then
EV%w_ix = neq+1
neq=neq+2
maxeq=maxeq+2
else
EV%w_ix=0
end if
!... Original:
! !Dark energy
! if (w_lam /= -1 .and. w_Perturb) then
! EV%w_ix = neq+1
! neq=neq+2
! maxeq=maxeq+2
! else
! EV%w_ix=0
! end if
! EFTCAMB MOD END
::::: subroutine CopyScalarVariableArray(y,yout, EV, EVout) :::::
Replace: After "yout(1:basic_num_eqns) = y(1:basic_num_eqns)"
Replace: After "yout(1:basic_num_eqns) = y(1:basic_num_eqns)"
! EFTCAMB MOD START
if ((w_lam /= -1.and.w_Perturb).or.&
&(CP%EFTflag/=0.and.EVout%EFTCAMBactive).or.&
&(CP%EFTflag/=0.and.EV%EFTCAMBactive)) then
yout(EVout%w_ix)=y(EV%w_ix)
yout(EVout%w_ix+1)=y(EV%w_ix+1)
end if
!... Original:
! if (w_lam /= -1 .and. w_Perturb) then
! yout(EVout%w_ix)=y(EV%w_ix)
! yout(EVout%w_ix+1)=y(EV%w_ix+1)
! end if
! EFTCAMB MOD END
if ((w_lam /= -1.and.w_Perturb).or.&
&(CP%EFTflag/=0.and.EVout%EFTCAMBactive).or.&
&(CP%EFTflag/=0.and.EV%EFTCAMBactive)) then
yout(EVout%w_ix)=y(EV%w_ix)
yout(EVout%w_ix+1)=y(EV%w_ix+1)
end if
!... Original:
! if (w_lam /= -1 .and. w_Perturb) then
! yout(EVout%w_ix)=y(EV%w_ix)
! yout(EVout%w_ix+1)=y(EV%w_ix+1)
! end if
! EFTCAMB MOD END
::::: subroutine GetNumEqns(EV) :::::
Add: Before "if (HighAccuracyDefault .and. CP%AccuratePolarization) then"
Add: Before "if (HighAccuracyDefault .and. CP%AccuratePolarization) then"
! EFTCAMB MOD START:
EV%EFTCAMBactive = .false.
! EFTCAMB MOD END
EV%EFTCAMBactive = .false.
! EFTCAMB MOD END
::::: subroutine output(EV,y, j,tau,sources) :::::
Add: After "real(dl) ISW"
Add: After "real(dl) ISW"
! EFTCAMB MOD START: definitions of EFTCAMB quantities
! background quantities. Designer approach.
real(dl) EFT_H0,Hdot, Hdotdot, adotdota
! 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 ! pi field equations
real(dl) EFTeomF, EFTeomN, EFTeomX, EFTeomY, EFTeomG, EFTeomU ! equations of motion
real(dl) EFTeomL, EFTeomM, EFTeomV, EFTeomNdot, EFTeomVdot, EFTeomXdot
real(dl) pidotdot
! EFT quantities relevant for sources:
real(dl) EFTISW, EFTLensing, EFTsigmadot, polterdot
! mu and gamma
real(dl) EFT_Psi, EFT_Phi, EFT_function_mu, EFT_function_gamma, EFT_xi
! Debug quantities
real(dl) temp1, temp2, temp3,temp4, temp5, temp6, temp7
! EFTCAMB MOD END
! background quantities. Designer approach.
real(dl) EFT_H0,Hdot, Hdotdot, adotdota
! 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 ! pi field equations
real(dl) EFTeomF, EFTeomN, EFTeomX, EFTeomY, EFTeomG, EFTeomU ! equations of motion
real(dl) EFTeomL, EFTeomM, EFTeomV, EFTeomNdot, EFTeomVdot, EFTeomXdot
real(dl) pidotdot
! EFT quantities relevant for sources:
real(dl) EFTISW, EFTLensing, EFTsigmadot, polterdot
! mu and gamma
real(dl) EFT_Psi, EFT_Phi, EFT_function_mu, EFT_function_gamma, EFT_xi
! Debug quantities
real(dl) temp1, temp2, temp3,temp4, temp5, temp6, temp7
! EFTCAMB MOD END
Replace: Before "grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhov_t"
! EFTCAMB MOD START: Computation of DE background density.
if (CP%EFTflag==0) then
grhov_t=grhov*a**(-1-3*w_lam)
else if (CP%EFTflag/=0) then
grhov_t=grhov*EFTw(a,2)
EFT_H0 = (CP%h0/c_EFT)*1000._dl
end if
!... Original:
! grhov_t=grhov*a**(-1-3*w_lam)
! EFTCAMB MOD END
if (CP%EFTflag==0) then
grhov_t=grhov*a**(-1-3*w_lam)
else if (CP%EFTflag/=0) then
grhov_t=grhov*EFTw(a,2)
EFT_H0 = (CP%h0/c_EFT)*1000._dl
end if
!... Original:
! grhov_t=grhov*a**(-1-3*w_lam)
! EFTCAMB MOD END
Replace: After "grho=grhob_t+grhoc_t+grhor_t+grhog_t+grhov_t"
! EFTCAMB MOD START: computation of gpres.
if (CP%EFTflag==0) then ! Normal CAMB code
gpres=(grhog_t+grhor_t)/3+grhov_t*w_lam
else if (CP%EFTflag/=0) then
gpres=(grhog_t+grhor_t)/3._dl +EFTw(a,0)*grhov_t
end if
!... Original:
! gpres=(grhog_t+grhor_t)/3+grhov_t*w_lam
! EFTCAMB MOD END
if (CP%EFTflag==0) then ! Normal CAMB code
gpres=(grhog_t+grhor_t)/3+grhov_t*w_lam
else if (CP%EFTflag/=0) then
gpres=(grhog_t+grhor_t)/3._dl +EFTw(a,0)*grhov_t
end if
!... Original:
! gpres=(grhog_t+grhor_t)/3+grhov_t*w_lam
! EFTCAMB MOD END
Replace: After "call MassiveNuVarsOut(EV,y,yprime,a,grho,gpres,dgrho,dgq,dgpi, dgpi_diff,pidot_sum)"
! EFTCAMB MOD START: initialization of DE perturbations.
if (CP%EFTflag==0) then
if (w_lam /= -1 .and. w_Perturb) then
clxq=y(EV%w_ix)
vq=y(EV%w_ix+1)
dgrho=dgrho + clxq*grhov_t
dgq = dgq + vq*grhov_t*(1+w_lam)
end if
else if (CP%EFTflag/=0.and.EV%EFTCAMBactive) then
clxq=y(EV%w_ix)
vq=y(EV%w_ix+1)
pidotdot = yprime(EV%w_ix+1)
end if
!... Original:
! if (w_lam /= -1 .and. w_Perturb) then
! clxq=y(EV%w_ix)
! vq=y(EV%w_ix+1)
! dgrho=dgrho + clxq*grhov_t
! dgq = dgq + vq*grhov_t*(1+w_lam)
! end if
! EFTCAMB MOD END
if (CP%EFTflag==0) then
if (w_lam /= -1 .and. w_Perturb) then
clxq=y(EV%w_ix)
vq=y(EV%w_ix+1)
dgrho=dgrho + clxq*grhov_t
dgq = dgq + vq*grhov_t*(1+w_lam)
end if
else if (CP%EFTflag/=0.and.EV%EFTCAMBactive) then
clxq=y(EV%w_ix)
vq=y(EV%w_ix+1)
pidotdot = yprime(EV%w_ix+1)
end if
!... Original:
! if (w_lam /= -1 .and. w_Perturb) then
! clxq=y(EV%w_ix)
! vq=y(EV%w_ix+1)
! dgrho=dgrho + clxq*grhov_t
! dgq = dgq + vq*grhov_t*(1+w_lam)
! end if
! EFTCAMB MOD END
Replace: After "adotoa=sqrt((grho+grhok)/3)" Before "if (EV%no_nu_multpoles) then"
! EFTCAMB MOD START: Computation of background quantities:
! 1) FRW background quantities.
if (CP%EFTflag/=0.and.a>=EFTturnonpi) then
adotdota=(adotoa*adotoa-gpres)/2._dl
Hdot =adotdota-adotoa**2._dl
Hdotdot = 2._dl*adotoa*Hdot &
& + 0.5_dl*adotoa*(grhob_t + grhoc_t + 8._dl*(grhog_t+grhor_t)/3._dl)&
& + 0.5_dl*adotoa*grhov_t*((1._dl+EFTw(a,0))*(1._dl+3._dl*EFTw(a,0))-a*EFTw(a,1))
end if
! 2) Compute and write once for all the values of the EFT functions.
if (CP%EFTflag/=0 .and. a>=EFTturnonpi*0.1_dl) then
EFTOmegaV = EFTOmega(a,0)
EFTOmegaP = EFTOmega(a,1)
EFTOmegaPP = EFTOmega(a,2)
EFTOmegaPPP = EFTOmega(a,3)
EFTAlpha1V = EFTAlpha1(a,0)
EFTAlpha1P = EFTAlpha1(a,1)
EFTAlpha2V = EFTAlpha2(a,0)
EFTAlpha2P = EFTAlpha2(a,1)
EFTAlpha3V = EFTAlpha3(a,0)
EFTAlpha3P = EFTAlpha3(a,1)
EFTAlpha4V = EFTAlpha4(a,0)
EFTAlpha4P = EFTAlpha4(a,1)
EFTAlpha4PP = EFTAlpha4(a,2)
EFTAlpha5V = EFTAlpha5(a,0)
EFTAlpha5P = EFTAlpha5(a,1)
EFTAlpha6V = EFTAlpha6(a,0)
EFTAlpha6P = EFTAlpha6(a,1)
end if
! 3) Computation of EFT background quantities.
if (CP%EFTflag/=0.and.a>=EFTturnonpi) then
if (CP%EFTflag==1) then
!8*pi*G*c*a^2
EFTc = (adotoa*adotoa - Hdot)*(EFTOmegaV + a*EFTOmegaP*0.5_dl) &
& - 0.5_dl*a2*adotoa*adotoa*EFTOmegaPP&
& + 0.5_dl*grhov_t*(1._dl+EFTw(a,0))
!8*pi*G*Lambda*a^2
EFTLambda = +EFTw(a,0)*grhov_t &
&-EFTOmegaV*(2._dl*Hdot+adotoa**2._dl) &
&-a*EFTOmegaP*(2._dl*adotoa**2._dl + Hdot) &
&-a2*adotoa**2._dl*EFTOmegaPP
!EFT C DOT: 8*pi*G*cdot*a^2
EFTcdot = -EFTOmegaV*(Hdotdot-4._dl*adotoa*Hdot+2._dl*adotoa*adotoa*adotoa) &
& + 0.5_dl*a*EFTOmegaP*(-Hdotdot+adotoa*Hdot+adotoa*adotoa*adotoa)&
& +0.5_dl*a2*adotoa*EFTOmegaPP*(adotoa*adotoa-3._dl*Hdot)&
& -0.5_dl*a*a2*adotoa*adotoa*adotoa*EFTOmegaPPP&
& +0.5_dl*adotoa*grhov_t*(-3._dl*(1._dl+EFTw(a,0))**2 + a*EFTw(a,1))
!EFT LAMBDA DOT: 8*pi*G*Ldot*a^2
EFTLambdadot = -2._dl*EFTOmegaV*(Hdotdot-adotoa*Hdot-adotoa*adotoa*adotoa)&
& - a*EFTOmegaP*(4._dl*adotoa*Hdot+Hdotdot)&
& -a2*EFTOmegaPP*adotoa*(3._dl*Hdot+2._dl*adotoa*adotoa)&
& -a*a2*EFTOmegaPPP*adotoa*adotoa*adotoa&
& +grhov_t*adotoa*(a*EFTw(a,1)-3._dl*EFTw(a,0)*(1._dl+EFTw(a,0)))
else if (CP%EFTflag==2) then
EFTc = EFTcTemp(a,0)
EFTcdot = EFTcTemp(a,1)
EFTLambda = EFTLambdaTemp(a,0)
EFTLambdadot = EFTLambdaTemp(a,1)
end if
EFTgrhoq = 2._dl*EFTc -EFTLambda -3._dl*a*adotoa*adotoa*EFTOmegaP
EFTgpresq = EFTLambda + a2*adotoa*adotoa*EFTOmegaPP&
& +a*EFTOmegaP*(Hdot + 2._dl*adotoa*adotoa)
EFTgrhodotq = -3._dl*adotoa*(EFTgrhoq+EFTgpresq) + 3._dl*a*adotoa**3._dl*EFTOmegaP
EFTgpresdotq = EFTLambdadot +a2*a*adotoa**3*EFTOmegaPPP + 3._dl*a2*adotoa*Hdot*EFTOmegaPP&
& +a*EFTOmegaP*Hdotdot +3._dl*a*adotoa*Hdot*EFTOmegaP +2._dl*a2*adotoa**3*EFTOmegaPP&
& -2._dl*a*adotoa**3*EFTOmegaP
end if
! Debug code. Print every background EFT quantity.
if (a>= EFTturnonpi*0.1_dl.and.DebugEFT) then
write (1,'(25e15.5)') a, tau, EFTOmegaV, EFTOmegaP, EFTOmegaPP, EFTOmegaPPP,&
& EFTAlpha1V, EFTAlpha1P, EFTAlpha2V, EFTAlpha2P, EFTAlpha3V, EFTAlpha3P,&
& EFTAlpha4V, EFTAlpha4P, EFTAlpha4PP, EFTAlpha5V, EFTAlpha5P, EFTAlpha6V, EFTAlpha6P
write (2,'(12e15.6)') a, tau, EFTc, EFTLambda, EFTcdot, EFTLambdadot, EFTgrhoq, EFTgpresq, EFTgrhodotq, EFTgpresdotq
write (3,'(8e15.5)') a, tau, adotoa, Hdot, Hdotdot, grhov_t
end if
! EFTCAMB MOD END
! EFTCAMB MOD START: computation of Einstein equations factors.
if (CP%EFTflag/=0.and.a>=EFTturnonpi) then
EFTeomF = 1.5_dl/(k*(1._dl+EFTOmegaV))*((EFTgrhoQ+EFTgpresQ)*clxq& ! Background op
& +adotoa*adotoa*a*EFTOmegaP*clxq+a*adotoa*EFTOmegaP*vq)&
& +1.5_dl*a*EFT_H0*EFTAlpha2V**3*(vq+adotoa*clxq)/(k*(1._dl+EFTOmegaV))& ! Alpha2
& -1.5_dl*EFTAlpha3V**2/(1._dl+EFTOmegaV)*(k-3._dl*(Hdot-adotoa**2)/(k))*clxq& ! ALpha3
& +1.5_dl*EFTAlpha4V**2/(1._dl+EFTOmegaV)*k*clxq& ! Alpha4
& -1.5_dl*EFTAlpha4V**2/(1._dl+EFTOmegaV)*(Hdot-adotoa**2)/k*clxq
EFTeomN = -a*adotoa*EFTOmegaP/(1._dl+EFTOmegaV)*k*clxq& ! Background op
& +2._dl*adotoa/(1._dl+EFTOmegaV)*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P)*k*clxq& !Alpha4
& +EFTAlpha4V**2/(1._dl+EFTOmegaV)*k*vq&
& +2._dl*EFTAlpha5V**2*k*(vq+adotoa*clxq)/(1._dl+EFTOmegaV) ! Alpha5
EFTeomNdot = -a*Hdot*EFTOmegaP/(1._dl+EFTOmegaV)*k*clxq& ! Backgrpund op
& -a*adotoa*EFTOmegaP/(1._dl+EFTOmegaV)*k*vq&
& -a*adotoa**2/(1._dl+EFTOmegaV)*(EFTOmegaP+a*EFTOmegaPP&
& -a*EFTOmegaP**2/(1._dl+EFTOmegaV))*k*clxq&
& +EFTAlpha4V**2*k*Pidotdot/(1._dl+EFTOmegaV)+ a*adotoa*k*vq/(1._dl+EFTOmegaV)*&! Alpha4
&(2._dl*EFTAlpha4V*EFTAlpha4P-EFTAlpha4V**2*EFTOmegaP/(1._dl+EFTOmegaV))&
& +2._dl*k/(1._dl+EFTOmegaV)*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P)*(Hdot*clxq+adotoa*vq)&
& +2._dl*a*adotoa**2*k*clxq/(1._dl+EFTOmegaV)*(a*EFTAlpha4V*EFTAlpha4PP+a*EFTAlpha4P**2&
& +3._dl*EFTAlpha4V*EFTAlpha4P-EFTOmegaP/(1._dl+EFTOmegaV)*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P))&
& +2._dl*EFTAlpha5V**2*k/(1._dl+EFTOmegaV)*(Pidotdot + adotoa*vq + Hdot*clxq)& ! Alpha5
& +2._dl*a*k*adotoa/(1._dl+EFTOmegaV)*(vq+adotoa*clxq)*(2._dl*EFTAlpha5V*EFTAlpha5P&
& -EFTAlpha5V**2*EFTOmegaP/(1._dl+EFTOmegaV))
EFTeomX = 1._dl& ! Backgrond op
& -EFTAlpha4V**2/(1._dl+EFTOmegaV) ! Alpha4
EFTeomXdot = -a*adotoa/(1._dl+EFTOmegaV)*(2._dl*EFTAlpha4V*EFTAlpha4P& ! Alpha4
& -EFTAlpha4V**2*EFTOmegaP/(1._dl+EFTOmegaV))
EFTeomY = +0.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)& ! Background op
& +1.5_dl/(1._dl+EFTOmegaV)*(EFTAlpha3V**2+a*EFTAlpha3V*EFTAlpha3P)& ! Alpha3
& +0.5_dl*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P)/(1._dl+EFTOmegaV) ! Alpha4
EFTeomG = +1._dl + 0.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)& ! Background op
& +0.5_dl*a*EFT_H0*EFTAlpha2V**3/adotoa/(1._dl+EFTOmegaV)& ! Alpha2
& +1.5_dl*EFTAlpha3V**2/(1._dl+EFTOmegaV)& ! Alpha3
& +EFTAlpha4V**2/(1._dl+EFTOmegaV) ! Alpha4
EFTeomU = 1._dl& ! Background op
& +1.5_dl*EFTAlpha3V**2/(1._dl+EFTOmegaV)& ! Alpha3
& +0.5_dl*EFTAlpha4V**2/(1._dl+EFTOmegaV) ! Alpha4
EFTeomL = -1.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)*(3._dl*adotoa*adotoa-Hdot)*clxq& ! Background op
& -1.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)*adotoa*vq&
& -0.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)*k2*clxq&
& +0.5_dl*clxq/(adotoa*(1._dl+EFTOmegaV))*EFTgrhodotQ&
& +(vq+adotoa*clxq)/(adotoa*(1+EFTOmegaV))*EFTc&
& +2._dl*a2*EFT_H0**2*EFTAlpha1V**4/adotoa/(1._dl+EFTOmegaV)*(vq+adotoa*clxq)& ! Alpha1
& +1.5_dl*a*EFT_H0*EFTAlpha2V**3/(1._dl+EFTOmegaV)*& ! Alpha2
&((Hdot/adotoa-2._dl*adotoa-k2/(3._dl*adotoa))*clxq-vq)&
& -1.5_dl*EFTAlpha3V**2/(1._dl+EFTOmegaV)*(k2-3._dl*(Hdot-adotoa**2))*clxq& ! Alpha3
& +3._dl*EFTAlpha4V**2/(1._dl+EFTOmegaV)*(Hdot-adotoa**2-k2/3._dl)*clxq& ! Alpha4
& +4._dl*EFTAlpha6V**2*k2/adotoa*(vq+adotoa*clxq)/(1._dl+EFTOmegaV) ! Alpha6
EFTeomM = EFTgpresdotQ*clxq + (EFTgrhoQ+ EFTgpresQ)*(vq+adotoa*clxq)& ! Background op
& +a2*adotoa**2*EFTOmegaPP*vq +a2*adotoa**3*EFTOmegaPP*clxq&
& +a*adotoa*EFTOmegaP*(pidotdot + (Hdot/adotoa+4._dl*adotoa)*vq&
& + (2._dl*Hdot+6._dl*adotoa**2 +2._dl/3._dl*k2)*clxq)&
& +a*EFT_H0*EFTAlpha2V**2*(EFTAlpha2V*pidotdot& ! Alpha2
& +(4._dl*EFTAlpha2V+3._dl*a*EFTAlpha2P)*adotoa*vq +(3._dl*adotoa**2*EFTAlpha2V&
& +Hdot*EFTAlpha2V +3._dl*a*adotoa**2*EFTAlpha2P)*clxq)&
& +EFTAlpha3V**2*(3._dl*adotoa**2-3._dl*Hdot+k2)*vq& ! Alpha3
& +EFTAlpha3V**2*(6._dl*adotoa**3-3._dl*Hdotdot)*clxq&
& +2._dl*adotoa*k2*clxq*(EFTAlpha3V**2+a*EFTAlpha3V*EFTAlpha3P)&
& -6._dl*a*adotoa*(Hdot-adotoa**2)*EFTAlpha3V*EFTAlpha3P*clxq&
& -EFTAlpha4V**2*(Hdot-adotoa**2-k2/3._dl)*vq& ! Alpha4
& -2._dl*adotoa*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P)*(Hdot-adotoa**2-k2/3._dl)*clxq&
& -EFTAlpha4V**2*(Hdotdot-2*adotoa*Hdot)*clxq&
& -4._dl*EFTAlpha5V**2*k2*(vq+adotoa*clxq)/3._dl ! Alpha5
EFTeomV = +0.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)& ! Background op
& -(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P)/(1._dl+EFTOmegaV) ! Alpha4
EFTeomVdot = +0.5_dl*a*adotoa/(1._dl+EFTOmegaV)*(EFTOmegaP+a*EFTOmegaPP& ! Background op
& -a*EFTOmegaP**2/(1._dl+EFTOmegaV))&
& -a*adotoa/(1._dl+EFTOmegaV)*(a*EFTAlpha4V*EFTAlpha4PP+a*EFTAlpha4P**2& ! Alpha4
& +3._dl*EFTAlpha4V*EFTAlpha4P-EFTOmegaP/(1._dl+EFTOmegaV)*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P))
end if
! Debug code. Prints every EFT quantity defined so far.
if (a>= EFTturnonpi .and. DebugEFTO) then
write (6,'(15E15.5)') a, tau, EFTeomF, EFTeomN, EFTeomX, EFTeomY, EFTeomG, EFTeomU,&
& EFTeomL, EFTeomM, EFTeomNdot, EFTeomVdot, EFTeomXdot
end if
! EFTCAMB MOD END
! EFTCAMB MOD START: compute z,dz before loading radiation and photon
if (CP%EFTflag==0.or.a
z=(0.5_dl*dgrho/k + etak)/adotoa
dz= -adotoa*z - 0.5_dl*dgrho/k
else if (CP%EFTflag/=0.and.a>=EFTturnonpi) then
! RSA approximation.
z= 1._dl/EFTeomG*(etak/adotoa + 0.5_dl*dgrho/(k*adotoa*(1._dl+EFTOmegaV)) + EFTeomL/(k*EFT_H0))
dz= 1._dl/EFTeomU*(-2._dl*adotoa*z*(1._dl+EFTeomY-0.5_dl*EFTeomG) &
& - 0.5_dl*dgrho/k/(1._dl+EFTOmegaV) -adotoa*EFTeomL/(k*EFT_H0) -1.5_dl/k/(1._dl+EFTOmegaV)*EFTeomM/EFT_H0)
end if
!equations of motion 1.
if (EV%no_nu_multpoles) then
clxr=-4*dz/k
qr=-4._dl/3*z
pir=0
pirdot=0
else
clxr=y(EV%r_ix)
qr =y(EV%r_ix+1)
pir =y(EV%r_ix+2)
pirdot=yprime(EV%r_ix+2)
end if
if (EV%no_phot_multpoles) then
clxg=-4*dz/k -4/k*opac(j)*(vb+z)
qg=-4._dl/3*z
pig=0
pigdot=0
octg=0
octgprime=0
qgdot = -4*dz/3
else
if (EV%TightCoupling) then
pig = EV%pig
!pigdot=EV%pigdot
if (second_order_tightcoupling) then
octg = (3._dl/7._dl)*pig*(EV%k_buf/opac(j))
ypol(2) = EV%pig/4 + EV%pigdot*(1._dl/opac(j))*(-5._dl/8._dl)
ypol(3) = (3._dl/7._dl)*(EV%k_buf/opac(j))*ypol(2)
else
ypol(2) = EV%pig/4
octg=0
end if
octgprime=0
else
pig =y(EV%g_ix+2)
pigdot=yprime(EV%g_ix+2)
octg=y(EV%g_ix+3)
octgprime=yprime(EV%g_ix+3)
end if
clxg=y(EV%g_ix)
qg =y(EV%g_ix+1)
qgdot =yprime(EV%g_ix+1)
end if
!... Original:
! if (EV%no_nu_multpoles) then
! z=(0.5_dl*dgrho/k + etak)/adotoa
! dz= -adotoa*z - 0.5_dl*dgrho/k
! clxr=-4*dz/k
! qr=-4._dl/3*z
! pir=0
! pirdot=0
! else
! clxr=y(EV%r_ix)
! qr =y(EV%r_ix+1)
! pir =y(EV%r_ix+2)
! pirdot=yprime(EV%r_ix+2)
! end if
! if (EV%no_phot_multpoles) then
! z=(0.5_dl*dgrho/k + etak)/adotoa
! dz= -adotoa*z - 0.5_dl*dgrho/k
! clxg=-4*dz/k -4/k*opac(j)*(vb+z)
! qg=-4._dl/3*z
! pig=0
! pigdot=0
! octg=0
! octgprime=0
! qgdot = -4*dz/3
! else
! if (EV%TightCoupling) then
! pig = EV%pig
! !pigdot=EV%pigdot
! if (second_order_tightcoupling) then
! octg = (3._dl/7._dl)*pig*(EV%k_buf/opac(j))
! ypol(2) = EV%pig/4 + EV%pigdot*(1._dl/opac(j))*(-5._dl/8._dl)
! ypol(3) = (3._dl/7._dl)*(EV%k_buf/opac(j))*ypol(2)
! else
! ypol(2) = EV%pig/4
! octg=0
! end if
! octgprime=0
! else
! pig =y(EV%g_ix+2)
! pigdot=yprime(EV%g_ix+2)
! octg=y(EV%g_ix+3)
! octgprime=yprime(EV%g_ix+3)
! end if
! clxg=y(EV%g_ix)
! qg =y(EV%g_ix+1)
! qgdot =yprime(EV%g_ix+1)
! end if
! EFTCAMB MOD END
! 1) FRW background quantities.
if (CP%EFTflag/=0.and.a>=EFTturnonpi) then
adotdota=(adotoa*adotoa-gpres)/2._dl
Hdot =adotdota-adotoa**2._dl
Hdotdot = 2._dl*adotoa*Hdot &
& + 0.5_dl*adotoa*(grhob_t + grhoc_t + 8._dl*(grhog_t+grhor_t)/3._dl)&
& + 0.5_dl*adotoa*grhov_t*((1._dl+EFTw(a,0))*(1._dl+3._dl*EFTw(a,0))-a*EFTw(a,1))
end if
! 2) Compute and write once for all the values of the EFT functions.
if (CP%EFTflag/=0 .and. a>=EFTturnonpi*0.1_dl) then
EFTOmegaV = EFTOmega(a,0)
EFTOmegaP = EFTOmega(a,1)
EFTOmegaPP = EFTOmega(a,2)
EFTOmegaPPP = EFTOmega(a,3)
EFTAlpha1V = EFTAlpha1(a,0)
EFTAlpha1P = EFTAlpha1(a,1)
EFTAlpha2V = EFTAlpha2(a,0)
EFTAlpha2P = EFTAlpha2(a,1)
EFTAlpha3V = EFTAlpha3(a,0)
EFTAlpha3P = EFTAlpha3(a,1)
EFTAlpha4V = EFTAlpha4(a,0)
EFTAlpha4P = EFTAlpha4(a,1)
EFTAlpha4PP = EFTAlpha4(a,2)
EFTAlpha5V = EFTAlpha5(a,0)
EFTAlpha5P = EFTAlpha5(a,1)
EFTAlpha6V = EFTAlpha6(a,0)
EFTAlpha6P = EFTAlpha6(a,1)
end if
! 3) Computation of EFT background quantities.
if (CP%EFTflag/=0.and.a>=EFTturnonpi) then
if (CP%EFTflag==1) then
!8*pi*G*c*a^2
EFTc = (adotoa*adotoa - Hdot)*(EFTOmegaV + a*EFTOmegaP*0.5_dl) &
& - 0.5_dl*a2*adotoa*adotoa*EFTOmegaPP&
& + 0.5_dl*grhov_t*(1._dl+EFTw(a,0))
!8*pi*G*Lambda*a^2
EFTLambda = +EFTw(a,0)*grhov_t &
&-EFTOmegaV*(2._dl*Hdot+adotoa**2._dl) &
&-a*EFTOmegaP*(2._dl*adotoa**2._dl + Hdot) &
&-a2*adotoa**2._dl*EFTOmegaPP
!EFT C DOT: 8*pi*G*cdot*a^2
EFTcdot = -EFTOmegaV*(Hdotdot-4._dl*adotoa*Hdot+2._dl*adotoa*adotoa*adotoa) &
& + 0.5_dl*a*EFTOmegaP*(-Hdotdot+adotoa*Hdot+adotoa*adotoa*adotoa)&
& +0.5_dl*a2*adotoa*EFTOmegaPP*(adotoa*adotoa-3._dl*Hdot)&
& -0.5_dl*a*a2*adotoa*adotoa*adotoa*EFTOmegaPPP&
& +0.5_dl*adotoa*grhov_t*(-3._dl*(1._dl+EFTw(a,0))**2 + a*EFTw(a,1))
!EFT LAMBDA DOT: 8*pi*G*Ldot*a^2
EFTLambdadot = -2._dl*EFTOmegaV*(Hdotdot-adotoa*Hdot-adotoa*adotoa*adotoa)&
& - a*EFTOmegaP*(4._dl*adotoa*Hdot+Hdotdot)&
& -a2*EFTOmegaPP*adotoa*(3._dl*Hdot+2._dl*adotoa*adotoa)&
& -a*a2*EFTOmegaPPP*adotoa*adotoa*adotoa&
& +grhov_t*adotoa*(a*EFTw(a,1)-3._dl*EFTw(a,0)*(1._dl+EFTw(a,0)))
else if (CP%EFTflag==2) then
EFTc = EFTcTemp(a,0)
EFTcdot = EFTcTemp(a,1)
EFTLambda = EFTLambdaTemp(a,0)
EFTLambdadot = EFTLambdaTemp(a,1)
end if
EFTgrhoq = 2._dl*EFTc -EFTLambda -3._dl*a*adotoa*adotoa*EFTOmegaP
EFTgpresq = EFTLambda + a2*adotoa*adotoa*EFTOmegaPP&
& +a*EFTOmegaP*(Hdot + 2._dl*adotoa*adotoa)
EFTgrhodotq = -3._dl*adotoa*(EFTgrhoq+EFTgpresq) + 3._dl*a*adotoa**3._dl*EFTOmegaP
EFTgpresdotq = EFTLambdadot +a2*a*adotoa**3*EFTOmegaPPP + 3._dl*a2*adotoa*Hdot*EFTOmegaPP&
& +a*EFTOmegaP*Hdotdot +3._dl*a*adotoa*Hdot*EFTOmegaP +2._dl*a2*adotoa**3*EFTOmegaPP&
& -2._dl*a*adotoa**3*EFTOmegaP
end if
! Debug code. Print every background EFT quantity.
if (a>= EFTturnonpi*0.1_dl.and.DebugEFT) then
write (1,'(25e15.5)') a, tau, EFTOmegaV, EFTOmegaP, EFTOmegaPP, EFTOmegaPPP,&
& EFTAlpha1V, EFTAlpha1P, EFTAlpha2V, EFTAlpha2P, EFTAlpha3V, EFTAlpha3P,&
& EFTAlpha4V, EFTAlpha4P, EFTAlpha4PP, EFTAlpha5V, EFTAlpha5P, EFTAlpha6V, EFTAlpha6P
write (2,'(12e15.6)') a, tau, EFTc, EFTLambda, EFTcdot, EFTLambdadot, EFTgrhoq, EFTgpresq, EFTgrhodotq, EFTgpresdotq
write (3,'(8e15.5)') a, tau, adotoa, Hdot, Hdotdot, grhov_t
end if
! EFTCAMB MOD END
! EFTCAMB MOD START: computation of Einstein equations factors.
if (CP%EFTflag/=0.and.a>=EFTturnonpi) then
EFTeomF = 1.5_dl/(k*(1._dl+EFTOmegaV))*((EFTgrhoQ+EFTgpresQ)*clxq& ! Background op
& +adotoa*adotoa*a*EFTOmegaP*clxq+a*adotoa*EFTOmegaP*vq)&
& +1.5_dl*a*EFT_H0*EFTAlpha2V**3*(vq+adotoa*clxq)/(k*(1._dl+EFTOmegaV))& ! Alpha2
& -1.5_dl*EFTAlpha3V**2/(1._dl+EFTOmegaV)*(k-3._dl*(Hdot-adotoa**2)/(k))*clxq& ! ALpha3
& +1.5_dl*EFTAlpha4V**2/(1._dl+EFTOmegaV)*k*clxq& ! Alpha4
& -1.5_dl*EFTAlpha4V**2/(1._dl+EFTOmegaV)*(Hdot-adotoa**2)/k*clxq
EFTeomN = -a*adotoa*EFTOmegaP/(1._dl+EFTOmegaV)*k*clxq& ! Background op
& +2._dl*adotoa/(1._dl+EFTOmegaV)*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P)*k*clxq& !Alpha4
& +EFTAlpha4V**2/(1._dl+EFTOmegaV)*k*vq&
& +2._dl*EFTAlpha5V**2*k*(vq+adotoa*clxq)/(1._dl+EFTOmegaV) ! Alpha5
EFTeomNdot = -a*Hdot*EFTOmegaP/(1._dl+EFTOmegaV)*k*clxq& ! Backgrpund op
& -a*adotoa*EFTOmegaP/(1._dl+EFTOmegaV)*k*vq&
& -a*adotoa**2/(1._dl+EFTOmegaV)*(EFTOmegaP+a*EFTOmegaPP&
& -a*EFTOmegaP**2/(1._dl+EFTOmegaV))*k*clxq&
& +EFTAlpha4V**2*k*Pidotdot/(1._dl+EFTOmegaV)+ a*adotoa*k*vq/(1._dl+EFTOmegaV)*&! Alpha4
&(2._dl*EFTAlpha4V*EFTAlpha4P-EFTAlpha4V**2*EFTOmegaP/(1._dl+EFTOmegaV))&
& +2._dl*k/(1._dl+EFTOmegaV)*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P)*(Hdot*clxq+adotoa*vq)&
& +2._dl*a*adotoa**2*k*clxq/(1._dl+EFTOmegaV)*(a*EFTAlpha4V*EFTAlpha4PP+a*EFTAlpha4P**2&
& +3._dl*EFTAlpha4V*EFTAlpha4P-EFTOmegaP/(1._dl+EFTOmegaV)*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P))&
& +2._dl*EFTAlpha5V**2*k/(1._dl+EFTOmegaV)*(Pidotdot + adotoa*vq + Hdot*clxq)& ! Alpha5
& +2._dl*a*k*adotoa/(1._dl+EFTOmegaV)*(vq+adotoa*clxq)*(2._dl*EFTAlpha5V*EFTAlpha5P&
& -EFTAlpha5V**2*EFTOmegaP/(1._dl+EFTOmegaV))
EFTeomX = 1._dl& ! Backgrond op
& -EFTAlpha4V**2/(1._dl+EFTOmegaV) ! Alpha4
EFTeomXdot = -a*adotoa/(1._dl+EFTOmegaV)*(2._dl*EFTAlpha4V*EFTAlpha4P& ! Alpha4
& -EFTAlpha4V**2*EFTOmegaP/(1._dl+EFTOmegaV))
EFTeomY = +0.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)& ! Background op
& +1.5_dl/(1._dl+EFTOmegaV)*(EFTAlpha3V**2+a*EFTAlpha3V*EFTAlpha3P)& ! Alpha3
& +0.5_dl*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P)/(1._dl+EFTOmegaV) ! Alpha4
EFTeomG = +1._dl + 0.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)& ! Background op
& +0.5_dl*a*EFT_H0*EFTAlpha2V**3/adotoa/(1._dl+EFTOmegaV)& ! Alpha2
& +1.5_dl*EFTAlpha3V**2/(1._dl+EFTOmegaV)& ! Alpha3
& +EFTAlpha4V**2/(1._dl+EFTOmegaV) ! Alpha4
EFTeomU = 1._dl& ! Background op
& +1.5_dl*EFTAlpha3V**2/(1._dl+EFTOmegaV)& ! Alpha3
& +0.5_dl*EFTAlpha4V**2/(1._dl+EFTOmegaV) ! Alpha4
EFTeomL = -1.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)*(3._dl*adotoa*adotoa-Hdot)*clxq& ! Background op
& -1.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)*adotoa*vq&
& -0.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)*k2*clxq&
& +0.5_dl*clxq/(adotoa*(1._dl+EFTOmegaV))*EFTgrhodotQ&
& +(vq+adotoa*clxq)/(adotoa*(1+EFTOmegaV))*EFTc&
& +2._dl*a2*EFT_H0**2*EFTAlpha1V**4/adotoa/(1._dl+EFTOmegaV)*(vq+adotoa*clxq)& ! Alpha1
& +1.5_dl*a*EFT_H0*EFTAlpha2V**3/(1._dl+EFTOmegaV)*& ! Alpha2
&((Hdot/adotoa-2._dl*adotoa-k2/(3._dl*adotoa))*clxq-vq)&
& -1.5_dl*EFTAlpha3V**2/(1._dl+EFTOmegaV)*(k2-3._dl*(Hdot-adotoa**2))*clxq& ! Alpha3
& +3._dl*EFTAlpha4V**2/(1._dl+EFTOmegaV)*(Hdot-adotoa**2-k2/3._dl)*clxq& ! Alpha4
& +4._dl*EFTAlpha6V**2*k2/adotoa*(vq+adotoa*clxq)/(1._dl+EFTOmegaV) ! Alpha6
EFTeomM = EFTgpresdotQ*clxq + (EFTgrhoQ+ EFTgpresQ)*(vq+adotoa*clxq)& ! Background op
& +a2*adotoa**2*EFTOmegaPP*vq +a2*adotoa**3*EFTOmegaPP*clxq&
& +a*adotoa*EFTOmegaP*(pidotdot + (Hdot/adotoa+4._dl*adotoa)*vq&
& + (2._dl*Hdot+6._dl*adotoa**2 +2._dl/3._dl*k2)*clxq)&
& +a*EFT_H0*EFTAlpha2V**2*(EFTAlpha2V*pidotdot& ! Alpha2
& +(4._dl*EFTAlpha2V+3._dl*a*EFTAlpha2P)*adotoa*vq +(3._dl*adotoa**2*EFTAlpha2V&
& +Hdot*EFTAlpha2V +3._dl*a*adotoa**2*EFTAlpha2P)*clxq)&
& +EFTAlpha3V**2*(3._dl*adotoa**2-3._dl*Hdot+k2)*vq& ! Alpha3
& +EFTAlpha3V**2*(6._dl*adotoa**3-3._dl*Hdotdot)*clxq&
& +2._dl*adotoa*k2*clxq*(EFTAlpha3V**2+a*EFTAlpha3V*EFTAlpha3P)&
& -6._dl*a*adotoa*(Hdot-adotoa**2)*EFTAlpha3V*EFTAlpha3P*clxq&
& -EFTAlpha4V**2*(Hdot-adotoa**2-k2/3._dl)*vq& ! Alpha4
& -2._dl*adotoa*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P)*(Hdot-adotoa**2-k2/3._dl)*clxq&
& -EFTAlpha4V**2*(Hdotdot-2*adotoa*Hdot)*clxq&
& -4._dl*EFTAlpha5V**2*k2*(vq+adotoa*clxq)/3._dl ! Alpha5
EFTeomV = +0.5_dl*a*EFTOmegaP/(1._dl+EFTOmegaV)& ! Background op
& -(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P)/(1._dl+EFTOmegaV) ! Alpha4
EFTeomVdot = +0.5_dl*a*adotoa/(1._dl+EFTOmegaV)*(EFTOmegaP+a*EFTOmegaPP& ! Background op
& -a*EFTOmegaP**2/(1._dl+EFTOmegaV))&
& -a*adotoa/(1._dl+EFTOmegaV)*(a*EFTAlpha4V*EFTAlpha4PP+a*EFTAlpha4P**2& ! Alpha4
& +3._dl*EFTAlpha4V*EFTAlpha4P-EFTOmegaP/(1._dl+EFTOmegaV)*(EFTAlpha4V**2+a*EFTAlpha4V*EFTAlpha4P))
end if
! Debug code. Prints every EFT quantity defined so far.
if (a>= EFTturnonpi .and. DebugEFTO) then
write (6,'(15E15.5)') a, tau, EFTeomF, EFTeomN, EFTeomX, EFTeomY, EFTeomG, EFTeomU,&
& EFTeomL, EFTeomM, EFTeomNdot, EFTeomVdot, EFTeomXdot
end if
! EFTCAMB MOD END
! EFTCAMB MOD START: compute z,dz before loading radiation and photon
if (CP%EFTflag==0.or.a
dz= -adotoa*z - 0.5_dl*dgrho/k
else if (CP%EFTflag/=0.and.a>=EFTturnonpi) then
! RSA approximation.
z= 1._dl/EFTeomG*(etak/adotoa + 0.5_dl*dgrho/(k*adotoa*(1._dl+EFTOmegaV)) + EFTeomL/(k*EFT_H0))
dz= 1._dl/EFTeomU*(-2._dl*adotoa*z*(1._dl+EFTeomY-0.5_dl*EFTeomG) &
& - 0.5_dl*dgrho/k/(1._dl+EFTOmegaV) -adotoa*EFTeomL/(k*EFT_H0) -1.5_dl/k/(1._dl+EFTOmegaV)*EFTeomM/EFT_H0)
end if
!equations of motion 1.
if (EV%no_nu_multpoles) then
clxr=-4*dz/k
qr=-4._dl/3*z
pir=0
pirdot=0
else
clxr=y(EV%r_ix)
qr =y(EV%r_ix+1)
pir =y(EV%r_ix+2)
pirdot=yprime(EV%r_ix+2)
end if
if (EV%no_phot_multpoles) then
clxg=-4*dz/k -4/k*opac(j)*(vb+z)
qg=-4._dl/3*z
pig=0
pigdot=0
octg=0
octgprime=0
qgdot = -4*dz/3
else
if (EV%TightCoupling) then
pig = EV%pig
!pigdot=EV%pigdot
if (second_order_tightcoupling) then
octg = (3._dl/7._dl)*pig*(EV%k_buf/opac(j))
ypol(2) = EV%pig/4 + EV%pigdot*(1._dl/opac(j))*(-5._dl/8._dl)
ypol(3) = (3._dl/7._dl)*(EV%k_buf/opac(j))*ypol(2)
else
ypol(2) = EV%pig/4
octg=0
end if
octgprime=0
else
pig =y(EV%g_ix+2)
pigdot=yprime(EV%g_ix+2)
octg=y(EV%g_ix+3)
octgprime=yprime(EV%g_ix+3)
end if
clxg=y(EV%g_ix)
qg =y(EV%g_ix+1)
qgdot =yprime(EV%g_ix+1)
end if
! if (EV%no_nu_multpoles) then
! z=(0.5_dl*dgrho/k + etak)/adotoa
! dz= -adotoa*z - 0.5_dl*dgrho/k
! clxr=-4*dz/k
! qr=-4._dl/3*z
! pir=0
! pirdot=0
! else
! clxr=y(EV%r_ix)
! qr =y(EV%r_ix+1)
! pir =y(EV%r_ix+2)
! pirdot=yprime(EV%r_ix+2)
! end if
! if (EV%no_phot_multpoles) then
! z=(0.5_dl*dgrho/k + etak)/adotoa
! dz= -adotoa*z - 0.5_dl*dgrho/k
! clxg=-4*dz/k -4/k*opac(j)*(vb+z)
! qg=-4._dl/3*z
! pig=0
! pigdot=0
! octg=0
! octgprime=0
! qgdot = -4*dz/3
! else
! if (EV%TightCoupling) then
! pig = EV%pig
! !pigdot=EV%pigdot
! if (second_order_tightcoupling) then
! octg = (3._dl/7._dl)*pig*(EV%k_buf/opac(j))
! ypol(2) = EV%pig/4 + EV%pigdot*(1._dl/opac(j))*(-5._dl/8._dl)
! ypol(3) = (3._dl/7._dl)*(EV%k_buf/opac(j))*ypol(2)
! else
! ypol(2) = EV%pig/4
! octg=0
! end if
! octgprime=0
! else
! pig =y(EV%g_ix+2)
! pigdot=yprime(EV%g_ix+2)
! octg=y(EV%g_ix+3)
! octgprime=yprime(EV%g_ix+3)
! end if
! clxg=y(EV%g_ix)
! qg =y(EV%g_ix+1)
! qgdot =yprime(EV%g_ix+1)
! end if
! EFTCAMB MOD END
Replace: After "dgpi = dgpi + grhor_t*pir + grhog_t*pig"
! EFTCAMB MOD START: equation of motion 3.
! Get sigma (shear) and z from the constraints
! have to get z from eta for numerical stability
if (CP%EFTflag==0 .or. a
z=(0.5_dl*dgrho/k + etak)/adotoa
sigma=(z+1.5_dl*dgq/k2)/EV%Kf(1)
else if (CP%EFTflag/=0 .and. a>=EFTturnonpi) then
z= 1._dl/EFTeomG*(etak/adotoa + 0.5_dl*dgrho/(k*adotoa*(1._dl+EFTOmegaV)) + EFTeomL/(k*EFT_H0))
dz=1._dl/EFTeomU*(-2*adotoa*(1._dl+EFTeomY)*z +etak -0.5_dl/k/(1._dl+EFTOmegaV)*(grhog_t*clxg+grhor_t*clxr)&
& -1.5_dl/k/(1._dl+EFTOmegaV)*EFTeomM/EFT_H0)
sigma= 1._dl/EFTeomX*(z*EFTeomU +1.5_dl*dgq/(k2*(1._dl+EFTOmegaV)) +EFTeomF/EFT_H0)
end if
!... Original:
! ! Get sigma (shear) and z from the constraints
! ! have to get z from eta for numerical stability
! z=(0.5_dl*dgrho/k + etak)/adotoa
! sigma=(z+1.5_dl*dgq/k2)/EV%Kf(1)
! EFTCAMB MOD END
! Get sigma (shear) and z from the constraints
! have to get z from eta for numerical stability
if (CP%EFTflag==0 .or. a
sigma=(z+1.5_dl*dgq/k2)/EV%Kf(1)
else if (CP%EFTflag/=0 .and. a>=EFTturnonpi) then
z= 1._dl/EFTeomG*(etak/adotoa + 0.5_dl*dgrho/(k*adotoa*(1._dl+EFTOmegaV)) + EFTeomL/(k*EFT_H0))
dz=1._dl/EFTeomU*(-2*adotoa*(1._dl+EFTeomY)*z +etak -0.5_dl/k/(1._dl+EFTOmegaV)*(grhog_t*clxg+grhor_t*clxr)&
& -1.5_dl/k/(1._dl+EFTOmegaV)*EFTeomM/EFT_H0)
sigma= 1._dl/EFTeomX*(z*EFTeomU +1.5_dl*dgq/(k2*(1._dl+EFTOmegaV)) +EFTeomF/EFT_H0)
end if
! ! Get sigma (shear) and z from the constraints
! ! have to get z from eta for numerical stability
! z=(0.5_dl*dgrho/k + etak)/adotoa
! sigma=(z+1.5_dl*dgq/k2)/EV%Kf(1)
! EFTCAMB MOD END
Replace: After "diff_rhopi = pidot_sum - (4*dgpi+ dgpi_diff )*adotoa" Before "! Doppler term"
! EFTCAMB MOD START: computation of quantities used to construct the sources.
if (CP%EFTflag/=0 .and. a>=EFTturnonpi) then
EFTsigmadot= 1._dl/EFTeomX*(-2._dl*adotoa*(1._dl+EFTeomV)*sigma +etak&
& -1._dl/k*dgpi/(1._dl+EFTOmegaV) +EFTeomN/EFT_H0)
EFTISW = 1._dl/EFTeomX*(-2._dl*(1._dl+EFTeomV)*(Hdot*sigma+ adotoa*EFTsigmadot)&
& -2._dl*adotoa*sigma*EFTeomVdot + 0.5_dl*dgq*(1._dl+EFTeomX)/(1._dl+EFTOmegaV)&
& +a*adotoa*EFTOmegaP/k*dgpi/(1._dl+EFTOmegaV)**2&
& -1._dl/(1._dl+EFTOmegaV)*(2._dl*adotoa*dgpi/k+diff_rhopi/k)&
& +(1._dl+EFTeomX)*k2/3._dl/EFT_H0*EFTeomF-EFTeomXdot*EFTsigmadot+EFTeomNdot/EFT_H0)
EFTLensing = 1._dl/EFTeomX*(-2._dl*adotoa*(1._dl+EFTeomV)*sigma +(1._dl+EFTeomX)*etak&
& -1._dl/k*dgpi/(1._dl+EFTOmegaV) +EFTeomN/EFT_H0)
end if
! EFTCAMB MOD END
! EFTCAMB MOD START: TT source
if (CP%EFTflag==0.or.a
!Maple's fortran output - see scal_eqs.map
!2phi' term (\phi' + \psi' in Newtonian gauge)
ISW = (4.D0/3.D0*k*EV%Kf(1)*sigma+(-2.D0/3.D0*sigma-2.D0/3.D0*etak/adotoa)*k &
-diff_rhopi/k**2-1.D0/adotoa*dgrho/3.D0+(3.D0*gpres+5.D0*grho)*sigma/k/3.D0 &
-2.D0/k*adotoa/EV%Kf(1)*etak)*expmmu(j)
!e.g. to get only late-time ISW
! if (1/a-1 < 30) ISW=0
!The rest, note y(9)->octg, yprime(9)->octgprime (octopoles)
sources(1)= ISW + ((-9.D0/160.D0*pig-27.D0/80.D0*ypol(2))/k**2*opac(j)+(11.D0/10.D0*sigma- &
3.D0/8.D0*EV%Kf(2)*ypol(3)+vb-9.D0/80.D0*EV%Kf(2)*octg+3.D0/40.D0*qg)/k-(- &
180.D0*ypolprime(2)-30.D0*pigdot)/k**2/160.D0)*dvis(j)+(-(9.D0*pigdot+ &
54.D0*ypolprime(2))/k**2*opac(j)/160.D0+pig/16.D0+clxg/4.D0+3.D0/8.D0*ypol(2)+(- &
21.D0/5.D0*adotoa*sigma-3.D0/8.D0*EV%Kf(2)*ypolprime(3)+vbdot+3.D0/40.D0*qgdot- &
9.D0/80.D0*EV%Kf(2)*octgprime)/k+(-9.D0/160.D0*dopac(j)*pig-21.D0/10.D0*dgpi-27.D0/ &
80.D0*dopac(j)*ypol(2))/k**2)*vis(j)+(3.D0/16.D0*ddvis(j)*pig+9.D0/ &
8.D0*ddvis(j)*ypol(2))/k**2+21.D0/10.D0/k/EV%Kf(1)*vis(j)*etak
else if (CP%EFTflag/=0 .and. a>=EFTturnonpi) then
polterdot = 9._dl/15._dl*ypolprime(2) + 0.1_dl*pigdot
ISW = expmmu(j)*(EFTISW)/k
sources(1) = ISW&
& +vis(j)* (clxg/4.D0 + polter/1.6d0 + vbdot/k -9.D0*(polterdot)/k2*opac(j)/16.D0 -9.D0/16.D0*dopac(j)*polter/k2&
& + 21.D0/10.D0*EFTsigmadot/k + 3.D0/40.D0*qgdot/k &
& +(-3.D0/8.D0*EV%Kf(2)*ypolprime(3) - 9.D0/80.D0*EV%Kf(2)*octgprime)/k)&
& +((-9.D0/160.D0*pig-27.D0/80.D0*ypol(2))/k**2*opac(j)+(11.D0/10.D0*sigma- &
& 3.D0/8.D0*EV%Kf(2)*ypol(3)+vb-9.D0/80.D0*EV%Kf(2)*octg+3.D0/40.D0*qg)/k-(- &
& 180.D0*ypolprime(2)-30.D0*pigdot)/k**2/160.D0)*dvis(j)&
& +(3.D0/16.D0*ddvis(j)*pig+9.D0/8.D0*ddvis(j)*ypol(2))/k**2
end if
! Debug code. Print quantities used to compute the sources.
if (a>= EFTturnonpi .and. DebugEFTO) then
write (7,'(10e15.5)') a, tau, z, sigma, EFTsigmadot, EFTISW, EFTLensing
end if
!... Original:
! !Maple's fortran output - see scal_eqs.map
! !2phi' term (\phi' + \psi' in Newtonian gauge)
! ISW = (4.D0/3.D0*k*EV%Kf(1)*sigma+(-2.D0/3.D0*sigma-2.D0/3.D0*etak/adotoa)*k &
! -diff_rhopi/k**2-1.D0/adotoa*dgrho/3.D0+(3.D0*gpres+5.D0*grho)*sigma/k/3.D0 &
! -2.D0/k*adotoa/EV%Kf(1)*etak)*expmmu(j)
! !e.g. to get only late-time ISW
! ! if (1/a-1 < 30) ISW=0
! !The rest, note y(9)->octg, yprime(9)->octgprime (octopoles)
! sources(1)= ISW + ((-9.D0/160.D0*pig-27.D0/80.D0*ypol(2))/k**2*opac(j)+ &
! (11.D0/10.D0*sigma- 3.D0/8.D0*EV%Kf(2)*ypol(3)+vb-9.D0/80.D0*EV%Kf(2)*octg+3.D0/40.D0*qg)/k- &
! (-180.D0*ypolprime(2)-30.D0*pigdot)/k**2/160.D0)*dvis(j) + &
! (-(9.D0*pigdot+ 54.D0*ypolprime(2))/k**2*opac(j)/160.D0+pig/16.D0+clxg/4.D0+3.D0/8.D0*ypol(2) + &
! (-21.D0/5.D0*adotoa*sigma-3.D0/8.D0*EV%Kf(2)*ypolprime(3) + &
! vbdot+3.D0/40.D0*qgdot- 9.D0/80.D0*EV%Kf(2)*octgprime)/k + &
! (-9.D0/160.D0*dopac(j)*pig-21.D0/10.D0*dgpi-27.D0/80.D0*dopac(j)*ypol(2))/k**2)*vis(j) + &
! (3.D0/16.D0*ddvis(j)*pig+9.D0/8.D0*ddvis(j)*ypol(2))/k**2+21.D0/10.D0/k/EV%Kf(1)*vis(j)*etak
! EFTCAMB MOD END
if (CP%EFTflag/=0 .and. a>=EFTturnonpi) then
EFTsigmadot= 1._dl/EFTeomX*(-2._dl*adotoa*(1._dl+EFTeomV)*sigma +etak&
& -1._dl/k*dgpi/(1._dl+EFTOmegaV) +EFTeomN/EFT_H0)
EFTISW = 1._dl/EFTeomX*(-2._dl*(1._dl+EFTeomV)*(Hdot*sigma+ adotoa*EFTsigmadot)&
& -2._dl*adotoa*sigma*EFTeomVdot + 0.5_dl*dgq*(1._dl+EFTeomX)/(1._dl+EFTOmegaV)&
& +a*adotoa*EFTOmegaP/k*dgpi/(1._dl+EFTOmegaV)**2&
& -1._dl/(1._dl+EFTOmegaV)*(2._dl*adotoa*dgpi/k+diff_rhopi/k)&
& +(1._dl+EFTeomX)*k2/3._dl/EFT_H0*EFTeomF-EFTeomXdot*EFTsigmadot+EFTeomNdot/EFT_H0)
EFTLensing = 1._dl/EFTeomX*(-2._dl*adotoa*(1._dl+EFTeomV)*sigma +(1._dl+EFTeomX)*etak&
& -1._dl/k*dgpi/(1._dl+EFTOmegaV) +EFTeomN/EFT_H0)
end if
! EFTCAMB MOD END
! EFTCAMB MOD START: TT source
if (CP%EFTflag==0.or.a
!2phi' term (\phi' + \psi' in Newtonian gauge)
ISW = (4.D0/3.D0*k*EV%Kf(1)*sigma+(-2.D0/3.D0*sigma-2.D0/3.D0*etak/adotoa)*k &
-diff_rhopi/k**2-1.D0/adotoa*dgrho/3.D0+(3.D0*gpres+5.D0*grho)*sigma/k/3.D0 &
-2.D0/k*adotoa/EV%Kf(1)*etak)*expmmu(j)
!e.g. to get only late-time ISW
! if (1/a-1 < 30) ISW=0
!The rest, note y(9)->octg, yprime(9)->octgprime (octopoles)
sources(1)= ISW + ((-9.D0/160.D0*pig-27.D0/80.D0*ypol(2))/k**2*opac(j)+(11.D0/10.D0*sigma- &
3.D0/8.D0*EV%Kf(2)*ypol(3)+vb-9.D0/80.D0*EV%Kf(2)*octg+3.D0/40.D0*qg)/k-(- &
180.D0*ypolprime(2)-30.D0*pigdot)/k**2/160.D0)*dvis(j)+(-(9.D0*pigdot+ &
54.D0*ypolprime(2))/k**2*opac(j)/160.D0+pig/16.D0+clxg/4.D0+3.D0/8.D0*ypol(2)+(- &
21.D0/5.D0*adotoa*sigma-3.D0/8.D0*EV%Kf(2)*ypolprime(3)+vbdot+3.D0/40.D0*qgdot- &
9.D0/80.D0*EV%Kf(2)*octgprime)/k+(-9.D0/160.D0*dopac(j)*pig-21.D0/10.D0*dgpi-27.D0/ &
80.D0*dopac(j)*ypol(2))/k**2)*vis(j)+(3.D0/16.D0*ddvis(j)*pig+9.D0/ &
8.D0*ddvis(j)*ypol(2))/k**2+21.D0/10.D0/k/EV%Kf(1)*vis(j)*etak
else if (CP%EFTflag/=0 .and. a>=EFTturnonpi) then
polterdot = 9._dl/15._dl*ypolprime(2) + 0.1_dl*pigdot
ISW = expmmu(j)*(EFTISW)/k
sources(1) = ISW&
& +vis(j)* (clxg/4.D0 + polter/1.6d0 + vbdot/k -9.D0*(polterdot)/k2*opac(j)/16.D0 -9.D0/16.D0*dopac(j)*polter/k2&
& + 21.D0/10.D0*EFTsigmadot/k + 3.D0/40.D0*qgdot/k &
& +(-3.D0/8.D0*EV%Kf(2)*ypolprime(3) - 9.D0/80.D0*EV%Kf(2)*octgprime)/k)&
& +((-9.D0/160.D0*pig-27.D0/80.D0*ypol(2))/k**2*opac(j)+(11.D0/10.D0*sigma- &
& 3.D0/8.D0*EV%Kf(2)*ypol(3)+vb-9.D0/80.D0*EV%Kf(2)*octg+3.D0/40.D0*qg)/k-(- &
& 180.D0*ypolprime(2)-30.D0*pigdot)/k**2/160.D0)*dvis(j)&
& +(3.D0/16.D0*ddvis(j)*pig+9.D0/8.D0*ddvis(j)*ypol(2))/k**2
end if
! Debug code. Print quantities used to compute the sources.
if (a>= EFTturnonpi .and. DebugEFTO) then
write (7,'(10e15.5)') a, tau, z, sigma, EFTsigmadot, EFTISW, EFTLensing
end if
! !Maple's fortran output - see scal_eqs.map
! !2phi' term (\phi' + \psi' in Newtonian gauge)
! ISW = (4.D0/3.D0*k*EV%Kf(1)*sigma+(-2.D0/3.D0*sigma-2.D0/3.D0*etak/adotoa)*k &
! -diff_rhopi/k**2-1.D0/adotoa*dgrho/3.D0+(3.D0*gpres+5.D0*grho)*sigma/k/3.D0 &
! -2.D0/k*adotoa/EV%Kf(1)*etak)*expmmu(j)
! !e.g. to get only late-time ISW
! ! if (1/a-1 < 30) ISW=0
! !The rest, note y(9)->octg, yprime(9)->octgprime (octopoles)
! sources(1)= ISW + ((-9.D0/160.D0*pig-27.D0/80.D0*ypol(2))/k**2*opac(j)+ &
! (11.D0/10.D0*sigma- 3.D0/8.D0*EV%Kf(2)*ypol(3)+vb-9.D0/80.D0*EV%Kf(2)*octg+3.D0/40.D0*qg)/k- &
! (-180.D0*ypolprime(2)-30.D0*pigdot)/k**2/160.D0)*dvis(j) + &
! (-(9.D0*pigdot+ 54.D0*ypolprime(2))/k**2*opac(j)/160.D0+pig/16.D0+clxg/4.D0+3.D0/8.D0*ypol(2) + &
! (-21.D0/5.D0*adotoa*sigma-3.D0/8.D0*EV%Kf(2)*ypolprime(3) + &
! vbdot+3.D0/40.D0*qgdot- 9.D0/80.D0*EV%Kf(2)*octgprime)/k + &
! (-9.D0/160.D0*dopac(j)*pig-21.D0/10.D0*dgpi-27.D0/80.D0*dopac(j)*ypol(2))/k**2)*vis(j) + &
! (3.D0/16.D0*ddvis(j)*pig+9.D0/8.D0*ddvis(j)*ypol(2))/k**2+21.D0/10.D0/k/EV%Kf(1)*vis(j)*etak
! EFTCAMB MOD END
Replace: After "sources(2)=0"
! EFTCAMB MOD START: Lensing source
if (CTransScal%NumSources > 2) then
!Get lensing sources
!Can modify this here if you want to get power spectra for other tracer
if (tau>taurend .and. CP%tau0-tau > 0.1_dl) then
if (CP%EFTflag==0 .or. a
!phi_lens = Phi - 1/2 kappa (a/k)^2 sum_i rho_i pi_i
!Neglect pi contributions because not accurate at late time anyway
phi = -(dgrho +3*dgq*adotoa/k)/(k2*EV%Kf(1)*2)
! - (grhor_t*pir + grhog_t*pig+ pinu*gpnu_t)/k2
sources(3) = -2*phi*f_K(tau-tau_maxvis)/(f_K(CP%tau0-tau_maxvis)*f_K(CP%tau0-tau))
!sources(3) = -2*phi*(tau-tau_maxvis)/((CP%tau0-tau_maxvis)*(CP%tau0-tau))
!We include the lensing factor of two here
else if (CP%EFTflag/=0 .and. a>=EFTturnonpi) then
sources(3) = -(EFTLensing)/k*f_K(tau-tau_maxvis)/(f_K(CP%tau0-tau_maxvis)*f_K(CP%tau0-tau))
end if
! Debug code:
if (a>EFTturnonpi.and.DebugEFTO) then
write (8,'(6e15.5)') a, tau, sources(1), sources(3)
end if
else
sources(3) = 0
end if
end if
!... Original:
! if (CTransScal%NumSources > 2) then
! !Get lensing sources
! !Can modify this here if you want to get power spectra for other tracer
! if (tau>taurend .and. CP%tau0-tau > 0.1_dl) then
! !phi_lens = Phi - 1/2 kappa (a/k)^2 sum_i rho_i pi_i
! !Neglect pi contributions because not accurate at late time anyway
! phi = -(dgrho +3*dgq*adotoa/k)/(k2*EV%Kf(1)*2)
! ! - dgpi/k2/2
! sources(3) = -2*phi*f_K(tau-tau_maxvis)/(f_K(CP%tau0-tau_maxvis)*f_K(CP%tau0-tau))
! ! sources(3) = -2*phi*(tau-tau_maxvis)/((CP%tau0-tau_maxvis)*(CP%tau0-tau))
! !We include the lensing factor of two here
! else
! sources(3) = 0
! end if
! end if
! EFTCAMB MOD END
! EFTCAMB MOD START: computation of mu and gamma. These are not used but it's nice to have the possibility to plot them...
if (CP%EFTflag/=0.and.a>=EFTturnonpi) then
EFT_Psi = EFTsigmadot/k + adotoa*sigma/k
EFT_Phi = etak/k - adotoa*sigma/k
EFT_function_mu = -2._dl*k*(EFTsigmadot+adotoa*sigma)/(dgrho)
EFT_function_gamma = (etak - adotoa*sigma)/(EFTsigmadot + adotoa*sigma)
EFT_xi = vq/(adotoa*clxq)
! This is used to print mu and gamma for different k values.
! If you don't want to print all the other quantities then
! turn off DenugEFTO and turn the other to true.
if (DebugEFTO.or..false.) then
write (9,'(10e15.5)') a, tau, 1._dl/a-1._dl, k, EFT_function_mu, EFT_function_gamma, EFT_xi
end if
end if
! EFTCAMB MOD END
if (CTransScal%NumSources > 2) then
!Get lensing sources
!Can modify this here if you want to get power spectra for other tracer
if (tau>taurend .and. CP%tau0-tau > 0.1_dl) then
if (CP%EFTflag==0 .or. a
!Neglect pi contributions because not accurate at late time anyway
phi = -(dgrho +3*dgq*adotoa/k)/(k2*EV%Kf(1)*2)
! - (grhor_t*pir + grhog_t*pig+ pinu*gpnu_t)/k2
sources(3) = -2*phi*f_K(tau-tau_maxvis)/(f_K(CP%tau0-tau_maxvis)*f_K(CP%tau0-tau))
!sources(3) = -2*phi*(tau-tau_maxvis)/((CP%tau0-tau_maxvis)*(CP%tau0-tau))
!We include the lensing factor of two here
else if (CP%EFTflag/=0 .and. a>=EFTturnonpi) then
sources(3) = -(EFTLensing)/k*f_K(tau-tau_maxvis)/(f_K(CP%tau0-tau_maxvis)*f_K(CP%tau0-tau))
end if
! Debug code:
if (a>EFTturnonpi.and.DebugEFTO) then
write (8,'(6e15.5)') a, tau, sources(1), sources(3)
end if
else
sources(3) = 0
end if
end if
! if (CTransScal%NumSources > 2) then
! !Get lensing sources
! !Can modify this here if you want to get power spectra for other tracer
! if (tau>taurend .and. CP%tau0-tau > 0.1_dl) then
! !phi_lens = Phi - 1/2 kappa (a/k)^2 sum_i rho_i pi_i
! !Neglect pi contributions because not accurate at late time anyway
! phi = -(dgrho +3*dgq*adotoa/k)/(k2*EV%Kf(1)*2)
! ! - dgpi/k2/2
! sources(3) = -2*phi*f_K(tau-tau_maxvis)/(f_K(CP%tau0-tau_maxvis)*f_K(CP%tau0-tau))
! ! sources(3) = -2*phi*(tau-tau_maxvis)/((CP%tau0-tau_maxvis)*(CP%tau0-tau))
! !We include the lensing factor of two here
! else
! sources(3) = 0
! end if
! end if
! EFTCAMB MOD END
! EFTCAMB MOD START: computation of mu and gamma. These are not used but it's nice to have the possibility to plot them...
if (CP%EFTflag/=0.and.a>=EFTturnonpi) then
EFT_Psi = EFTsigmadot/k + adotoa*sigma/k
EFT_Phi = etak/k - adotoa*sigma/k
EFT_function_mu = -2._dl*k*(EFTsigmadot+adotoa*sigma)/(dgrho)
EFT_function_gamma = (etak - adotoa*sigma)/(EFTsigmadot + adotoa*sigma)
EFT_xi = vq/(adotoa*clxq)
! This is used to print mu and gamma for different k values.
! If you don't want to print all the other quantities then
! turn off DenugEFTO and turn the other to true.
if (DebugEFTO.or..false.) then
write (9,'(10e15.5)') a, tau, 1._dl/a-1._dl, k, EFT_function_mu, EFT_function_gamma, EFT_xi
end if
end if
! EFTCAMB MOD END
::::: subroutine initial(EV,y, tau) ::::::
Add: After "EV%TightSwitchoffTime = min(tight_tau,Thermo_OpacityToTime(EV%k_buf/ep))"
Add: After "EV%TightSwitchoffTime = min(tight_tau,Thermo_OpacityToTime(EV%k_buf/ep))"
! EFTCAMB MOD START: compute the time at which turn on EFT.
EV%EFTturnOnTime = DeltaTime(0._dl,EFTturnonpi)
! EFTCAMB MOD END
EV%EFTturnOnTime = DeltaTime(0._dl,EFTturnonpi)
! EFTCAMB MOD END
Replace: After "y(EV%g_ix+1)=InitVec(i_qg)"
! EFTCAMB MOD START: initial values.
if (w_lam /= -1 .and. w_Perturb) then
y(EV%w_ix) = InitVec(i_clxq)
y(EV%w_ix+1) = InitVec(i_vq)
end if
!... Original:
! if (w_lam /= -1 .and. w_Perturb) then
! y(EV%w_ix) = InitVec(i_clxq)
! y(EV%w_ix+1) = InitVec(i_vq)
! end if
! EFTCAMB MOD END
if (w_lam /= -1 .and. w_Perturb) then
y(EV%w_ix) = InitVec(i_clxq)
y(EV%w_ix+1) = InitVec(i_vq)
end if
!... Original:
! if (w_lam /= -1 .and. w_Perturb) then
! y(EV%w_ix) = InitVec(i_clxq)
! y(EV%w_ix+1) = InitVec(i_vq)
! end if
! EFTCAMB MOD END
::::: subroutine derivs(EV,n,tau,ay,ayprime) ::::::
Add: After "real(dl) cothxor !1/tau in flat case"
Add: After "real(dl) cothxor !1/tau in flat case"
! EFTCAMB MOD START: definition of variables
! background quantities. Designer approach.
real(dl) EFT_H0,Hdot, Hdotdot
! storage for EFT functions.
real(dl) EFTOmegaV, EFTOmegaP, EFTOmegaPP, EFTOmegaPPP
real(dl) EFTAlpha1V, EFTAlpha1P, EFTAlpha2V, EFTAlpha2P, EFTAlpha3V, EFTAlpha3P
real(dl) EFTAlpha4V, EFTAlpha4P, EFTAlpha4PP, EFTAlpha5V, EFTAlpha5P, EFTAlpha6V, EFTAlpha6P
! Background quantities.
real(dl) EFTc, EFTLambda, EFTcdot, EFTLambdadot
real(dl) EFTgrhoq, EFTgpresq, EFTgrhodotq, EFTgpresdotq
! perturbations quantities.
real(dl) EFTpiA, EFTpiB, EFTpiC, EFTpiD, EFTpiE ! pi field equations
real(dl) EFTeomF, EFTeomN, EFTeomX, EFTeomY, EFTeomG, EFTeomU ! equations of motion
real(dl) EFTeomL, EFTeomM, EFTeomV, EFTeomNdot, EFTeomVdot, EFTeomXdot
real(dl) pidotdot
! Debug quantities
real(dl) temp1, temp2, temp3,temp4, temp5, temp6, temp7
! EFTCAMB MOD END
! background quantities. Designer approach.
real(dl) EFT_H0,Hdot, Hdotdot
! storage for EFT functions.
real(dl) EFTOmegaV, EFTOmegaP, EFTOmegaPP, EFTOmegaPPP
real(dl) EFTAlpha1V, EFTAlpha1P, EFTAlpha2V, EFTAlpha2P, EFTAlpha3V, EFTAlpha3P
real(dl) EFTAlpha4V, EFTAlpha4P, EFTAlpha4PP, EFTAlpha5V, EFTAlpha5P, EFTAlpha6V, EFTAlpha6P
! Background quantities.
real(dl) EFTc, EFTLambda, EFTcdot, EFTLambdadot
real(dl) EFTgrhoq, EFTgpresq, EFTgrhodotq, EFTgpresdotq
! perturbations quantities.
real(dl) EFTpiA, EFTpiB, EFTpiC, EFTpiD, EFTpiE ! pi field equations
real(dl) EFTeomF, EFTeomN, EFTeomX, EFTeomY, EFTeomG, EFTeomU ! equations of motion
real(dl) EFTeomL, EFTeomM, EFTeomV, EFTeomNdot, EFTeomVdot, EFTeomXdot
real(dl) pidotdot
! Debug quantities
real(dl) temp1, temp2, temp3,temp4, temp5, temp6, temp7
! EFTCAMB MOD END
Replace: Before "! Get sound speed and ionisation fraction."
! 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,2)
EFT_H0 = (CP%h0/c_EFT)*1000._dl
end if
!... Original:
! 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,2)
EFT_H0 = (CP%h0/c_EFT)*1000._dl
end if
!... Original:
! if (w_lam==-1._dl) then
! grhov_t=grhov*a2
! else
! grhov_t=grhov*a**(-1-3*w_lam)
! end if
! EFTCAMB MOD END
-1-