From 056d116db845bc7162e148cb5a9849360d949f0d Mon Sep 17 00:00:00 2001
From: Aaron LI <aly@aaronly.me>
Date: Thu, 1 Nov 2018 11:51:45 +0800
Subject: clusters/halo: Support to check whether halo is genuine

Add the 'is_genuine()' method to check whether the radio halo is genuine
by comparing the flux density to the fiducial value.  Two new configuration
options (fiducial_freq and fiducial_factor) are added to specify the
frequency and limit whether the radio halo is regarded as genuine.
---
 fg21sim/configs/config.spec            |  7 +++++++
 fg21sim/extragalactic/clusters/halo.py | 33 +++++++++++++++++++++++++++++++++
 2 files changed, 40 insertions(+)

diff --git a/fg21sim/configs/config.spec b/fg21sim/configs/config.spec
index 0d7e7aa..fefb8b8 100644
--- a/fg21sim/configs/config.spec
+++ b/fg21sim/configs/config.spec
@@ -440,3 +440,10 @@ stream = option("stderr", "stdout", "", default="stderr")
   # the acceleration is turned off and only leaves energy loss mechanisms.
   # Unit: [Gyr]
   time_init = float(default=1.0, min=0)
+
+  # The frequency and factor used to determine whether the radio halo is
+  # genuine, i.e., the flux density at ``fiducial_freq`` is at least
+  # ``fiducial_factor`` times the fiducial value.
+  # Unit: [MHz]
+  fiducial_freq = float(default=150.0)
+  fiducial_factor = float(default=3.0, min=1)
diff --git a/fg21sim/extragalactic/clusters/halo.py b/fg21sim/extragalactic/clusters/halo.py
index ed6b4d3..0db0b0d 100644
--- a/fg21sim/extragalactic/clusters/halo.py
+++ b/fg21sim/extragalactic/clusters/halo.py
@@ -63,6 +63,7 @@ from scipy import integrate
 
 from . import helper
 from .solver import FokkerPlanckSolver
+from .emission import HaloEmission
 from ...share import CONFIGS, COSMO
 from ...utils.units import (Units as AU,
                             UnitConversions as AUC,
@@ -164,6 +165,8 @@ class RadioHalo:
         self.time_step = configs.getn(comp+"/time_step")
         self.time_init = configs.getn(comp+"/time_init")
         self.injection_index = configs.getn(comp+"/injection_index")
+        self.fiducial_freq = configs.getn(comp+"/fiducial_freq")
+        self.fiducial_factor = configs.getn(comp+"/fiducial_factor")
 
     def _set_solver(self):
         self.fpsolver = FokkerPlanckSolver(
@@ -478,12 +481,42 @@ class RadioHalo:
             self._acceleration_disabled = True
             self.fpsolver.tstep = self.time_step * 2  # To save time
 
+        logger.debug("Calculating the %s electron spectrum ..." %
+                     ("[fiducial]" if fiducial else ""))
         n_e = self.fpsolver.solve(u0=n0_e, tstart=tstart, tstop=tstop)
         self._acceleration_disabled = False
         self.fpsolver.tstep = self.time_step
 
         return n_e
 
+    def is_genuine(self, n_e):
+        """
+        Check whether the radio halo is genuine/observable by comparing the
+        radio flux density to the fiducial value, which is calculated from
+        the fiducial electron spectrum derived with turbulent acceleration
+        turned off.
+
+        Parameters
+        ----------
+        n_e : float 1D `~numpy.ndarray`
+            The finally derived electron spectrum.
+            Unit: [cm^-3]
+
+        Returns
+        -------
+        genuine : bool
+        """
+        haloem = HaloEmission(gamma=self.gamma, n_e=n_e,
+                              B=self.B_obs, radius=self.radius,
+                              redshift=self.z_obs)
+        flux = haloem.calc_flux(self.fiducial_freq)
+
+        ne_fiducial = self.calc_electron_spectrum(fiducial=True)
+        haloem.n_e = ne_fiducial
+        flux_fiducial = haloem.calc_flux(self.fiducial_freq)
+
+        return flux >= flux_fiducial * self.fiducial_factor
+
     def fp_injection(self, gamma, t=None):
         """
         Electron injection (rate) term for the Fokker-Planck equation.
-- 
cgit v1.2.2