[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 


         MG_mu = 4.d0/3.d0-1.d0/3.d0*(mubar*a)**2/(k2+(mubar*a)**2)

         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




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

Add the Hu-Sawicki model to MGCAMB