import numpy as np
import pyccl as ccl
from scipy import interpolate, integrate
from .Base import BaseBFGProfiles, hyper_params
from .Schneider19 import Gas, DarkMatterBaryon, TwoHalo
from ..utils.constants import (G, Msun_to_Kg, Mpc_to_m, Pth_to_Pe, m_p, m_to_cm, kb_cgs, sigma_T_cgs, c_cgs, m_e_cgs, X)
from ..utils.Tabulate import _set_parameter, _get_parameter
from ..utils.Xray import EmissivityTable
from .Schneider19 import model_params as S19_mp
from .Arico20 import model_params as A20_mp
from .Mead20 import model_params as M20_mp
model_params = list({*S19_mp, *A20_mp, *M20_mp})
#Technically P(r -> infty) is zero, but we may need finite
#value for numerical reasons (interpolator). This is a
#computatational constant.
Pressure_at_infinity = 1e-200
__all__ = ['Pressure', 'NonThermalFrac', 'NonThermalFracGreen20',
'Temperature', 'ThermalSZ', 'ElectronPressure', 'GasNumberDensity',
'Emissivity', 'Metallicity', 'XrayCounts', 'XraySkyCounts']
class BaseThermodynamicProfile(BaseBFGProfiles):
model_param_names = model_params
@property
def model_params(self):
"""
Returns a dictionary containing all model parameters and their current values.
Returns
-------
params : dict
Dictionary of model parameters.
"""
if hasattr(self, 'prof4params'):
params = {k:v for k,v in vars(self.prof4params).items() if k in self.model_param_names}
else:
params = {k:v for k,v in vars(self).items() if k in self.model_param_names}
return params
@property
def hyper_params(self):
"""
Returns a dictionary containing all hyper parameters oof the calculation and their current values.
Returns
-------
params : dict
Dictionary of hyper parameters.
"""
if hasattr(self, 'prof4params'):
params = {k:v for k,v in vars(self.prof4params).items() if k in self.hyper_param_names}
else:
params = {k:v for k,v in vars(self).items() if k in self.hyper_param_names}
params['c_M_relation'] = self._c_M_relation #Swap this one specifically
params['use_fftlog_projection'] = self._use_fftlog_projection #This one isn't saved normally so do it here
return params
[docs]class Pressure(BaseThermodynamicProfile):
"""
Class for computing the gas pressure profile in halos.
This class extends `SchneiderProfiles` to compute the gas pressure profile within halos
under the assumption of hydrostatic equilibrium. The gas pressure is derived using a
total mass profile and a gas density profile. We define a pressure gradient from the
assumption of hydrostatic equilibrium, and integrate to obtain the pressure.
This gives only the *total gas pressure*. If you want the electron pressure
see `ElectronPressure`, and if you want to only the thermal/non-thermal
pressure see `NonThermalFrac`.
Inherits from
-------------
SchneiderProfiles : Base class for halo profiles.
Parameters
----------
gas : Gas, optional
An instance of the `Gas` class defining the gas density profile. If not provided,
a default `Gas` object is created using `kwargs`.
darkmatterbaryon : DarkMatterBaryon, optional
An instance of the `DarkMatterBaryon` class defining the combined dark matter
and baryonic mass profile. If not provided, a default `DarkMatterBaryon` object
is created using `kwargs`.
nonthermal_model : object, optional
An instance defining a model for non-thermal pressure contributions. Default is None.
**kwargs
Additional keyword arguments passed to initialize the `Gas`, `DarkMatterBaryon`,
and other parameters from `SchneiderProfiles`.
Notes
-----
- This class calculates the pressure assuming hydrostatic equilibrium, which gives:
.. math::
\\frac{dP}{dr} = -\\frac{GM(<r)\\rho_{\\text{gas}}(r)}{r^2}
where:
- \( G \) is the gravitational constant.
- \( M(<r) \) is the cumulative mass within radius \( r \).
- \( \\rho_{\\text{gas}}(r) \) is the gas density profile.
- The gas pressure \( P(r) \) is then obtained by integrating \( dP/dr \):
.. math::
P(r) = \\int_r^{\\infty} -\\frac{GM(<r')\\rho_{\\text{gas}}(r')}{r'^2} r' d\\ln r'
- The integration is performed numerically, and the result is converted to CGS units
for practical applications.
- An exponential cutoff is applied to the profile to prevent numerical overflow at large radii.
Methods
-------
_real(cosmo, r, M, a)
Computes the gas pressure profile based on the given cosmology, radii, mass,
scale factor, and mass definition.
"""
model_param_names = model_params
def __init__(self, gas = None, darkmatterbaryon = None, **kwargs):
self.Gas = gas
self.DarkMatterBaryon = darkmatterbaryon
#The subtraction in DMB case is so we only have the 1halo term
if self.Gas is None: self.Gas = Gas(**kwargs)
if self.DarkMatterBaryon is None: self.DarkMatterBaryon = DarkMatterBaryon(**kwargs) - TwoHalo(**kwargs)
#Now make sure the cutoff is sufficiently high
#We don't want small cutoff of 1-halo term when computing the TRUE pressure profile.
#The cutoff is reapplied to the derives pressure profiles in _real()
_set_parameter(self.Gas, 'cutoff', 1000)
_set_parameter(self.DarkMatterBaryon, 'cutoff', 1000)
self.prof4params = self.Gas
super().__init__(**kwargs)
[docs] def _real(self, cosmo, r, M, a):
"""
Computes the gas pressure profile using hydrostatic equilibrium.
This method calculates the gas pressure profile for a specified cosmology, radii,
halo mass, and scale factor. The pressure is computed by integrating the pressure
gradient derived from the hydrostatic equilibrium condition. The final profile is
in units of comoving volume. Use a factor of 1/a^3 (not 1/a^4) to convert to physical pressure.
Parameters
----------
cosmo : object
A CCL cosmology instance containing the cosmological parameters used for calculations.
r : array_like
Radii at which to evaluate the pressure profile, in comoving Mpc.
M : float or array_like
Halo mass or array of halo masses, in solar masses.
a : float
Scale factor, related to redshift by `a = 1 / (1 + z)`.
Returns
-------
prof : ndarray
Pressure profile corresponding to the input radii and halo mass, in CGS units.
Notes
-----
- The pressure gradient is calculated using the formula:
.. math::
\\frac{dP}{dr} = -\\frac{GM(<r)\\rho_{\\text{gas}}(r)}{r^2}
where:
- \( G \) is the gravitational constant.
- \( M(<r) \) is the cumulative mass within radius \( r \), calculated from the
total density profile.
- \( \\rho_{\\text{gas}}(r) \) is the gas density profile.
- The pressure profile \( P(r) \) is obtained by numerically integrating the pressure gradient:
.. math::
P(r) = \\int_r^{\\infty} -\\frac{GM(<r')\\rho_{\\text{gas}}(r')}{r'^2} r' d\\ln r'
- The integration is performed from large radii towards smaller radii to satisfy the boundary
condition \( P(r \\to \\infty) = 0 \).
- The profile is converted to CGS units using appropriate conversion factors.
Examples
--------
Compute the gas pressure profile for a given cosmology and halo:
>>> pressure_model = Pressure(gas=gas_profile, darkmatterbaryon=dmb_profile, cosmo=my_cosmology)
>>> r = np.logspace(-2, 1, 50) # Radii in comoving Mpc
>>> M = 1e14 # Halo mass in solar masses
>>> a = 0.5 # Scale factor corresponding to redshift z
>>> pressure_profile = pressure_model._real(my_cosmology, r, M, a)
"""
r_use = np.atleast_1d(r)
M_use = np.atleast_1d(M)
z = 1/a - 1
R = self.mass_def.get_radius(cosmo, M_use, a)/a #in comoving Mpc
r_integral = np.geomspace(self.r_min_int, self.r_max_int, self.r_steps)
rho_total = self.DarkMatterBaryon.real(cosmo, r_integral, M_use, a)
rho_gas = self.Gas.real(cosmo, r_integral, M_use, a)
#Integrate total density profile to get the cumulative mass distribution
dlnr = np.log(r_integral[1]) - np.log(r_integral[0])
dV = 4 * np.pi * r_integral**3 * dlnr
M_total = integrate.cumulative_simpson(dV * rho_total, axis = -1, initial = 0) + dV[0] * rho_total[:, [0]]
#Assuming hydrostatic equilibrium to get dP/dr = -G*M(<r)*rho(r)/r^2
dP_dr = - G * M_total * rho_gas / r_integral**2
#Make it have the right shape that ccl expects (size(M), size(r))
if len(dP_dr.shape) < 2:
dP_dr = dP_dr[np.newaxis, :]
#integrate to get actual pressure, P(r). Boundary condition is P(r -> infty) = 0.
#So we start from the boundary and integrate inwards. We reverse array once to
#flip the integral direction, and flip it second time so P(r) goes from r = 0 to r = infty
#We use trapezoid rule here because simpson was causing odd oscillatory errors because
#Some profiles have sharp transitions in their pressure/gas profiles (eg. Mead)
intgr = (dP_dr * r_integral)[:, ::-1] * dlnr
prof = -np.array([integrate.cumulative_trapezoid(intgr[i], initial = 0)[::-1] + intgr[i, 0] for i in range(intgr.shape[0])])
prof = interpolate.PchipInterpolator(np.log(r_integral), np.log(prof + Pressure_at_infinity), axis = 1, extrapolate = False)
prof = np.exp(prof(np.log(r_use))) - Pressure_at_infinity
prof = np.where(np.isfinite(prof), prof, 0) #Get rid of pesky NaN and inf values if they exist! They break CCL spline interpolator
#Convert to CGS. Using only one factor of Mpc_to_m is correct here!
#The factor of 1/a is so the temperature piece of Pressure = Temp x density
#is always in physical units. So the comoving units are just in the density.
prof = prof * (Msun_to_Kg * 1e3) / (Mpc_to_m * 1e2)
prof = prof / a
#Now do cutoff
arg = (r_use[None, :] - self.cutoff)
arg = np.where(arg > 30, np.inf, arg) #This is to prevent an overflow in the exponential
kfac = 1/( 1 + np.exp(2*arg) ) #Extra exponential cutoff
prof = prof * kfac
#Handle dimensions so input dimensions are mirrored in the output
if np.ndim(r) == 0: prof = np.squeeze(prof, axis=-1)
if np.ndim(M) == 0: prof = np.squeeze(prof, axis=0)
return prof
[docs]class NonThermalFrac(BaseThermodynamicProfile):
"""
Class for computing the non-thermal pressure fraction profile in halos.
This class extends `SchneiderProfiles` to compute the fraction of pressure that is
non-thermal within halos. The non-thermal fraction is modelled in a redshift, and radius
dependent manner, following Equations 15/16 in `Pandey et. al 2025 <https://arxiv.org/pdf/2401.18072>`_.
This can model be applied to any profiles using simple multiplication of the initialized classes,
`ThermalPressure = Pressure(**kwargs) * (1 - NonThermalFrac(**kwargs))`
The `real()` method of `ThermalPressure` will then account for non-thermal pressure
effects as well.
Inherits from
-------------
SchneiderProfiles : Base class for halo profiles.
Parameters
----------
alpha_nt : float
Normalization factor for the non-thermal pressure fraction.
nu_nt : float
Parameter controlling the redshift dependence of the non-thermal fraction.
gamma_nt : float
Parameter controlling the radial dependence of the non-thermal fraction.
**kwargs
Additional keyword arguments passed to initialize other parameters from `SchneiderProfiles`.
Notes
-----
The `NonThermalFrac` class is used to compute the non-thermal pressure fraction,
which represents the fraction of total pressure that is non-thermal (e.g., due to
turbulence or magnetic fields) as a function of radius and redshift. This fraction
is used to modify the total pressure profile, accounting for contributions that are
not purely thermal.
The non-thermal pressure fraction \( f_{\\text{nt}}(r, z) \) is calculated using:
.. math::
f_{\\text{nt}}(r, z) = \\alpha_{\\text{nt}} \\times f_z \\times \\left( \\frac{r}{R} \\right)^{\\gamma_{\\text{nt}}}
where:
- \( \\alpha_{\\text{nt}} \) is the normalization factor.
- \( f_z \) is the redshift-dependent factor, defined as:
.. math::
f_z = \\min\\left[(1 + z)^{\\nu_{\\text{nt}}}, \\left(f_{\\text{max}} - 1\\right) \\tanh\\left(\\nu_{\\text{nt}} z\\right) + 1\\right]
- \( R \) is the halo radius based on the mass definition.
- \( \\gamma_{\\text{nt}} \) controls the radial dependence.
"""
def __init__(self, alpha_nt, nu_nt, gamma_nt, **kwargs):
super().__init__(**kwargs)
self.alpha_nt = alpha_nt
self.nu_nt = nu_nt
self.gamma_nt = gamma_nt
def _real(self, cosmo, r, M, a):
r_use = np.atleast_1d(r)
M_use = np.atleast_1d(M)
z = 1/a - 1
R = self.mass_def.get_radius(cosmo, M_use, a)/a #in comoving Mpc
f_max = 6**-self.gamma_nt/self.alpha_nt
f_z = np.min([(1 + z)**self.nu_nt, (f_max - 1)*np.tanh(self.nu_nt * z) + 1])
f_nt = self.alpha_nt * f_z * (r_use/R[:, None])**self.gamma_nt
f_nt = np.clip(f_nt, 0, 1) #Enforce 0 < f_nt < 1
prof = f_nt #Rename just for consistency sake
#Handle dimensions so input dimensions are mirrored in the output
if np.ndim(r) == 0: prof = np.squeeze(prof, axis=-1)
if np.ndim(M) == 0: prof = np.squeeze(prof, axis=0)
return prof
[docs]class NonThermalFracGreen20(BaseThermodynamicProfile):
"""
Class for computing the non-thermal pressure fraction profile using the Green et al. (2020) model.
Notes
-----
The model is based on parameters calibrated to simulations and is specifically defined
with respect to \( R_{200m} \), the radius within which the mean density is 200 times
the mean matter density of the universe.
The non-thermal pressure fraction \( f_{\\text{nt}}(r) \) is calculated using:
.. math::
f_{\\text{nt}}(r) = 1 - a \\left(1 + \\exp\\left(-\\left(\\frac{x}{b}\\right)^c\\right)\\right)
\\left(\\frac{\\nu_M}{4.1}\\right)^{\\frac{d}{1 + \\left(\\frac{x}{e}\\right)^f}}
where:
- \( x = \\frac{r}{R_{200m}} \)
- \( \\nu_M = \\frac{1.686}{\\sigma(M_{200m})} \) is the peak height parameter.
- \( a, b, c, d, e, f \) are model parameters calibrated to fit simulation data.
There are no free parameters in this model; it is completely specified by the halo mass and redshift.
"""
def _real(self, cosmo, r, M, a):
r_use = np.atleast_1d(r)
M_use = np.atleast_1d(M)
z = 1/a - 1
R = self.mass_def.get_radius(cosmo, M_use, a)/a #in comoving Mpc
#They define the model with R200m, so gotta use that redefinition here.
mdef = ccl.halos.massdef.MassDef(200, 'matter')
cnvrt = ccl.halos.mass_translator(mass_in = self.mass_def, mass_out = mdef, concentration = self.mass_def.concentration)
M200m = cnvrt(cosmo, M_use, a)
R200m = mdef.get_radius(cosmo, M_use, a)/a #in comoving distance
x = r_use/R200m[:, None]
nu_M = 1.686/ccl.sigmaM(cosmo, M200m, a)
nu_M = nu_M[:, None]
A, b, c, d, e, f = 0.495, 0.719, 1.417,-0.166, 0.265, -2.116
nth = 1 - A * (1 + np.exp(-(x/b)**c)) * (nu_M/4.1)**(d/(1 + (x/e)**f))
nth = np.clip(nth, 0, 1)
prof = nth #Rename just for consistency sake
#Handle dimensions so input dimensions are mirrored in the output
if np.ndim(r) == 0: prof = np.squeeze(prof, axis=-1)
if np.ndim(M) == 0: prof = np.squeeze(prof, axis=0)
return prof
[docs]class ElectronPressure(Pressure):
"""
Class for computing the electron pressure profile in halos.
This class extends the `Pressure` class to compute the electron pressure
profile from the total gas pressure. The conversion factor is
\( P_{\\text{e}} = P_{\\text{th}} \\times P_{\\text{th-to-Pe}} \), where
\( P_{\\text{th-to-Pe}} = (4 - 2Y)/(8 - 5Y), with Y = 0.24\).
Inherits from
-------------
Pressure : Base class for computing gas pressure profiles in halos.
Notes
-----
The `ElectronPressure` class is used to compute the electron pressure profile
within halos, which is relevant for understanding the thermal Sunyaev-Zel'dovich
(tSZ) effect. The conversion is done assuming pressure equilibrium between electrons
and protons in a system dominated by hydrogen and helium species (hence the use of Y).
"""
def _real(self, cosmo, r, M, a):
prof = Pth_to_Pe * super()._real(cosmo, r, M, a)
return prof
[docs]class GasNumberDensity(BaseThermodynamicProfile):
"""
Class for computing the gas number density profile in halos.
This class extends `SchneiderProfiles` to compute the gas number density profile
within halos. The number density is derived from the gas density profile by dividing
by the mean molecular weight and the mass of the proton, and then converting to
proper CGS units.
Inherits from
-------------
SchneiderProfiles : Base class for halo profiles.
Parameters
----------
gas : Gas, optional
An instance of the `Gas` class defining the gas density profile. If not provided,
a default `Gas` object is created using `kwargs`.
mean_molecular_weight : float, optional
Mean molecular weight of the gas. Default is 1.15, which is typical for ionized
hydrogen with a small fraction of helium.
**kwargs
Additional keyword arguments passed to initialize the `Gas` profile and other
parameters from `SchneiderProfiles`.
Notes
-----
The `GasNumberDensity` class is used to compute the number density of gas particles
in halos, which is relevant for understanding the baryonic content of halos and for
modeling various astrophysical processes, such as cooling and star formation.
The gas number density \( n_{\\text{gas}} \) is calculated by dividing the gas density
profile \( \\rho_{\\text{gas}} \) by the mean molecular weight and the mass of the proton:
.. math::
n_{\\text{gas}}(r) = \\frac{\\rho_{\\text{gas}}(r)}{\\mu \\cdot m_p}
where:
- \( \\mu \) is the mean molecular weight of the gas.
- \( m_p \) is the mass of the proton.
- The result is converted to the proper units (number per cubic centimeter).
"""
def __init__(self, gas = None, **kwargs):
self.Gas = gas
if self.Gas is None: self.Gas = Gas(**kwargs)
super().__init__(**kwargs)
assert 'mean_molecular_weight' in kwargs, "Please provide a mean_molecular_weight input"
self.mean_molecular_weight = kwargs['mean_molecular_weight']
#Convert from Msun --> n_proton
#Then convert from 1/Mpc^3 to 1/cm^3
#The projected profile will just be in units of Mpc/cm^3
self.factor = 1 / (self.mean_molecular_weight * m_p) / (Mpc_to_m * m_to_cm)**3
self.prof4params = self.Gas
#Need to explicitly call real and projected routines here.
#We call "_real" since if we define "real" over "_real" then CCL will complain
def _real(self, cosmo, r, M, a): return self.Gas._real(cosmo, r, M, a) * self.factor
[docs] def projected(self, cosmo, r, M, a): return self.Gas.projected(cosmo, r, M, a) * self.factor
[docs]class Temperature(BaseThermodynamicProfile):
"""
Class for computing the temperature profile in halos.
The temperature is derived from the thermal pressure and the number density profiles,
of a species using the ideal gas law. The temperature profile is important for understanding
the thermal state of the intracluster medium and its impact on various astrophysical processes.
For this model to be correct, the input pressure must be the *thermal pressure*, i.e. the
non-thermal pressure must have already been accounted for in the model passed to this class.
Parameters
----------
pressure : Pressure, optional
An instance of the `Pressure` class defining the thermal gas pressure profile.
If non-thermal pressure is relevant for your problem, it must be included in this
profile; see `Pressure` or `NonThermalFrac` for more details.
If this parameter is not provided, a default `Pressure` object is created using `kwargs`.
gasnumberdensity : GasNumberDensity, optional
An instance of the `GasNumberDensity` class defining the gas number density profile.
If not provided, a default `GasNumberDensity` object is created using `kwargs`.
**kwargs
Additional keyword arguments passed to initialize the `Pressure`, `GasNumberDensity`,
and other parameters from `SchneiderProfiles`.
Notes
-----
The `Temperature` class computes the temperature profile of the gas in halos by dividing
the gas pressure by the gas number density and the Boltzmann constant. This calculation
assumes the ideal gas law, which relates pressure, number density, and temperature.
The gas temperature \( T \) is calculated using:
.. math::
T(r) = \\frac{P}(r)}{n(r) \\cdot k_B}
where:
- \( P(r) \) is the Thermal pressure profile of a species.
- \( n(r) \) is the number density profile of a species.
- \( k_B \) is the Boltzmann constant (in eV).
"""
def __init__(self, thermalpressure = None, gasnumberdensity = None, **kwargs):
self.Pressure = thermalpressure
self.GasNumberDensity = gasnumberdensity
if self.Pressure is None: self.Pressure = Pressure(**kwargs) * (1 - NonThermalFrac(**kwargs))
if self.GasNumberDensity is None: self.GasNumberDensity = GasNumberDensity(**kwargs)
super().__init__(**kwargs)
if hasattr(self.Pressure, 'prof4params'):
self.prof4params = self.Pressure.prof4params
elif hasattr(self.GasNumberDensity, 'prof4params'):
self.prof4params = self.GasNumberDensity.prof4params
else:
self.prof4params = self
def _real(self, cosmo, r, M, a):
P = self.Pressure.real(cosmo, r, M, a)
n = self.GasNumberDensity.real(cosmo, r, M, a)
#We'll have instances of n == 0, which isn't a problem so let's ignore
#warnings of divide errors, because we know they happen here.
#Instead we will fix them by replacing the temperature with 0s,
#since there is no gas in those regions to use anyway.
with np.errstate(divide = 'ignore', invalid = 'ignore'):
prof = P/(n * kb_cgs)
prof = np.where(n == 0, 0, prof)
return prof
[docs] def projected(self, cosmo, r, M, a):
"""
Compute the projected temperature profile along the line of sight.
This method calculates the "average temperature" along the line of sight,
which is a physically meaningful quantity for comparing with observations
such as X-ray or Sunyaev-Zel'dovich measurements. It differs from the
"integrated temperature," which lacks physical relevance in most astrophysical contexts.
Parameters
----------
cosmo : Cosmology
The cosmology object containing cosmological parameters.
r : array_like
The projected radial distances at which to compute the profile, in units of Mpc/h.
M : float
The halo mass, in units of solar masses.
a : float
The scale factor of the Universe.
Returns
-------
prof : array_like
The projected average temperature profile, in units of eV.
Notes
-----
The projected temperature is computed using the ideal gas law:
.. math::
T_{\text{proj}}(r) = \\frac{P_{\text{proj}}(r)}{n_{\text{proj}}(r) \\cdot k_B}
where:
- \( P_{\text{proj}}(r) \) is the projected thermal pressure profile.
- \( n_{\text{proj}}(r) \) is the projected number density profile.
- \( k_B \) is the Boltzmann constant (in eV).
Regions with zero gas density (\( n_{\text{proj}}(r) = 0 \)) are assigned a temperature of 0
to avoid division errors, as these regions lack gas to support a meaningful temperature.
"""
P = self.Pressure.projected(cosmo, r, M, a)
n = self.GasNumberDensity.projected(cosmo, r, M, a)
#We'll have instances of n == 0, which isn't a problem so let's ignore
#warnings of divide errors, because we know they happen here.
#Instead we will fix them by replacing the temperature with 0s,
#since there is no gas in those regions to use anyway.
with np.errstate(divide = 'ignore', invalid = 'ignore'):
prof = P/(n * kb_cgs)
prof = np.where(n == 0, 0, prof)
return prof
[docs]class ThermalSZ(BaseThermodynamicProfile):
"""
Class for computing the thermal Sunyaev-Zel'dovich (tSZ) effect profile in halos.
This class extends `SchneiderProfiles` to compute the tSZ effect, which is caused
by the inverse Compton scattering of cosmic microwave background (CMB) photons
off hot electrons in the intracluster medium of galaxy clusters. The tSZ effect
is represented by the Compton-y parameter, which is proportional to the line-of-sight
integral of the electron pressure.
In practice, this scale uses the `projected` method of the input `pressure` object.
It accounts for the right units, to provide a dimensionless compton-y parameter.
Inherits from
-------------
SchneiderProfiles : Base class for halo profiles.
Parameters
----------
pressure : Pressure, optional
An instance of the `Pressure` class defining the thermal gas pressure profile.
If not provided, a default `Pressure` object is created using `kwargs`.
**kwargs
Additional keyword arguments passed to initialize the `Pressure` profile and other
parameters from `SchneiderProfiles`.
Notes
-----
The `ThermalSZ` class computes the tSZ effect by calculating the projected electron
pressure profile along the line of sight.
- The tSZ effect is computed by projecting the electron pressure along the line of sight:
.. math::
y(r) = \\frac{\\sigma_T}{m_e c^2} \\int P_{\\text{e}}(r') \\, dr'
where \( P_{\\text{e}}(r') \) is the electron pressure profile.
Methods
-------
Pgas_to_Pe(cosmo, r, M, a)
Returns the conversion factor from gas pressure to electron pressure.
projected(cosmo, r, M, a)
Computes the projected tSZ profile (Compton-y parameter) based on the given
cosmology, radii, mass, scale factor, and mass definition.
real(cosmo, r, M, a)
This computes the real pressure profile scaled by the same constant as the
projected one. Generally not to be used unless you know what units you
want.
"""
def __init__(self, pressure = None, **kwargs):
self.Pressure = pressure
if self.Pressure is None: self.Pressure = Pressure(**kwargs)
super().__init__(**kwargs)
if hasattr(self.Pressure, 'prof4params'):
self.prof4params = self.Pressure.prof4params
else:
self.prof4params = self
[docs] def Pgas_to_Pe(self, cosmo, r, M, a):
"""
Returns the precomputed conversion factor from gas pressure to electron pressure.
Can be redefined by user if they wish to use a different (mass-dependent) value
for this quantity.
"""
return Pth_to_Pe
def _real(self, cosmo, r, M, a):
r_use = np.atleast_1d(r)
M_use = np.atleast_1d(M)
z = 1/a - 1
R = self.mass_def.get_radius(cosmo, M_use, a)/a #in comoving Mpc
#Now a series of units changes to the projected profile.
prof = self.Pressure.real(cosmo, r_use, M_use, a) #generate profile
prof = prof * (Mpc_to_m * 1e2) #Line-of-sight integral is done in Mpc, we want cm
prof = prof * sigma_T_cgs/(m_e_cgs*c_cgs**2) #Convert to SZ (dimensionless units)
prof = prof * self.Pgas_to_Pe(cosmo, r_use, M_use, a) #Then convert from gas pressure to electron pressure
return prof
def _projected(self, cosmo, r, M, a):
#This is the 1/(1 + z) that shows up in the tracer kernel
return super()._projected(cosmo, r, M, a) * a
[docs]class Emissivity(BaseThermodynamicProfile):
"""
Class representing the X-ray emissivity profile of halo gas.
This class computes the radial emissivity profile by combining a gas
temperature profile, a metallicity profile, and an external emissivity
lookup function or table. The emissivity is evaluated as a function of
the local gas temperature and metallicity at each radius.
The temperature and metallicity profiles can either be supplied directly
or automatically instantiated using the `Temperature` and `Metallicity`
classes with the provided keyword arguments.
See `BaseThermodynamicProfile` for more docstring details.
Notes
-----
The emissivity is computed by evaluating an external emissivity function
or table:
.. math::
\\epsilon(r) = \\Lambda(T(r), Z(r), a)
where:
- :math:`T(r)` is the gas temperature profile.
- :math:`Z(r)` is the gas metallicity profile.
- :math:`a` is the cosmological scale factor.
- :math:`\\Lambda` is the emissivity lookup function provided through
``EmissivityTable``.
The emissivity table is expected to accept temperature, metallicity, and
scale factor as inputs and return the corresponding emissivity. This
output emissivity is intended to be in photon counts, though there is
no problem if the user supplies any quantity they want.
If the temperature profile defines the attribute ``prof4params`` (used
internally for parameter handling in some profile implementations), this
class forwards that attribute to ensure consistent parameter handling.
Parameters
----------
EmissivityTable : callable
Function or lookup table returning emissivity as a function of
temperature, metallicity, and scale factor. It must have the form
``EmissivityTable(T, Z, a)``
where ``T`` and ``Z`` are arrays matching the radial grid.
temperature : BaseThermodynamicProfile, optional
Temperature profile object used to compute :math:`T(r)`. If not
provided, a default `Temperature` profile is created using the
provided ``kwargs``.
metallicity : BaseThermodynamicProfile, optional
Metallicity profile object used to compute :math:`Z(r)`. If not
provided, a default `Metallicity` profile is created using the
provided ``kwargs``.
**kwargs
Additional keyword arguments passed to the internally constructed
`Temperature` and `Metallicity` profiles when they are not supplied.
"""
def __init__(self, EmissivityTable, temperature = None, metallicity = None, **kwargs):
self.Temperature = temperature if temperature is not None else Temperature(**kwargs)
self.Metallicity = metallicity if metallicity is not None else Metallicity(**kwargs)
self.EmissivityTable = EmissivityTable
super().__init__()
if hasattr(self.Temperature, 'prof4params'):
self.prof4params = self.Temperature.prof4params
else:
self.prof4params = self
def _real(self, cosmo, r, M, a):
r_use = np.atleast_1d(r)
M_use = np.atleast_1d(M)
z = 1/a - 1
R = self.mass_def.get_radius(cosmo, M_use, a)/a #in comoving Mpc
T = self.Temperature.real(cosmo, r_use, M, a)
Z = self.Metallicity.real(cosmo, r_use, M, a)
E = self.EmissivityTable(T, Z, a)
prof = E #Just renaming for simplicity
return prof
[docs]class XrayCounts(BaseThermodynamicProfile):
"""
Class representing the radial X-ray count-rate profile of halo gas.
This class computes the X-ray counts profile by combining an emissivity
profile with electron and hydrogen number density profiles. The resulting
profile is proportional to the local X-ray emission rate per unit volume,
assuming the standard dependence on the product of electron and hydrogen
number densities.
The emissivity profile must be supplied explicitly. If the electron or
hydrogen number density profiles are not provided, they are constructed
automatically using `GasNumberDensity`. By default, the hydrogen number
density is initialized with a mean molecular weight appropriate for hydrogen
of X = 0.76. The user can alter this by supplying their own profile.
See `BaseThermodynamicProfile` for more docstring details.
Notes
-----
The X-ray counts profile is computed as
.. math::
C(r) = n_e(r)\,n_{\\rm H}(r)\,\\Lambda(r),
where:
- :math:`n_e(r)` is the electron number density profile.
- :math:`n_{\\rm H}(r)` is the hydrogen number density profile.
- :math:`\\Lambda(r)` is the emissivity profile returned by ``self.Emissivity``.
In this implementation, the emissivity profile is produced by an
`Emissivity` object, which itself will depend on gas temperature,
metallicity, and scale factor.
Parameters
----------
emissivity : BaseThermodynamicProfile
Emissivity profile object used to compute :math:`\\Lambda(r)`.
electronnumberdensity : BaseThermodynamicProfile, optional
Electron number density profile object used to compute :math:`n_e(r)`.
If not provided, a default `GasNumberDensity` profile is constructed
from ``kwargs``.
hydrogennumberdensity : BaseThermodynamicProfile, optional
Hydrogen number density profile object used to compute
:math:`n_{\\rm H}(r)`. If not provided, a default `GasNumberDensity`
profile is constructed from ``kwargs``, with
``mean_molecular_weight = 1/X``.
**kwargs
Additional keyword arguments passed to internally constructed
`GasNumberDensity` profiles.
"""
def __init__(self, emissivity, electronnumberdensity = None, hydrogennumberdensity = None, **kwargs):
self.Emissivity = emissivity
self.ElectronNumberDensity = electronnumberdensity
self.HydrogenNumberDensity = hydrogennumberdensity
np_kwargs = {k:v for k,v in kwargs.items() if k != 'mean_molecular_weight'} #Drop this so we can replace it after (if user doesn't supply it)
if self.ElectronNumberDensity is None: self.ElectronNumberDensity = GasNumberDensity(**kwargs)
if self.HydrogenNumberDensity is None: self.HydrogenNumberDensity = GasNumberDensity(**np_kwargs, mean_molecular_weight = 1/X)
super().__init__()
if hasattr(self.ElectronNumberDensity, 'prof4params'):
self.prof4params = self.ElectronNumberDensity.prof4params
elif hasattr(self.HydrogenNumberDensity, 'prof4params'):
self.prof4params = self.HydrogenNumberDensity.prof4params
else:
self.prof4params = self
def _real(self, cosmo, r, M, a):
r_use = np.atleast_1d(r)
M_use = np.atleast_1d(M)
z = 1/a - 1
R = self.mass_def.get_radius(cosmo, M_use, a)/a #in comoving Mpc
ne = self.ElectronNumberDensity.real(cosmo, r_use, M, a)
nH = self.HydrogenNumberDensity.real(cosmo, r_use, M, a)
J = self.Emissivity.real(cosmo, r_use, M, a)
prof = ne * nH * J
return prof
[docs]class XraySkyCounts(BaseThermodynamicProfile):
"""
Class representing the observable X-ray sky-counts profile of halo gas.
This class converts the intrinsic X-ray volume emissivity/count profile into
a profile more directly suited for observed X-ray surface brightness or sky
counts. It wraps an `XrayCounts` profile and applies the appropriate geometric
and cosmological conversion factors so that the resulting quantity corresponds
to an observed line-of-sight sky signal rather than a purely local volumetric
emission rate.
See `BaseThermodynamicProfile` for more docstring details.
Notes
-----
The underlying `XrayCounts` profile is first evaluated as a local volumetric
quantity proportional to
.. math::
C(r) = n_e(r)\,n_{\\rm H}(r)\,\\Lambda(r),
where :math:`n_e` is the electron number density, :math:`n_{\\rm H}` is the
hydrogen number density, and :math:`\\Lambda` is the emissivity profile.
This class then applies a sequence of conversion factors to obtain a profile
appropriate for observable sky counts:
.. math::
C_{\\rm sky}(r)
=
C(r)
\\times (\\mathrm{Mpc \\to cm})
\\times a^4
\\times \\frac{1}{4\\pi}.
These factors account for:
- **Line-of-sight unit conversion**: the projected integral is performed in
comoving Mpc, but observational X-ray quantities are typically expressed
using cgs length units, so the profile is multiplied by the conversion from
Mpc to cm.
- **Cosmological surface-brightness dimming**: the observed signal is reduced
by a factor of :math:`(1+z)^{-4} = a^4`.
- **Solid-angle conversion**: the factor :math:`1/(4\pi)` converts the
isotropically emitted volumetric signal into a per-steradian sky quantity.
The resulting profile is therefore more appropriate for comparison with
observed X-ray surface brightness or sky-count measurements than the raw
`XrayCounts` profile.
Parameters
----------
xraycounts : BaseThermodynamicProfile, optional
Intrinsic X-ray counts profile to be converted into sky counts.
**kwargs
Additional keyword arguments passed to `BaseThermodynamicProfile` and,
when needed, to the internally constructed `XrayCounts` profile.
"""
def __init__(self, xraycounts = None, **kwargs):
self.XrayCounts = xraycounts
super().__init__(**kwargs)
if hasattr(self.XrayCounts, 'prof4params'):
self.prof4params = self.XrayCounts.prof4params
else:
self.prof4params = self
def _real(self, cosmo, r, M, a):
r_use = np.atleast_1d(r)
M_use = np.atleast_1d(M)
z = 1/a - 1
R = self.mass_def.get_radius(cosmo, M_use, a)/a #in comoving Mpc
#Now a series of units changes to the projected profile.
prof = self.XrayCounts.real(cosmo, r_use, M_use, a) #generate profile
prof = prof * (Mpc_to_m * m_to_cm) #Line-of-sight integral is done in Mpc, we want cm
prof = prof * a**3 #Cosmic dimming causes a 1/(1 + z)^3 factor (we use counts, not energy, so one factor is missing)
prof = prof * 1/(4*np.pi) #Converting 1/cm^3 into 1/steradians
return prof
def _projected(self, cosmo, r, M, a):
#This is the 1/(1 + z) that shows up in the tracer kernel
return super()._projected(cosmo, r, M, a) * a