fix phase fn calc for small x

Author: scottprahl

CHANGELOG.rst

@@ -1,6 +1,10 @@
 Changelog
 =========
 
+v2.5.0 (8/4/2023)
+-------------------
+*   fix scattering function for very small spheres
+
 v2.4.0 (6/10/2023)
 -------------------
 *   add mie_phase_matrix() to calculate scattering (Mueller) matrix

miepython/__init__.py

@@ -42,7 +42,7 @@
     miepython.mie_phase_matrix(m, x, mu)
 
 """
-__version__ = '2.4.0'
+__version__ = '2.5.0'
 __author__ = 'Scott Prahl'
 __email__ = 'scott.prahl@oit.edu'
 __copyright__ = 'Copyright 2017-23, Scott Prahl'

miepython/miepython.py

@@ -476,8 +476,8 @@ def _small_mie_S1_S2(m, x, mu):
     return [S1, S2]
 
 
-@njit((complex128[:], complex128[:], float64, int32), cache=True)
-def normalization_factor(a, b, x, norm_int):
+@njit((complex128, float64, int32), cache=True)
+def normalization_factor(m, x, norm_int):
     """
     Figure out scattering function normalization.
 
@@ -492,7 +492,7 @@ def normalization_factor(a, b, x, norm_int):
     """
     # Qsca normalization
     if norm_int == 3:
-        return np.sqrt(np.pi * x**2)
+        return x * np.sqrt(np.pi)
 
     # Bohren Normalization
     if norm_int == 5:
@@ -503,26 +503,23 @@ def normalization_factor(a, b, x, norm_int):
         return 1
 
     # calculate qsca and qext
-    n = np.arange(1, len(a) + 1)
-    cn = 2.0 * n + 1.0
-    qext = 2 * np.sum(cn * (a.real + b.real)) / x**2
-    qsca = 2 * np.sum(cn * (np.abs(a)**2 + np.abs(b)**2)) / x**2
+    qext, qsca, _, _ = _mie_scalar(m, x)
 
     # albedo Normalization
     if norm_int == 0:
-        return np.sqrt(np.pi * x**2 * qext)
+        return x * np.sqrt(np.pi * qext)
 
     # Unity normalization
     if norm_int == 1:
-        return np.sqrt(qsca * np.pi * x**2)
+        return x * np.sqrt(qsca * np.pi)
 
     # 4pi Normalization
     if norm_int == 2:
-        return np.sqrt(qsca * x**2 / 4)
+        return x * np.sqrt(qsca / 4)
 
     # Qext Normalization
     if norm_int == 4:  # 4pi
-        return np.sqrt(qsca * np.pi * x**2 / qext)
+        return x * np.sqrt(qsca * np.pi / qext)
 
     raise ValueError("norm-int must be in the range 0..6")
 
@@ -609,7 +606,7 @@ def _mie_S1_S2(m, x, mu, norm_int):
             pi_nm1 = ((2 * n + 1) * mu[k] * pi_nm1 - (n + 1) * pi_nm2) / n
             pi_nm2 = temp
 
-    normalization = normalization_factor(a, b, x, norm_int)
+    normalization = normalization_factor(m, x, norm_int)
 
     S1 /= normalization
     S2 /= normalization

miepython/miepython_nojit.py

@@ -465,13 +465,12 @@ def _small_mie_S1_S2(m, x, mu):
     return [S1, S2]
 
 
-def normalization_factor(a, b, x, norm_str):
+def normalization_factor(m, x, norm_str):
     """
     Figure out scattering function normalization.
 
     Args:
-        a: complex array of An coefficients
-        b: complex array of Bn coefficients
+        m: complex index of refraction of sphere
         x: dimensionless sphere size
         norm_str: string describing type of normalization
 
@@ -486,26 +485,22 @@ def normalization_factor(a, b, x, norm_str):
     if norm in ['wiscombe']:
         return 1
 
-    n = np.arange(1, len(a) + 1)
-    cn = 2.0 * n + 1.0
-    qext = 2 * np.sum(cn * (a.real + b.real)) / x**2
+    if norm in ['qsca', 'scattering_efficiency']:
+        return x * np.sqrt(np.pi)
 
-    if norm in ['a', 'albedo']:
-        return np.sqrt(np.pi * x**2 * qext)
+    qext, qsca, _, _ = _mie_scalar(m, x)
 
-    qsca = 2 * np.sum(cn * (np.abs(a)**2 + np.abs(b)**2)) / x**2
+    if norm in ['a', 'albedo']:
+        return x * np.sqrt(np.pi * qext)
 
     if norm in ['1', 'one', 'unity']:
-        return np.sqrt(qsca * np.pi * x**2)
+        return x * np.sqrt(qsca * np.pi)
 
     if norm in ['four_pi', '4pi']:
-        return np.sqrt(qsca * x**2 / 4)
-
-    if norm in ['qsca', 'scattering_efficiency']:
-        return np.sqrt(np.pi * x**2)
+        return x * np.sqrt(qsca / 4)
 
     if norm in ['qext', 'extinction_efficiency']:
-        return np.sqrt(qsca * np.pi * x**2 / qext)
+        return x * np.sqrt(qsca * np.pi / qext)
 
     raise ValueError("normalization must be one of 'albedo' (default), 'one'"
                      "'4pi', 'qext', 'qsca', 'bohren', or 'wiscombe'")
@@ -549,7 +544,7 @@ def mie_S1_S2(m, x, mu, norm='albedo'):
             pi_nm1 = ((2 * n + 1) * mu[k] * pi_nm1 - (n + 1) * pi_nm2) / n
             pi_nm2 = temp
 
-    normalization = normalization_factor(a, b, x, norm)
+    normalization = normalization_factor(m, x, norm)
 
     S1 /= normalization
     S2 /= normalization

tests/test_jit.py

@@ -467,6 +467,12 @@ def test_i_par_i_per_01(self):
             total = (total1 + total2) / 2
             self.assertAlmostEqual(total / expected[i], 1, delta=1e-3)
 
+    def test_molecular_hydrogen(self):
+        m = 1.00013626
+        x = 0.0006403246172921872
+        mu = np.linspace(-1, 1, 100)
+        ph = miepython.i_unpolarized(m, x, mu)
+        self.assertAlmostEqual(ph[1], 0.1169791, delta=1e-5)
 
 class MiePhaseMatrix(unittest.TestCase):
     def test_mie_phase_matrix_basic(self):

tests/test_nojit.py

@@ -467,6 +467,12 @@ def test_i_par_i_per_01(self):
             total = (total1 + total2) / 2
             self.assertAlmostEqual(total / expected[i], 1, delta=1e-3)
 
+    def test_molecular_hydrogen(self):
+        m = 1.00013626
+        x = 0.0006403246172921872
+        mu = np.linspace(-1, 1, 100)
+        ph = miepython.i_unpolarized(m, x, mu)
+        self.assertAlmostEqual(ph[1], 0.1169791, delta=1e-5)
 
 class MiePhaseMatrix(unittest.TestCase):
     def test_mie_phase_matrix_basic(self):