From aadfca56f00233e40ff09f997a3f331c0ee570ad Mon Sep 17 00:00:00 2001
From: John Ragland <john.ragland@whoi.edu>
Date: Thu, 21 May 2026 11:55:54 -0400
Subject: [PATCH] fix vertical ray bug

---
 src/pygenray/integration_processes.py |  7 ++++-
 tests/test_physics.py                 | 40 +++++++++++++++++++++++++++
 2 files changed, 46 insertions(+), 1 deletion(-)

diff --git a/src/pygenray/integration_processes.py b/src/pygenray/integration_processes.py
index 6ab7163..30876fd 100644
--- a/src/pygenray/integration_processes.py
+++ b/src/pygenray/integration_processes.py
@@ -85,7 +85,12 @@ def derivsrd(
     cp = bilinear_interp(x,z,rin,zin,cpin)
 
     # calculate derivatives
-    fact=1/np.sqrt(1-(c**2)*(pz**2))
+    # clamp to avoid division by zero when a Runge-Kutta intermediate stage
+    # lands on an exactly vertical ray before the vertical_ray event can fire
+    arg = 1.0 - (c**2) * (pz**2)
+    if arg <= 0.0:
+        arg = 1e-30
+    fact=1/np.sqrt(arg)
     dydx = np.array([
         fact/c,
         c*pz*fact,
diff --git a/tests/test_physics.py b/tests/test_physics.py
index 62cae0d..2f38ee8 100644
--- a/tests/test_physics.py
+++ b/tests/test_physics.py
@@ -322,3 +322,43 @@ def test_regression(self, request):
         np.testing.assert_array_equal(rf.n_surfs, ref['n_surfs'])
 
 
+# ---------------------------------------------------------------------------
+# F.  Near-vertical ray — regression for ZeroDivisionError in derivsrd
+# ---------------------------------------------------------------------------
+
+class TestNearVerticalRay:
+    """
+    Shooting a ray at or near 90° causes pz = sin(θ)/c ≈ 1/c, making
+    1 - c²·pz² → 0 in derivsrd. The vertical_ray event terminates
+    integration, but scipy evaluates derivsrd at RK intermediate stages
+    before checking events.  These tests verify no ZeroDivisionError is raised.
+    """
+
+    def test_near_vertical_no_crash(self):
+        """A ray launched at 89.9° (user: -89.9) must not raise ZeroDivisionError."""
+        env = _const_c_env(c0=1500.0, z_max=5000.0, r_max=100e3, bathy_depth=4500.0)
+        try:
+            ray = shoot_ray(200.0, 0.0, -89.9, 10e3, 50, env,
+                            rtol=1e-6, flatearth=False, debug=False)
+        except ZeroDivisionError:
+            pytest.fail("ZeroDivisionError raised for near-vertical ray (θ ≈ 89.9°)")
+
+    def test_exactly_vertical_no_crash(self):
+        """A ray launched at exactly 90° must not raise ZeroDivisionError."""
+        env = _const_c_env(c0=1500.0, z_max=5000.0, r_max=100e3, bathy_depth=4500.0)
+        try:
+            ray = shoot_ray(200.0, 0.0, -90.0, 10e3, 50, env,
+                            rtol=1e-6, flatearth=False, debug=False)
+        except ZeroDivisionError:
+            pytest.fail("ZeroDivisionError raised for exactly vertical ray (θ = 90°)")
+
+    @pytest.mark.parametrize("angle", [-85.0, -87.0, -89.0, -89.9, -90.0])
+    def test_steep_rays_no_crash(self, angle):
+        """Steep rays near-vertical must not crash regardless of launch angle."""
+        env = _munk_env(r_max=50e3)
+        try:
+            shoot_ray(1000.0, 0.0, angle, 10e3, 50, env,
+                      rtol=1e-6, flatearth=False, debug=False)
+        except ZeroDivisionError:
+            pytest.fail(f"ZeroDivisionError raised for steep ray at θ = {angle}°")
+