PROGRAM MeanY
  USE cosmo
  USE linearpk
  USE sigma
  USE growth
  USE angular_distance
! This code calculates the mean Compton Y from halos.
! The methodology is described in Klaus, Komatsu & Sunyaev, arXiv:1509.05134
! The output, "redshift_dydln1pz.txt", contains:
! 1st: redshift
! 2nd: dY/dln(1+z)
! July 8, 2015: E.Komatsu
  IMPLICIT none
  double precision :: z1=0.005d0,z2=6d0 ! redshift limits
  character(len=128) :: filename
  integer :: npk,i,j,n=1001 ! 201 is good enough for l>100
  double precision :: Y,dYdlnx,lnx,dlnx,lnx1,lnx2
  external dYdlnx
! Specify cosmological parameters
! The data type has been defined in MODULE cosmo.
  ! Planck 2015 Parameters (with CMB lensing)
  !!$ om0=0.308d0
  !!$ ode0=0.692d0
  !!$ h0=0.6781d0
  !!$ obh2=0.02226d0
  !!$ w=-1d0
  ! WMAP9 Parameters
  om0=0.282d0
  ode0=0.718d0
  h0=0.697d0
  obh2=0.02240d0
  w=-1d0
  ! Magneticum Parameters
  !!$ om0=0.272d0
  !!$ ode0=0.728d0
  !!$ h0=0.704d0
  !!$ obh2=0.0456*h0**2d0
  !!$ w=-1d0
! read in and tabulate P(k)
  ! Planck 2015 Parameters (with CMB lensing)
  !!$ filename='planck15_matterpower.dat'
  !!$ npk=707 ! # of lines in the file
  ! WMAP9 Parameters
  filename='wmap9_matterpower.dat'
  npk=706 ! # of lines in the file
  ! Magneticum Parameters
  !!$ filename='wavenumber_pk_at_z=0_magneticum.txt'
  !!$ npk=896 ! # of lines in the file
  CALL open_linearpk(filename,npk)
! fit sigma^2(R) to Chebyshev polynomials
  CALL compute_sigma2
  CALL close_linearpk
! compute the growth factor
  CALL setup_growth
! compute and tabulate da(z) 
  CALL setup_da
! compute the mean Compton Y
  open(1,file='redshift_dydln1pz_yz.txt')
  lnx1=dlog(1d0+z1)
  lnx2=dlog(1d0+z2)
  dlnx=(lnx2-lnx1)/dble(n-1)
  Y=0d0
  Y=dYdlnx(lnx1)*(0.5d0*dlnx)
  Y=Y+dYdlnx(lnx2)*(0.5d0*dlnx)
  do i=2,n-1
     lnx=lnx1+dble(i-1)*dlnx
     Y=Y+dYdlnx(lnx)*dlnx
     write(1,'(3E15.5)')dexp(lnx)-1d0,dYdlnx(lnx),Y
  enddo
  print*,'total mean Y=',Y
  close(1)
END PROGRAM MeanY
!--------------------------------------------------------------
double precision function dYdlnx(lnx)
  ! dYdlnx = (1+z) (dV/dz) int dlnM dn/dlnM Y
  ! lnx=ln(1+z)
  IMPLICIT none
  double precision, intent(IN) :: lnx
  double precision :: M1=5d5,M2=5d15 ! mass limits, h^-1 Msun
  double precision :: integrand,lnM1,lnM2,lnM,z
  integer :: i
  external integrand
  z=exp(lnx)-1d0
  lnM1=dlog(M1) ! minimum mass, h^-1 Msun
  lnM2=dlog(M2) ! maximum mass, h^-1 Msun
  CALL qgaus1(integrand,lnM1,lnM2,dYdlnx,z) ! Gaussian quadrature
  dYdlnx=dYdlnx*(1d0+z)
  return
end function dYdlnx