Fix negative-focal-length lens reflection; use ideal ray-slope thin-lens law (#102)

src/optiverse/raytracing/elements/lens.py

@@ -1,11 +1,9 @@
 """
 Lens element implementation.
 
-Implements an ideal thin lens using the exact (non-paraxial) deflection formula.
+Implements an ideal thin lens using the exact ray-slope (f·tanθ) deflection law.
 """
 
-import math
-
 import numpy as np
 
 from ...core.raytracing_math import normalize
@@ -15,13 +13,21 @@
 
 class LensElement(IOpticalElement):
     """
-    Ideal thin lens element (non-paraxial).
+    Ideal thin lens element.
+
+    Uses the ray-slope deflection law:  tan(θ_out) = tan(θ_in) − y/f
 
-    Uses the exact deflection formula: θ_out = θ_in - arctan(y/f)
-    which produces perfect focusing at all ray heights, not just near
-    the optical axis. The paraxial approximation θ_out = θ_in - y/f
-    is the small-angle linearization that introduces spurious aberration
-    at large impact parameters.
+    This is the exact transformation for an ideal (distortion-free) thin lens:
+    a bundle of parallel rays incident at any angle θ_in converges to a single
+    point in the focal plane at height f·tan(θ_in), for any ray height y. The
+    earlier θ_out = θ_in − arctan(y/f) form only focuses on-axis bundles and
+    smears off-axis ones; the plain-angle θ_out = θ_in − y/f form is the
+    small-angle linearization that aberrates at large ray heights.
+
+    Building the exit direction as (normal + slope·tangent) also keeps its
+    forward component positive, so the lens is always transmissive — it never
+    reflects, for either sign of f. This matches the web engine's thin-lens
+    deflection (GeometricEngine.thinLens / optiverse_engine _thin_lens).
     """
 
     def __init__(self, p1: np.ndarray, p2: np.ndarray, efl_mm: float):
@@ -45,12 +51,15 @@ def interact(
         self, ray: RayState, hit_point: np.ndarray, normal: np.ndarray, tangent: np.ndarray
     ) -> list[RayState]:
         """
-        Deflect ray using ideal (non-paraxial) thin lens equation.
+        Deflect ray using the ideal thin lens ray-slope law.
 
         Physics:
-        - Exact formula: θ_out = θ_in - arctan(y/f)
-        - This ensures a collimated beam at any height y converges
-          exactly to the focal point, without spherical-like aberration.
+        - Exact law: tan(θ_out) = tan(θ_in) − y/f, i.e. the exit ray slope
+          (measured from the lens normal) is reduced by y/f. A parallel bundle
+          at any angle focuses to a single point in the focal plane (f·tanθ).
+        - Constructing the exit direction as (normal + slope·tangent) keeps the
+          forward component positive, so the lens never reflects: a converging
+          lens focuses, a diverging (f < 0) lens diverges — never a mirror.
         - Polarization unchanged through ideal lens
         """
         # Ensure normal points in ray propagation direction
@@ -62,22 +71,18 @@ def interact(
         y = float(np.dot(hit_point - center, tangent))
 
         # Decompose ray direction into normal and tangent components
-        a_n = float(np.dot(ray.direction, normal))
-        a_t = float(np.dot(ray.direction, tangent))
-
-        # Compute incident angle
-        theta_in = math.atan2(a_t, a_n)
-
-        # Apply ideal thin lens equation: θ_out = θ_in - arctan(y/f)
-        # The arctan form is the exact angle subtended by height y at
-        # distance f, giving perfect focusing at all ray heights.
-        if abs(self.efl_mm) > 1e-12:
-            theta_out = theta_in - math.atan2(y, self.efl_mm)
+        a_n = float(np.dot(ray.direction, normal))  # cos(θ_in) ≥ 0 (normal flipped forward)
+        a_t = float(np.dot(ray.direction, tangent))  # sin(θ_in); a_t / a_n = tan(θ_in)
+
+        # Apply the ideal thin lens: tan(θ_out) = tan(θ_in) − y/f. The exit
+        # direction (normal + slope·tangent) always points forward, so the lens
+        # stays transmissive for either sign of f. Skip when the focal length is
+        # infinite or the ray grazes the lens (tan θ_in undefined) → pass through.
+        if abs(self.efl_mm) > 1e-12 and a_n >= 1e-9:
+            slope = a_t / a_n - y / self.efl_mm
+            direction_out = normalize(normal + slope * tangent)
         else:
-            theta_out = theta_in  # Infinite focal length = no deflection
-
-        # Reconstruct direction from angle
-        direction_out = normalize(math.cos(theta_out) * normal + math.sin(theta_out) * tangent)
+            direction_out = normalize(ray.direction)
 
         # Polarization unchanged through ideal lens
         EPS_ADV = 1e-3
@@ -120,14 +125,16 @@ def transform_q(
             center = 0.5 * (self.p1 + self.p2)
             y = float(np.dot(hit_point - center, tvec))
             a_t = float(np.dot(v, tvec))
-            theta_in = math.atan2(a_t, a_n)
-            theta_out = theta_in - math.atan2(y, f)
-            f_local = f * (1.0 + (y / f) ** 2)
-            direction_out = math.cos(theta_out) * n + math.sin(theta_out) * tvec
+            # Same ray-slope deflection as interact(), so the Gaussian ABCD uses
+            # the matching exit direction (mirrors web GeometricEngine.lensQ):
+            # tan(θ_out) = tan(θ_in) − y/f, direction = normal + slope·tangent.
+            slope = a_t / cos_theta_in - y / f
+            direction_out = n + slope * tvec
             norm_out = float(np.linalg.norm(direction_out))
             if norm_out > 1e-12:
                 direction_out = direction_out / norm_out
             cos_theta_out = max(abs(float(np.dot(direction_out, n))), 1e-12)
+            f_local = f * (1.0 + (y / f) ** 2)
             A = cos_theta_out / cos_theta_in
             C = -1.0 / (f_local * cos_theta_in)
             return apply_abcd(q, A, 0.0, C, 1.0)

tests/raytracing/test_optical_elements.py

@@ -234,6 +234,86 @@ def test_lens_off_axis_ray_deflected(self):
         assert refracted.direction[0] > 0  # Still going forward
         assert refracted.direction[1] < 0  # Deflected downward
 
+    def test_negative_focal_length_diverges_forward(self):
+        """Regression (#102): a negative-focal-length lens must stay transmissive.
+
+        A diverging lens should pass the ray forward and bend it *away* from the
+        optical axis. The ray-slope law tan(θ_out) = tan(θ_in) − y/f builds the
+        exit direction as (normal + slope·tangent), whose forward component is
+        always positive, so the lens can never reflect for either sign of f.
+        """
+        from optiverse.raytracing.elements import LensElement
+        from optiverse.raytracing.ray import Polarization, RayState
+
+        # Vertical diverging lens at x=0, f=-50mm
+        lens = LensElement(p1=np.array([0.0, -15.0]), p2=np.array([0.0, 15.0]), efl_mm=-50.0)
+
+        # Off-axis ray at y=10mm, parallel to axis
+        ray = RayState(
+            position=np.array([-10.0, 10.0]),
+            direction=np.array([1.0, 0.0]),
+            intensity=1.0,
+            polarization=Polarization.horizontal(),
+            wavelength_nm=633.0,
+            path=[],
+            events=0,
+        )
+
+        hit_point = np.array([0.0, 10.0])
+        normal = np.array([1.0, 0.0])
+        tangent = np.array([0.0, 1.0])
+
+        refracted = lens.interact(ray, hit_point, normal, tangent)[0]
+
+        # Must keep going forward (transmit), NOT reflect backward like a mirror.
+        assert refracted.direction[0] > 0
+        # Must diverge: an above-axis ray bends further away from the axis (upward).
+        assert refracted.direction[1] > 0
+        # Virtual extension crosses the axis on the incident side at distance |f|:
+        # the outgoing slope magnitude equals |y / f| = 10 / 50 = 0.2.
+        slope = refracted.direction[1] / refracted.direction[0]
+        assert slope == pytest.approx(10.0 / 50.0, rel=1e-6)
+
+    def test_oblique_parallel_bundle_focuses_to_point(self):
+        """Ideal-lens law: a parallel bundle at ANY angle focuses to one point.
+
+        With tan(θ_out) = tan(θ_in) − y/f, rays of a common input angle but
+        different heights y all cross the back focal plane (x = f) at the same
+        height f·tan(θ_in). This f·tanθ focusing is what makes the ray-slope law
+        exact; the older θ_out = θ_in − arctan(y/f) form smeared off-axis bundles.
+        """
+        from optiverse.raytracing.elements import LensElement
+        from optiverse.raytracing.ray import Polarization, RayState
+
+        f = 100.0
+        lens = LensElement(p1=np.array([0.0, -40.0]), p2=np.array([0.0, 40.0]), efl_mm=f)
+        normal = np.array([1.0, 0.0])
+        tangent = np.array([0.0, 1.0])
+
+        theta_in = 0.30  # oblique bundle (~17°)
+        in_dir = np.array([math.cos(theta_in), math.sin(theta_in)])
+
+        focal_plane_heights = []
+        for y in (-20.0, -10.0, 0.0, 10.0, 20.0):
+            ray = RayState(
+                position=np.array([-10.0, y]),
+                direction=in_dir,
+                intensity=1.0,
+                polarization=Polarization.horizontal(),
+                wavelength_nm=633.0,
+                path=[],
+                events=0,
+            )
+            hit_point = np.array([0.0, y])
+            out = lens.interact(ray, hit_point, normal, tangent)[0]
+            # Propagate the exit ray from the lens (x=0) to the focal plane x=f.
+            t = f / out.direction[0]
+            focal_plane_heights.append(hit_point[1] + t * out.direction[1])
+
+        target = f * math.tan(theta_in)
+        for h in focal_plane_heights:
+            assert h == pytest.approx(target, abs=1e-6)
+
 
 class TestRefractiveElement:
     """Test RefractiveElement implementation"""