*** R.m.s. mass fluctuation, sigma(R) ***
August 23, 2008: E.Komatsu

An r.m.s. mass fluctuation, sigma(R), where R is a top-hat smoothing
scale in units of h^-1 Mpc (NOTE h^-1!) is 

sigma^2(R) = int_0^infinity dlnk k^3 P(k)/(2pi^2) W^2(kR)

where P(k) is the matter power spectrum in units of h^3 Mpc^-3 (NOTE h^3!),
and W(kR) is a window function given by

W(x) = (3/x)*(sin(x)/x^2-cos(x)/x)

Here we provide a simple subroutine, sigma.f90, for computing the r.m.s.
mass fluctuation, sigma(R), as well as its logarithmic derivative,
dln(sigma)/dlnR. R is always measured in units of h^-1 Mpc.

To compute sigma(R), 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
(5a) 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!
(5b) To compute dln(sigma^2)/dlnRh, use a single-precision function 
CHEBEV(lnR1,lnR2,cder,ndim,lnRh), i.e., replace "c" in (4a) with "cder".
(6a) To convert ln(sigma^2) to sigma, use exp(0.5*ln(sigma^2)).
(6b) To convert dln(sigma^2)/dlnRh to dlnsigma/dlnRh, use 0.5*dln(sigma^2)/dlnRh.

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

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

 USE linearpk
 USE sigma
 REAL :: chebev
 REAL :: lnsigma2,lnRh,Rh,dlnsigma2dlnRh
 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)
    print*,'sigma(R) at R=',Rh,' h^-1 Mpc is',exp(0.5*lnsigma2)
    dlnsigma2dlnRh = CHEBEV(lnR1,lnR2,cder,ndim,lnRh)
    print*,'dln(sigma)/dlnR at R=',Rh,' h^-1 Mpc is',0.5*dlnsigma2dlnRh
 endif
 end

Another program, "compute_sigma.f90", for generating sigma(R) and 
dln(sigma)/dlnR for many R (h^-1 Mpc) is

PROGRAM Compute_Sigma
  USE linearpk
  USE sigma
  ! A sample program for computing sigma(Rh) and dln(sigma)/dlnRh
  ! August 23, 2008: E.Komatsu
  IMPLICIT none
  real :: chebev
  real :: lnsigma2,lnRh,Rh,dlnsigma2dlnRh
  character(len=128) :: filename
  integer :: n
! read in and tabulate P(k)
  filename='wmap5baosn_max_likelihood_matterpower.dat'
  n=896 ! # of lines in the file
  CALL open_linearpk(filename,n)
! fit sigma^2(R) to Chebyshev polynomials
  CALL compute_sigma2
 CALL close_linearpk
! now output sigma(R) and dln(sigma)/dlnR as a function of R [h^-1 Mpc]...'
  open(1,file='Rh_sigma.txt')
  open(2,file='Rh_dlnsigmadlnRh.txt')
  do lnRh=lnR1,lnR2,0.01
     lnsigma2 = CHEBEV(lnR1,lnR2,c,ndim,lnRh)
     Rh=exp(lnRh) ! h^-1 Mpc
     write(1,'(2E13.5)') Rh,exp(0.5*lnsigma2)
     dlnsigma2dlnRh = CHEBEV(lnR1,lnR2,cder,ndim,lnRh)
     write(2,'(2E13.5)') Rh,0.5*dlnsigma2dlnRh
  enddo
  close(1)
  close(2)
END PROGRAM Compute_Sigma

An extended version, "compute_sigma2.f90", computes sigma(R) and dln(sigma)/dlnR at 
some redshift. The input power spectrum is given at z=30. Note that dln(sigma)/dlnR is
independent of z, by definition.

PROGRAM Compute_Sigma
  USE cosmo
  USE linearpk
  USE sigma
  USE growth
  ! A sample program for computing sigma(Rh) and dln(sigma)/dlnRh
  ! at some redshift. Note that dln(sigma)/dlnRh is independent of z.
  ! August 24, 2008: E.Komatsu
  IMPLICIT none
  real :: chebev
  real :: lnsigma2,lnRh,Rh,dlnsigma2dlnRh
  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
! 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 output redshift
  print*,'output redshift?'
  read*,zout
! compute the growth factor
  CALL setup_growth
  scalefactor=g(zout)/g(zin)*(1d0+zin)/(1d0+zout)
  print*,'omega matter=',om0
  print*,'omega de=',ode0
  print*,'w=',w
! now output sigma(R) and dln(sigma)/dlnR as a function of R [h^-1 Mpc]...'
  open(1,file='Rh_sigma.txt')
  open(2,file='Rh_dlnsigmadlnRh.txt')
  do lnRh=lnR1,lnR2,0.01
     lnsigma2 = CHEBEV(lnR1,lnR2,c,ndim,lnRh)
     Rh=exp(lnRh) ! h^-1 Mpc
     write(1,'(2E13.5)') Rh,exp(0.5*lnsigma2)*scalefactor
     dlnsigma2dlnRh = CHEBEV(lnR1,lnR2,cder,ndim,lnRh)
     write(2,'(2E13.5)') Rh,0.5*dlnsigma2dlnRh
  enddo
  close(1)
  close(2)
END PROGRAM Compute_Sigma

A program for computing P(k) and dlnP(k)/dlnk at a given wavenumber, 
"sample_linearpk.f90", is

 USE linearpk
 double precision :: linear_pk,dlnpkdlnk,k_ov_h
 character(len=128) :: filename
 integer :: n
 external linear_pk,dlnpkdlnk
 filename='wmap5baosn_max_likelihood_matterpower.dat'
 n=896 ! # of lines in the file
 CALL open_linearpk(filename,n)
 k_ov_h=1d0
 print*,'P(k) at k=',k_ov_h,' h Mpc^-1 is',linear_pk(k_ov_h),' h^3 Mpc^-3'
 print*,'dlnP(k)/dlnk at k=',k_ov_h,' h Mpc^-1 is',dlnpkdlnk(k_ov_h)
 CALL close_linearpk
 end

A program for linearly evolving the input power spectrum to another redshift,
"evolve_linearpk.f90", is

PROGRAM Evolve_LinearPk
 USE cosmo
 USE growth
 ! A program for evolving the linear power spectrum from one redshift to another
 ! August 24, 2008: E.Komatsu
 IMPLICIT none
 double precision :: linear_pk,k_ov_h,g,zin,zout,scalefactor
 character(len=128) :: infile,outfile
 integer :: i,n
 external g
 ! input power spectrum data
 infile='wmap5baosn_max_likelihood_matterpower.dat'
 n=896 ! # of lines in the file
 zin=0d0 ! redshift of the data in the file
! Specify three cosmological parameters
! The data type has been defined in MODULE cosmo.
 ode0=0.723d0
 om0=0.277d0
 w=-1d0
 CALL setup_growth ! tabulate the growth factor
 ! output power spectrum data
 outfile='evolved_linearpk.txt'
 print*,'output redshift?'
 read*,zout
 ! compute the power ratio
 scalefactor=(g(zout)/g(zin)*(1d0+zin)/(1d0+zout))**2d0
 print*,'omega matter=',om0
 print*,'omega de=',ode0
 print*,'w=',w
 print*,'input redshift=',zin
 print*,'output redshift=',zout
 print*,'P(k,out)/P(k,in)=',scalefactor
 open(1,file=infile,status='old')
 open(2,file=outfile,status='unknown')
 do i=1,n
    read(1,*)k_ov_h,linear_pk
    write(2,'(2E16.5)')k_ov_h,linear_pk*scalefactor
 enddo 
 close(1)
 close(2)
END PROGRAM Evolve_LinearPk