!----------------------------------------------------------------------------- ! ! Halo multiplicity function in double precision, mf(lnnu,z), ! where lnnu=ln(nu) and nu is called the "threshold," ! given by nu = [delta_c/sigma(R,z)]^2. Here, delta_c=1.6865, sigma(R,z) ! is the r.m.s. mass fluctuation within a top-hat smoothing of scale R ! at a redshift z. ! ! The functional form of the mass function is from Tinker et al., ApJ, ! 688, 709 (2008). Note that "mf(lnnu,z)" here corresponds to f(sigma)/2 ! in this reference, where sigma=sqrt(1.6865/nu). ! ! ** This version uses the parameters and the redshift evolution given by ! ** Eq.(4) of Bocquet et al., arXiv:1502.07357 for Delta=200m and "DMOnly" ! ! IMPORTANT: This mass function is NOT normalized, i.e., ! ! int_{-infinity}^{+infinity} dln(nu) mf(lnnu,z) = 1 ! ! is NOT satisfied! Also, we use the fitting function for Delta=200m. ! (That is, the mass is defined within a spherical region whose mean ! density is equal to 200 times the MEAN MASS density of the universe.) ! ! The halo mass function can be computed in the following way. ! ! (1) First, the comoving number density of halos per ln(nu) is related ! to the mass function, dn/dM, as ! ! dM M dn/dM = rho_m dln(nu) mf(lnnu,z) ! ! where rho_m = (average mass density of the universe today) ! ! Dividing both sides by dM, one finds ! ! dn/dlnM = rho_m dln(nu)/dlnM mf(lnnu,z)/M ! ! (2) nu is more directly related to R, as nu=[delta_c/sigma(R,z)]^2. ! So, it is more convenient to write dn/dM as ! ! dn/dlnM = dlnR/dlnM dn/dlnR ! ! and compute dn/dlnR first, via ! ! dn/dlnR = rho_m dln(nu)/dlnR mf(lnnu,z)/M(R) ! ! Now, we can use M(R)=(4pi/3)rho_m R^3 to obtain ! ! dn/dlnR = (3/4pi) dln(nu)/dlnR mf(lnnu,z)/R^3 ! ! with ln(nu) and dln(nu)/dlnR related to lnR via ! - ln(nu) = 2*ln(1.6865)-ln[sigma^2(R,z)] ! - dln(nu)/dlnR = -dln[sigma^2(R)]/dlnR (which is indep. of z) ! ! (3) Once dn/dlnR is obtained as a function of R, one can ! compute dn/dlnM using dlnR/dlnM=1/3, and obtain ! ! dn/dlnM = (1/3) dn/dlnR ! ! with M(R)=(4pi/3)rho_m R^3, ! and rho_m=2.775d11 (omega_matter h^2) M_solar Mpc^-3. ! ! ! << sample code >> ! ! USE linearpk ! USE sigma ! double precision :: deltac=1.6865d0,mf,dndlnRh,dndlnMh ! double precision :: lnnu,dlnnudlnRh,Mh ! 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 get_linearpk(filename,n) ! CALL compute_sigma2 ! Rh=8. ! h^-1 Mpc ! lnRh = alog(Rh) ! if(lnRh>lnR2.or.lnRh