Add scattering after correcting the Rayleigh absorption sign

CHANGELOG.rst

@@ -3,6 +3,15 @@ Changelog
 
 unreleased
 -------------------
+*   fix ``rayleigh.efficiencies_mx`` subtracting scattering from extinction.  With
+    the m = n - ik convention the absorption term is negative, so
+    ``abs(absorption + qsca)`` gave absorption minus scattering, and qext for
+    m = 1.5 - 1e-6j at x = 0.01 came out 21% low (thanks @antonsoo).  The x**4 term
+    now follows Bohren & Huffman's ``Re(ratio**2)`` rather than ``|ratio|**2``,
+    which agrees with Mie to 2e-5 instead of 6e-4 at m = 1.5 - 1j, x = 0.1.  ``m``
+    is conjugated first, as ``core.py`` does.  The Rayleigh ``qext`` also sets the
+    ``albedo`` and ``qext`` normalizations, so ``rayleigh.S1_S2`` and the
+    intensities change too
 *   support Python 3.11 through 3.15.  ``requires-python`` is now ``>=3.11``, 3.10 is
     dropped from the classifiers, the black target versions and the CI matrix, and
     3.15 is added to all three

miepython/rayleigh.py

@@ -32,10 +32,19 @@ def efficiencies_mx(m, x):
         qback: the backscatter efficiency
         g: the average cosine of the scattering phase function
     """
+    # ensure imaginary part of refractive index is negative, as core.py does
+    if np.isscalar(m):
+        m = np.conj(m) if np.imag(m) > 0 else m
+    else:
+        m = np.where(np.imag(m) > 0, np.conj(m), m)
+
     ratio = (m**2 - 1) / (m**2 + 2)
     qsca = 8 / 3 * x**4 * abs(ratio) ** 2
-    qext = 4 * x * ratio * (1 + x**2 / 15 * ratio * (m**4 + 27 * m**2 + 38) / (2 * m**2 + 3))
-    qext = abs(qext.imag + qsca)
+    # Bohren & Huffman write m = n + ik, so with m = n - ik the absorption term
+    # changes sign.  The x**4 term is Re(ratio**2), not |ratio|**2: the two only
+    # agree for a non-absorbing sphere, and the difference matters near resonance.
+    qabs = -4 * x * (ratio * (1 + x**2 / 15 * ratio * (m**4 + 27 * m**2 + 38) / (2 * m**2 + 3))).imag
+    qext = qabs + 8 / 3 * x**4 * (ratio**2).real
     qback = 4 * x**4 * abs(ratio) ** 2
     g = 0
     return qext, qsca, qback, g

tests/test_rayleigh.py

@@ -37,6 +37,39 @@ def test_efficiencies(m, case):
     assert np.isclose(g_ray, g_mie, atol=2e-5), f"g mismatch for {case}"  # g should be 0 for Rayleigh
 
 
+@pytest.mark.parametrize("index", [1.33, 1.5])
+@pytest.mark.parametrize("absorption", [1e-7, 1e-6])
+@pytest.mark.parametrize("size", [0.01, 0.05])
+def test_weak_absorption_is_positive(index, absorption, size):
+    """Extinction must include positive absorption under either index convention."""
+    m = index - 1j * absorption
+    negative = rayleigh.efficiencies_mx(m, size)
+    positive = rayleigh.efficiencies_mx(np.conj(m), size)
+    np.testing.assert_allclose(negative, positive, rtol=1e-13)
+    qabs = negative[0] - negative[1]
+    assert qabs > 0
+    # Leading Rayleigh absorption: 4*x*Im((m**2-1)/(m**2+2)), with
+    # a positive imaginary index. Finite-size corrections are O(x**2).
+    polarizability = (np.conj(m) ** 2 - 1) / (np.conj(m) ** 2 + 2)
+    np.testing.assert_allclose(qabs, 4 * size * polarizability.imag, rtol=0.003)
+
+
+@pytest.mark.parametrize(
+    "m, size, rtol",
+    [
+        (1.5 - 1e-6j, 0.01, 1e-5),  # the sign error made qext 21% low here
+        (1.5 - 0.1j, 0.05, 1e-5),
+        (1.5 - 1j, 0.1, 1e-4),  # |ratio|**2 in the x**4 term is 6e-4 high here
+        (0.05 - 1.43j, 0.05, 2e-3),  # silver-like, near the Frohlich resonance
+    ],
+)
+def test_extinction_matches_mie_for_absorbing_spheres(m, size, rtol):
+    """The 1% tolerance of test_efficiencies hides errors this size, so compare tightly."""
+    qext_ray = rayleigh.efficiencies_mx(m, size)[0]
+    qext_mie = mie.efficiencies_mx(m, size)[0]
+    assert qext_ray == pytest.approx(qext_mie, rel=rtol)
+
+
 @pytest.mark.parametrize("m, case", test_cases)
 def test_S1_S2(m, case):
     """Compare scattering amplitudes S1 and S2 to Mie theory."""