"""Cosmology-related functions for halo mass function and bias calculations."""
import numpy as np
from colossus.lss import mass_function, bias
from .config import HMF_NUM_POINTS
[docs]
def get_halo_mass_function(z, M_min=6, M_max=15, num_points=None,
mdef='vir', model='watson13'):
"""Compute the halo mass function at redshift z.
Parameters
----------
z : float
Redshift
M_min : float, optional
Minimum log10(M/M_sun) in physical units. Default is 6.
M_max : float, optional
Maximum log10(M/M_sun) in physical units. Default is 15.
num_points : int, optional
Number of mass bins. Default is from config.
mdef : str, optional
Mass definition: 'vir', 'fof', '200m', etc. Default is 'vir'.
model : str, optional
HMF model: 'watson13', 'tinker08', 'reed07', etc. Default is 'watson13'.
Returns
-------
log10_mass : ndarray
Array of log10(M/M_sun) in physical units
hmf : ndarray
Halo mass function dn/dln(M) in physical units of Mpc^-3
Notes
-----
Colossus expects mass in M_sun/h and returns dn/dlnM in (h/Mpc)^3.
To convert between physical and comoving units:
- Mass: M [M_sun/h] = M [M_sun] * h (factor out h by multiplying)
- HMF: dn/dlnM [Mpc^-3] = dn/dlnM [(h/Mpc)^3] * h^3 (evaluate h)
"""
if num_points is None:
num_points = HMF_NUM_POINTS
h = 0.7 # H0/100
# Create mass array in physical units [M_sun]
log10_mass_physical = np.linspace(M_min, M_max, num_points)
mass_physical = 10**log10_mass_physical
# Convert to comoving units for Colossus [M_sun/h]
# To factor out h: M [M_sun/h] = M [M_sun] * h
mass_comoving = mass_physical * h
# Get HMF from Colossus in comoving units [(h/Mpc)^3]
hmf_comoving = mass_function.massFunction(
mass_comoving, z,
mdef=mdef,
model=model,
q_out='dndlnM'
)
# Convert HMF to physical units [Mpc^-3]
# To evaluate h: dn/dM [Mpc^-3] = dn/dlnM [(h/Mpc)^3] * h^3 * ln(10)
hmf_physical = hmf_comoving * np.log(10) * h**3
return log10_mass_physical, hmf_physical
[docs]
def get_halo_bias(M, z, mdef='vir', model='tinker10'):
"""Compute halo bias as a function of mass and redshift.
Parameters
----------
M : float or array_like
Halo mass in M_sun (physical units)
z : float
Redshift
mdef : str, optional
Mass definition. Default is 'vir'.
model : str, optional
Bias model. Default is 'tinker10'.
Returns
-------
b : float or ndarray
Halo bias (dimensionless)
Notes
-----
Colossus expects mass in M_sun/h.
To factor out h: M [M_sun/h] = M [M_sun] * h
"""
h = 0.7 # H0/100
M_comoving = np.atleast_1d(M) * h # Convert to M_sun/h
return bias.haloBias(M_comoving, model=model, z=z, mdef=mdef)
[docs]
class HaloMassFunction:
"""Class for computing and caching halo mass functions.
This class provides a convenient interface for HMF calculations
with optional caching of results.
Attributes
----------
z : float
Redshift
M_min : float
Minimum log10(M/M_sun)
M_max : float
Maximum log10(M/M_sun)
num_points : int
Number of mass bins
mdef : str
Mass definition
model : str
HMF model
"""
def __init__(self, z, M_min=6, M_max=15, num_points=None,
mdef='vir', model='watson13'):
"""Initialize HMF calculator.
Parameters
----------
z : float
Redshift
M_min : float, optional
Minimum log10(M/M_sun)
M_max : float, optional
Maximum log10(M/M_sun)
num_points : int, optional
Number of mass bins
mdef : str, optional
Mass definition
model : str, optional
HMF model
"""
self.z = z
self.M_min = M_min
self.M_max = M_max
self.num_points = num_points if num_points else HMF_NUM_POINTS
self.mdef = mdef
self.model = model
# Cache
self._log10_mass = None
self._hmf = None
self._bias = None
@property
def log10_mass(self):
"""Get log10 mass array."""
if self._log10_mass is None:
self._compute()
return self._log10_mass
@property
def mass(self):
"""Get linear mass array."""
return 10**self.log10_mass
@property
def hmf(self):
"""Get halo mass function."""
if self._hmf is None:
self._compute()
return self._hmf
@property
def halo_bias(self):
"""Get halo bias."""
if self._bias is None:
self._bias = get_halo_bias(self.mass, self.z,
mdef=self.mdef, model='tinker10')
return self._bias
def _compute(self):
"""Compute HMF and cache results."""
self._log10_mass, self._hmf = get_halo_mass_function(
self.z, self.M_min, self.M_max, self.num_points,
self.mdef, self.model
)
def __repr__(self):
return (f"HaloMassFunction(z={self.z}, M_min={self.M_min}, "
f"M_max={self.M_max}, model='{self.model}')")