clusters: Implement calc_{power,flux,brightness_mean} in helper

Help development & debug ...

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

2 files changed, 116 insertions, 40 deletions

diff --git a/fg21sim/extragalactic/clusters/halo.py b/fg21sim/extragalactic/clusters/halo.py
index d11004f..f056735 100644
--- a/fg21sim/extragalactic/clusters/halo.py
+++ b/fg21sim/extragalactic/clusters/halo.py
@@ -48,7 +48,6 @@ from .emission import SynchrotronEmission
 from ...share import CONFIGS, COSMO
 from ...utils.units import (Units as AU,
                             UnitConversions as AUC)
-from ...utils.convert import Fnu_to_Tb_fast
 
 
 logger = logging.getLogger(__name__)
@@ -285,15 +284,9 @@ class RadioHalo:
 
     def calc_power(self, emissivity):
         """
-        Calculate the synchrotron power (i.e., power *emitted* per
+        Calculate the halo synchrotron power (i.e., power *emitted* per
         unit frequency) from emissivity.
 
-        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
         ----------
         emissivity : float, or 1D `~numpy.ndarray`
@@ -306,10 +299,7 @@ class RadioHalo:
             The calculated synchrotron power w.r.t. each input emissivity.
             Unit: [W/Hz]
         """
-        emissivity = np.asarray(emissivity)
-        power = emissivity * self.volume  # [erg/s/Hz]
-        power *= 1e-7  # [erg/s/Hz] -> [W/Hz]
-        return power
+        return helper.calc_power(emissivity, volume=self.volume)
 
     def calc_flux(self, emissivity):
         """
@@ -329,29 +319,21 @@ class RadioHalo:
             Unit: [Jy] = 1e-23 [erg/s/cm^2/Hz] = 1e-26 [W/m^2/Hz]
         """
         power = self.calc_power(emissivity)  # [W/Hz]
-        DL = COSMO.DL(self.z_obs) * AUC.Mpc2m  # [m]
-        flux = 1e26 * power / (4*np.pi * DL*DL)  # [Jy]
-        return flux
+        return helper.calc_flux(power, z=self.z_obs)
 
-    def calc_brightness_mean(self, emissivity, frequencies, pixelsize=None):
+    def calc_brightness_mean(self, emissivity, frequency, pixelsize=None):
         """
         Calculate the mean surface brightness (power observed per unit
-        frequency and per unit solid angle) at the specified frequencies
-        from emissivity.
-
-        NOTE
-        ----
-        If the solid angle that the radio halo extends is smaller than
-        the specified pixel area, then the halo is assumed to cover a
-        solid angle of 1 pixel.
+        frequency and per unit solid angle) expressed in *brightness
+        temperature* at the specified frequencies from emissivity.
 
         Parameters
         ----------
         emissivity : float, or 1D `~numpy.ndarray`
             The synchrotron emissivity at multiple frequencies.
             Unit: [erg/s/cm^3/Hz]
-        frequencies : float, or 1D `~numpy.ndarray`
-            The frequencies where to calculate the synchrotron emissivity.
+        frequency : float, or 1D `~numpy.ndarray`
+            The frequencies where the synchrotron emissivity is calculated.
             Unit: [MHz]
         pixelsize : float, optional
             The pixel size of the output simulated sky image.
@@ -364,21 +346,9 @@ class RadioHalo:
             Unit: [K] <-> [Jy/pixel]
         """
         omega = np.pi * self.angular_radius**2  # [arcsec^2]
-        if pixelsize:
-            pixelarea = pixelsize ** 2  # [arcsec^2]
-            if omega < pixelarea:
-                omega = pixelarea
-                logger.warning("Halo sky coverage (%.2f [arcsec^2])" % omega +
-                               " < pixel area (%.2f [arcsec^2])" % pixelarea)
-
         flux = self.calc_flux(emissivity)
-        Tb = [Fnu_to_Tb_fast(Fnu, omega, freq)
-              for Fnu, freq in zip(np.array(flux, ndmin=1),
-                                   np.array(frequencies, ndmin=1))]
-        if len(Tb) == 1:
-            return Tb[0]
-        else:
-            return np.array(Tb)
+        return helper.calc_brightness_mean(flux, frequency=frequency,
+                                           omega=omega, pixelsize=pixelsize)
 
     def fp_injection(self, gamma, t=None):
         """
diff --git a/fg21sim/extragalactic/clusters/helper.py b/fg21sim/extragalactic/clusters/helper.py
index 0d66dee..4b7185f 100644
--- a/fg21sim/extragalactic/clusters/helper.py
+++ b/fg21sim/extragalactic/clusters/helper.py
@@ -27,6 +27,7 @@ References
    http://adsabs.harvard.edu/abs/2014MNRAS.438..124Z
 """
 
+import logging
 
 import numpy as np
 from scipy import integrate
@@ -35,6 +36,10 @@ from ...share import CONFIGS, COSMO
 from ...utils.units import (Units as AU,
                             Constants as AC,
                             UnitConversions as AUC)
+from ...utils.convert import Fnu_to_Tb_fast
+
+
+logger = logging.getLogger(__name__)
 
 
 def radius_virial(mass, z=0.0):
@@ -293,3 +298,104 @@ def magnetic_field(mass):
     M_mean = 1.6e15  # [Msun]
     B = b_mean * (mass/M_mean) ** b_index
     return B
+
+
+def calc_power(emissivity, volume):
+    """
+    Calculate the synchrotron power (i.e., power *emitted* per unit
+    frequency) from emissivity, which assumed to be uniform within
+    the 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
+    ----------
+    emissivity : float, or 1D `~numpy.ndarray`
+        The synchrotron emissivity at multiple frequencies.
+        Unit: [erg/s/cm^3/Hz]
+    volume : float
+        The volume of the radio halo
+        Unit: [kpc^3]
+
+    Returns
+    -------
+    power : float, or 1D `~numpy.ndarray`
+        The calculated synchrotron power w.r.t. each input emissivity.
+        Unit: [W/Hz]
+    """
+    emissivity = np.asarray(emissivity)
+    power = emissivity * (volume * AUC.kpc2cm**3)  # [erg/s/Hz]
+    power *= 1e-7  # [erg/s/Hz] -> [W/Hz]
+    return power
+
+
+def calc_flux(power, z):
+    """
+    Calculate the synchrotron flux density (i.e., power *observed*
+    per unit frequency) from radio power at a certain redshift (i.e.,
+    distance).
+
+    Parameters
+    ----------
+    power : float, or 1D `~numpy.ndarray`
+        The synchrotron power at multiple frequencies.
+        Unit: [W/Hz]
+
+    Returns
+    -------
+    flux : float, or 1D `~numpy.ndarray`
+        The calculated synchrotron flux w.r.t. each input power.
+        Unit: [Jy] = 1e-23 [erg/s/cm^2/Hz] = 1e-26 [W/m^2/Hz]
+    """
+    DL = COSMO.DL(z) * AUC.Mpc2m  # [m]
+    flux = 1e26 * power / (4*np.pi * DL*DL)  # [Jy]
+    return flux
+
+
+def calc_brightness_mean(flux, frequency, omega, 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 from flux.
+
+    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
+    ----------
+    flux : float, or 1D `~numpy.ndarray`
+        The synchrotron flux densities at multiple frequencies.
+        Unit: [Jy]
+    frequency : float, or 1D `~numpy.ndarray`
+        The frequencies where the above flux calculated.
+        Unit: [MHz]
+    omega : float
+        The sky coverage (angular size) of the object.
+        Unit: [arcsec^2]
+    pixelsize : float, optional
+        The pixel size of the output simulated sky image.
+        Unit: [arcsec]
+
+    Returns
+    -------
+    Tb : float, or 1D `~numpy.ndarray`
+        The mean surface brightness at each frequency.
+        Unit: [K] <-> [Jy/pixel]
+    """
+    if pixelsize and (omega < pixelsize**2):
+        omega = pixelsize ** 2  # [arcsec^2]
+        logger.warning("Object sky coverage < 1 pixel; force to be 1 pixel")
+
+    Tb = [Fnu_to_Tb_fast(Fnu, omega, freq)
+          for Fnu, freq in zip(np.array(flux, ndmin=1),
+                               np.array(frequency, ndmin=1))]
+    if len(Tb) == 1:
+        return Tb[0]
+    else:
+        return np.array(Tb)