Nikosarcevic ccl v3 paper plots

examples/ccl_v3_paper/classes/___init___.py

@@ -0,0 +1,6 @@
+from .gas_halo_profile import GasHaloProfile
+from .mass_fraction import MassFraction
+from .power_spectra_calculator import PowerSpectraCalculator
+from .presets import Presets
+from .profile_interpolation import ProfileInterpolation
+from .stellar_halo_profile import StellarHaloProfile

examples/ccl_v3_paper/classes/baryon_halo_model.py

@@ -0,0 +1,41 @@
+from .gas_halo_profile import GasHaloProfile
+from .mass_fraction import MassFraction
+from .stellar_halo_profile import StellarHaloProfile
+from .power_spectra_calculator import PowerSpectraCalculator
+from .presets import Presets
+from .profile_interpolation import ProfileInterpolation
+
+
+class BaryonHaloModel:
+    """
+    Class to encapsulate the baryon halo model calculations.
+    Methods are used to calculate the gas halo profile, the mass fraction, the stellar halo profile,
+    the power spectra, and the profile interpolation.
+    Parameters
+    ----------
+    cosmology : object
+        A pyCCL Cosmology object.
+    k_array : array_like
+        Array of wavenumbers.
+    scale_factor : float
+        Scale factor.
+    halo_mass_definition : ccl.halos.MassDef, optional.
+    """
+
+    def __init__(self, cosmology, k_array, scale_factor, halo_mass_definition=None):
+        # Initialize the Presets with cosmology settings and optional halo mass definition
+        self.presets = Presets(cosmology, k_array, scale_factor, halo_mass_definition)
+
+        # Initialize other classes that depend on the presets
+        self.gas_halo_profile = GasHaloProfile(self.presets)
+        self.mass_fraction = MassFraction(self.presets)
+        self.stellar_halo_profile = StellarHaloProfile(self.presets)
+        self.power_spectra_calculator = PowerSpectraCalculator(self.presets)
+        self.profile_interpolation = ProfileInterpolation(self.presets)
+
+    def __getattr__(self, name):
+        for component in (self.gas_halo_profile, self.mass_fraction, self.stellar_halo_profile,
+                          self.power_spectra_calculator, self.profile_interpolation, self.presets):
+            if hasattr(component, name):
+                return getattr(component, name)
+        raise AttributeError(f"'{type(self).__name__}' object has no attribute '{name}'")

examples/ccl_v3_paper/classes/gas_halo_profile.py

@@ -0,0 +1,180 @@
+from .mass_fraction import MassFraction
+import numpy as np
+from .presets import Presets
+import pyccl as ccl
+
+
+class GasHaloProfile(ccl.halos.HaloProfile):
+    """
+    Class to calculate the gas profile of a halo using the formula:
+    ρ_g(r) = ρ_g,0 / [(1 + u) ^ β * (1 + v^2) ^ ((7 - β) / 2)], where u = r / r_co and v = r / r_ej.
+    For more information, see Fedeli (2014), arXiv:1401.2997.
+    Methods
+    -------
+    _rs(cosmo, M, a)
+        Calculate the scaled radius of the halo.
+    _rint(Rd, r)
+        Calculate the radial integrand used for the normalization of the gas profile.
+    _norm(M, Rd)
+        Calculate the normalization constant for the gas profile density.
+    _real(cosmo, r, M, a)
+        Calculate the gas profile density at given radii for specified halo masses.
+    Attributes
+    ----------
+    cosmology : Cosmology object
+        Cosmology object from pyCCL used for calculations.
+    scale_factor : float
+        Scale factor at which to evaluate the profile.
+    halo_mass_definition : MassDef object
+        Mass definition used for the halo profile calculations.
+    parameters : dict
+        Dictionary containing the parameters for the gas profile calculation.
+    mass_frac : MassFraction object
+        MassFraction object used for gas mass fraction calculations.
+    """
+
+    def __init__(self, presets):
+        """
+        Initialize the GasHaloProfile class with the provided presets.
+        Parameters
+        ----------
+        presets : Presets instance
+            Presets object containing the cosmology and other parameters.
+        """
+        if not isinstance(presets, Presets):
+            raise TypeError("Expected a Presets instance.")
+        super().__init__(mass_def=presets.halo_mass_definition)
+        # Set the instance attributes based on the provided presets
+        self.cosmology = presets.cosmology
+        self.scale_factor = presets.scale_factor
+        self.halo_mass_definition = presets.halo_mass_definition
+        self.parameters = presets.parameters
+        self.mass_frac = MassFraction(presets)
+
+    def _rs(self, cosmo=None, M=None, a=None):
+        """
+        Calculate the scaled radius of the halo, defaulting to instance attributes if
+        certain parameters are not provided.
+
+        Parameters:
+        cosmo : Cosmology object, optional
+            The cosmology used to calculate the radius. Defaults to self.cosmology if None.
+        M : float or array_like
+            Mass of the halo. This parameter is required.
+        a : float, optional
+            Scale factor. Defaults to self.scale_factor if None.
+
+        Returns:
+        float or array_like
+            The scaled radius of the halo.
+        """
+        if cosmo is None:
+            cosmo = self.cosmology
+        if a is None:
+            a = self.scale_factor
+        if M is None:
+            raise ValueError("Mass 'M' must be provided for radius calculation.")
+
+        radius = self.halo_mass_definition.get_radius(cosmo, M, a) / a
+        return radius
+
+    def _rint(self, Rd, r):
+        """
+        Calculate the radial integrand used for the normalization of the gas profile, based on the formula:
+        ρ_g(r) = ρ_g,0 / [(1 + u) ^ β * (1 + v^2) ^ ((7 - β) / 2)],
+        where u = r / r_co and v = r / r_ej.
+
+        Parameters:
+        Rd : float
+            Characteristic radius of the halo, typically derived from halo properties.
+        r : float or array_like
+            Radial distance at which to evaluate the integrand.
+
+        Returns:
+        float or array_like
+            Value of the integrand at radius r.
+        """
+        beta = self.parameters['gas_profile']['beta']
+        r_co = 0.1 * Rd  # Define core radius as 10% of the characteristic radius
+        r_ej = 4.5 * Rd  # Define envelope radius as 450% of the characteristic radius
+
+        # Define u and v according to the provided expressions
+        u = r / r_co
+        v = r / r_ej
+
+        # Calculate the integrand for the gas profile using the defined variables
+        radial_distance = r ** 2 / ((1 + u) ** beta * (1 + v ** 2) ** ((7 - beta) / 2))
+        return radial_distance
+
+    def _norm(self, M, Rd):
+        """
+        Calculate the normalization constant for the gas profile density using trapezoidal integration.
+
+        Parameters:
+        M : float or array_like
+            Mass of the halo.
+        Rd : float or array_like
+            Characteristic radii of the halos.
+
+        Returns:
+        numpy.ndarray
+            Normalization constants for each halo mass.
+        """
+        f_gas = self.mass_frac.gas_mass_fraction(M)
+        r_grid = np.linspace(1E-3, 50, 500)  # Create a radial grid for integration
+        r_integrands = np.array([self._rint(rd, r_grid) for rd in Rd])  # Compute integrands for each Rd
+
+        # Perform trapezoidal integration over the radial grid for each characteristic radius
+        integrals = np.trapz(r_integrands, r_grid, axis=1)
+
+        norm = f_gas * M / (4 * np.pi * integrals)
+        return norm
+
+    def _real(self, cosmo=None, r=None, M=None, a=None):
+        """
+        Calculate the gas profile density at given radii for specified halo masses.
+
+        Parameters:
+        cosmo : Cosmology object, optional
+            The cosmology used for calculations. Defaults to self.cosmology if None.
+        r : float or array_like, optional
+            Radial distances to evaluate the profile. Must be provided.
+        M : float or array_like, optional
+            Mass of the halo. Must be provided.
+        a : float, optional
+            Scale factor. Defaults to self.scale_factor if None.
+
+        Returns:
+        numpy.ndarray
+            Gas profile densities for each (r, M) combination.
+        """
+        if cosmo is None:
+            cosmo = self.cosmology
+        if a is None:
+            a = self.scale_factor
+        if r is None or M is None:
+            raise ValueError("Both radial distances 'r' and masses 'M' must be provided.")
+
+        r_use = np.atleast_1d(r)
+        M_use = np.atleast_1d(M)
+
+        r_delta = self._rs(cosmo, M_use, a)
+        r_co = 0.1 * r_delta
+        r_ej = 4.5 * r_delta
+
+        beta = self.parameters['gas_profile']['beta']
+        norm = self._norm(M_use, r_delta)
+
+        # Define u and v for readability
+        u = r_use[None, :] / r_co[:, None]
+        v = r_use[None, :] / r_ej[:, None]
+
+        # Compute the profile using defined u and v
+        prof = norm[:, None] / ((1 + u) ** beta * (1 + v ** 2) ** ((7 - beta) / 2))
+
+        # Reshape output to match the dimensions of the inputs
+        if np.isscalar(r):
+            prof = np.squeeze(prof, axis=-1)
+        if np.isscalar(M):
+            prof = np.squeeze(prof, axis=0)
+        return prof

examples/ccl_v3_paper/classes/mass_fraction.py

@@ -0,0 +1,101 @@
+import numpy as np
+from .presets import Presets
+from scipy import integrate
+from scipy.special import erf
+
+
+class MassFraction:
+
+    def __init__(self, presets):
+        if not isinstance(presets, Presets):
+            raise TypeError("Expected a Presets instance.")
+
+        self.cosmology = presets.cosmology
+        self.scale_factor = presets.scale_factor
+        self.halo_mass_definition = presets.halo_mass_definition
+        self.parameters = presets.parameters
+        self.densities = presets.densities
+        self.mass_ranges = presets.mass_ranges
+        self.halo_model_quantities = presets.halo_model_quantities
+
+    def gas_mass_fraction(self, mass):
+        m_0 = self.parameters['gas_mass_fraction']['m_0']
+        sigma = self.parameters['gas_mass_fraction']['sigma']
+        omega_b = self.cosmology['Omega_b']
+        omega_m = self.cosmology['Omega_m']
+
+        omega_ratio = omega_b / omega_m
+        mass_ratio = np.atleast_1d(mass) / m_0
+        gas_mass_frac = np.zeros_like(mass_ratio)
+
+        # Apply the original conditional logic
+        valid_indices = mass_ratio >= 1  # Only apply erf calculation to valid mass ratios
+        gas_mass_frac[valid_indices] = omega_ratio * erf((np.log10(mass_ratio[valid_indices])) / sigma)
+
+        # Return scalar if input was scalar
+        if np.isscalar(mass):
+            return gas_mass_frac[0]
+        return gas_mass_frac
+
+    def stellar_mass_fraction(self, mass):
+        # Retrieve constants and necessary data
+        rho_star = self.densities['stars']
+        m_0 = self.parameters['stellar_mass_fraction']['m_0']
+        sigma = self.parameters['stellar_mass_fraction']['sigma']
+        mass_min = self.mass_ranges['stars']['min']
+        mass_max = self.mass_ranges['stars']['max']
+
+        # Use the halo mass function from the pre-defined quantities in the class
+        halo_mass_function = self.halo_model_quantities['halo_mass_function']
+
+        # Define the integrand function using the halo mass function
+        def smf_integrand(m):
+            mass_dex = (m * np.log(10))
+            mf = halo_mass_function(self.cosmology, m, self.scale_factor) / mass_dex
+            return m * np.exp(-(np.log10(m / m_0)) ** 2 / (2 * sigma ** 2)) * mf
+
+        # Integration using scipy.integrate.quad for normalization
+        integral, error = integrate.quad(smf_integrand, mass_min, mass_max, epsabs=0, epsrel=1E-3, limit=5000)
+
+        # Compute A
+        A = rho_star / integral
+
+        # Compute f_star for the provided mass array
+        star_mass_frac = A * np.exp(-(np.log10(mass / m_0)) ** 2 / (2 * sigma ** 2))
+
+        return star_mass_frac
+
+    def dark_matter_mass_fraction(self, mass):
+        omega_b = self.cosmology['Omega_b']
+        omega_m = self.cosmology['Omega_m']
+
+        omega_ratio = omega_b / omega_m
+        # My mass frac is 1 - omega_ratio
+        mass_ratio = np.atleast_1d(mass)
+        dark_matter_frac = np.ones_like(mass_ratio) - omega_ratio
+
+        # Return scalar if input was scalar
+        if np.isscalar(mass):
+            return dark_matter_frac[0]
+        return dark_matter_frac
+
+    def mass_fraction_dict(self, mass):
+        """
+        Compute the mass fraction of each component for a given mass.
+        Parameters
+        ----------
+            mass : float or array_like
+            The mass of the halo
+
+        Returns
+        -------
+            A dictionary containing the mass fractions of each component.
+            Components are: stars, gas, dark_matter.
+        """
+        mass_frac_dict = {
+            "dark_matter": self.dark_matter_mass_fraction(mass),
+            "gas": self.gas_mass_fraction(mass),
+            "stars": self.stellar_mass_fraction(mass),
+        }
+
+        return mass_frac_dict