From 415f05c3fe3f5b484068f70eded92ab7ef121bc8 Mon Sep 17 00:00:00 2001 From: Anton Soloviev Date: Tue, 27 Oct 2026 16:45:44 -0700 Subject: [PATCH 1/2] Add scattering after correcting Rayleigh absorption sign --- miepython/rayleigh.py | 2 +- tests/test_rayleigh.py | 17 +++++++++++++++++ 2 files changed, 18 insertions(+), 1 deletion(-) diff --git a/miepython/rayleigh.py b/miepython/rayleigh.py index 86d7035..34ebab3 100644 --- a/miepython/rayleigh.py +++ b/miepython/rayleigh.py @@ -35,7 +35,7 @@ def efficiencies_mx(m, x): 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) + qext = abs(qext.imag) + qsca qback = 4 * x**4 * abs(ratio) ** 2 g = 0 return qext, qsca, qback, g diff --git a/tests/test_rayleigh.py b/tests/test_rayleigh.py index fb8127b..3b8753b 100644 --- a/tests/test_rayleigh.py +++ b/tests/test_rayleigh.py @@ -37,6 +37,23 @@ 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, case", test_cases) def test_S1_S2(m, case): """Compare scattering amplitudes S1 and S2 to Mie theory.""" From 12e63820ab54f3295bec349a24789f125a9bbee9 Mon Sep 17 00:00:00 2001 From: Scott Prahl Date: Wed, 28 Oct 2026 10:19:54 -0700 Subject: [PATCH 2/2] Use Re(ratio**2) and conjugate m in Rayleigh extinction Follow-up to the Rayleigh absorption sign fix: - conjugate m up front, as core.py does, so absorption is simply -Im(...) instead of abs() hiding the sign - use Bohren & Huffman's Re(ratio**2) for the x**4 extinction term rather than |ratio|**2; at m = 1.5 - 1j, x = 0.1 the error against Mie falls from 6e-4 to 2e-5, and near resonance (m = 0.05 - 1.43j, x = 0.05) from 0.43% to 0.12% - compare qext tightly against mie.efficiencies_mx for weak, moderate, strong and near-resonant absorption; the 1% tolerance of test_efficiencies hid both the sign bug and the |ratio|**2 error - changelog entry crediting @antonsoo Co-Authored-By: Claude Opus 5.5 (1M context) --- CHANGELOG.rst | 9 +++++++++ miepython/rayleigh.py | 13 +++++++++++-- tests/test_rayleigh.py | 16 ++++++++++++++++ 3 files changed, 36 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 97d5a4e..79b8fc5 100644 --- a/CHANGELOG.rst +++ b/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 diff --git a/miepython/rayleigh.py b/miepython/rayleigh.py index 34ebab3..123531d 100644 --- a/miepython/rayleigh.py +++ b/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 diff --git a/tests/test_rayleigh.py b/tests/test_rayleigh.py index 3b8753b..0e2ef7c 100644 --- a/tests/test_rayleigh.py +++ b/tests/test_rayleigh.py @@ -54,6 +54,22 @@ def test_weak_absorption_is_positive(index, absorption, size): 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."""