*** Non-linear matter density power spectrum from 3rd-order perturbation thoery ***
August 23, 2008: E.Komatsu

Non-linear matter power spectrum computed from the 3rd-order cosmological
perturbation theory is given by

P_dd(k) = P_11(k) + P_22(k) + 2P_13(k)

where P_11(k) is the linear power spectrum, supplied by the user.
Ref. e.g., Eq.(11)-(14) of Jeong & Komatsu, ApJ, 651, 619 (2006)

*** 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.)
 
One can scale this result to an arbitrary redshift, z, by rescaling 
the linear growth factor as

P_dd(k,z) = [D(z)/D(z*)]^2 P_11(k) + [D(z)/D(z*)]^4 [P_22(k) + 2P_13(k)]

where D(z)/D(z*) is the linear growth factor relative to the input linear
power spectrum computed at z=z*.

Here we provide a simple subroutine, pkd.f90, for computing the 3rd-order
perturbation theory spectra, P_11, P_22, and P_13.

To compute the non-linear P(k), it is necessary to use the input linear power spectrum.
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".

(1) Put "USE linearpk" at the beginning of the program
(2) CALL open_linearpk(filename,n) where "filename" contains the name of power spectrum
file (character), and n is the number of lines in the file (integer)
(3) CALL pkd(k_ov_h,p11,p22,p13), where all the arguments are in double precision.
 - k_ov_h: k in units of h Mpc^-1
 - p11: P_11(k) in units of h^-3 Mpc^3
 - p22: P_22(k) in units of h^-3 Mpc^3
 - p13: P_13(k) in units of h^-3 Mpc^3

- To compile and use the sample programs (given below), edit Makefile
and simply "./make"
- It will generate executables called "sample" and "compute_pkd"

=========================================================================
A simple program to use "pkd.f90" (also included as "sample.f90" in
this directory):

 USE linearpk
 double precision :: k_ov_h,p11,p22,p13
 character(len=128) :: filename
 integer :: n
 filename='wmap5baosn_max_likelihood_matterpower.dat'
 n=896 ! # of lines in the file
 CALL open_linearpk(filename,n)
 k_ov_h=0.01d0 ! h Mpc^-1
 CALL pkd(k_ov_h,p11,p22,p13)
 print*,'P_11(k) at k=',k_ov_h,' h Mpc^-1 is',p11,' h^-3 Mpc^3'
 print*,'P_22(k) at k=',k_ov_h,' h Mpc^-1 is',p22,' h^-3 Mpc^3'
 print*,'P_13(k) at k=',k_ov_h,' h Mpc^-1 is',p13,' h^-3 Mpc^3'
 print*,'** Total power spectrum is P_11 + P_22 + 2*P_13. **'
 CALL close_linearpk
 end

A program for generating P(k) as a function of k, for a given redshift,
is provided as "compute_pkd.f90". The input linear spectrum is 
"wmap5baosn_max_likelihood_matterpower_at_z=30.dat".  It looks like this:

PROGRAM Compute_PkD
  USE cosmo
  USE growth
  USE linearpk
  ! A sample program for computing the non-linear power spectrum
  ! of density fluctuations.
  ! - k is in units of h Mpc^-1
  ! - P(k) is in units of h^-3 Mpc^3
  ! August 23, 2008: E.Komatsu
  IMPLICIT none
  integer :: j
  double precision :: g,z,D
  double precision :: ak,p11,p22,p13,dlnk,lnk,zin
  external g
  character(len=128) :: filename
  integer :: n
! Specify three cosmological parameters in double precision
! The data type has been defined in MODULE cosmo.
  ode0=0.723d0
  om0=0.277d0
  w=-1d0
! tabulate g(z) by calling "setup_growth"
  call setup_growth
! read in and tabulate P(k)
  filename='wmap5baosn_max_likelihood_matterpower_at_z=30.dat'
  n=896 ! # of lines in the file
  zin=30d0
  CALL open_linearpk(filename,n)
! ask for redshift
  print*,'redshift?'
  read*,z
  D=g(z)/g(zin)*(1d0+zin)/(1d0+z) ! growth relative to z=zin
! now output P_dd(k,z) as a function of k
  open(1,file='wavenumber_pkd_p11.txt')
  ak=4d-3
  dlnk=1d-2
  do while (ak<0.45d0)
     CALL pkd(ak,p11,p22,p13)
     write(1,'(3E18.8)')ak,D**2d0*p11+D**4d0*(p22+2d0*p13),D**2d0*p11
     lnk=dlog(ak)+dlnk
     ak=dexp(lnk)
  enddo
  dlnk=1d-1
  do while (ak<15d0)
     CALL pkd(ak,p11,p22,p13)
     write(1,'(3E18.8)')ak,D**2d0*p11+D**4d0*(p22+2d0*p13),D**2d0*p11
     lnk=dlog(ak)+dlnk
     ak=dexp(lnk)
  enddo
  close(1)
  CALL close_linearpk
END PROGRAM Compute_PkD