[A] Changes in module mgvariables in equations_ppf.f90
Add the following new variables,
real(dl) :: fR0, fRn ! the |f_{R0}| and n parameters in the Hu-Sawicki model
[B1] Changes in subroutine derivs in equations_ppf.f90
[1] Add the following just after ‘use mgvariables’, which is at the very beginning of this subroutine.
use constants, only : c
[2] In the parameter definition part (also at the beginning of this subroutine), add the following,
real(dl) mubar, Rbar, R0bar, Rdot, mudot
[3] Just after ‘if (tempmodel /= 0) then’, change the next line to,
if (model==1 .or.model==4 .or.model==5.or.model==6 .or. model==8) then ! model 8 is the Hu-Sawicki model (model 7 is reserved for DGP)
[4] Just after model 6 block [if (model==6) then ...... end if], add the following piece of code for Hu-Sawicki,
if (model==8) then ! Hu-Sawicki
omm=CP%omegab+CP%omegac
Rbar=3.d0*((CP%H0/c*1.d3)**2)*(omm*(a**(-3))+4.d0*(1.d0-omm))
R0bar=3.d0*((CP%H0/c*1.d3)**2)*(omm+4.d0*(1.d0-omm))
mubar=1.d0/3.d0/(1.d0+fRn)
mubar=mubar*Rbar/fR0
mubar=mubar*((Rbar/R0bar)**(1.d0+fRn))
mubar=sqrt(mubar)
MG_mu = 4.d0/3.d0-1.d0/3.d0*(mubar*a)**2/(k2+(mubar*a)**2)
Rdot=-9.d0*((CP%H0/c*1.d3)**2)*omm*adotoa*a**(-3)
mudot=1.d0/3.d0/(1.d0+fRn)/fR0/((R0bar)**(1.d0+fRn))
mudot=sqrt(mudot)
mudot=mudot*(1.d0+fRn/2.d0)*(Rbar**(fRn/2.d0))*Rdot
MG_mudot = -2.d0/3.d0*k2/((k2+(mubar*a)**2)**2)
MG_mudot = MG_mudot*(mubar*mudot*(a**2)+(mubar*a)**2*adotoa)
MG_gamma = 1.d0-2.d0*k2/(3.d0*(mubar*a)**2+4.d0*k2)
MG_gammadot = 12.d0*k2*(mubar*a)**2
MG_gammadot = MG_gammadot*(adotoa+mudot/mubar)
MG_gammadot = MG_gammadot/(3.d0*(mubar*a)**2+4.d0*k2)**2
end if
[B2] Do the same changes in subroutine output in equations_ppf.f90
[C] Changes in inidriver.F90
Add the Hu-Sawicki model just after the model 6 block [if (model==6) then ...... end if]
else if (model ==8) then ! Hu-Sawicki
fR0 = Ini_Read_Double('|fR0|')
fRn = Ini_Read_Double('fRn')
write(*,*) '|fR0|=',fR0, 'fRn=',fRn
[D] Add the following stuff to params.ini
Add the Hu-Sawicki model
#model= 8 : Hu-Sawicki
and the related parameters
|fR0|=1E-4
fRn=1.0
[E] Recompile the code and run the code using the following settings (the Planck parameters) in params.ini
model = 8
GRtrans= 0.1
use_physical = T
ombh2 = 0.022161
omch2 = 0.11889
omnuh2 = 0
omk = 0
hubble = 67.77
initial_power_num = 1
pivot_scalar = 0.05
pivot_tensor = 0.05
scalar_amp(1) = 2.21381e-9
scalar_spectral_index(1) = 0.9611
scalar_nrun(1) = 0
transfer_high_precision = T
transfer_kmax = 20
transfer_k_per_logint = 0
transfer_num_redshifts = 1
transfer_interp_matterpower = T
transfer_redshift(1) = 0
If everything works, you should be able to reproduce my figure at the top of this page. Note that it is slightly different from Fig 4 in the Hu-Sawicki paper
because different cosmological parameters used.