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