[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.