From 200ac8cdc3661af0306b85d7589a13759a66ee41 Mon Sep 17 00:00:00 2001
From: Aaron LI <aaronly.me@outlook.com>
Date: Mon, 24 Oct 2016 22:21:31 +0800
Subject: convert.py: Optimize "Fnu_to_Tb()" and "Sb_to_Tb()"

* Optimize the "Fnu_to_Tb()" and "Sb_to_Tb()" functions by explicitly
  calculating the conversions, avoiding the slow `astropy.unit`
  conversions (which is rather slow).  The new fast functions are named
  as "Fnu_to_Tb_fast()" and "Sb_to_Tb_fast()".
* Optimize the new functions using `numba.jit` further.
* Add `numba` to the dependencies.
---
 fg21sim/utils/convert.py | 65 ++++++++++++++++++++++++++++++++++++++++++++++++
 requirements.txt         |  1 +
 setup.py                 |  1 +
 3 files changed, 67 insertions(+)

diff --git a/fg21sim/utils/convert.py b/fg21sim/utils/convert.py
index 1212780..fc9a228 100644
--- a/fg21sim/utils/convert.py
+++ b/fg21sim/utils/convert.py
@@ -5,7 +5,9 @@
 Utilities for conversion among common astronomical quantities.
 """
 
+import numpy as np
 import astropy.units as au
+import numba
 
 
 def Fnu_to_Tb(Fnu, omega, freq):
@@ -61,3 +63,66 @@ def Sb_to_Tb(Sb, freq):
     omega = 1.0 * au.sr
     Fnu = Sb * omega
     return Fnu_to_Tb(Fnu, omega, freq)
+
+
+@numba.jit(nopython=True)
+def Sb_to_Tb_fast(Sb, freq):
+    """Convert surface brightness to brightness temperature, using the
+    Rayleigh-Jeans law, in the Rayleigh-Jeans limit.
+
+    This function does the calculations explicitly, and does NOT rely
+    on the `astropy.units`, therefore it is faster.  However, the input
+    parameters must be in right units.
+
+        Tb = Sb * c^2 / (2 * k_B * nu^2)
+
+    where `SB` is the surface brightness density measured at a certain
+    frequency (unit: [ Jy/rad^2 ] = [ erg/s/cm^2/Hz/rad^2 ]).
+
+    Parameters
+    ----------
+    Sb : float
+        Input surface brightness, unit [ Jy/deg^2 ]
+    freq : float
+        Frequency where the flux density measured, unit [ MHz ]
+
+    Returns
+    -------
+    Tb : float
+        Calculated brightness temperature, unit [ K ]
+    """
+    # NOTE: `radian` is dimensionless
+    rad2_to_deg2 = np.rad2deg(1.0) * np.rad2deg(1.0)
+    Sb_rad2 = Sb * rad2_to_deg2  # unit: [ Jy/rad^2 ] -> [ Jy ]
+    c = 29979245800.0  # speed of light, [ cm/s ]
+    k_B = 1.3806488e-16  # Boltzmann constant, [ erg/K ]
+    coef = 1e-35  # take care the unit conversions
+    Tb = coef * (Sb_rad2 * c*c) / (2 * k_B * freq*freq)  # unit: [ K ]
+    return Tb
+
+
+@numba.jit(nopython=True)
+def Fnu_to_Tb_fast(Fnu, omega, freq):
+    """Convert flux density to brightness temperature, using the
+    Rayleigh-Jeans law, in the Rayleigh-Jeans limit.
+
+    This function does the calculations explicitly, and does NOT rely
+    on the `astropy.units`, therefore it is faster.  However, the input
+    parameters must be in right units.
+
+    Parameters
+    ----------
+    Fnu : float
+        Input flux density, unit [ Jy ]
+    omega : float
+        Source angular size/area, unit [ deg^2 ]
+    freq : float
+        Frequency where the flux density measured, unit [ MHz ]
+
+    Returns
+    -------
+    Tb : float
+        Calculated brightness temperature, unit [ K ]
+    """
+    Sb = Fnu / omega
+    return Sb_to_Tb_fast(Sb, freq)
diff --git a/requirements.txt b/requirements.txt
index 5e9a12c..60199f0 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,4 +1,5 @@
 numpy
+numba
 scipy
 pandas
 astropy
diff --git a/setup.py b/setup.py
index b7725e3..1fce5c6 100755
--- a/setup.py
+++ b/setup.py
@@ -80,6 +80,7 @@ setup(
     ],
     install_requires=[
         "numpy",
+        "numba",
         "scipy",
         "pandas",
         "astropy",
-- 
cgit v1.2.2