From c0b81153caea9b4299118f25dbc5856b2cfe6fde Mon Sep 17 00:00:00 2001
From: Aaron LI <aly@aaronly.me>
Date: Sat, 12 Aug 2017 23:33:35 +0800
Subject: clusters/main.py: Implement "_draw_halos()" method

New functions "halo_rprofile()" and "draw_halo()" added to helper.py

Signed-off-by: Aaron LI <aly@aaronly.me>
---
 fg21sim/extragalactic/clusters/helper.py | 62 ++++++++++++++++++++++++++++++++
 fg21sim/extragalactic/clusters/main.py   | 29 +++++++++++++--
 2 files changed, 89 insertions(+), 2 deletions(-)

diff --git a/fg21sim/extragalactic/clusters/helper.py b/fg21sim/extragalactic/clusters/helper.py
index 4b7185f..becb608 100644
--- a/fg21sim/extragalactic/clusters/helper.py
+++ b/fg21sim/extragalactic/clusters/helper.py
@@ -22,6 +22,10 @@ References
    Cassano et al. 2012, A&A, 548, A100
    http://adsabs.harvard.edu/abs/2012A%26A...548A.100C
 
+.. [murgia2009]
+   Murgia et al. 2009, A&A, 499, 679
+   http://adsabs.harvard.edu/abs/2009A%26A...499..679M
+
 .. [zandanel2014]
    Zandanel, Pfrommer & Prada 2014, MNRAS, 438, 124
    http://adsabs.harvard.edu/abs/2014MNRAS.438..124Z
@@ -37,6 +41,8 @@ from ...utils.units import (Units as AU,
                             Constants as AC,
                             UnitConversions as AUC)
 from ...utils.convert import Fnu_to_Tb_fast
+from ...utils.draw import circle
+from ...utils.transform import circle2ellipse
 
 
 logger = logging.getLogger(__name__)
@@ -399,3 +405,59 @@ def calc_brightness_mean(flux, frequency, omega, pixelsize=None):
         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.
+
+    NOTE
+    ----
+    The exponential radial profile is adopted for the radio halos:
+        I(r) = I0 * exp(-r/re)
+    with the e-folding radius ``re ~ R_halo / 3``.
+
+    Parameters
+    ----------
+    re : float
+        The e-folding radius in unit of pixels.
+    num_re : float, optional
+        The times of ``re`` to determine the maximum radius.
+        Default: 5, i.e., rmax = 5 * re
+    I0 : float
+        The intensity/brightness at the center (i.e., r=0)
+        Default: 1.0
+
+    Returns
+    -------
+    rprofile : 1D `~numpy.ndarray`
+        The values along the radial pixels (0, 1, 2, ...)
+
+    References: Ref.[murgia2009],Eq.(1)
+    """
+    rmax = round(re * num_re)
+    r = np.arange(rmax+1)
+    rprofile = I0 * np.exp(-r/re)
+    return rprofile
+
+
+def draw_halo(rprofile, felong, rotation=0.0):
+    """
+    Draw the template image of one halo, which is used to simulate
+    the image at requested frequencies by adjusting the brightness
+    values.
+
+    Parameters
+    ----------
+    rprofile
+
+    Returns
+    -------
+    image : 2D `~numpy.ndarray`
+        2D array of the drawn halo template image.
+        The image is normalized to have sum of 1.
+    """
+    image = circle(rprofile=rprofile)
+    image = circle2ellipse(image, bfraction=felong, rotation=rotation)
+    image /= image.sum()
+    return image
diff --git a/fg21sim/extragalactic/clusters/main.py b/fg21sim/extragalactic/clusters/main.py
index 959e877..2d8f4d5 100644
--- a/fg21sim/extragalactic/clusters/main.py
+++ b/fg21sim/extragalactic/clusters/main.py
@@ -30,6 +30,7 @@ from ...share import CONFIGS, COSMO
 from ...utils.io import dataframe_to_csv, pickle_dump
 from ...utils.ds import dictlist_to_dataframe
 from ...sky import get_sky
+from . import helper
 
 
 logger = logging.getLogger(__name__)
@@ -294,7 +295,30 @@ class GalaxyClusters:
         logger.info("Done halos data conversion.")
 
     def _draw_halos(self):
-        pass
+        """
+        Draw the template images for each halo, and cache them for
+        simulating the superimposed halos at requested frequencies.
+
+        NOTE
+        ----
+        The drawn template images are append to the dictionaries of
+        the corresponding halo within the ``self.halos``.
+        The templates are normalized to have sum of 1.
+        """
+        num = len(self.halos)
+        logger.info("Draw template images for %d halos ..." % num)
+        self.halos = []
+        i = 0
+        for hdict in self.halos:
+            i += 1
+            if i % 50 == 0:
+                logger.info("[%d/%d] %.1f%% ..." % (i, num, 100*i/num))
+            theta_e = hdict["angular_radius"] / self.sky.pixelsize
+            rprofile = helper.halo_rprofile(re=theta_e)
+            template = helper.draw_halo(rprofile, felong=hdict["felong"],
+                                        rotation=hdict["rotation"])
+            hdict["template"] = template
+        logger.info("Done drawn halo template images.")
 
     def preprocess(self):
         """
@@ -314,7 +338,8 @@ class GalaxyClusters:
         self._process_catalog()
         self._simulate_mergers()
         self._simulate_halos()
-        # TODO ???
+        self._draw_halos()
+
         self._preprocessed = True
 
     def postprocess(self):
-- 
cgit v1.2.2