*** Sheth&Tormen's Halo Bias ***
December 5, 2010: E.Komatsu

Ref. Eqs[68-70] of Cooray & Sheth, Phys. Rept., 372, 1 (2002) with typos corrected.

The linear bias parameter, b1, is given by 

b1(nu) = 1d0+ep1+E1,

where

ep1 = (q*nu-1d0)/delc,
E1  = 2d0*p/delc/(1d0+(q*nu)**p),

and

q=0.707, p=0.3, delc=1.6865.

Non-linear bias parameters, b2 and b3, are given in the above reference.

Here, the argument "nu" is related to the smoothing radius, R, via

nu = [delc/sigma(R,z)]^2.

Finally, R is related to the mass, M, via

M = (4pi/3)rho_matter R^3,

where

rho_matter = 2.775e11*omega_matter [h^2 Msun/Mpc^3]

To compute sigma(R,z), 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) Put "USE sigma" at the beginning of the program
(3) 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)
(4) CALL compute_sigma2
(5) To compute ln(sigma^2), use a single-precision function CHEBEV(lnR1,lnR2,c,ndim,lnRh),
 where lnR1, lnR2, c and ndim have been pre-computed already by the step (3). 
The only argument you need to enter is the logarithm of the scale at which you want to 
compute ln(sigma^2), lnRh, which is in single precision. Note again that Rh is in units 
of h^-1 Mpc!
(6) To convert ln(sigma^2) to ln(nu), use 2*ln(deltac)-ln(sigma^2), where deltac=1.6865
(7) Compute b1(nu), b2(nu), and b3(nu)

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

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

 USE linearpk
 USE sigma
 double precision :: deltac=1.6865d0,b1,b2,b3
 double precision :: lnnu,Mh
 REAL :: chebev
 REAL :: lnsigma2,lnRh,Rh
 character(len=128) :: filename
 integer :: n
 filename='wmap5baosn_max_likelihood_matterpower.dat'
 n=896 ! # of lines in the file
 CALL open_linearpk(filename,n)
 CALL compute_sigma2
 CALL close_linearpk
 Rh=8. ! h^-1 Mpc
 lnRh = alog(Rh)
 if(lnRh>lnR2.or.lnRh<lnR1)then
    print*,'lnR=',lnRh,' is out of range. It must be between',lnR1,' and',lnR2
    stop
 else
    lnsigma2 = CHEBEV(lnR1,lnR2,c,ndim,lnRh)          ! ln(sigma^2)
    lnnu = 2d0*dlog(deltac)-dble(lnsigma2) ! ln(nu)
    Mh = (4d0*3.1415926535d0/3d0)*2.775d11*Rh**3d0 ! in units of omega_matter h^-1 M_solar
    print*,'b1 at M=',Mh,' omega_matter h^-1 M_solar is ',b1(dexp(lnnu))
    print*,'b2 at M=',Mh,' omega_matter h^-1 M_solar is ',b2(dexp(lnnu))
    print*,'b3 at M=',Mh,' omega_matter h^-1 M_solar is ',b3(dexp(lnnu))
 endif
 end


Another program, "compute_bias.f90", for generating b1, b2, and b3  
for a given mass (h^-1 Msun) and redshift is

PROGRAM Compute_Bias
  USE cosmo
  USE linearpk
  USE sigma
  USE growth
  ! A sample program for computing b1, b2, and b3 using
  ! Sheth&Tormen's formula for a given mass and redshift.
  ! Units: Mh in h^-1 M_solar
  ! December 5, 2010: E.Komatsu
  IMPLICIT none
  double precision :: deltac=1.6865d0,b1,b2,b3
  double precision :: lnnu,Mh
  real :: chebev
  real :: lnsigma2,lnRh,Rh
  double precision :: zin,zout,g,scalefactor
  character(len=128) :: filename
  integer :: n
  external g
! Specify three cosmological parameters
! The data type has been defined in MODULE cosmo.
  ode0=0.723d0
  om0=0.277d0
  w=-1d0
  print*,'omega matter=',om0
  print*,'omega de=',ode0
  print*,'w=',w
! read in and tabulate P(k)
  filename='wmap5baosn_max_likelihood_matterpower_at_z=30.dat'
  n=896 ! # of lines in the file
  CALL open_linearpk(filename,n)
  zin=30d0
! fit sigma^2(R) to Chebyshev polynomials
  CALL compute_sigma2
  CALL close_linearpk
! ask for the mass and redshift
  print*,'mass [h^-1 Msun] and redshift?'
  read*,Mh,zout
! compute the growth factor
  CALL setup_growth
  scalefactor=g(zout)/g(zin)*(1d0+zin)/(1d0+zout)
  Rh=((3./4./3.14159)*Mh/(2.775d11*om0))**(1./3.) 
  lnsigma2 = CHEBEV(lnR1,lnR2,c,ndim,alog(Rh)) ! ln(sigma^2)
  lnnu = 2d0*dlog(deltac)-dble(lnsigma2)-2d0*dlog(scalefactor) ! ln(nu)
  print*,'b1,b2,b3=',b1(dexp(lnnu)),b2(dexp(lnnu)),b3(dexp(lnnu))
END PROGRAM Compute_Bias