*** Power Spectrum of the Sunyaev-Zel'dovich Effect ***
This version uses Planck 2013's pressure profile with a hydrostatic mass bias
correction applied, and reproduces the result presented in Dolag, Komatsu & 
Sunyaev (2015; DKS).
March 6, 2015: E.Komatsu

Ref: Komatsu & Seljak, MNRAS, 336, 1256 (2002)
     Planck Collaboration, A&A, 550, A131 (2013)
     Dolag, Komatsu & Sunyaev, arXiv:1509.05134

====
Acknowledgment(Nov 29, 2016): The code now uses an efficient integrator for 
oscillating integrals to compute the 2d Fourier transform (in "integrand.f90").
This integration routine was developed by Prof. Takuya Ooura at Kyoto University,
which is available publicly at http://www.kurims.kyoto-u.ac.jp/~ooura/intde.html
EK thanks Prof. Issha Kayo (Tokyo Univ. of Technology) for pointing me to this
routine, and Dr. Ryu Makiya (MPA) for implementing it in this code.
====

This code calculates the power spectrum of the Sunyaev-Zel'dovich (SZ)
effect:

Cl = int dz dV/dz int dlnM dn/dlnM Tsz^2

where dV/dz is the comoving volume element, dn/dlnM is the halo mass 
function, and Tsz is the 2d Fourier transform of the SZ effect *at 150GHz*, 
i.e., Tsz=-0.953652*y*2.725e6 [uK]. To calculate Cl at
different frequencies, you will need to multiply it by [gnu/(-0.953652)]^2, 
where gnu is the appropriate spectrum function of the SZ effect.

** This code is different from the original version used for 
Komatsu & Seljak (2002) in three respects:
1. Mass function: this code uses Bocquet et al. (2015); KS02 used 
Jenkins et al. (2001). 
2. Concentration: this code uses Duffy et al. (2008); KS02 used
Seljak (2000).
3. Pressure profile: this code uses Planck 2013; KS02 used
Komatsu&Seljak (2001).


To compute the SZ power spectrum, it is necessary to use the linear 
matter power spectrum. We provide the sample data, 
"planck15_matterpower.dat," which was generated using CAMB 
code for the "TT+lowP+lensing" parameters given in Table 4 of 
Planck 2015 (arXiv:1502.01589). The input file 
for CAMB is also provided (planck15_params.ini). NOTE THAT 
THIS POWER SPECTRUM IS COMPUTED AT Z=0.

We also provide the matter power spectrum used by Magneticum Pathfinder Simulation
in "wavenumber_pk_at_z=0_magneticum.txt". This power spectrum was generated by
Eisentein&Hu's fitting function (which can be computed by one of the codes in CRL)

<< NOTE ON NFW CONCENTRATION PARAMETER >>

The default concentration parameter is Duffy et al. (2008):
cvir=7.85d0*(mvir/2d12)**(-0.081)/(1d0+z)**0.71

You may also use the concentration parameter of Seljak (2000):
cvir=10d0*(mvir/3.42d12)**(-0.2)/(1d0+z)

which was originally used by Komatsu&Seljak (2002).

<< NOTE ON MASS FUNCTION >>

The default mass function is Eq.(4) of Bocquet et al. (2015) with 
Delta=200m and "Hydro".

Three more mass functions are provided: 
- Sheth&Tormen (mf_shethtormen.f90)
- Jenkins et al. (mf_jenkins.f90; which was originally used by Komatsu&Seljak 2002).
- Bocquet et al. with Delta=200m and "DM Only"
- Tinker et al. (2008)
- Tinker et al. (2010)

To use these mass functions, change the following:

1. In Makefile, change "mf_magneticum.o" to: 
"mf_jenkins.o", 
"mf_shethtormen.o", 
"mf_magneticum_dm.o", 
"mf_tinker_redshift.o", or 
"mf_tinker10.o" 

2. For Jenkins et al. and Sherth&Tormen, change "mf(lnnu,z)" to "mf(lnnu)" 
in integrand.f90 [because Jenkins et al.'s and Sheth&Tormen's functions 
do not have explicit redshift dependence].

Specifically: 
dndlnRh=(3d0/4d0/pi)*dlnnudlnRh*mf(lnnu,z)/Rh**3d0

to 

dndlnRh=(3d0/4d0/pi)*dlnnudlnRh*mf(lnnu)/Rh**3d0

in integrand.f90.

3.  Also, since Jenkins et al. use a different mass definition, you need to 
comment out

  ! Sheth&Tormen and Tinker et al.'s mass functions are given for the 
  ! overdensity mass M200d (with respect to the mean mass density rather than 
  ! the critical density).
  CALL mvir2mdel(mvir,rs,cvir,200d0*omega*rhoc,m200d)
  Rh=(3d0*m200d/4d0/pi/om0/2.775d11)**(1d0/3d0) ! h^-1 Mpc

and use

  ! Alternatively, one may wish to use Jenkins et al.'s mass function,
  ! which is given for the overdensity mass M180d (with respect to the mean 
  ! mass density rather than the critical density).
  CALL mvir2mdel(mvir,rs,cvir,180d0*omega*rhoc,m180d)

in "integrand.f90".



- To compile and use the program, edit Makefile and simply "make"
- It will generate an executable called "szfast"
- Running szfast will generate the data file called "multipole_szpower_150ghz.txt", 
which contains:

1st: multipole
2nd: l(l+1)Cl/twopi [uK^2]

For your convenience, the result from the default options (Duffy et al.'s 
concentration parameter and Bocquet et al.'s mass function) is provided as
"multipole_szpower_150ghz_dks.txt". This is what is presented in Dolag, 
Komatsu & Sunyaev (2015).