Source code for rabacus.cosmology.mass_function.mass_function

""" 
A halo mass function module.  See the following reference for a discussion, 
http://adsabs.harvard.edu/abs/2007ApJ...671.1160L
"""

import sys

import numpy as np
#import matplotlib.pyplot as plt
#import matplotlib as mpl



#import ..cosmology
#import bbks86
#from ....constants import physical
#import idle_hands.constants.physical_constants_units as pcu

import scipy.integrate
import scipy.interpolate

__all__ = ['MassFunction']




[docs]class MassFunction: """ A mass function class. During initialization the normalization of the power spectrum is set to match the sigma8 from cosmo. Args: `cosmo` (:class:`~rabacus.cosmology.general.Cosmology`): an instance of the cosmology class. `tf` (class): an instance of a transfer function class. For example, :class:`~rabacus.cosmology.transfer_functions.bbks86.TransferBBKS`. """ def __init__(self, cosmo, tf): """ cosmo is an instance of Cosmology and tf is an instance of TransferFunction """ self.cosmo = cosmo self.tf = tf self.cu = cosmo.cu self.power_norm = 1.0 sig2 = self.sigma2_R( 8.0 * self.cu.Mpc / self.cu.h, z=0.0 ) self.power_norm = self.cosmo.sigma8**2 / sig2 self.delta_c0 = 1.68647 self.map_sig2_R()
[docs] def calc_mf( self, z, fit='Warren06' ): """ Calculate mass function from high to low mass. The only redshift dependence is in f(sigma) via a rescaling of sigma """ N = self.Mw_arr.size self.dndM = np.zeros( N ) * (self.cu.h / self.cu.Mpc)**3 self.dndlogM = np.zeros( N ) * (self.cu.h / self.cu.Mpc)**3 rho_mean = self.cosmo.rho_crit0 * self.cosmo.OmegaM delta_c = self.delta_c0 / self.D1( z ) for i in reversed(np.arange(N)): # go from high mass to low mass sigma = self.sig_arr[i] mass = self.Mw_arr[i] logM = np.log10( np.array( mass ) ) f = self.mult_func( sigma, z, fit=fit ) dlnisigdM = self.lnisig_of_Mw.derivatives(mass)[1] dlnisigdlogM = self.lnisig_of_logMw.derivatives(logM)[1] # http://adsabs.harvard.edu/abs/2008ApJ...688..709T Eq. 2 dndM = (rho_mean / mass) * f * dlnisigdM dndlogM = (rho_mean / mass) * f * dlnisigdlogM #print 'M,dndM,dndlogM: ', mass, dndM, dndlogM self.dndM[i] = dndM self.dndlogM[i] = dndlogM print 'attribute dndM now contains dn/dM at z = ', str(z) print 'attribute dndlogM now contains dn/dlogM at z = ', str(z)
[docs] def mult_func( self, sigma_in, z, fit='Warren06' ): r""" The multiplicity function :math:`f(\sigma)`. This function determines the shape of the mass function given the variation of :math:`\sigma` with scale. The variable `fit` determines the form of the multiplicity function. A convenient variable is :math:`\nu = \delta_{c,0} / \sigma`. The following values for `fit` lead to the following multiplicity funcitons, - ``PS``: `Press & Schechter 74 <http://adsabs.harvard.edu/abs/1974ApJ...187..425P>`_ .. math:: f = \sqrt{ \frac{2}{\pi} } \, \nu \exp( -\nu^2/2 ) - ``ST``: `Sheth & Tormen 99 <http://adsabs.harvard.edu/abs/1999MNRAS.308..119S>`_ .. math:: f = A \sqrt{ \frac{2a}{\pi} } [1 + (a \nu^2)^{-p}] \nu \exp(-a \nu^2/2) A = 0.3222, \, a = 0.75, \, p = 0.3 - ``Jenkins01``: `Jenkins 01 <http://adsabs.harvard.edu/abs/2001MNRAS.321..372J>`_ .. math:: f = 0.315 \exp( -| \ln \sigma^{-1} + 0.61 |^{3.8} ) - ``Warren06``: `Warren 06 <http://adsabs.harvard.edu/abs/2006ApJ...646..881W>`_ .. math:: f = A ( \sigma^{-a} + b ) \exp(-c/\sigma^2) A=0.7234, \, a=1.625, \, b=0.2538, \, c=1.1982 - ``Tinker08``: `Tinker 08 <http://adsabs.harvard.edu/abs/2008ApJ...688..709T>`_ .. math:: f = A \left[ \left( \frac{\sigma}{b} \right)^{-a} + 1 \right] \exp(-c/\sigma^2) \Delta=300 A = 0.1 \, \log \Delta - 0.05 a = 1.43 + ( \log \Delta - 2.3 )^{1.5} b = 1.00 + ( \log \Delta - 1.6 )^{-1.5} c = 1.20 + ( \log \Delta - 2.35 )^{1.6} """ sigma = sigma_in * self.D1(z) nu = self.delta_c0 / sigma if fit == 'PS': # http://adsabs.harvard.edu/abs/2007ApJ...671.1160L Eq. 9 f = np.sqrt(2/np.pi) * nu * np.exp(-0.5*nu*nu) elif fit == 'ST': # http://adsabs.harvard.edu/abs/2007ApJ...671.1160L Eq. 11 A = 0.3222 a = 0.75 # a = 0.707 p = 0.3 t1 = A * np.sqrt( 2*a/np.pi ) t2 = 1 + (a*nu*nu)**(-p) t3 = nu * np.exp( -0.5*a*nu*nu ) f = t1 * t2 * t3 elif fit == "Jenkins01": # http://adsabs.harvard.edu/abs/2007ApJ...671.1160L Eq. 12 t1 = np.abs( np.log(1.0/sigma) + 0.61 ) f = 0.315 * np.exp( -t1**(3.8) ) elif fit == "Warren06": # http://adsabs.harvard.edu/abs/2007ApJ...671.1160L Eq. 13 A=0.7234 a=1.625 b=0.2538 c=1.1982 t1 = sigma**(-a) + b f = A * t1 * np.exp(1.0)**( -c / (sigma*sigma) ) elif fit == "Tinker08": # http://adsabs.harvard.edu/abs/2008ApJ...688..709T Eq. 3 Delta = 300.0 logD = np.log10( Delta ) if Delta < 1600.0: A = 0.1 * logD - 0.05 else: A = 0.26 a = 1.43 + ( logD - 2.3 )**(1.5) b = 1.00 + ( logD - 1.6 )**(-1.5) c = 1.20 + ( logD - 2.35 )**(1.6) t1 = (sigma/b)**(-a) + 1 f = A * t1 * np.exp(1.0)**( -c / (sigma*sigma) ) else: print 'fit not recognized' sys.exit(1) return f
[docs] def map_sig2_R( self, nbins=200 ): """ Map out the relationship between sigma^2 and R. In this routine we map out the relationship at z=0 and assume that different redshifts can be accomodated through a simple scaling with the growth function D1(z). """ z = 0.0 Rw_min = 8.0e-4 * self.cu.Mpc / self.cu.h Rw_max = 8.0e1 * self.cu.Mpc / self.cu.h logRw_min = np.log10( np.array(Rw_min) ) logRw_max = np.log10( np.array(Rw_max) ) rho_mean = self.cosmo.rho_crit0 * self.cosmo.OmegaM V_min = (4.0/3.0) * np.pi * Rw_min**3 V_max = (4.0/3.0) * np.pi * Rw_max**3 M_min = rho_mean * V_min M_max = rho_mean * V_max print 'min/max Rw correspond to ... ' print 'Rw min/max: ', Rw_min, Rw_max print 'M min/max: ', M_min, M_max print 'calculating relationship between R and sigma^2' print 'this may take a minute or so ... ' N = nbins dlogRw = ( logRw_max - logRw_min ) / N self.logRw_arr = np.linspace( logRw_min, logRw_max, N ) self.Rw_arr = 10**self.logRw_arr * self.cu.Mpc / self.cu.h self.Vw_arr = (4.0/3.0) * np.pi * self.Rw_arr**3 self.Mw_arr = rho_mean * self.Vw_arr self.logVw_arr = np.log10( np.array( self.Vw_arr ) ) self.logMw_arr = np.log10( np.array( self.Mw_arr ) ) self.logsig_arr = np.zeros( N ) self.sig_arr = np.zeros( N ) self.isig_arr = np.zeros( N ) self.lnisig_arr = np.zeros( N ) for i in range(N): Rw = self.Rw_arr[i] sig2 = self.sigma2_R( Rw, z ) sig = np.sqrt(sig2) self.logsig_arr[i] = np.log10( sig ) self.sig_arr[i] = sig #print 'i,Rw,sig2: ', i,Rw,sig self.isig_arr = 1.0 / self.sig_arr self.lnisig_arr = np.log( self.isig_arr ) # create spline representations #---------------------------------------------------------- fit = scipy.interpolate.UnivariateSpline( self.logsig_arr, self.logRw_arr, s=0) self.logRw_of_logsig = fit fit = scipy.interpolate.UnivariateSpline( self.logRw_arr, self.logsig_arr, s=0) self.logsig_of_logRw = fit fit = scipy.interpolate.UnivariateSpline( self.Mw_arr, self.lnisig_arr, s=0) self.lnisig_of_Mw = fit fit = scipy.interpolate.UnivariateSpline( self.logMw_arr, self.lnisig_arr, s=0) self.lnisig_of_logMw = fit fit = scipy.interpolate.UnivariateSpline( self.Mw_arr, self.sig_arr, s=0) self.sig_of_Mw = fit
[docs] def delta_cz( self, z ): """ Redshift dependent critical collapse value Args: `z` (real): redshift .. math:: \delta_c(z) = \delta_{c,0} / D_1(z) """ dz = self.delta_c0 / self.D1( z ) return dz
[docs] def Pk(self, k, z): """ Power spectrum from `tf` normalized to match sigma8 from `cosmo`. """ if hasattr(k,'units'): k.units = 'hh/Mpc' else: raise ValueError, '\n Input variable k must have units \n' Pk = self.tf.Pk_cdm(k) return Pk * self.power_norm * self.D1(z)**2
[docs] def Delta2k(self, k, z): r""" Dimensionless power spectrum, Args: `k` (real or array): wavenumber. `z` (real): redshift .. math:: \Delta^2(k,z) = \frac{k^3}{2 \pi^2} P(k,z) """ if hasattr(k,'units'): k.units = 'hh/Mpc' else: raise ValueError, '\n Input variable k must have units \n' kk = np.array(k) Dk = kk * kk * kk / (2.0 * np.pi * np.pi) * self.Pk(k,z) return Dk
[docs] def D1(self,z): """ Linear growth function from `cosmo` """ D = self.cosmo.D1z( z ) return D
[docs] def WkR(self, k, R): r""" The fourier transform of a real space spherical top hat filter. Args: `k` (real or array): wavenumber. `R` (real): filter scale. .. math:: W = \frac{ 3 j_1(x) }{ x }, \, x = k R \\ j_1 = [\sin(x) - x \cos(x)] \, x^{-2} """ if hasattr(k,'units'): k.units = 'hh/Mpc' else: raise ValueError, '\n Input variable k must have units \n' if hasattr(R,'units'): R.units = 'Mpc/hh' else: raise ValueError, '\n Input variable R must have units \n' y = k * R yi = 1.0 / y yi2 = yi * yi yi3 = yi2 * yi t1 = np.sin(y) * yi3 t2 = np.cos(y) * yi2 W = 3.0 * (t1-t2) return W
[docs] def W2kR(self, k,R): r""" The square of the fourier transform of a real space spherical top hat filter. Args: `k` (real or array): wavenumber. `R` (real): filter scale. .. math:: W^2 = \left[ \frac{ 3 j_1(x) }{ x } \right]^2, \, x = k R \\ j_1 = [\sin(x) - x \cos(x)] \, x^{-2} """ if hasattr(k,'units'): k.units = 'hh/Mpc' else: raise ValueError, '\n Input variable k must have units \n' if hasattr(R,'units'): R.units = 'Mpc/hh' else: raise ValueError, '\n Input variable R must have units \n' W = self.WkR(k,R) W2 = W * W return W2
[docs] def W2lnkR(self, lnk, R): r""" The square of the fourier transform of a real space spherical top hat filter as a function of the natural log of `k`. Args: `lnk` (real or array): natural log of wavenumber, ln( k [h/Mpc] ). `R` (real): filter scale. """ if hasattr(R,'units'): R.units = 'Mpc/hh' else: raise ValueError, '\n Input variable R must have units \n' kk = np.exp(1.0)**lnk / self.cu.Mpc / self.cu.h W2 = self.W2kR( kk, R ) return W2
[docs] def dsig2_dk( self, k, R, z ): r""" Integrand for calculation of sigma^2. Inputs k and R must have units. Args: `k` (real or array): wavenumber. `R` (real): filter scale. `z` (real): redshift .. math:: \frac{d\sigma^2}{dk} = \frac{ \Delta^2(k,z) W^2(k,R) }{ k } """ if hasattr(k,'units'): k.units = 'hh/Mpc' else: raise ValueError, '\n Input variable k must have units \n' if hasattr(R,'units'): R.units = 'Mpc/hh' else: raise ValueError, '\n Input variable R must have units \n' kk = np.array(k) dsdk = self.Delta2k(k,z) * self.W2kR(k,R) / kk return dsdk
[docs] def dsig2_dlnk( self, lnk, R, z ): r""" Input is ln(k) which is unitless but k must have units of h/Mpc Args: `lnk` (real or array): wavenumber, ln( k [h/Mpc] ). `R` (real): filter scale. `z` (real): redshift .. math:: \frac{d\sigma^2}{d \ln k} = \Delta^2(k,z) W^2(k,R) """ if hasattr(R,'units'): R.units = 'Mpc/hh' else: raise ValueError, '\n Input variable R must have units \n' kk = np.exp(1.0)**lnk / self.cu.Mpc / self.cu.h dsdlnk = self.Delta2k(kk,z) * self.W2kR(kk,R) return dsdlnk
[docs] def dsig2_dlogk( self, logk, R, z ): r""" Input is log10(k) which is unitless but k must have units of h/Mpc Args: `logk` (real or array): wavenumber, log10( k [h/Mpc] ) `R` (real): filter scale. `z` (real): redshift .. math:: \frac{d\sigma^2}{d \log k} = \ln(10) \Delta^2(k,z) W^2(k,R) """ if hasattr(R,'units'): R.units = 'Mpc/hh' else: raise ValueError, '\n Input variable R must have units \n' kk = 10**logk / self.cu.Mpc / self.cu.h dsdlogk = np.log(10.) * self.Delta2k(kk,z) * self.W2kR(kk,R) return dsdlogk
[docs] def sigma2_R(self, R, z, method='romberg'): r""" Variance of density field smoothed on scale R Args: `R` (real): filter scale. `z` (real): redshift .. math:: \sigma^2(R) = \int \Delta^2(k,z) W^2(k,R) d \ln k """ if hasattr(R,'units'): R.units = 'Mpc/hh' else: raise ValueError, '\n Input variable R must have units \n' kw = np.array( 1.0/R ) lnkw = np.log( kw ) lnk_min = lnkw - 6.0 lnk_max = lnkw + 6.0 if method == 'romberg': sig2 = scipy.integrate.romberg( self.dsig2_dlnk, lnk_min, lnk_max, args=(R,z), tol=1.0e-12, # rtol=1.0e-12, divmax=20, vec_func=True ) elif method == 'quad': sig2 = scipy.integrate.quad( self.dsig2_dlnk, lnk_min, lnk_max, args=(R,z), epsabs=1.0e-8, epsrel=1.0e-8 )[0] return sig2
[docs] def sigma2_V(self, V, z): r""" Variance of density field smoothed on scale V Args: `V` (real): filter scale. `z` (real): redshift .. math:: R = \left( \frac{3 V}{4 \pi} \right)^{1/3}, \\ \sigma^2(R) = \int \Delta^2(k,z) W^2(k,R) d \ln k """ if hasattr(V,'units'): V.units = 'Mpc^3/hh^3' else: raise ValueError, '\n Input variable V must have units \n' R = ( (3 * V) / (4 * np.pi) )**(1./3) total = self.sigma2_R(R,z) return total
[docs] def sigma2_M(self, M, z): r""" Variance of density field smoothed on scale V Args: `V` (real): filter scale. `z` (real): redshift .. math:: V = \frac{M}{\rho_{c,0} \, \Omega_M}, \, R = \left( \frac{3 V}{4 \pi} \right)^{1/3}, \\ \sigma^2(R) = \int \Delta^2(k,z) W^2(k,R) d \ln k """ if hasattr(M,'units'): M.units = 'Msun/hh' else: raise ValueError, '\n Input variable M must have units \n' rho_mean = self.cosmo.rho_crit0 * self.cosmo.OmegaM vol = M / rho_mean R = ( (3 * vol) / (4 * np.pi) )**(1./3) total = self.sigma2_R( R, z ) return total
[docs] def n_eff_approx( self, vol, dlogM, z ): dlnM = dlogM * np.log(10.) mass = np.array(vol) * 10**(dlogM/2.) * self.cu.Msunh sig2 = self.sigma2_M( mass, z ) sig = np.sqrt(sig2) t1 = 1.0 / sig mass = np.array(vol) * 10**(-dlogM/2.) * self.cu.Msunh sig2 = self.sigma2_M( mass, z ) sig = np.sqrt(sig2) t2 = 1.0 / sig dlnisig = np.abs( np.log(t1) - np.log(t2) ) neff = 6.0 * (dlnisig / dlnM) - 3.0 return neff
[docs] def frac_mass( self, sig, delta_c ): """ fraction of mass in halos per unit dln(1/sigma) """ nu = delta_c / sig nu_prime = np.sqrt(0.707) * nu ln_isig = np.log( 1.0/sig ) ln_gauss1 = np.exp( -(ln_isig-0.4)**2 / (2.0 * 0.6**2) ) t1 = 0.3222 * np.sqrt(2.0/np.pi) * nu_prime t2 = np.exp( -1.08 * nu_prime**2 / 2.0 ) t3 = (1.0 + 1.0/nu_prime**0.6 + 0.2*ln_gauss1) frac = t1 * t2 * t3 return frac # def plot_W2kR( self ): # """ Plots the window function for R=8 Mpc/h """ # Rw = 8.0 * self.CU.Mpc / self.CU.h # kw = np.array( 1.0/Rw ) # lnkw = np.log( kw ) # logkw = np.log10(kw) # N = 5000 # lnk_min = lnkw - 2.0 # lnk_max = lnkw + 4.0 # kk_min = np.exp(1.0)**lnk_min # kk_max = np.exp(1.0)**lnk_max # logk_min = np.log10( kk_min ) # logk_max = np.log10( kk_max ) # lnk = np.linspace( lnk_min, lnk_max, N ) # kk = np.exp(1.0)**lnk / self.CU.Mpc / self.CU.h # logk = np.log10( np.array(kk) ) # fig = plt.figure( figsize=(10,8) ) # fig.subplots_adjust( top=0.95, right = 0.95, # bottom=0.12, left=0.14 ) # ax = fig.add_subplot(111) # WW = self.W2lnkR( lnk, Rw ) # ax.plot( logk, np.log10(WW), lw=2.0, color='blue' ) # ax.plot( [logkw,logkw], [-10,10] ) # ax.set_xlabel( r'$\log(k \, [\rm{h/Mpc}])$', fontsize=26 ) # txt = r'$ \log \, [W_R(k)] $' # ax.set_ylabel( txt, fontsize=26 ) # ax.set_xlim( (logk_min,logk_max) ) # ax.set_ylim( (-6.0,0.5) ) # def plot_dsig2_dk( self ): # """ Plot the sigma^2(R) integrand = Delta^2(k) [W(k,R)]^2 for several # values of R """ # mpl.rcParams['xtick.labelsize'] = 20 # mpl.rcParams['ytick.labelsize'] = 20 # z = 0.0 # N = 100000 # R_arr = [800.0, 80.0, 8.0, 0.8, 0.08] # c_arr = ['red', 'orange', 'green', 'blue', 'purple'] # fig = plt.figure( figsize=(10,8) ) # fig.subplots_adjust( top=0.86, right = 0.95, # bottom=0.12, left=0.14 ) # ax = fig.add_subplot(111) # first plot the unitless power spectrum #-------------------------------------------------- # plt_lnk_min = -8.0 # plt_lnk_max = 6.0 # plt_kk_min = np.exp(1.0)**plt_lnk_min # plt_kk_max = np.exp(1.0)**plt_lnk_max # plt_logk_min = np.log10( plt_kk_min ) # plt_logk_max = np.log10( plt_kk_max ) # logk = np.linspace( plt_logk_min, plt_logk_max, N ) # kk = 10**logk / self.CU.Mpc / self.CU.h # lnk = np.log( np.array(kk) ) # D2k = self.Delta2k(kk, z) # ax.plot( logk, np.log10(D2k), lw=3.0, color='black' ) # for i in range(len(R_arr)): # color = c_arr[i] # Rw = R_arr[i] * self.CU.Mpc / self.CU.h # kw = np.array( 1.0/Rw ) # logkw = np.log10( kw ) # lnkw = np.log( kw ) # lnk_min = lnkw - 4.0 # lnk_max = lnkw + 8.0 # kk_min = np.exp(1.0)**lnk_min # kk_max = np.exp(1.0)**lnk_max # logk_min = np.log10( kk_min ) # logk_max = np.log10( kk_max ) # logk = np.linspace( logk_min, logk_max, N ) # kk = 10**logk # lnk = np.log( np.array(kk) ) # dsig2_dlogk = self.dsig2_dlogk( logk, Rw, z ) # dsig2_dlnk = self.dsig2_dlnk( lnk, Rw, z ) # ax.plot( logk, np.log10(dsig2_dlnk), lw=1.5, color=color, # alpha=0.5) # xx = (logkw - plt_logk_min) / (plt_logk_max - plt_logk_min) # ax.arrow( xx, 0.98, 0.0, -0.05, # transform=ax.transAxes, color=color ) # txt = str(np.array(Rw)) + r' $ \rm {Mpc/h}$' # ax.text( 0.02, 0.82-0.08*i, txt, # horizontalalignment='left', # verticalalignment='center', # transform=ax.transAxes, color=color, # fontsize=20) # ax.grid() # ax2 = ax.twiny() # ax2.plot( [1,1], [1,1], alpha=0.0 ) # xlo = np.log10( 1.0/np.exp(1.0)**plt_lnk_max ) # xhi = np.log10( 1.0/np.exp(1.0)**plt_lnk_min ) # ax2.set_xlim( xhi, xlo ) # ax.set_xlim( (plt_logk_min, plt_logk_max) ) # ax.set_ylim( (-7.4,1.7) ) # ax.set_xlabel( r'$\log(k \, [\rm{h/Mpc}])$', fontsize=26 ) # ax2.text( 0.5, 1.10, r'$\log(R \, [\rm{Mpc/h}])$', # horizontalalignment='center', # verticalalignment='center', # transform=ax2.transAxes, color='black', # fontsize=26) # txt = r'$ \log \, [\Delta^2(k) W_R(k)] = \log \,[$' # txt = txt + r'$ \frac{d \sigma^2_R}{ d \, \ln k }] $' # ax.set_ylabel( txt, fontsize=26 )