*** Mean Compton Y with Planck 2013 Pressure Profile (with HSE bias) ***
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).
July 8, 2015: E.Komatsu

Ref: Barbosa, Bartlett, Blanchard & Oukbir, A&A, 314, 13 (1996)
     Planck Collaboration, A&A, 550, A131 (2013)
     Dolag, Komatsu & Sunyaev, arXiv:1509.05134

This code calculates the mean Compton Y:

<Y> = int dz dV/dz int dlnM dn/dlnM Y

where dV/dz is the comoving volume element, dn/dlnM is the halo mass 
function, and Y is the monopole of the 2d Fourier transform of the Compton Y profile.

To compute the mean Y, 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).
This power spectrum and the associated cosmological parameters need to be chosen
in "meany.f90" to reproduce the dashed line in Fig.7 of DKS. Finally, we also provide
the matter power spectrum of WMAP9 in "wmap9_matterpower.dat". (The CAMB parameter 
file is "wmap9_params.ini".) This parameter was used by Hill et al., arXiv:1507.01583

<< 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.(8)-(12) of Tinker et al. (2010) for Delta=200.

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 "Hydro"
- Bocquet et al. with Delta=200m and "DM Only"
- Tinker et al. (2008)

To use these mass functions, change the following:

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

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 "meany"
- Running meany will generate the data file called "redshift_dydln1pz.txt", 
which contains:

1st: redshift
2nd: dY/dln(1+z)

and also return a value of the mean Y on the command line.

For your convenience, the result from the default options (Duffy et al.'s 
concentration parameter and Tinker et al. (2010)'s mass function) is provided as
"redshift_dydln1pz_dks.txt". The mean Y from this configuration is <Y>=1.31928e-6