clusters/halo.py: Significant cleanups

Several methods/functions have been migrated into "helper.py", while
other methods/functions are obsolete.

Signed-off-by: Aaron LI <aly@aaronly.me>

1 files changed, 2 insertions, 551 deletions

diff --git a/fg21sim/extragalactic/clusters/halo.py b/fg21sim/extragalactic/clusters/halo.py
index c78c2c3..958dd6c 100644
--- a/fg21sim/extragalactic/clusters/halo.py
+++ b/fg21sim/extragalactic/clusters/halo.py
@@ -20,11 +20,7 @@ References
 import logging
 
 import numpy as np
-import scipy.interpolate
-import scipy.integrate
-import scipy.optimize
 
-from .formation import ClusterFormation
 from .solver import FokkerPlanckSolver
 from ...utils import cosmo
 from ...utils.units import (Units as AU,
@@ -54,101 +50,10 @@ class RadioHalo:
         Unit: [Msun]
     z0 : float
         Redshift from where to simulate former merging history.
-    configs : `ConfigManager`
-        A `ConfigManager` instance containing default and user configurations.
-        For more details, see the example configuration specifications.
-
-    Attributes
-    ----------
-    mec : float
-        Unit for electron momentum (p): mec = m_e * c, p = gamma * mec,
-        therefore value of p is the Lorentz factor.
-    mtree : `~MergerTree`
-        Merging history of this cluster.
     """
-    # Unit for electron momentum (p), thus its value is the Lorentz factor
-    mec = AU.mec  # [g cm / s]
-    # Merger tree (i.e., merging history) of this cluster
-    mtree = None
-
-    def __init__(self, M0, z0, configs):
-        self.M0 = M0  # [Msun]
+    def __init__(self, M0, z0):
+        self.M0 = M0
         self.z0 = z0
-        self.configs = configs
-        self._set_configs()
-
-    def _set_configs(self):
-        """
-        Set up the necessary class attributes according to the configs.
-        """
-        comp = "extragalactic/halos"
-        self.zmax = self.configs.getn(comp+"/zmax")
-        self.zbinsize = self.configs.getn(comp+"/zbinsize")
-        # Mass threshold of the sub-cluster for a significant merger
-        self.merger_mass_th = self.configs.getn(comp+"/merger_mass_th")
-        self.radius = self.configs.getn(comp+"/radius")
-        self.b_mean = self.configs.getn(comp+"/b_mean")
-        self.b_index = self.configs.getn(comp+"/b_index")
-        self.eta_t = self.configs.getn(comp+"/eta_t")
-        self.eta_e = self.configs.getn(comp+"/eta_e")
-        self.pmin = self.configs.getn(comp+"/pmin")
-        self.pmax = self.configs.getn(comp+"/pmax")
-        self.pgrid_num = self.configs.getn(comp+"/pgrid_num")
-        self.buffer_np = self.configs.getn(comp+"/buffer_np")
-        self.time_step = self.configs.getn(comp+"/time_step")
-        self.injection_index = self.configs.getn(comp+"/injection_index")
-        logger.info("Loaded and set up configurations")
-
-    @property
-    def zgrid(self):
-        """
-        Redshift grid between 0 and ``self.zmax``.
-        """
-        return np.arange(start=0.0, stop=self.zmax+self.zbinsize,
-                         step=self.zbinsize)
-
-    @property
-    def pgrid(self):
-        """
-        Electron momentum values of the adopted logarithmic grid
-        used for solving the Fokker-Planck equation.
-        """
-        grid = np.logspace(np.log10(self.pmin), np.log10(self.pmax),
-                           num=self.pgrid_num)
-        return grid
-
-    @property
-    def magnetic_field(self):
-        """
-        Calculate the mean magnetic field of this cluster from the
-        scaling relation between magnetic field and cluster mass.
-
-        Unit: [uG]
-
-        Reference: Ref.[3],Eq.(1)
-        """
-        M_mean = 1.6e15  # [Msun]
-        B = self.b_mean * (self.M0/M_mean) ** self.b_index
-        return B
-
-    def simulate_mergertree(self):
-        """
-        Simulate the merging history of the cluster using the extended
-        Press-Schechter formalism.
-
-        Attributes
-        ----------
-        mtree : `~MergerTree`
-            Generated merger tree of this cluster.
-
-        Returns
-        -------
-        mtree : `~MergerTree`
-            Generated merger tree of this cluster.
-        """
-        self.formation = ClusterFormation(self.M0, self.z0, self.configs)
-        self.mtree = self.formation.simulate_mergertree()
-        return self.mtree
 
     def calc_electron_spectrum(self, zbegin=None, zend=None, n0_e=None):
         """
@@ -204,261 +109,6 @@ class RadioHalo:
         n_e = fpsolver.solve(u0=n0_e, tstart=tstart, tstop=tstop)
         return (p, n_e)
 
-    def kT_mass(self, mass):
-        """
-        Estimate the cluster ICM temperature from its mass by assuming
-        an (observed) temperature-mass relation.
-
-        TODO: upgrade this M-T relation.
-
-        Parameters
-        ----------
-        mass : float
-            Mass (unit: Msun) of the cluster
-
-        Returns
-        -------
-        kT : float
-            Temperature of the ICM (unit: keV)
-
-        References
-        ----------
-        [1] Nevalainen et al. 2000, ApJ, 532, 694
-            Ettori et al, 2013, Space Science Review, 177, 119-154
-            NOTE: H0 = 50 * h50 [km/s/Mpc]
-        """
-        kT = 10 * (mass/1.23e15) ** (1/1.79)  # [keV]
-        return kT
-
-    def _radius_virial(self, mass, z=0.0):
-        """
-        Calculate the virial radius of a cluster.
-
-        Parameters
-        ----------
-        mass : float
-            Mass (unit: Msun) of the cluster
-        z : float
-            Redshift
-
-        Returns
-        -------
-        Rvir : float
-            Virial radius (unit: kpc) of the cluster at given redshift
-        """
-        Dc = cosmo.overdensity_virial(z)
-        rho = cosmo.rho_crit(z)  # [g/cm^3]
-        R_vir = (3*mass*AUC.Msun2g / (4*np.pi * Dc * rho))**(1/3)  # [cm]
-        R_vir *= AUC.cm2kpc  # [kpc]
-        return R_vir
-
-    def _radius_stripping(self, mass, M_main, z):
-        """
-        Calculate the stripping radius of the sub-cluster at which
-        equipartition between static and ram pressure is established,
-        and the stripping is efficient outside this stripping radius.
-
-        Note that the value of the stripping radius obtained would
-        give the *mean value* of the actual stripping radius during
-        a merger.
-
-        Parameters
-        ----------
-        mass : float
-            The mass (unit: Msun) of the sub-cluster.
-        M_main : float
-            The mass (unit: Msun) of the main cluster.
-        z : float
-            Redshift
-
-        Returns
-        -------
-        rs : float
-            The stripping radius of the sub-cluster.
-            Unit: kpc
-
-        References
-        ----------
-        [1] Cassano & Brunetti 2005, MNRAS, 357, 1313
-            http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C
-            Eq.(11)
-        """
-        vi = self._velocity_impact(M_main, mass, z) * 1e5  # [cm/s]
-        kT = self.kT_mass(mass) * AUC.keV2erg  # [erg]
-        coef = kT / (AC.mu * AC.u * vi**2)  # dimensionless
-        rho_avg = self._density_average(M_main, z)  # [g/cm^3]
-
-        def equation(r):
-            return coef * self.density_profile(r, mass, z) / rho_avg - 1
-
-        r_vir = self._radius_virial(mass, z)  # [kpc]
-        rs = scipy.optimize.brentq(equation, a=0, b=r_vir)  # [kpc]
-        return rs
-
-    def _density_average(self, mass, z=0.0):
-        """
-        Average density of the cluster ICM.
-
-        Returns
-        -------
-        rho : float
-            Average ICM density (unit: g/cm^3)
-
-        References
-        ----------
-        [1] Cassano & Brunetti 2005, MNRAS, 357, 1313
-            http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C
-            Eq.(12)
-        """
-        f_baryon = cosmo.Ob0 / cosmo.Om0
-        Rv = self._radius_virial(mass, z) * AUC.kpc2cm  # [cm]
-        V = (4*np.pi / 3) * Rv**3  # [cm^3]
-        rho = f_baryon * mass*AUC.Msun2g / V  # [g/cm^3]
-        return rho
-
-    def density_profile(self, r, mass, z):
-        """
-        ICM (baryon) density profile, assuming the beta model.
-
-        Parameters
-        ----------
-        r : float
-            Radius (unit: kpc) where to calculate the density
-        mass : float
-            Cluster mass (unit: Msun)
-        z : float
-            Redshift
-
-        Returns
-        -------
-        rho_r : float
-            Density at the specified radius (unit: g/cm^3)
-
-        References
-        ----------
-        [1] Cassano & Brunetti 2005, MNRAS, 357, 1313
-            http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C
-            Eq.(13)
-        """
-        f_baryon = cosmo.Ob0 / cosmo.Om0
-        M_ICM = mass * f_baryon * AUC.Msun2g  # [g]
-        r *= AUC.kpc2cm  # [cm]
-        Rv = self._radius_virial(mass, z) * AUC.kpc2cm  # [cm]
-        rc = self._beta_rc(Rv)
-        beta = self._beta_beta()
-        norm = self._beta_norm(M_ICM, beta, rc, Rv)  # [g/cm^3]
-        rho_r = norm * (1 + (r/rc)**2) ** (-3*beta/2)  # [g/cm^3]
-        return rho_r
-
-    @staticmethod
-    def _beta_rc(r_vir):
-        """
-        Core radius of the beta model for the ICM density profile.
-
-        TODO: upgrade this!
-        """
-        return 0.1*r_vir
-
-    @staticmethod
-    def _beta_beta():
-        """
-        Beta value of the beta model for the ICM density profile.
-
-        TODO: upgrade this!
-        """
-        return 0.8
-
-    @staticmethod
-    def _beta_norm(mass, beta, rc, r_vir):
-        """
-        Calculate the normalization of the beta model for the ICM
-        density profile.
-
-        Parameters
-        ----------
-        mass : float
-            The mass (unit: g) of ICM
-        beta : float
-            Beta value of the assumed beta profile
-        rc : float
-            Core radius (unit: cm) of the assumed beta profile
-        r_vir : float
-            The virial radius (unit: cm) of the cluster
-
-        Returns
-        -------
-        norm : float
-            Normalization of the beta model (unit: g/cm^3)
-
-        References
-        ----------
-        [1] Cassano & Brunetti 2005, MNRAS, 357, 1313
-            http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C
-            Eq.(14)
-        """
-        integration = scipy.integrate.quad(
-            lambda r: r*r * (1+(r/rc)**2) ** (-3*beta/2),
-            0, r_vir)[0]
-        norm = mass / (4*np.pi * integration)  # [g/cm^3]
-        return norm
-
-    def _velocity_impact(self, M_main, M_sub, z=0.0):
-        """
-        Calculate the relative impact velocity between the two merging
-        clusters when they are at a distance of virial radius.
-
-        Parameters
-        ----------
-        M_main : float
-            Mass of the main cluster (unit: Msun)
-        M_sub : float
-            Mass of the sub cluster (unit: Msun)
-        z : float
-            Redshift
-
-        Returns
-        -------
-        vi : float
-            Relative impact velocity (unit: km/s)
-
-        References
-        ----------
-        [1] Cassano & Brunetti 2005, MNRAS, 357, 1313
-            http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C
-            Eq.(9)
-        """
-        eta_v = 4 * (1 + M_main/M_sub) ** (1/3)
-        R_vir = self._radius_virial(M_main, z) * AUC.kpc2cm  # [cm]
-        vi = np.sqrt(2*AC.G * (1-1/eta_v) *
-                     (M_main+M_sub)*AUC.Msun2g / R_vir)  # [cm/s]
-        vi /= AUC.km2cm  # [km/s]
-        return vi
-
-    def _time_crossing(self, M_main, M_sub, z):
-        """
-        Calculate the crossing time of the sub-cluster during a merger.
-
-        Parameters
-        ----------
-        M_main : float
-            Mass of the main cluster (unit: Msun)
-        M_sub : float
-            Mass of the sub cluster (unit: Msun)
-        z : float
-            Redshift where the merger occurs.
-
-        Returns
-        -------
-        time : float
-            Crossing time (unit: Gyr)
-        """
-        R_vir = self._radius_virial(M_main, z)  # [kpc]
-        vi = self._velocity_impact(M_main, M_sub, z)  # [km/s]
-        # Unit conversion coefficient: [s kpc/km] => [Gyr]
-        uconv = AUC.kpc2km * AUC.s2Gyr
-        time = uconv * R_vir / vi  # [Gyr]
-        return time
-
     def _z_end(self, z_begin, time):
         """
         Calculate the ending redshift from ``z_begin`` after elapsing
@@ -479,165 +129,6 @@ class RadioHalo:
             z_end = cosmo.redshift(t_end)
         return z_end
 
-    @property
-    def merger_events(self):
-        """
-        Trace only the main cluster, and filter out the significant
-        merger events.
-
-        Returns
-        -------
-        mevents : list[dict]
-            List of dictionaries that records all the merger events
-            of the main cluster.
-            NOTE:
-            The merger events are ordered by increasing redshifts.
-        """
-        events = []
-        tree = self.mtree
-        while tree:
-            if (tree.main and tree.sub and
-                    tree.sub.data["mass"] >= self.merger_mass_th and
-                    tree.main.data["z"] <= self.zmax):
-                events.append({
-                    "z": tree.main.data["z"],
-                    "age": tree.main.data["age"],
-                    "M_main": tree.main.data["mass"],
-                    "M_sub": tree.sub.data["mass"],
-                })
-            tree = tree.main
-        return events
-
-    def _coef_acceleration(self, z):
-        """
-        Calculate the electron-acceleration coefficient at arbitrary
-        redshift.
-
-        Parameters
-        ----------
-        z : float
-            Redshift where to calculate the acceleration coefficient.
-
-        Returns
-        -------
-        chi : float
-            The calculated electron-acceleration coefficient.
-            (unit: Gyr^-1)
-
-        XXX/NOTE
-        --------
-        This coefficient may be very small and even zero, then the
-        diffusion coefficient of the Fokker-Planck equation is thus
-        very small and even zero, which cause problems for calculating
-        some quantities (e.g., w(x), C(x)) and wrong/invalid results.
-        To avoid these problems, force the minimal value of this
-        coefficient to be 1/(10*t0), which t0 is the present-day age
-        of the universe.
-        """
-        zgrid, chigrid = self._chi_data
-        # XXX: force a minimal value instead of zero or too small
-        chi_min = 1 / (10 * cosmo.age0)
-
-        try:
-            zi = np.where(z < zgrid)[0][0]
-            slope = (chigrid[zi]-chigrid[zi-1]) / (zgrid[zi]-zgrid[zi-1])
-            chi = (z-zgrid[zi]) * slope + chigrid[zi]
-        except IndexError:
-            # Given z >= zgrid[-1]
-            chi = 0.0
-
-        if chi > chi_min:
-            return chi
-        else:
-            return chi_min
-
-    @property
-    def _chi_data(self):
-        """
-        Returns
-        -------
-        zgrid : 1D `~numpy.ndarray`
-            Redshift positions where the chi values are calculated.
-            Same as the attribute ``self.zgrid``
-        chigrid : 1D `~numpy.ndarray`
-            Values of acceleration coefficients at each redshift.
-        """
-        if hasattr(self, "_chi_data_"):
-            zgrid, chigrid = self._chi_data_
-        else:
-            # Order the merger events by decreasing redshifts
-            mevents = list(reversed(self.merger_events))
-            chi_z = [self._chi_at_zidx(zidx, mevents)
-                     for zidx in range(len(mevents))]
-            zgrid = self.zgrid
-            chigrid = np.zeros(zgrid.shape)
-            for (chi_, zbegin, zend) in chi_z:
-                # NOTE: zbegin > zend
-                mask = (zgrid <= zbegin) & (zgrid >= zend)
-                chigrid[mask] += chi_
-            self._chi_data_ = (zgrid, chigrid)
-        return (zgrid, chigrid)
-
-    def _chi_at_zidx(self, zidx, mevents):
-        """
-        Calculate electron-acceleration coefficient at the specified
-        merger event which is specified with a redshift index.
-
-        Parameters
-        ----------
-        zidx : int
-            Index of the redshift where to calculate the coefficient.
-        mevents : list[dict]
-            List of dictionaries that records all the merger events
-            of the main cluster.
-            NOTE:
-            The merger events should be ordered by increasing time
-            (or decreasing redshifts).
-
-        Returns
-        -------
-        chi : float
-            The calculated electron-acceleration coefficient.
-            Unit: [Gyr^-1]
-        zbegin : float
-            The redshift when the merger begins
-        zend : float
-            The redshift when the merger ends
-            NOTE: zbegin > zend
-
-        References: Ref.[1],Eq.(40)
-        """
-        redshifts = np.array([ev["z"] for ev in mevents])
-        zbegin = mevents[zidx]["z"]
-        M_main = mevents[zidx]["M_main"]
-        M_sub = mevents[zidx]["M_sub"]
-        t_crossing = self._time_crossing(M_main, M_sub, zbegin)  # [Gyr]
-        zend = self._z_end(zbegin, t_crossing)
-        try:
-            zend_idx = np.where(redshifts < zend)[0][0]
-        except IndexError:
-            # Specified redshift already the last/smallest one
-            zend_idx = zidx + 1
-        #
-        coef = 2.23e-16 * self.eta_t / (self.radius/500)**3  # [s^-1]
-        coef *= AUC.Gyr2s  # [Gyr^-1]
-        chi = 0.0
-        for ev in mevents[zidx:zend_idx]:
-            M_main = ev["M_main"]
-            M_sub = ev["M_sub"]
-            z = ev["z"]
-            R_vir = self._radius_virial(M_main, z)
-            rs = self._radius_stripping(M_sub, M_main, z)
-            kT = self.kT_mass(M_main)
-            term1 = ((M_main+M_sub)/2e15 * (2.6e3/R_vir)) ** (3/2)
-            term2 = (rs/500)**2 / np.sqrt(kT/7)
-            if rs <= self.radius:
-                term3 = 1.0
-            else:
-                term3 = (self.radius/rs) ** 2
-            chi += coef * term1 * term2 * term3
-        return (chi, zbegin, zend)
-
     def fp_injection(self, p, t=None):
         """
         Electron injection term for the Fokker-Planck equation.
@@ -776,47 +267,7 @@ class RadioHalo:
                 ((self.magnetic_field/3.2)**2 + (1+z)**4))
         return dpdt
 
-    @property
-    def e_thermal(self):
-        """
-        Calculate the thermal energy density of the ICM.
-
-        Returns
-        -------
-        e_th : float
-            Energy density of the ICM (unit: erg/cm^3)
         """
-        mass = self.M0
-        f_baryon = cosmo.Ob0 / cosmo.Om0
-        kT = self.kT_mass(mass)  # [keV]
-        N = mass * AUC.Msun2g * f_baryon / (AC.mu * AC.u)
-        E_th = kT*AUC.keV2erg * N  # [erg]
-        Rv = self._radius_virial(mass, self.z0) * AUC.kpc2cm  # [cm]
-        V = (4*np.pi / 3) * Rv**3  # [cm^3]
-        e_th = E_th / V  # [erg/cm^3]
-        return e_th
-
-    def _n_thermal(self, mass, z=0.0):
-        """
-        Calculate the number density of the ICM thermal plasma.
 
-        Parameters
         ----------
-        mass : float
-            Mass of the cluster at the redshift
-            Unit: [Msun]
-        z : float
-            Redshift
-
-        Returns
-        -------
-        n_th : float
-            Number density of the ICM
-            Unit: [cm^-3]
         """
-        f_baryon = cosmo.Ob0 / cosmo.Om0
-        N = mass * AUC.Msun2g * f_baryon / (AC.mu * AC.u)
-        Rv = self._radius_virial(mass, z) * AUC.kpc2cm  # [cm]
-        V = (4*np.pi / 3) * Rv**3  # [cm^3]
-        n_th = N / V  # [cm^-3]
-        return n_th