From 226af8a558ea70e5d9d3fa39855e3c3d00ec57e3 Mon Sep 17 00:00:00 2001
From: Aaron LI <aly@aaronly.me>
Date: Sat, 29 Jul 2017 15:18:45 +0800
Subject: galactic/freefree: Add config "electron_temperature"

Update method "_calc_ratio_a()" to "_calc_factor_a()"

Signed-off-by: Aaron LI <aly@aaronly.me>
---
 fg21sim/configs/10-galactic.conf.spec |  5 +++++
 fg21sim/galactic/freefree.py          | 33 +++++++++++++++++++++++----------
 2 files changed, 28 insertions(+), 10 deletions(-)

diff --git a/fg21sim/configs/10-galactic.conf.spec b/fg21sim/configs/10-galactic.conf.spec
index 1ef15eb..b082a5d 100644
--- a/fg21sim/configs/10-galactic.conf.spec
+++ b/fg21sim/configs/10-galactic.conf.spec
@@ -62,6 +62,11 @@
   # The unit of the above dust map (e.g., "MJy/sr")
   dustmap_unit = string(default=None)
 
+  # The electron temperature assumed for the ionized interstellar medium
+  # that generating H{\alpha} emission.
+  # Unit: [K]
+  electron_temperature = float(default=7000.0, min=1000)
+
   # Filename prefix for this component
   prefix = string(default="gfree")
   # Whether save this component to disk
diff --git a/fg21sim/galactic/freefree.py b/fg21sim/galactic/freefree.py
index 41e84e6..a4d8ad6 100644
--- a/fg21sim/galactic/freefree.py
+++ b/fg21sim/galactic/freefree.py
@@ -75,6 +75,7 @@ class FreeFree:
         self.dustmap_path = self.configs.get_path(comp+"/dustmap")
         self.dustmap_unit = au.Unit(
             self.configs.getn(comp+"/dustmap_unit"))
+        self.Te = self.configs.getn(comp+"/electron_temperature")  # [K]
         self.prefix = self.configs.getn(comp+"/prefix")
         self.save = self.configs.getn(comp+"/save")
         self.output_dir = self.configs.get_path(comp+"/output_dir")
@@ -143,15 +144,27 @@ class FreeFree:
         self._dust_corrected = True
         logger.info("Done dust absorption correction")
 
-    def _calc_ratio_a(self, Te, nu_GHz):
-        """Calculate the ratio factor a(T, nu), which will be used to
-        transform the Halpha emission (Rayleigh) to free-free emission
-        brightness temperature (mK).
+    def _calc_factor_a(self, nu):
+        """
+        Calculate the ratio factor a(Te, ν), which will be used to
+        convert the Halpha emission [Rayleigh] to free-free emission
+        brightness temperature [K].
+
+        Parameters
+        ----------
+        nu : float
+            The frequency where to calculate the factor a(nu).
+            Unit: [MHz]
+
+        Returns
+        -------
+        a : float
+            The factor for Hα to free-free conversion.
 
-        References: [Dickinson2003], Eq.(8)
+        References: [dickinson2003],Eq.(8)
         """
-        term1 = 0.366 * nu_GHz**0.1 * Te**(-0.15)
-        term2 = np.log(4.995e-2 / nu_GHz) + 1.5*np.log(Te)
+        term1 = 0.183 * nu**0.1 * self.Te**(-0.15)
+        term2 = 3.91 - np.log(nu) + 1.5*np.log(self.Te)
         a = term1 * term2
         return a
 
@@ -256,12 +269,12 @@ class FreeFree:
         logger.info("Simulating {name} map at {freq} ({unit}) ...".format(
             name=self.name, freq=frequency, unit=self.freq_unit))
         # Assumed electron temperature [ K ]
-        Te = 7000.0
-        T4 = Te / 1e4
+        T4 = self.Te / 1e4
         nu = frequency * self.freq_unit.to(au.GHz)  # frequency [ GHz ]
-        ratio_a = self._calc_ratio_a(Te, nu)
         # NOTE: ignored the "10^3" term in the referred equation
         ratio_mK_R = (8.396 * ratio_a * nu**(-2.1) *
+        factor_a = self._calc_factor_a(frequency)
+        ratio_mK_R = (8.396 * factor_a * nu**(-2.1) *
                       T4**0.667 * 10**(0.029/T4) * (1+0.08))
         # Use "Kelvin" as the brightness temperature unit
         ratio_K_R = ratio_mK_R * au.mK.to(au.K)
-- 
cgit v1.2.2