clusters/halo: split emissivity/power/flux calculations into HaloEmission

Add emission.HaloEmission class to calculate halo emissivity, power, flux,
brightness etc.

2 files changed, 200 insertions, 175 deletions

diff --git a/fg21sim/extragalactic/clusters/emission.py b/fg21sim/extragalactic/clusters/emission.py
index ddecab3..a798470 100644
--- a/fg21sim/extragalactic/clusters/emission.py
+++ b/fg21sim/extragalactic/clusters/emission.py
@@ -31,7 +31,11 @@ import numpy as np
 import scipy.special
 from scipy import integrate, interpolate
 
-from ...utils.units import (Units as AU, Constants as AC)
+from ...share import COSMO
+from ...utils.convert import Fnu_to_Tb
+from ...utils.units import (Units as AU,
+                            UnitConversions as AUC,
+                            Constants as AC)
 
 
 logger = logging.getLogger(__name__)
@@ -223,3 +227,198 @@ class SynchrotronEmission:
             return syncem[0]
         else:
             return syncem
+
+
+class HaloEmission:
+    """
+    Calculate the synchrotron emission of a (giant) radio halo.
+
+    Parameters
+    ----------
+    gamma : 1D `~numpy.ndarray`
+        The Lorentz factors γ of the electron spectrum.
+    n_e : 1D `~numpy.ndarray`
+        The electron spectrum (w.r.t. Lorentz factors γ).
+        Unit: [cm^-3]
+    B : float
+        The magnetic field strength.
+        Unit: [uG]
+    radius : float, optional
+        The radio halo radius.
+        Required to calculate the power.
+        Unit: [kpc]
+    redshift : float, optional
+        The redshift to the radio halo.
+        Required to calculate the flux, which also requires ``radius``.
+    """
+    def __init__(self, gamma, n_e, B, radius=None, redshift=None):
+        self.gamma = np.asarray(gamma)
+        self.n_e = np.asarray(n_e)
+        self.B = B
+        self.radius = radius
+        self.redshift = redshift
+
+    @property
+    def angular_radius(self):
+        """
+        The angular radius of the radio halo.
+        Unit: [arcsec]
+        """
+        if self.redshift is None:
+            raise RuntimeError("parameter 'redshift' is required")
+        if self.radius is None:
+            raise RuntimeError("parameter 'radius' is required")
+
+        DA = COSMO.DA(self.redshift) * 1e3  # [Mpc] -> [kpc]
+        theta = self.radius / DA  # [rad]
+        return theta * AUC.rad2arcsec
+
+    @property
+    def volume(self):
+        """
+        The halo volume.
+        Unit: [kpc^3]
+        """
+        if self.radius is None:
+            raise RuntimeError("parameter 'radius' is required")
+
+        return (4*np.pi/3) * self.radius**3
+
+    def calc_emissivity(self, frequencies):
+        """
+        Calculate the synchrotron emissivity for the derived electron
+        spectrum.
+
+        Parameters
+        ----------
+        frequencies : float, or 1D `~numpy.ndarray`
+            The frequencies where to calculate the synchrotron emissivity.
+            Unit: [MHz]
+
+        Returns
+        -------
+        emissivity : float, or 1D `~numpy.ndarray`
+            The calculated synchrotron emissivity at each specified
+            frequency.
+            Unit: [erg/s/cm^3/Hz]
+        """
+        syncem = SynchrotronEmission(gamma=self.gamma, n_e=self.n_e, B=self.B)
+        emissivity = syncem.emissivity(frequencies)
+        return emissivity
+
+    def calc_power(self, frequencies, emissivity=None):
+        """
+        Calculate the halo synchrotron power (i.e., power *emitted* per
+        unit frequency) by assuming the emissivity is uniform throughout
+        the halo volume.
+
+        NOTE
+        ----
+        The calculated power (a.k.a. spectral luminosity) is in units of
+        [W/Hz] which is common in radio astronomy, instead of [erg/s/Hz].
+            1 [W] = 1e7 [erg/s]
+
+        Parameters
+        ----------
+        frequencies : float, or 1D `~numpy.ndarray`
+            The frequencies where to calculate the synchrotron power.
+            Unit: [MHz]
+        emissivity : float, or 1D `~numpy.ndarray`, optional
+            The synchrotron emissivity at the input frequencies.
+            If not provided, then invoke above ``calc_emissivity()``
+            method to calculate them.
+            Unit: [erg/s/cm^3/Hz]
+
+        Returns
+        -------
+        power : float, or 1D `~numpy.ndarray`
+            The calculated synchrotron power at each input frequency.
+            Unit: [W/Hz]
+        """
+        frequencies = np.asarray(frequencies)
+        if emissivity is None:
+            emissivity = self.calc_emissivity(frequencies=frequencies)
+        else:
+            emissivity = np.asarray(emissivity)
+        power = emissivity * (self.volume * AUC.kpc2cm**3)  # [erg/s/Hz]
+        power *= 1e-7  # [erg/s/Hz] -> [W/Hz]
+        return power
+
+    def calc_flux(self, frequencies):
+        """
+        Calculate the synchrotron flux density (i.e., power *observed*
+        per unit frequency) of the halo, with k-correction considered.
+
+        NOTE
+        ----
+        The *k-correction* must be applied to the flux density (Sν) or
+        specific luminosity (Lν) because the redshifted object is emitting
+        flux in a different band than that in which you are observing.
+        And the k-correction depends on the spectrum of the object in
+        question.  For any other spectrum (i.e., vLv != const.), the flux
+        density Sv is related to the specific luminosity Lv by:
+            Sv = (1+z) L_v(1+z) / (4π DL^2),
+        where
+        * L_v(1+z) is the specific luminosity emitting at frequency v(1+z),
+        * DL is the luminosity distance to the object at redshift z.
+
+        Reference: Ref.[hogg1999],Eq.(22)
+
+        Returns
+        -------
+        flux : float, or 1D `~numpy.ndarray`
+            The calculated flux density w.r.t. each input frequency.
+            Unit: [Jy] = 1e-23 [erg/s/cm^2/Hz] = 1e-26 [W/m^2/Hz]
+        """
+        if self.redshift is None:
+            raise RuntimeError("parameter 'redshift' is required")
+
+        freqz = np.asarray(frequencies) * (1+self.redshift)
+        power = self.calc_power(freqz)  # [W/Hz]
+        DL = COSMO.DL(self.redshift) * AUC.Mpc2m  # [m]
+        flux = 1e26 * (1+self.redshift) * power / (4*np.pi * DL*DL)  # [Jy]
+        return flux
+
+    def calc_brightness_mean(self, frequencies, flux=None, pixelsize=None):
+        """
+        Calculate the mean surface brightness (power observed per unit
+        frequency and per unit solid angle) expressed in *brightness
+        temperature* at the specified frequencies.
+
+        NOTE
+        ----
+        If the solid angle that the object extends is smaller than the
+        specified pixel area, then is is assumed to have size of 1 pixel.
+
+        Parameters
+        ----------
+        frequencies : float, or 1D `~numpy.ndarray`
+            The frequencies where to calculate the mean brightness temperature
+            Unit: [MHz]
+        flux : float, or 1D `~numpy.ndarray`, optional
+            The flux density w.r.t. each input frequency.
+            Unit: [Jy]
+        pixelsize : float, optional
+            The pixel size of the output simulated sky image.
+            If not provided, then invoke above ``calc_flux()`` method to
+            calculate them.
+            Unit: [arcsec]
+
+        Returns
+        -------
+        Tb : float, or 1D `~numpy.ndarray`
+            The mean brightness temperature at each frequency.
+            Unit: [K] <-> [Jy/pixel]
+        """
+        frequencies = np.asarray(frequencies)
+        if flux is None:
+            flux = self.calc_flux(frequencies=frequencies)  # [Jy]
+        else:
+            flux = np.asarray(flux)
+        omega = np.pi * self.angular_radius**2  # [arcsec^2]
+        if pixelsize and (omega < pixelsize**2):
+            omega = pixelsize ** 2  # [arcsec^2]
+            logger.warning("Halo size < 1 pixel; force to be 1 pixel!")
+
+        Tb = Fnu_to_Tb(flux, omega, frequencies)  # [K]
+        return Tb
diff --git a/fg21sim/extragalactic/clusters/halo.py b/fg21sim/extragalactic/clusters/halo.py
index 4e5f30f..29262b8 100644
--- a/fg21sim/extragalactic/clusters/halo.py
+++ b/fg21sim/extragalactic/clusters/halo.py
@@ -62,12 +62,10 @@ import numpy as np
 
 from . import helper
 from .solver import FokkerPlanckSolver
-from .emission import SynchrotronEmission
 from ...share import CONFIGS, COSMO
 from ...utils.units import (Units as AU,
                             UnitConversions as AUC,
                             Constants as AC)
-from ...utils.convert import Fnu_to_Tb
 
 
 logger = logging.getLogger(__name__)
@@ -494,178 +492,6 @@ class RadioHalo:
         else:
             raise ValueError("given electron spectrum has wrong shape!")
 
-    def calc_emissivity(self, frequencies, n_e=None, gamma=None, B=None):
-        """
-        Calculate the synchrotron emissivity for the derived electron
-        spectrum.
-
-        Parameters
-        ----------
-        frequencies : float, or 1D `~numpy.ndarray`
-            The frequencies where to calculate the synchrotron emissivity.
-            Unit: [MHz]
-        n_e : 1D `~numpy.ndarray`, optional
-            The electron spectrum (w.r.t. Lorentz factors γ).
-            If not provided, then use the cached ``self.electron_spec``
-            that was solved at above.
-            Unit: [cm^-3]
-        gamma : 1D `~numpy.ndarray`, optional
-            The Lorentz factors γ of the electron spectrum.
-            If not provided, then use ``self.gamma``.
-        B : float, optional
-            The magnetic field strength.
-            If not provided, then use ``self.B_obs``.
-            Unit: [uG]
-
-        Returns
-        -------
-        emissivity : float, or 1D `~numpy.ndarray`
-            The calculated synchrotron emissivity at each specified
-            frequency.
-            Unit: [erg/s/cm^3/Hz]
-        """
-        if n_e is None:
-            n_e = self.electron_spec
-        if gamma is None:
-            gamma = self.gamma
-        if B is None:
-            B = self.B_obs
-        syncem = SynchrotronEmission(gamma=gamma, n_e=n_e, B=B)
-        emissivity = syncem.emissivity(frequencies)
-        return emissivity
-
-    def calc_power(self, frequencies, emissivity=None, **kwargs):
-        """
-        Calculate the halo synchrotron power (i.e., power *emitted* per
-        unit frequency) by assuming the emissivity is uniform throughout
-        the halo volume.
-
-        NOTE
-        ----
-        The calculated power (a.k.a. spectral luminosity) is in units of
-        [W/Hz] which is common in radio astronomy, instead of [erg/s/Hz].
-            1 [W] = 1e7 [erg/s]
-
-        Parameters
-        ----------
-        frequencies : float, or 1D `~numpy.ndarray`
-            The frequencies where to calculate the synchrotron power.
-            Unit: [MHz]
-        emissivity : float, or 1D `~numpy.ndarray`, optional
-            The synchrotron emissivity at the input frequencies.
-            If not provided, then invoke above ``calc_emissivity()``
-            method to calculate them.
-            Unit: [erg/s/cm^3/Hz]
-        **kwargs : optional arguments, i.e., ``n_e``, ``gamma``, and ``B``.
-
-        Returns
-        -------
-        power : float, or 1D `~numpy.ndarray`
-            The calculated synchrotron power at each input frequency.
-            Unit: [W/Hz]
-        """
-        frequencies = np.asarray(frequencies)
-        if emissivity is None:
-            emissivity = self.calc_emissivity(frequencies=frequencies,
-                                              **kwargs)
-        else:
-            emissivity = np.asarray(emissivity)
-            if emissivity.shape != frequencies.shape:
-                raise ValueError("input 'frequencies' and 'emissivity' "
-                                 "do not match")
-        power = emissivity * (self.volume * AUC.kpc2cm**3)  # [erg/s/Hz]
-        power *= 1e-7  # [erg/s/Hz] -> [W/Hz]
-        return power
-
-    def calc_flux(self, frequencies, **kwargs):
-        """
-        Calculate the synchrotron flux density (i.e., power *observed*
-        per unit frequency) of the halo, with k-correction considered.
-
-        NOTE
-        ----
-        The *k-correction* must be applied to the flux density (Sν) or
-        specific luminosity (Lν) because the redshifted object is emitting
-        flux in a different band than that in which you are observing.
-        And the k-correction depends on the spectrum of the object in
-        question.  For any other spectrum (i.e., vLv != const.), the flux
-        density Sv is related to the specific luminosity Lv by:
-            Sv = (1+z) L_v(1+z) / (4π DL^2),
-        where
-        * L_v(1+z) is the specific luminosity emitting at frequency v(1+z),
-        * DL is the luminosity distance to the object at redshift z.
-
-        Reference: Ref.[hogg1999],Eq.(22)
-
-        Parameters
-        ----------
-        frequencies : float, or 1D `~numpy.ndarray`
-            The frequencies where to calculate the flux density.
-            Unit: [MHz]
-        **kwargs : optional arguments, i.e., ``n_e``, ``gamma``, and ``B``.
-
-        Returns
-        -------
-        flux : float, or 1D `~numpy.ndarray`
-            The calculated flux density w.r.t. each input frequency.
-            Unit: [Jy] = 1e-23 [erg/s/cm^2/Hz] = 1e-26 [W/m^2/Hz]
-        """
-        z = self.z_obs
-        freqz = np.asarray(frequencies) * (1+z)
-        power = self.calc_power(freqz, **kwargs)  # [W/Hz]
-        DL = COSMO.DL(self.z_obs) * AUC.Mpc2m  # [m]
-        flux = 1e26 * (1+z) * power / (4*np.pi * DL*DL)  # [Jy]
-        return flux
-
-    def calc_brightness_mean(self, frequencies, flux=None, pixelsize=None,
-                             **kwargs):
-        """
-        Calculate the mean surface brightness (power observed per unit
-        frequency and per unit solid angle) expressed in *brightness
-        temperature* at the specified frequencies.
-
-        NOTE
-        ----
-        If the solid angle that the object extends is smaller than the
-        specified pixel area, then is is assumed to have size of 1 pixel.
-
-        Parameters
-        ----------
-        frequencies : float, or 1D `~numpy.ndarray`
-            The frequencies where to calculate the mean brightness temperature
-            Unit: [MHz]
-        flux : float, or 1D `~numpy.ndarray`, optional
-            The flux density w.r.t. each input frequency.
-            Unit: [Jy]
-        pixelsize : float, optional
-            The pixel size of the output simulated sky image.
-            If not provided, then invoke above ``calc_flux()`` method to
-            calculate them.
-            Unit: [arcsec]
-        **kwargs : optional arguments, i.e., ``n_e``, ``gamma``, and ``B``.
-
-        Returns
-        -------
-        Tb : float, or 1D `~numpy.ndarray`
-            The mean brightness temperature at each frequency.
-            Unit: [K] <-> [Jy/pixel]
-        """
-        frequencies = np.asarray(frequencies)
-        if flux is None:
-            flux = self.calc_flux(frequencies=frequencies, **kwargs)  # [Jy]
-        else:
-            flux = np.asarray(flux)
-            if flux.shape != frequencies.shape:
-                raise ValueError("input 'frequencies' and 'flux' do not match")
-
-        omega = np.pi * self.angular_radius**2  # [arcsec^2]
-        if pixelsize and (omega < pixelsize**2):
-            omega = pixelsize ** 2  # [arcsec^2]
-            logger.warning("Halo size < 1 pixel; force to be 1 pixel!")
-
-        Tb = Fnu_to_Tb(flux, omega, frequencies)  # [K]
-        return Tb
-
     def fp_injection(self, gamma, t=None):
         """
         Electron injection (rate) term for the Fokker-Planck equation.