extragalactic/clusters.py: Fix multiple bugs (Test OK)

* Fix several bugs, e.g., typo, wrong quantity operation, etc.
* Add new method "_calc_specindex()", also allow further improvement
* Improve docstrings, comments, and log messages.

1 files changed, 45 insertions, 10 deletions

diff --git a/fg21sim/extragalactic/clusters.py b/fg21sim/extragalactic/clusters.py
index b20d261..8249e3a 100644
--- a/fg21sim/extragalactic/clusters.py
+++ b/fg21sim/extragalactic/clusters.py
@@ -175,7 +175,7 @@ class GalaxyClusters:
         logger.info("Adopted dimensionless Hubble parameter: {0}".format(h))
         # Cluster masses, unit: solMass (NOTE: h dependence)
         self.catalog["mass"] = (self.catalog["m"] *
-                                self.catalog_prop["m_particles"].value / h)
+                                self.catalog_prop["m_particle"].value / h)
         self.units["mass"] = au.solMass
         logger.info("Catalog: calculated cluster masses")
         # Cluster distances from the observer, unit: Mpc
@@ -199,7 +199,8 @@ class GalaxyClusters:
         fullsky = 4*np.pi * au.sr
         factor = float(fullsky / self.catalog_prop["coverage"])
         n0_cluster = len(self.catalog)
-        logger.info("Halo fraction in clusters: %f" % self.halo_fraction)
+        logger.info("Radio halo fraction in clusters: {0}".format(
+            self.halo_fraction))
         # Total number of clusters on the full sky
         N_cluster = int(n0_cluster * factor * self.halo_fraction)
         logger.info("Total number of clusters on the full sky: {0:,}".format(
@@ -281,10 +282,10 @@ class GalaxyClusters:
         Attributes
         ----------
         catalog["r_vir"] : 1D `~numpy.ndarray`
-           The virial radii (unit: Mpc) calculated from the cluster masses
+            The virial radii (unit: Mpc) calculated from the cluster masses
         catalog["size_major"], catalog["size_minor] : 1D `~numpy.ndarray`
-           The major and minor axes (unit: degree) of the clusters calculated
-           from the above virial radii and the random eccentricities.
+            The major and minor axes (unit: degree) of the clusters calculated
+            from the above virial radii and the random eccentricities.
 
         NOTE
         ----
@@ -308,7 +309,7 @@ class GalaxyClusters:
         # Calculate (elliptical) angular sizes, i.e., major and minor axes
         logger.info("Calculating the elliptical angular sizes ...")
         distance = self.catalog["distance"].data * self.units["distance"]
-        theta2 = float(2 * r_vir / distance)  # [ rad ]
+        theta2 = (r_vir / distance).decompose().value  # [ rad ]
         size_major = theta2 / (1 - self.catalog["eccentricity"]**2) ** 0.25
         size_minor = theta2 * (1 - self.catalog["eccentricity"]**2) ** 0.25
         self.catalog["size_major"] = size_major * au.rad.to(au.deg)
@@ -325,6 +326,16 @@ class GalaxyClusters:
         Then, derive the radio luminosity by employing the scaling
         relation between X-ray and radio luminosity.
 
+        Attributes
+        ----------
+        catalog["luminosity"] : 1D `~numpy.ndarray`
+            The luminosity density (at 1.4 GHz) of each cluster.
+        catalog_prop["luminosity_freq"] : `~astropy.units.Quantity`
+            The frequency (as an ``astropy`` quantity) where the above
+            luminosity derived.
+        units["luminosity"] : `~astropy.units.Unit`
+            The unit used by the above luminosity.
+
         XXX/TODO
         --------
         The scaling relations used here may be outdated, and some of the
@@ -367,7 +378,9 @@ class GalaxyClusters:
         # Hubble parameter conversion factor
         h_conv1 = (h_our / h_others) ** (b_X-2)
         # X-ray luminosity (0.1-2.4 keV) [ erg/s ]
-        L_X = (a_X * 1e45 * float(mass_RB / (1e15*au.solMass))**b_X) * h_conv1
+        L_X = ((a_X * 1e45 *
+                (mass_RB / (1e15*au.solMass)).decompose().value ** b_X) *
+               h_conv1)
         # Calculate the radio luminosity from X-ray luminosity
         a_r = 2.78
         b_r = 1.94
@@ -380,6 +393,24 @@ class GalaxyClusters:
         self.units["luminosity"] = au.W / au.Hz
         logger.info("Done Calculate the radio luminosity")
 
+    def _calc_specindex(self):
+        """Calculate the radio spectral indexes for each cluster.
+
+        Attributes
+        ----------
+        catalog["specindex"] : 1D `~numpy.ndarray`
+            The radio spectral index of each cluster.
+
+        XXX/TODO
+        --------
+        Currently, a common/uniform spectral index (1.2) is assumed for all
+        clusters, which may be improved by investigating more recent results.
+        """
+        specindex = 1.2
+        logger.info("Use common spectral index for all clusters: "
+                    "{0}".format(specindex))
+        self.catalog["specindex"] = specindex
+
     def _calc_Tb(self, luminosity, distance, specindex, frequency, size):
         """Calculate the brightness temperature at requested frequency
         by assuming a power-law spectral shape.
@@ -453,6 +484,9 @@ class GalaxyClusters:
         logger.info("Simulating HEALPix templates for each cluster ...")
         templates = []
         resolution = self.resolution.to(au.deg).value
+        # Make sure the index is reset, therefore, the *row indexes* can be
+        # simply used to identify the corresponding template image.
+        self.catalog.reset_index(inplace=True)
         # XXX/TODO: be parallel
         for row in self.catalog.itertuples():
             # TODO: progress bar
@@ -492,8 +526,8 @@ class GalaxyClusters:
         --------
         `self._simulate_template()` for more detailed description.
         """
-        name = data.name
-        hpidx, hpval = self.templates[name]
+        index = data.Index
+        hpidx, hpval = self.templates[index]
         # Calculate the brightness temperature
         luminosity = data.luminosity
         distance = data.distance
@@ -525,7 +559,7 @@ class GalaxyClusters:
         history) for the simulated products.
         """
         header = fits.Header()
-        header["COMP"] = ("extragalactic clusters of galaxies (radio halos)",
+        header["COMP"] = ("Extragalactic clusters of galaxies",
                           "Emission component")
         header["UNIT"] = ("Kelvin", "Map unit")
         header["CREATOR"] = (__name__, "File creator")
@@ -586,6 +620,7 @@ class GalaxyClusters:
         #
         self._calc_sizes()
         self._calc_luminosity()
+        self._calc_specindex()
         #
         self._simulate_templates()
         #