*** Non-linear Matter Spectra from Peacock&Dodds (1996); Ma et al. (1999); and Halofit [Smith et al. (2003)] ***
January 1, 2012: E.Komatsu

This program computes non-linear matter power spectra, P(k,z), as a function of wavenumbers k [in units of h/Mpc] and redshifts z, using 3 different empirical prescriptions in the literature:

1. Peacock and Dodds, MNRAS, 280, L19 (1996), which uses the scaling argument developed by HKLM, calibrated to N-body simulations.
2. Ma et al., ApJL, 521, L1 (1999), which extends Peacock&Dodds by including effects of dynamical dark energy models.
3. Halofit [Smith et al., MNRAS, 341, 1311 (2003)], which uses a halo model calibrated to N-body simulations.

*** P(k) is in units of h^-3 Mpc^3, and k is in units of h Mpc^-1 ***
(This is the same as CAMB definition.)

In order to calculate non-linear power spectra, we need the input linear power spectrum. For this purpose we provide the sample data, "wmap5baosn_max_likelihood_matterpower.dat," which was generated using CAMB code for the maximum likelihood parameters given in Table I of Komatsu et al.(2008) [WMAP 5-year interpretation paper] with "WMAP5+BAO+SN". The input file for CAMB is also provided (wmap5baosn_max_likelihood_params.ini). NOTE THAT THIS POWER SPECTRUM IS COMPUTED AT Z=0.

Another power spectrum, evolved back to z=30, is provided as "wmap5baosn_max_likelihood_matterpower_at_z=30.dat".

In addition, we provide the linear power spectrum WITHOUT BARYONIC OSCILLATION, computed from Eisenstein&Hu's smooth transfer function. This was computed by "compute_pk_nowiggle," which is available also on CRL. The data are called "wavenumber_pknowiggle.txt" for z=0 and "wavenumber_pknowiggle_at_z=30.txt" for z=30.

Finally, we provide the non-linear power spectrum computed from the 3rd-order perturbation theory (3PT) for z=0, which can be compared with the outputs of this program. It is called "wavenumber_pkd_p11.txt" and the format is: (1st) wavenumber [h/Mpc]; (2nd) non-linear P(k) [Mpc^3/h^3]; (3rd) linear P(k )[Mpc^3/h^3]. These data were generated by "compute_pkd," which is available also on CRL.

- To compile and use the program, edit Makefile and simply "make"
- It will generate an executable called "compare_pk"
- Running it will generate 3 output files:
  - wavenumber_pkpd96.txt: Peacock&Dodds (1996)
  - wavenumber_pkma99.txt: Ma et al. (1999)
  - wavenumber_pkhalofit.txt: Halofit [Smith et al. (2003)]

The format is
(1st) wavenumber [h Mpc^-1]
(2nd) non-linear P(k) [h^-3 Mpc^3]

******************************************************************
Note on halofit.f [available on http://www.roe.ac.uk/~jap/haloes/]
******************************************************************
You can reproduce the results from this program (with the nowiggle linear power spectrum) using the original routine, "halofit.f," provided by Smith et al. You need to modify "function p_cdm" as follows. Note that the cosmological parameters are hard-coded here, and thus the comparison can be made only for the WMAP 5-year WMAP+BAO+SN best-fit parameters.

With this modification, the input parameters must be
 input z=0 matter density (om_m0)
0.277
 input z=0 vacuum density (om_v0)
0.723

"sig8" and "gams" do not play any role, and thus you can put in any numbers.

c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c Bond & Efstathiou (1984) approximation to the linear CDM power spectrum 

      function p_cdm(rk,gams,sig8)
      implicit none 
      real*8 p_cdm,rk,gams,sig8,p_index,rkeff,q,q8,tk,tk8
      real*8 alpha,sound,shape,aq,omegamh2,fb,om0,h0
      real*8 ak,deltaR2
      deltaR2=2.46d-9
      h0=0.702d0
      om0=0.277d0
      omegamh2=om0*h0**2
      fb=0.0459d0/om0
      p_index=0.962d0
      ak=rk*h0 ! 1/Mpc. "rk" is in units of h/Mpc.
      alpha=1d0-0.328d0*dlog(431d0*omegamh2)*fb
     &+0.38d0*dlog(22.3d0*omegamh2)*fb**2d0 
      sound=44.5d0*dlog(9.83d0/omegamh2)
     &/dsqrt(1d0+10d0*(fb*omegamh2)**(3d0/4d0))
      shape=omegamh2*(alpha+(1d0-alpha)/(1d0+(0.43d0*ak*sound)**4d0)) 
      aq=ak*(2.725d0/2.7d0)**2d0/shape
      tk=dlog(2d0*dexp(1d0)+1.8d0*aq)/(dlog(2d0*dexp(1d0)+1.8d0*aq)
     &+(14.2d0+731d0/(1d0+62.5d0*aq))*aq*aq)
      p_cdm=deltaR2*(2d0*ak**2d0*2998d0**2d0/h0**2d0/5d0/om0)**2d0
     &*tk*tk*(ak/0.002d0)**(p_index-1d0)
     &*0.7646**2d0 ! value of D(z=0) normalized by (1+z)D(z)->1 during MD
c      p_index=1.
c      rkeff=0.172+0.011*log(gams/0.36)*log(gams/0.36)
c      q=1.e-20 + rk/gams
c      q8=1.e-20 + rkeff/gams
c      tk=1/(1+(6.4*q+(3.0*q)**1.5+(1.7*q)**2)**1.13)**(1/1.13)
c      tk8=1/(1+(6.4*q8+(3.0*q8)**1.5+(1.7*q8)**2)**1.13)**(1/1.13)
c      p_cdm=sig8*sig8*((q/q8)**(3.+p_index))*tk*tk/tk8/tk8
      return
      end