clusters/halo: rewrite power/flux/Tb calc with k-correction considered

3 files changed, 98 insertions, 130 deletions

diff --git a/fg21sim/extragalactic/clusters/halo.py b/fg21sim/extragalactic/clusters/halo.py
index 0f8752a..beaff9e 100644
--- a/fg21sim/extragalactic/clusters/halo.py
+++ b/fg21sim/extragalactic/clusters/halo.py
@@ -36,6 +36,10 @@ References
 .. [sarazin1999]
    Sarazin 1999, ApJ, 520, 529
    http://adsabs.harvard.edu/abs/1999ApJ...520..529S
+
+.. [hogg1999]
+   Hogg 1999, arXiv:astro-ph/9905116
+   http://adsabs.harvard.edu/abs/1999astro.ph..5116H
 """
 
 import logging
@@ -48,6 +52,7 @@ 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
 
 
 logger = logging.getLogger(__name__)
@@ -364,73 +369,137 @@ class RadioHalo:
         emissivity = syncem.emissivity(frequencies)
         return emissivity
 
-    def calc_power(self, emissivity):
+    def calc_power(self, frequencies, emissivity=None, **kwargs):
         """
         Calculate the halo synchrotron power (i.e., power *emitted* per
-        unit frequency) from emissivity.
+        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
         ----------
-        emissivity : float, or 1D `~numpy.ndarray`
-            The synchrotron emissivity at multiple frequencies.
+        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`` and ``gamma``
 
         Returns
         -------
         power : float, or 1D `~numpy.ndarray`
-            The calculated synchrotron power w.r.t. each input emissivity.
+            The calculated synchrotron power at each input frequency.
             Unit: [W/Hz]
         """
-        return helper.calc_power(emissivity, volume=self.volume)
+        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, emissivity):
+    def calc_flux(self, frequencies, **kwargs):
         """
         Calculate the synchrotron flux density (i.e., power *observed*
-        per unit frequency) from emissivity.
+        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
         ----------
-        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 flux density.
+            Unit: [MHz]
+        **kwargs : optional arguments, i.e., ``n_e`` and ``gamma``
 
         Returns
         -------
         flux : float, or 1D `~numpy.ndarray`
-            The calculated synchrotron flux w.r.t. each input emissivity.
+            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]
         """
-        power = self.calc_power(emissivity)  # [W/Hz]
-        return helper.calc_flux(power, z=self.z_obs)
+        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, emissivity, frequency, pixelsize=None):
+    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 from emissivity.
+        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
         ----------
-        emissivity : float, or 1D `~numpy.ndarray`
-            The synchrotron emissivity at multiple frequencies.
-            Unit: [erg/s/cm^3/Hz]
-        frequency : float, or 1D `~numpy.ndarray`
-            The frequencies where the synchrotron emissivity is calculated.
+        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`` and ``gamma``
 
         Returns
         -------
         Tb : float, or 1D `~numpy.ndarray`
-            The mean surface brightness at each frequency.
+            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]
-        flux = self.calc_flux(emissivity)
-        return helper.calc_brightness_mean(flux, frequency=frequency,
-                                           omega=omega, pixelsize=pixelsize)
+        if pixelsize and (omega < pixelsize**2):
+            omega = pixelsize ** 2  # [arcsec^2]
+            logger.warning("Object 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):
         """
diff --git a/fg21sim/extragalactic/clusters/helper.py b/fg21sim/extragalactic/clusters/helper.py
index e553da9..6b835a0 100644
--- a/fg21sim/extragalactic/clusters/helper.py
+++ b/fg21sim/extragalactic/clusters/helper.py
@@ -40,7 +40,6 @@ 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
 from ...utils.draw import circle
 from ...utils.transform import circle2ellipse
 
@@ -304,107 +303,6 @@ def magnetic_field(mass):
     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(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)
-
-
 def halo_rprofile(re, num_re=5, I0=1.0):
     """
     Generate the radial profile of a halo.
diff --git a/fg21sim/extragalactic/clusters/main.py b/fg21sim/extragalactic/clusters/main.py
index 70d5e6b..e6c6154 100644
--- a/fg21sim/extragalactic/clusters/main.py
+++ b/fg21sim/extragalactic/clusters/main.py
@@ -303,9 +303,10 @@ class GalaxyClusters:
             halo.set_electron_spectrum(hdict["n_e"])
 
             emissivity = halo.calc_emissivity(frequencies=self.frequencies)
-            power = halo.calc_power(emissivity)
-            flux = halo.calc_flux(emissivity)
-            Tb_mean = halo.calc_brightness_mean(emissivity, self.frequencies,
+            power = halo.calc_power(self.frequencies, emissivity=emissivity)
+            # k-correction considered
+            flux = halo.calc_flux(self.frequencies)
+            Tb_mean = halo.calc_brightness_mean(self.frequencies, flux=flux,
                                                 pixelsize=self.sky.pixelsize)
             # Update or add new items
             hdict["frequency"] = self.frequencies  # [MHz]