Normalize scattering CDFs before inverse angle sampling

miepython/monte_carlo.py

@@ -52,21 +52,25 @@ def cdf(m, x, num):
     sampled uniformly and therefore huge array sizes are needed to adequately
     sample highly anisotropic phase functions.
 
-    Since this is a cumulative distribution function, the maximum value
-    should be 1.
+    The cumulative trapezoidal integral starts at zero and is normalized to
+    end at one. Increase num to resolve a sharply peaked phase function.
 
     Args:
         m: the complex index of refraction of the sphere
         x: the size parameter of the sphere
-        num: length of desired CDF array
+        num: length of desired CDF array, at least two
 
     Returns:
         mu: array of cosines of angles
         cdf: array of cumulative distribution function values
     """
+    if num < 2:
+        raise ValueError("num must be at least 2")
     mu = np.linspace(-1, 1, num)
-    intensity_per_mu = mie.i_unpolarized(m, x, mu, norm="4pi") / num
-    cdf_values = np.cumsum(intensity_per_mu)
+    intensity_per_mu = mie.i_unpolarized(m, x, mu, norm="4pi") / 2
+    cdf_values = np.zeros(num)
+    cdf_values[1:] = np.cumsum((intensity_per_mu[1:] + intensity_per_mu[:-1]) * np.diff(mu) / 2)
+    cdf_values /= cdf_values[-1]
     return mu, cdf_values
 
 
@@ -89,12 +93,14 @@ def mu_with_uniform_cdf(m, x, num):
     Args:
         m: the complex index of refraction of the sphere
         x: the size parameter of the sphere
-        num: length of desired CDF array
+        num: length of desired CDF array, at least two
 
     Returns:
         mu: array of cosines of angles (irregularly spaced)
         cdf: array of cumulative distribution function values
     """
+    if num < 2:
+        raise ValueError("num must be at least 2")
     big_num = 2000  # large to work with x up to 10
     big_mu, big_cdf = cdf(m, x, big_num)
     mu = np.empty(num)

tests/test_monte_carlo.py

@@ -42,6 +42,12 @@ def draw(mu_table, count, seed):
 class TestCdf:
     """Test the cumulative distribution sampled uniformly in cos(theta)."""
 
+    def test_small_sphere_matches_rayleigh_cdf(self):
+        """The small-sphere phase density is 3/8 * (1 + mu**2)."""
+        mu, cdf = mc.cdf(M_WATER, 0.001, 100)
+        expected = 0.5 + 3 * mu / 8 + mu**3 / 8
+        np.testing.assert_allclose(cdf, expected, atol=2e-5)
+
     @pytest.mark.parametrize("num", [20, 50, 500])
     def test_shape_and_angles(self, num):
         """The angles span the full exit range at even spacing."""
@@ -55,27 +61,20 @@ def test_is_non_decreasing_and_starts_near_zero(self, num):
         """A cumulative distribution may never run backwards."""
         _, cdf = mc.cdf(M_WATER, X_DROPLET, num)
         assert np.all(np.diff(cdf) >= 0)
-        assert 0 < cdf[0] < 2.0 / num
-
-    def test_endpoint_approaches_one_from_above(self):
-        """The quadrature overshoots, and the error falls off like 1/num.
+        assert cdf[0] == 0.0
 
-        ``cdf`` sums ``i_unpolarized / num`` rather than multiplying by the true
-        spacing 2/(num-1), so the last value is not exactly 1: it is about
-        1 + 1.3/num for this sphere.  The docstring claims a maximum of 1, so the
-        overshoot is worth pinning rather than leaving as a surprise -- at the
-        num=20 the notebook plots, the CDF reaches 1.067.
-        """
-        errors = {}
-        for num in (50, 500, 5000):
-            _, cdf = mc.cdf(M_WATER, X_DROPLET, num)
-            errors[num] = cdf[-1] - 1.0
-            assert errors[num] > 0, "overshoot expected"
+    @pytest.mark.parametrize("num", [2, 20, 50, 500, 5000])
+    def test_endpoint_is_one(self, num):
+        """A CDF contains all the probability at its last angle."""
+        _, cdf = mc.cdf(M_WATER, X_DROPLET, num)
+        assert cdf[-1] == 1.0
 
-        # ten times as many points, a tenth of the error
-        assert errors[50] / errors[500] == pytest.approx(10.0, rel=0.1)
-        assert errors[500] / errors[5000] == pytest.approx(10.0, rel=0.1)
-        assert errors[5000] < 1e-3
+    @pytest.mark.parametrize("num", [0, 1])
+    @pytest.mark.parametrize("builder", [mc.cdf, mc.mu_with_uniform_cdf])
+    def test_at_least_two_points_are_required(self, num, builder):
+        """One angle cannot define an interval of scattering probability."""
+        with pytest.raises(ValueError, match="at least 2"):
+            builder(M_WATER, X_DROPLET, num)
 
     def test_more_forward_scattering_delays_the_rise(self):
         """A larger sphere throws light forward, so the cdf climbs later."""
@@ -88,6 +87,14 @@ def test_more_forward_scattering_delays_the_rise(self):
 class TestUniformCdfTable:
     """Test the inverse-transform table used for sampling."""
 
+    def test_rayleigh_inverse_angles_are_physical(self):
+        """Fine tables must not interpolate across the two CDF endpoints."""
+        mu, _ = mc.mu_with_uniform_cdf(M_WATER, 0.001, 4096)
+        assert np.all(mu >= -1)
+        assert np.all(mu <= 1)
+        assert np.all(np.diff(mu) >= 0)
+        assert np.all(np.isfinite(np.arccos(mu)))
+
     @pytest.mark.parametrize("num", [20, 200, 3000])
     def test_cdf_values_are_exactly_uniform(self, num):
         """That is the whole point: equal probability per interval."""