*** Non-linear galaxy power spectrum from 3rd-order perturbation thoery ***
November 5, 2008: E.Komatsu

Galaxy power spectrum with non-linear bias, computed using Patrick McDonald's
re-parametrized bias parameters. For a given input wavenumber, k, in units of
h Mpc^-1, This subroutine will return P_b2(k) and P_b22(k), which appear in
the galaxy power spectrum as

P_galaxy(k) = P0 + b1^2*[ P_dd(k) + b2*P_b2(k) + b2^2*P_b22(k) ]

where P_dd(k) is the non-linear matter density power spectrum computed
from the 3rd-order perturbation theory, and P_b2(k) and P_b22(k) are
the non-linear bias terms computed from the linear matter power spectrum,
supplied by the user. 

The bias parameters are:
- P0: constant term, in units of h^-3 Mpc^3
- b1: linear bias parameter
- b2: non-linear bias parameter

Ref. McDonald, PRD, 74, 103512
     [See Eq.(10)&(16) of Erratum, PRD, 74, 129901(E)]
     See also Eq.(2)&(3) of Jeong & Komatsu, arXiv:0805.2632 (ApJ in press)

*** 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.)

Both P_b2(k) and P_b22(k) will grow in time as [D(z)]^4 where D(z) is the 
linear growth factor, i.e.,

P_b2(k,z) = [D(z)/D(z*)]^4 P_b2(k)
P_b22(k,z)= [D(z)/D(z*)]^4 P_b22(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, nonlinearbias.f90, for computing 
P_b2(k) and P_b22(k) for a given k.

To compute P_b2&P_b22, 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".

We also provide the data for the non-linear power spectrum,
"wavenumber_pkd11_pkd22_pkd13_at_z=30.txt", computed from
the 3rd-order perturbation theory.

(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 nonlinearbias(k_ov_h,pb2,pb22), where all the arguments are in double precision.
 - k_ov_h: k in units of h Mpc^-1
 - pb2:  P_b2(k)  in units of h^-3 Mpc^3
 - pb22: P_b22(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_pkgalaxy"

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

 USE linearpk
 double precision :: k_ov_h,pb2,pb22
 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 nonlinearbias(k_ov_h,pb2,pb22)
 print*,'P_b2(k) at k=',k_ov_h,' h Mpc^-1 is',pb2,' h^-3 Mpc^3'
 print*,'P_b22(k) at k=',k_ov_h,' h Mpc^-1 is',pb22,' h^-3 Mpc^3'
 print*,'** Total power spectrum is b1^2*(P_dd + b2*P_b22 + b2^2*P_b22) **'
 CALL close_linearpk
 end

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

PROGRAM Compute_PkGalaxy
  USE cosmo
  USE growth
  USE linearpk
  USE nonlinearpk
  ! A sample program for computing the non-linear power spectrum
  ! of galaxies
  ! The output file contains
  ! - 1st column: k in units of h Mpc^-1
  ! - 2nd column: non-linear matter density P(k)
  ! - 3rd column: non-linear bias term 1: Pb2(k)
  ! - 4th column: non-linear bias term 2: Pb22(k)
  ! - 5th column: non-linear galaxy P(k)
  ! November 3, 2008: E.Komatsu
  IMPLICIT none
  integer :: j
  double precision :: g,z,D,P0,b1,b2
  double precision :: ak,nonlinear_pkd,pb2,pb22,dlnk,lnk,zin
  external g
  character(len=128) :: filename
  integer :: n
  external nonlinear_pkd
! Specify three cosmological parameters
! 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 linear 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
  ! read in and tabulate non-linear P(k) density
  filename='wavenumber_pkd11_pkd22_pkd13_at_z=30.txt'
  n=508    ! # of lines in the file
  zin=30d0 ! redshift of the input power spectrum
  CALL open_nonlinearpkd(filename,n,D)
! ask for three bias parameters
  print*,'P0[h^-3 Mpc^3],b1,b2?'
  read*,P0,b1,b2
! now output P_dd(k,z),p_b2(k,z),P_b22(k,z),and P_galaxy(k,z) as a function of k
  open(1,file='wavenumber_pkd_pb2_pb22_pkgalaxy.txt')
  ak=4d-3
  dlnk=1d-2
  do while (ak<0.45d0)
     CALL nonlinearbias(ak,pb2,pb22)
     write(1,'(5E18.8)')ak,nonlinear_pkd(ak),&
          D**4d0*pb2,D**4d0*pb22,&
          P0 + b1**2d0*( nonlinear_pkd(ak)+b2*D**4d0*pb2+b2**2d0*D**4d0*pb22 )
     lnk=dlog(ak)+dlnk
     ak=dexp(lnk)
  enddo
  dlnk=1d-1
  do while (ak<15d0)
     CALL nonlinearbias(ak,pb2,pb22)
     write(1,'(5E18.8)')ak,nonlinear_pkd(ak),&
          D**4d0*pb2,D**4d0*pb22,&
          P0 + b1**2d0*( nonlinear_pkd(ak)+b2*D**4d0*pb2+b2**2d0*D**4d0*pb22 )
     lnk=dlog(ak)+dlnk
     ak=dexp(lnk)
  enddo
  close(1)
  CALL close_linearpk
  CALL close_nonlinearpkd
END PROGRAM Compute_PkGalaxy