extragalactic/clusters.py: Use "Fnu_to_Tb_fast()"

1 files changed, 15 insertions, 13 deletions

diff --git a/fg21sim/extragalactic/clusters.py b/fg21sim/extragalactic/clusters.py
index 8249e3a..7da778e 100644
--- a/fg21sim/extragalactic/clusters.py
+++ b/fg21sim/extragalactic/clusters.py
@@ -23,7 +23,7 @@ import pandas as pd
 
 from ..utils import write_fits_healpix
 from ..utils.random import spherical_uniform
-from ..utils.convert import Fnu_to_Tb
+from ..utils.convert import Fnu_to_Tb_fast
 from ..utils.grid import make_grid_ellipse, map_grid_to_healpix
 
 
@@ -228,7 +228,6 @@ class GalaxyClusters:
         glat = 90.0 - np.degrees(theta)
         self.catalog["glon"] = glon
         self.catalog["glat"] = glat
-        self.units["coord"] = au.deg
         logger.info("Done add random positions for each cluster")
 
     def _add_random_eccentricity(self):
@@ -353,7 +352,7 @@ class GalaxyClusters:
           *cosmic mean mass density ρ_m = Ω_m * ρ_c,
           therefore, the mass needs following conversion (which is an
           approximation):
-              M_{R&B} ~= M * sqrt(OmegaM0)
+              M_{R&B} ≈ M * sqrt(OmegaM0)
         - The derived X-ray luminosity is for the 0.1-2.4 keV energy band.
         - The X-ray-radio luminosity scaling relation adopted here is
           derived at 1.4 GHz.
@@ -370,7 +369,7 @@ class GalaxyClusters:
         #
         logger.info("Calculating the radio luminosity (at 1.4 GHz) ...")
         # Calculate the X-ray luminosity from mass
-        # XXX: why the mass conversion ??
+        # NOTE: mass conversion (see also the above notes)
         mass_RB = (self.catalog["mass"].data * self.units["mass"] *
                    self.catalog_prop["omega_m"]**0.5)
         a_X = 0.449
@@ -389,7 +388,7 @@ class GalaxyClusters:
         # Radio luminosity density (at 1.4 GHz) [ W/Hz ]
         L_r = (a_r * 1e24 * (L_X / 1e45)**b_r) * h_conv2
         self.catalog["luminosity"] = L_r
-        self.catalog_prop["luminosity_freq"] = 1.4 * au.GHz
+        self.catalog_prop["luminosity_freq"] = 1400 * self.freq_unit
         self.units["luminosity"] = au.W / au.Hz
         logger.info("Done Calculate the radio luminosity")
 
@@ -428,10 +427,10 @@ class GalaxyClusters:
             The spectral index of the power-law spectrum.
             Note the *negative* sign in the formula.
         frequency : float
-            The frequency (unit: `self.freq_unit`) where the brightness
+            The frequency (unit: [ MHz ]) where the brightness
             temperature requested.
         size : 2-float tuple
-            The (major, minor) axes (unit: `self.units["size"]`).
+            The (major, minor) axes (unit: [ deg ]).
             The order of major and minor can be arbitrary.
 
         Returns
@@ -447,14 +446,15 @@ class GalaxyClusters:
         be calculated by extrapolating the spectrum, then convert the flux
         density to derive the brightness temperature.
         """
-        freq = frequency * self.freq_unit
-        freq_ref = self.catalog_prop["luminosity_freq"]
+        freq = frequency  # [ MHz ]
+        freq_ref = self.catalog_prop["luminosity_freq"].value
         luminosity = luminosity * self.units["luminosity"]
-        Lnu = luminosity * float(freq / freq_ref) ** (-specindex)
+        Lnu = luminosity * (freq / freq_ref) ** (-specindex)
         Fnu = Lnu / (4*np.pi * (distance*self.units["distance"])**2)
-        omega = size[0]*self.units["size"] * size[1]*self.units["size"]
-        Tb = Fnu_to_Tb(Fnu, omega, freq)
-        return Tb.value
+        Fnu_Jy = Fnu.to(au.Jy).value  # [ Jy ]
+        omega = size[0] * size[1]  # [ deg^2 ]
+        Tb = Fnu_to_Tb_fast(Fnu_Jy, omega, freq)
+        return Tb
 
     def _simulate_templates(self):
         """Simulate the template (HEALPix) images for each cluster, and
@@ -649,7 +649,9 @@ class GalaxyClusters:
         logger.info("Simulating {name} map at {freq} ({unit}) ...".format(
             name=self.name, freq=frequency, unit=self.freq_unit))
         hpmap_f = np.zeros(hp.nside2npix(self.nside))
+        # XXX/TODO: be parallel
         for row in self.catalog.itertuples():
+            # TODO: progress bar
             hpidx, hpval = self._simulate_single(row, frequency)
             hpmap_f[hpidx] += hpval
         #