Repository: HarrisonKramer/optiland

Issue:
Add an operand to control clearances during optimization
## Description

Implement a `clearance` operand that allows to control clearances during optimization.
Here's how it should work:

The function computes the signed perpendicular distance (in the YZ plane) from a line A to a point B. 
This is useful for creating clearance or interference constraints when you are optimizing
off-axis reflective systems, such as TMA configurations.

The line A is a reference ray (ray RA) at field FA after it leaves a specified surface SA, and the point B is the
intersection of any reference ray (ray RB) at field FB with a given surface SB.
The sign convention is such that for a ray (line A) propagating generally in the +Z direction (N direction cosine > 0),
the signed distance is positive if the “point” is on the +Y side of the line. 
Therefore the distance is negative in the following figure:

![Image](https://github.com/user-attachments/assets/a7e90165-aa63-44c1-8174-ea4175ad22e5)

## Requirements

- [ ] Identify or write a function to easily retreive coordinates of a designed ray on a surface
- [ ] Compute the equation of "Line A"
- [ ] Compute the distance "Line A to Point B"
- [ ] Create the operand that constraints this distance to a target
- [ ] Add documentation
- [ ] Add unit tests.

Hints from issue comments:

Solution patch:
diff --git a/optiland/optimization/operand/ray.py b/optiland/optimization/operand/ray.py
index adca3311..c0dc512f 100644
--- a/optiland/optimization/operand/ray.py
+++ b/optiland/optimization/operand/ray.py
@@ -323,3 +323,77 @@ class RayOperand:
         opd = wavefront_data.opd
         delta = (opd - be.mean(opd)) * weights
         return be.mean(be.abs(delta))
+
+    @staticmethod
+    def clearance(
+        optic,
+        line_ray_surface_idx,
+        line_ray_field_coords,
+        line_ray_pupil_coords,
+        point_ray_surface_idx,
+        point_ray_field_coords,
+        point_ray_pupil_coords,
+        wavelength,
+    ):
+        """Computes the signed perpendicular distance in the YZ plane from a
+        reference line (Line A) to a reference point (Point B).
+
+        Line A is defined by a ray (RA) traced at field FA, after it leaves
+        surface SA. Point B is the intersection of a ray (RB) traced at
+        field FB with surface SB.
+
+        This operand is useful for creating clearance or interference constraints,
+        particularly in off-axis reflective systems.
+
+        The sign convention is such that for Line A propagating generally in the
+        +Z direction (N direction cosine > 0), the signed distance is positive
+        if Point B is on the +Y side of Line A. If Line A propagates generally
+        in the -Z direction (N direction cosine < 0), this sign is flipped.
+
+        Args:
+            optic: The optical system model.
+            line_ray_surface_idx: The index of the surface (SA) from which
+                Line A originates (i.e., ray data is taken *after* this
+                surface).
+            line_ray_field_coords: A tuple (Hx, Hy) representing the
+                normalized field coordinates for the ray defining Line A (FA).
+            line_ray_pupil_coords: A tuple (Px, Py) representing the
+                normalized pupil coordinates for the ray defining Line A (FA).
+            point_ray_surface_idx: The index of the surface (SB) with which
+                the ray defining Point B intersects.
+            point_ray_field_coords: A tuple (Hx, Hy) representing the
+                normalized field coordinates for the ray defining Point B (FB).
+            point_ray_pupil_coords: A tuple (Px, Py) representing the
+                normalized pupil coordinates for the ray defining Point B (FB).
+            wavelength: The wavelength at which to trace the rays.
+
+        Returns:
+            float: The signed perpendicular distance in the YZ plane from
+                   Line A to Point B. Returns 0.0 if Line A has zero length
+                   in the YZ plane (i.e., mA and nA are both zero).
+        """
+        FA_Hx, FA_Hy = line_ray_field_coords
+        FA_Px, FA_Py = line_ray_pupil_coords
+        optic.trace_generic(FA_Hx, FA_Hy, FA_Px, FA_Py, wavelength)
+        yA = optic.surface_group.y[line_ray_surface_idx, 0]
+        zA = optic.surface_group.z[line_ray_surface_idx, 0]
+        mA = optic.surface_group.M[line_ray_surface_idx, 0]
+        nA = optic.surface_group.N[line_ray_surface_idx, 0]
+
+        FB_Hx, FB_Hy = point_ray_field_coords
+        FB_Px, FB_Py = point_ray_pupil_coords
+        optic.trace_generic(FB_Hx, FB_Hy, FB_Px, FB_Py, wavelength)
+        yB = optic.surface_group.y[point_ray_surface_idx, 0]
+        zB = optic.surface_group.z[point_ray_surface_idx, 0]
+
+        denominator = be.sqrt(mA**2 + nA**2)
+        epsilon = 1e-9
+
+        if be.abs(denominator) < epsilon:
+            d = 0.0
+        else:
+            numerator = nA * (yB - yA) - mA * (zB - zA)
+            d = numerator / denominator
+            if nA < 0:
+                d = -d
+        return d


Test patch:
diff --git a/tests/test_aberrations.py b/tests/test_aberrations.py
index e65f14e8..02f75da3 100644
--- a/tests/test_aberrations.py
+++ b/tests/test_aberrations.py
@@ -62,8 +62,16 @@ class TestDoubleGaussAberrations:
 
     def test_seidels(self, set_test_backend, double_gauss):
         S = double_gauss.aberrations.seidels()
-        assert_allclose(S, [-0.003929457875534847, 0.0003954597633218682, 0.0034239055031729947, 
-                            -0.016264753735226404, -0.046484107476755930])
+        assert_allclose(
+            S,
+            [
+                -0.003929457875534847,
+                0.0003954597633218682,
+                0.0034239055031729947,
+                -0.016264753735226404,
+                -0.046484107476755930,
+            ],
+        )
 
     def test_third_order(self, set_test_backend, double_gauss):
         data = double_gauss.aberrations.third_order()
diff --git a/tests/test_analysis.py b/tests/test_analysis.py
index 4f2f0741..bac072d0 100644
--- a/tests/test_analysis.py
+++ b/tests/test_analysis.py
@@ -744,7 +744,6 @@ def test_spotdiagram_invalid_coordinates(cooke_triplet):
 
 
 def test_generate_field_data_local(set_test_backend, cooke_triplet):
-
     spot = analysis.SpotDiagram(cooke_triplet, coordinates="local")
 
     # Pick the first field and wavelength
@@ -769,7 +768,6 @@ def test_generate_field_data_local(set_test_backend, cooke_triplet):
 
 
 def test_generate_field_data_global(set_test_backend, cooke_triplet):
-
     spot = analysis.SpotDiagram(cooke_triplet, coordinates="global")
 
     # Pick the first field and wavelength
@@ -1168,13 +1166,13 @@ class TestIncoherentIrradiance:
         # Invalid slice index for cross-x
         irr.view(cross_section=("cross-x", res_val[0] + 5))  # Index out of bounds
         mock_print.assert_any_call(
-            f"[IncoherentIrradiance] Warning: X-slice index {res_val[0]+5} is out of bounds for map shape {(res_val[0],res_val[1])}. Skipping plot."
+            f"[IncoherentIrradiance] Warning: X-slice index {res_val[0] + 5} is out of bounds for map shape {(res_val[0], res_val[1])}. Skipping plot."
         )
 
         # Invalid slice index for cross-y
         irr.view(cross_section=("cross-y", res_val[1] + 5))  # Index out of bounds
         mock_print.assert_any_call(
-            f"[IncoherentIrradiance] Warning: Y-slice index {res_val[1]+5} is out of bounds for map shape {(res_val[0],res_val[1])}. Skipping plot."
+            f"[IncoherentIrradiance] Warning: Y-slice index {res_val[1] + 5} is out of bounds for map shape {(res_val[0], res_val[1])}. Skipping plot."
         )
 
         # Invalid cross_section_info format (not tuple)
diff --git a/tests/test_aperture.py b/tests/test_aperture.py
index f2169acd..5c99e900 100644
--- a/tests/test_aperture.py
+++ b/tests/test_aperture.py
@@ -8,7 +8,12 @@ from .utils import assert_allclose
 
 @pytest.mark.parametrize(
     "ap_type, value",
-    [("EPD", 10), ("imageFNO", 3.2), ("objectNA", 0.265), ("float_by_stop_size", 1.142857)],
+    [
+        ("EPD", 10),
+        ("imageFNO", 3.2),
+        ("objectNA", 0.265),
+        ("float_by_stop_size", 1.142857),
+    ],
 )
 def test_aperture_generate(set_test_backend, ap_type, value):
     """Check instantiation of aperture"""
diff --git a/tests/test_coordinate_system.py b/tests/test_coordinate_system.py
index 09cc5b20..4a6059cf 100644
--- a/tests/test_coordinate_system.py
+++ b/tests/test_coordinate_system.py
@@ -46,12 +46,12 @@ def test_coordinate_system_localize(set_test_backend):
     cs = CoordinateSystem(10, 20, 30, 0.5, 0.6, 0.7, ref_cs)
     rays = RealRays(1, 2, 3, 0.1, 0.2, 0.3, 1, 1)
     cs.localize(rays)
-    assert_allclose(rays.x, -1.826215)       
-    assert_allclose(rays.y, -26.51361)      
-    assert_allclose(rays.z, -32.55361)      
-    assert_allclose(rays.L, -0.0122750)     
-    assert_allclose(rays.M, 0.2697627)      
-    assert_allclose(rays.N, 0.2589931)      
+    assert_allclose(rays.x, -1.826215)
+    assert_allclose(rays.y, -26.51361)
+    assert_allclose(rays.z, -32.55361)
+    assert_allclose(rays.L, -0.0122750)
+    assert_allclose(rays.M, 0.2697627)
+    assert_allclose(rays.N, 0.2589931)
 
 
 def test_coordinate_system_globalize(set_test_backend):
@@ -80,12 +80,12 @@ def test_coordinate_system_globalize(set_test_backend):
         1,
     )
     cs.globalize(rays)
-    assert_allclose(rays.x, 4.654692)        
-    assert_allclose(rays.y, -12.50279)      
-    assert_allclose(rays.z, 21.51751)       
-    assert_allclose(rays.L, 0.0866579)      
-    assert_allclose(rays.M, 0.3526931)      
-    assert_allclose(rays.N, 0.08998877)    
+    assert_allclose(rays.x, 4.654692)
+    assert_allclose(rays.y, -12.50279)
+    assert_allclose(rays.z, 21.51751)
+    assert_allclose(rays.L, 0.0866579)
+    assert_allclose(rays.M, 0.3526931)
+    assert_allclose(rays.N, 0.08998877)
 
 
 def test_coordinate_system_transform(set_test_backend):
diff --git a/tests/test_distribution.py b/tests/test_distribution.py
index b7ebac47..aae61e2c 100644
--- a/tests/test_distribution.py
+++ b/tests/test_distribution.py
@@ -85,15 +85,19 @@ def test_cross(set_test_backend, num_points):
     x_line_x_expected_full = be.linspace(-1, 1, num_points)
     x_line_y_expected_full = be.zeros(num_points)
 
-    if num_points % 2 == 1: # Odd number of points
+    if num_points % 2 == 1:  # Odd number of points
         # Remove the middle element from the x-axis line as it's the duplicated origin
         mid_idx = num_points // 2
-        x_line_x_to_concat = be.concatenate((x_line_x_expected_full[:mid_idx], x_line_x_expected_full[mid_idx+1:]))
-        x_line_y_to_concat = be.concatenate((x_line_y_expected_full[:mid_idx], x_line_y_expected_full[mid_idx+1:]))
-    else: # Even number of points (origin is not in the middle of linspace for an odd-length array)
+        x_line_x_to_concat = be.concatenate(
+            (x_line_x_expected_full[:mid_idx], x_line_x_expected_full[mid_idx + 1 :])
+        )
+        x_line_y_to_concat = be.concatenate(
+            (x_line_y_expected_full[:mid_idx], x_line_y_expected_full[mid_idx + 1 :])
+        )
+    else:  # Even number of points (origin is not in the middle of linspace for an odd-length array)
         x_line_x_to_concat = x_line_x_expected_full
         x_line_y_to_concat = x_line_y_expected_full
-    
+
     # Concatenate in the same order as in the implementation
     expected_x = be.concatenate((y_line_x_expected, x_line_x_to_concat))
     expected_y = be.concatenate((y_line_y_expected, x_line_y_to_concat))
diff --git a/tests/test_distribution_points.py b/tests/test_distribution_points.py
index 183b1eb4..be46a54a 100644
--- a/tests/test_distribution_points.py
+++ b/tests/test_distribution_points.py
@@ -10,7 +10,7 @@ from optiland.distribution import (
     RingDistribution,
     LineXDistribution,
     LineYDistribution,
-    GaussianQuadrature
+    GaussianQuadrature,
 )
 
 
@@ -27,15 +27,17 @@ def count_unique_points(x_coords, y_coords, tol=1e-9):
     # The number of decimal places can be derived from tol
     # e.g., tol=1e-9 means 9 decimal places
     decimals = int(-np.log10(tol))
-    
-    points = np.stack((np.round(x_coords_writable, decimals), np.round(y_coords_writable, decimals)), axis=-1)
+
+    points = np.stack(
+        (np.round(x_coords_writable, decimals), np.round(y_coords_writable, decimals)),
+        axis=-1,
+    )
     unique_points = np.unique(points, axis=0)
     return unique_points.shape[0]
 
 
 @pytest.mark.parametrize("backend", ["numpy", "torch"])
 class TestDistributionPoints:
-
     def test_cross_distribution_odd_points(self, backend):
         be.set_backend(backend)
         num_points = 3  # Odd number
@@ -43,8 +45,9 @@ class TestDistributionPoints:
         dist.generate_points(num_points=num_points)
         expected_total_points = 2 * num_points - 1  # Origin counted once
         unique_count = count_unique_points(dist.x, dist.y)
-        assert unique_count == expected_total_points, \
+        assert unique_count == expected_total_points, (
             f"CrossDistribution (odd) expected {expected_total_points} unique points, got {unique_count}"
+        )
 
     def test_cross_distribution_even_points(self, backend):
         be.set_backend(backend)
@@ -54,7 +57,7 @@ class TestDistributionPoints:
         # For even num_points > 0, origin (0,0) is NOT generated by EITHER linspace.
         # So total points = num_points (vertical) + num_points (horizontal) = 2 * num_points.
         expected_total_points = 2 * num_points
-        if num_points == 0: # Edge case
+        if num_points == 0:  # Edge case
             expected_total_points = 0
             # When num_points is 0, generate_points will result in empty x, y
             # y_line_x/y will be empty. x_line_x/y will be empty.
@@ -62,9 +65,9 @@ class TestDistributionPoints:
             # So dist.x/y will be concatenation of empty arrays, which are empty.
 
         unique_count = count_unique_points(dist.x, dist.y)
-        assert unique_count == expected_total_points, \
+        assert unique_count == expected_total_points, (
             f"CrossDistribution (even) expected {expected_total_points} unique points, got {unique_count}"
-
+        )
 
     @pytest.mark.parametrize("num_rings", [0, 1, 2, 3])
     def test_hexagonal_distribution(self, backend, num_rings):
@@ -76,26 +79,28 @@ class TestDistributionPoints:
         else:
             # Formula: 1 (center) + sum_{i=1 to num_rings} 6*i = 1 + 6 * num_rings * (num_rings+1) / 2
             expected_total_points = 1 + 3 * num_rings * (num_rings + 1)
-        
+
         unique_count = count_unique_points(dist.x, dist.y)
-        assert unique_count == expected_total_points, \
+        assert unique_count == expected_total_points, (
             f"HexagonalDistribution ({num_rings} rings) expected {expected_total_points} unique points, got {unique_count}"
+        )
 
     def test_random_distribution(self, backend):
         be.set_backend(backend)
         num_points = 100
-        dist = RandomDistribution(seed=42) # Seed for reproducibility
+        dist = RandomDistribution(seed=42)  # Seed for reproducibility
         dist.generate_points(num_points=num_points)
         # For random, exact count might be tricky if true duplicates are astronomically rare
         # but we check if the count matches num_points, assuming no collisions for typical N.
         # The test primarily ensures the generation logic doesn't inherently create duplicates.
         unique_count = count_unique_points(dist.x, dist.y)
-        assert unique_count == num_points, \
+        assert unique_count == num_points, (
             f"RandomDistribution expected {num_points} unique points, got {unique_count}"
+        )
 
     def test_uniform_distribution(self, backend):
         be.set_backend(backend)
-        num_grid_points = 5 # This means 5x5 grid initially for the underlying meshgrid
+        num_grid_points = 5  # This means 5x5 grid initially for the underlying meshgrid
         dist = UniformDistribution()
         dist.generate_points(num_points=num_grid_points)
         # The actual number of points is num_grid_points^2 filtered by r2 <= 1
@@ -104,11 +109,11 @@ class TestDistributionPoints:
         # So, number of generated points should equal number of unique points.
         num_generated_points = len(dist.x)
         unique_count = count_unique_points(dist.x, dist.y)
-        assert unique_count == num_generated_points, \
+        assert unique_count == num_generated_points, (
             f"UniformDistribution expected {num_generated_points} unique points, got {unique_count}"
+        )
         assert num_generated_points > 0, "UniformDistribution generated no points"
 
-
     def test_ring_distribution(self, backend):
         be.set_backend(backend)
         num_points = 10
@@ -116,45 +121,55 @@ class TestDistributionPoints:
         dist.generate_points(num_points=num_points)
         expected_total_points = num_points
         unique_count = count_unique_points(dist.x, dist.y)
-        assert unique_count == expected_total_points, \
+        assert unique_count == expected_total_points, (
             f"RingDistribution expected {expected_total_points} unique points, got {unique_count}"
+        )
 
     def test_line_x_distribution(self, backend):
         be.set_backend(backend)
         num_points = 5
-        dist = LineXDistribution() # positive_only=False by default
+        dist = LineXDistribution()  # positive_only=False by default
         dist.generate_points(num_points=num_points)
         expected_total_points = num_points
         unique_count = count_unique_points(dist.x, dist.y)
-        assert unique_count == expected_total_points, \
+        assert unique_count == expected_total_points, (
             f"LineXDistribution expected {expected_total_points} unique points, got {unique_count}"
+        )
 
     def test_line_y_distribution(self, backend):
         be.set_backend(backend)
         num_points = 5
-        dist = LineYDistribution() # positive_only=False by default
+        dist = LineYDistribution()  # positive_only=False by default
         dist.generate_points(num_points=num_points)
         expected_total_points = num_points
         unique_count = count_unique_points(dist.x, dist.y)
-        assert unique_count == expected_total_points, \
+        assert unique_count == expected_total_points, (
             f"LineYDistribution expected {expected_total_points} unique points, got {unique_count}"
+        )
 
-    @pytest.mark.parametrize("num_points", [1, 3, 5]) # Corresponds to radius_dict keys in GaussianQuadrature's _get_radius
+    @pytest.mark.parametrize(
+        "num_points", [1, 3, 5]
+    )  # Corresponds to radius_dict keys in GaussianQuadrature's _get_radius
     @pytest.mark.parametrize("is_symmetric", [True, False])
-    def test_gaussian_quadrature_distribution(self, backend, num_points, is_symmetric): # num_points here is num_rings for GaussianQuad
+    def test_gaussian_quadrature_distribution(
+        self, backend, num_points, is_symmetric
+    ):  # num_points here is num_rings for GaussianQuad
         be.set_backend(backend)
         dist = GaussianQuadrature(is_symmetric=is_symmetric)
-        dist.generate_points(num_rings=num_points) # num_points in test param is num_rings for method
-        
+        dist.generate_points(
+            num_rings=num_points
+        )  # num_points in test param is num_rings for method
+
         # Expected points calculation based on GaussianQuadrature logic
         # The number of points is determined by the definition of Gaussian Quadrature
         # and not straightforwardly by num_rings like other distributions.
         # For GQ, all generated points should inherently be unique by definition.
         num_generated_points = len(dist.x)
-        
+
         unique_count = count_unique_points(dist.x, dist.y)
-        assert unique_count == num_generated_points, \
-             f"GaussianQuadrature ({num_points} rings, symmetric={is_symmetric}) expected all {num_generated_points} generated points to be unique, but got {unique_count}"
+        assert unique_count == num_generated_points, (
+            f"GaussianQuadrature ({num_points} rings, symmetric={is_symmetric}) expected all {num_generated_points} generated points to be unique, but got {unique_count}"
+        )
         assert num_generated_points > 0, "GaussianQuadrature generated no points"
 
     def test_cross_distribution_zero_points(self, backend):
@@ -168,8 +183,9 @@ class TestDistributionPoints:
         # So, dist.x and dist.y will be empty.
         expected_total_points = 0
         unique_count = count_unique_points(dist.x, dist.y)
-        assert unique_count == expected_total_points, \
+        assert unique_count == expected_total_points, (
             f"CrossDistribution (zero points) expected {expected_total_points} unique points, got {unique_count}"
+        )
         assert len(dist.x) == 0
         assert len(dist.y) == 0
 
@@ -181,10 +197,11 @@ class TestDistributionPoints:
         # If num_rings < 0, range(num_rings) is empty.
         # It initializes x = be.zeros(1), y = be.zeros(1). Loop doesn't run. So 1 point (origin).
         dist.generate_points(num_rings=num_rings)
-        expected_total_points = 1 # Center point only
+        expected_total_points = 1  # Center point only
         unique_count = count_unique_points(dist.x, dist.y)
-        assert unique_count == expected_total_points, \
+        assert unique_count == expected_total_points, (
             f"HexagonalDistribution ({num_rings} rings) expected {expected_total_points} unique point (origin), got {unique_count}"
+        )
 
     def test_ring_distribution_zero_points(self, backend):
         be.set_backend(backend)
@@ -195,7 +212,8 @@ class TestDistributionPoints:
         # So cos([]) = [], sin([]) = []. dist.x/y will be empty.
         expected_total_points = 0
         unique_count = count_unique_points(dist.x, dist.y)
-        assert unique_count == expected_total_points, \
+        assert unique_count == expected_total_points, (
             f"RingDistribution (0 points) expected {expected_total_points} unique points, got {unique_count}"
+        )
         assert len(dist.x) == 0
         assert len(dist.y) == 0
diff --git a/tests/test_fileio.py b/tests/test_fileio.py
index 711f6628..1afc9ded 100644
--- a/tests/test_fileio.py
+++ b/tests/test_fileio.py
@@ -37,16 +37,21 @@ def zemax_file_reader(zemax_file):
 def zemax_file_formats_reader(zemax_file_formats):
     return ZemaxFileReader(zemax_file_formats)
 
+
 @pytest.fixture
 def zemax_file_fold_mirrors():
     current_dir = os.path.dirname(os.path.abspath(__file__))
-    filename = os.path.join(current_dir, "zemax_files/complicated_fold_mirrors_setup_v2.zmx")
+    filename = os.path.join(
+        current_dir, "zemax_files/complicated_fold_mirrors_setup_v2.zmx"
+    )
     return filename
 
+
 @pytest.fixture
 def zemax_fold_mirrors_reader(zemax_file_fold_mirrors):
     return ZemaxFileReader(zemax_file_fold_mirrors)
 
+
 def replace_line_in_zmx(zmx_file, line_prefix, new_line):
     temp_dir = tempfile.mkdtemp()
     temp_file = os.path.join(temp_dir, "temp.zmx")
@@ -99,7 +104,9 @@ class TestZemaxFileReader:
         assert zemax_file_reader.filename == "temp.zmx"
 
     @patch("builtins.open", new_callable=mock_open)
-    def test_configure_source_input_file(self, mock_open, set_test_backend, zemax_file_formats_reader):
+    def test_configure_source_input_file(
+        self, mock_open, set_test_backend, zemax_file_formats_reader
+    ):
         zemax_file_formats_reader.source = "local_file.zmx"
         zemax_file_formats_reader._configure_source_input()
 
@@ -167,7 +174,9 @@ class TestZemaxFileReader:
         zemax_file_reader._read_config_data(data)
         assert zemax_file_reader.data["fields"]["type"] == "object_height"
 
-    def test_read_field_type_parax_img_height(self, set_test_backend, zemax_file_reader):
+    def test_read_field_type_parax_img_height(
+        self, set_test_backend, zemax_file_reader
+    ):
         data = ["FTYP", "2", "0", "0", "0", "0", "0", "0"]
         zemax_file_reader._read_config_data(data)
         assert zemax_file_reader.data["fields"]["type"] == "paraxial_image_height"
@@ -331,36 +340,57 @@ class TestZemaxToOpticConverter:
         zemax_file_reader.data["fields"]["vignette_decenter_y"] = [0.5]
         zemax_file_reader.generate_lens()
 
-    def test_read_surface_coordinate_break(self, set_test_backend, zemax_fold_mirrors_reader):
+    def test_read_surface_coordinate_break(
+        self, set_test_backend, zemax_fold_mirrors_reader
+    ):
         data = ["TYPE", "COORDBRK"]
-        zemax_fold_mirrors_reader._current_surf_data = {} 
+        zemax_fold_mirrors_reader._current_surf_data = {}
         zemax_fold_mirrors_reader._read_surf_type(data)
-        assert zemax_fold_mirrors_reader._current_surf_data["type"] == "coordinate_break"
+        assert (
+            zemax_fold_mirrors_reader._current_surf_data["type"] == "coordinate_break"
+        )
 
     def test_read_glass_mirror(self, set_test_backend, zemax_fold_mirrors_reader):
-        zemax_fold_mirrors_reader._current_surf_data = {} 
-        data = ["GLAS", "MIRROR", "0", "0", "0", "0"] 
+        zemax_fold_mirrors_reader._current_surf_data = {}
+        data = ["GLAS", "MIRROR", "0", "0", "0", "0"]
         zemax_fold_mirrors_reader._read_glass(data)
         assert zemax_fold_mirrors_reader._current_surf_data["material"] == "mirror"
         assert "index" not in zemax_fold_mirrors_reader._current_surf_data
         assert "abbe" not in zemax_fold_mirrors_reader._current_surf_data
 
-    def test_read_coordinate_break_params(self, set_test_backend, zemax_fold_mirrors_reader):
+    def test_read_coordinate_break_params(
+        self, set_test_backend, zemax_fold_mirrors_reader
+    ):
         zemax_fold_mirrors_reader._current_surf_data = {"type": "coordinate_break"}
-        
-        params_to_test_values = [1.0, 2.0, 30.0, 5.0, -10.0, 0.0] 
-        expected_keys = ["param_0", "param_1", "param_2", "param_3", "param_4", "param_5"]
+
+        params_to_test_values = [1.0, 2.0, 30.0, 5.0, -10.0, 0.0]
+        expected_keys = [
+            "param_0",
+            "param_1",
+            "param_2",
+            "param_3",
+            "param_4",
+            "param_5",
+        ]
 
         for i, val in enumerate(params_to_test_values):
-            zemax_fold_mirrors_reader._read_surface_parameter(["PARM", str(i + 1), str(val)])
-        
+            zemax_fold_mirrors_reader._read_surface_parameter(
+                ["PARM", str(i + 1), str(val)]
+            )
+
         for i, key in enumerate(expected_keys):
-            assert zemax_fold_mirrors_reader._current_surf_data[key] == params_to_test_values[i]
+            assert (
+                zemax_fold_mirrors_reader._current_surf_data[key]
+                == params_to_test_values[i]
+            )
+
+    def test_generate_complicated_system(
+        self, set_test_backend, zemax_fold_mirrors_reader
+    ):
+        lens = zemax_fold_mirrors_reader.generate_lens()
+        assert lens is not None
+        assert isinstance(lens, Optic)
 
-    def test_generate_complicated_system(self, set_test_backend, zemax_fold_mirrors_reader):
-            lens = zemax_fold_mirrors_reader.generate_lens()
-            assert lens is not None
-            assert isinstance(lens, Optic)
 
 def test_save_load_json_obj(set_test_backend):
     mat = Material("SF11")
diff --git a/tests/test_geometries.py b/tests/test_geometries.py
index 4c33fb1d..02ab7918 100644
--- a/tests/test_geometries.py
+++ b/tests/test_geometries.py
@@ -356,11 +356,9 @@ class TestEvenAsphere:
 class TestPolynomialGeometry:
     def test_str(self, set_test_backend):
         cs = CoordinateSystem()
-        coefficients = be.array([
-            [0.0, 1e-2, -2e-3],
-            [0.1, 1e-2, -1e-3],
-            [0.2, 1e-2, 0.0]
-        ])
+        coefficients = be.array(
+            [[0.0, 1e-2, -2e-3], [0.1, 1e-2, -1e-3], [0.2, 1e-2, 0.0]]
+        )
         geometry = geometries.PolynomialGeometry(
             cs,
             radius=22.0,
@@ -371,11 +369,9 @@ class TestPolynomialGeometry:
 
     def test_sag(self, set_test_backend):
         cs = CoordinateSystem()
-        coefficients = be.array([
-            [0.0, 1e-2, -2e-3],
-            [0.1, 1e-2, -1e-3],
-            [0.2, 1e-2, 0.0]
-        ])
+        coefficients = be.array(
+            [[0.0, 1e-2, -2e-3], [0.1, 1e-2, -1e-3], [0.2, 1e-2, 0.0]]
+        )
         geometry = geometries.PolynomialGeometry(
             cs,
             radius=22.0,
@@ -400,11 +396,9 @@ class TestPolynomialGeometry:
 
     def test_distance(self, set_test_backend):
         cs = CoordinateSystem()
-        coefficients = be.array([
-            [0.0, 1e-2, 2e-3],
-            [0.1, -1e-2, 1e-3],
-            [0.2, 1e-2, 2e-4]
-        ])
+        coefficients = be.array(
+            [[0.0, 1e-2, 2e-3], [0.1, -1e-2, 1e-3], [0.2, 1e-2, 2e-4]]
+        )
         geometry = geometries.PolynomialGeometry(
             cs,
             radius=-26.0,
@@ -442,11 +436,9 @@ class TestPolynomialGeometry:
 
     def test_surface_normal(self, set_test_backend):
         cs = CoordinateSystem()
-        coefficients = be.array([
-            [0.0, 1e-2, 2e-3],
-            [0.1, -1e-2, 1e-3],
-            [0.2, 1e-2, 2e-4]
-        ])
+        coefficients = be.array(
+            [[0.0, 1e-2, 2e-3], [0.1, -1e-2, 1e-3], [0.2, 1e-2, 2e-4]]
+        )
         geometry = geometries.PolynomialGeometry(
             cs,
             radius=-26.0,
@@ -462,11 +454,9 @@ class TestPolynomialGeometry:
 
     def test_to_dict(self, set_test_backend):
         cs = CoordinateSystem()
-        coefficients = be.array([
-             [0.0, 1e-2, 2e-3],
-             [0.1, -1e-2, 1e-3],
-             [0.2, 1e-2, 2e-4]
-        ])
+        coefficients = be.array(
+            [[0.0, 1e-2, 2e-3], [0.1, -1e-2, 1e-3], [0.2, 1e-2, 2e-4]]
+        )
         geometry = geometries.PolynomialGeometry(
             cs,
             radius=-26.0,
@@ -487,11 +477,9 @@ class TestPolynomialGeometry:
 
     def test_from_dict(self, set_test_backend):
         cs = CoordinateSystem()
-        coefficients = be.array([
-            [0.0, 1e-2, 2e-3],
-            [0.1, -1e-2, 1e-3],
-            [0.2, 1e-2, 2e-4]
-        ])
+        coefficients = be.array(
+            [[0.0, 1e-2, 2e-3], [0.1, -1e-2, 1e-3], [0.2, 1e-2, 2e-4]]
+        )
         geometry = geometries.PolynomialGeometry(
             cs,
             radius=-26.0,
@@ -507,11 +495,9 @@ class TestPolynomialGeometry:
 class TestChebyshevGeometry:
     def test_str(self, set_test_backend):
         cs = CoordinateSystem()
-        coefficients = be.array([
-            [0.0, 1e-2, -2e-3],
-            [0.1, 1e-2, -1e-3],
-            [0.2, 1e-2, 0.0]
-        ])
+        coefficients = be.array(
+            [[0.0, 1e-2, -2e-3], [0.1, 1e-2, -1e-3], [0.2, 1e-2, 0.0]]
+        )
         geometry = geometries.ChebyshevPolynomialGeometry(
             cs,
             radius=22.0,
@@ -524,11 +510,9 @@ class TestChebyshevGeometry:
 
     def test_sag(self, set_test_backend):
         cs = CoordinateSystem()
-        coefficients = be.array([
-            [0.0, 1e-2, -2e-3],
-            [0.1, 1e-2, -1e-3],
-            [0.2, 1e-2, 0.0]
-        ])
+        coefficients = be.array(
+            [[0.0, 1e-2, -2e-3], [0.1, 1e-2, -1e-3], [0.2, 1e-2, 0.0]]
+        )
         geometry = geometries.ChebyshevPolynomialGeometry(
             cs,
             radius=22.0,
@@ -555,11 +539,9 @@ class TestChebyshevGeometry:
 
     def test_distance(self, set_test_backend):
         cs = CoordinateSystem()
-        coefficients = be.array([
-            [0.0, 1e-2, -2e-3],
-            [0.1, 1e-2, -1e-3],
-            [0.2, 1e-2, 0.0]
-        ])
+        coefficients = be.array(
+            [[0.0, 1e-2, -2e-3], [0.1, 1e-2, -1e-3], [0.2, 1e-2, 0.0]]
+        )
         geometry = geometries.ChebyshevPolynomialGeometry(
             cs,
             radius=-26.0,
@@ -599,11 +581,9 @@ class TestChebyshevGeometry:
 
     def test_surface_normal(self, set_test_backend):
         cs = CoordinateSystem()
-        coefficients = be.array([
-            [0.0, 1e-2, -2e-3],
-            [0.1, 1e-2, -1e-3],
-            [0.2, 1e-2, 0.0]
-        ])
+        coefficients = be.array(
+            [[0.0, 1e-2, -2e-3], [0.1, 1e-2, -1e-3], [0.2, 1e-2, 0.0]]
+        )
         geometry = geometries.ChebyshevPolynomialGeometry(
             cs,
             radius=-26.0,
@@ -621,11 +601,9 @@ class TestChebyshevGeometry:
 
     def test_invalid_input(self, set_test_backend):
         cs = CoordinateSystem()
-        coefficients = be.array([
-            [0.0, 1e-2, -2e-3],
-            [0.1, 1e-2, -1e-3],
-            [0.2, 1e-2, 0.0]
-        ])
+        coefficients = be.array(
+            [[0.0, 1e-2, -2e-3], [0.1, 1e-2, -1e-3], [0.2, 1e-2, 0.0]]
+        )
         geometry = geometries.ChebyshevPolynomialGeometry(
             cs,
             radius=-26.0,
@@ -640,11 +618,9 @@ class TestChebyshevGeometry:
 
     def test_to_dict(self, set_test_backend):
         cs = CoordinateSystem()
-        coefficients = be.array([
-            [0.0, 1e-2, -2e-3],
-            [0.1, 1e-2, -1e-3],
-            [0.2, 1e-2, 0.0]
-        ])
+        coefficients = be.array(
+            [[0.0, 1e-2, -2e-3], [0.1, 1e-2, -1e-3], [0.2, 1e-2, 0.0]]
+        )
         geometry = geometries.ChebyshevPolynomialGeometry(
             cs,
             radius=-26.0,
@@ -669,11 +645,9 @@ class TestChebyshevGeometry:
 
     def test_from_dict(self, set_test_backend):
         cs = CoordinateSystem()
-        coefficients = be.array([
-            [0.0, 1e-2, -2e-3],
-            [0.1, 1e-2, -1e-3],
-            [0.2, 1e-2, 0.0]
-        ])
+        coefficients = be.array(
+            [[0.0, 1e-2, -2e-3], [0.1, 1e-2, -1e-3], [0.2, 1e-2, 0.0]]
+        )
         geometry = geometries.ChebyshevPolynomialGeometry(
             cs,
             radius=-26.0,
@@ -770,11 +744,7 @@ class TestOddAsphere:
 class TestZernikeGeometry:
     def test_str(self, set_test_backend):
         cs = CoordinateSystem()
-        coefficients = be.array([
-            [0.0, 1e-2, -2e-3],
-            [0, 0, 0],
-            [0, 0, 0]
-        ])
+        coefficients = be.array([[0.0, 1e-2, -2e-3], [0, 0, 0], [0, 0, 0]])
         geometry = geometries.ZernikePolynomialGeometry(
             cs,
             radius=22.0,
@@ -791,57 +761,61 @@ def basic_toroid_geometry(set_test_backend):
     """Provides a basic ToroidalGeometry instance for testing"""
     cs = CoordinateSystem(x=0, y=0, z=0)
     radius_rotation = 100.0  # R (X-Z radius)
-    radius_yz = 50.0         # R_y (Y-Z radius)
-    conic = -0.5             # k_yz (YZ conic)
-    coeffs_poly_y = [1e-5]   
+    radius_yz = 50.0  # R_y (Y-Z radius)
+    conic = -0.5  # k_yz (YZ conic)
+    coeffs_poly_y = [1e-5]
     return geometries.ToroidalGeometry(
         coordinate_system=cs,
         radius_rotation=radius_rotation,
         radius_yz=radius_yz,
         conic=conic,
-        coeffs_poly_y=coeffs_poly_y
+        coeffs_poly_y=coeffs_poly_y,
     )
 
+
 @pytest.fixture
 def cylinder_x_geometry(set_test_backend):
     """Provides a cylindrical geometry (flat in X). R_rot = inf"""
     cs = CoordinateSystem(x=0, y=0, z=0)
-    radius_rotation = be.inf # Flat in X-Z plane
+    radius_rotation = be.inf  # Flat in X-Z plane
     radius_yz = -50.0
     conic = 0.0
     return geometries.ToroidalGeometry(
         coordinate_system=cs,
         radius_rotation=radius_rotation,
         radius_yz=radius_yz,
-        conic=conic
+        conic=conic,
     )
 
+
 @pytest.fixture
 def cylinder_y_geometry(set_test_backend):
     """Provides a cylindrical geometry (flat in Y). R_yz = inf"""
     cs = CoordinateSystem(x=0, y=0, z=0)
-    radius_rotation = 100.0 
-    radius_yz = be.inf      # Flat in Y-Z plane
+    radius_rotation = 100.0
+    radius_yz = be.inf  # Flat in Y-Z plane
     conic = 0.0
     return geometries.ToroidalGeometry(
         coordinate_system=cs,
         radius_rotation=radius_rotation,
         radius_yz=radius_yz,
-        conic=conic
+        conic=conic,
     )
 
+
 @pytest.fixture
 def toroid_no_x_curvature_coeffs(set_test_backend):
     """Toroidal surface with no curvature in x, conic_yz = -1.1, and coeffs."""
     cs = CoordinateSystem(x=0, y=0, z=0)
     return geometries.ToroidalGeometry(
         coordinate_system=cs,
-        radius_rotation=be.inf, 
+        radius_rotation=be.inf,
         radius_yz=50.0,
         conic=-1.1,
-        coeffs_poly_y=[1e-5, -2e-6]
+        coeffs_poly_y=[1e-5, -2e-6],
     )
 
+
 @pytest.fixture
 def toroid_no_y_curvature_coeffs(set_test_backend):
     """Toroidal surface with no curvature in y, conic_yz = -0.9, and coeffs.
@@ -851,13 +825,13 @@ def toroid_no_y_curvature_coeffs(set_test_backend):
     return geometries.ToroidalGeometry(
         coordinate_system=cs,
         radius_rotation=100.0,
-        radius_yz=be.inf, # No curvature in Y (base YZ curve is flat before polynomials)
-        conic=-0.9, # This k_yz will not be used if R_yz is inf for conic part
-        coeffs_poly_y=[1e-5, -2e-6]
+        radius_yz=be.inf,  # No curvature in Y (base YZ curve is flat before polynomials)
+        conic=-0.9,  # This k_yz will not be used if R_yz is inf for conic part
+        coeffs_poly_y=[1e-5, -2e-6],
     )
 
+
 class TestToroidalGeometry:
-    
     def test_toroidal_str(self, set_test_backend):
         """Test string representation."""
         cs = CoordinateSystem()
@@ -874,17 +848,19 @@ class TestToroidalGeometry:
         """Test sag at the vertex (0, 0). Should be 0 by definition."""
         x = be.array([0.0])
         y = be.array([0.0])
-        be.allclose(basic_toroid_geometry.sag(x, y), be.array([0.0]), rtol=1e-5, atol=1e-6)
+        be.allclose(
+            basic_toroid_geometry.sag(x, y), be.array([0.0]), rtol=1e-5, atol=1e-6
+        )
 
     def test_toroidal_normal_vertex(self, basic_toroid_geometry, set_test_backend):
-        """" Test the normal vector at the vertex (0, 0). Should be (0, 0, -1) """
+        """ " Test the normal vector at the vertex (0, 0). Should be (0, 0, -1)"""
         x = be.array([0.0])
         y = be.array([0.0])
         nx, ny, nz = basic_toroid_geometry._surface_normal(x, y)
         be.allclose(nx, be.array([0.0]), rtol=1e-5, atol=1e-6)
         be.allclose(ny, be.array([0.0]), rtol=1e-5, atol=1e-6)
         be.allclose(nz, be.array([-1.0]), rtol=1e-5, atol=1e-6)
-        
+
     def test_toroidal_sag_known_points(self, basic_toroid_geometry, set_test_backend):
         """Test sag at specific points"""
         x = be.array([0.0, 10.0, 0.0])
@@ -894,11 +870,13 @@ class TestToroidalGeometry:
         expected_z_10_0 = 0.50125628
         expected_z_0_5 = 0.25056330
 
-        expected_z = be.array([expected_z_0_10, expected_z_10_0, expected_z_0_5])  
+        expected_z = be.array([expected_z_0_10, expected_z_10_0, expected_z_0_5])
         calculated_z = basic_toroid_geometry.sag(x, y)
-        assert be.allclose(calculated_z, expected_z, rtol=1e-5, atol=1e-6)    
-    
-    def test_toroidal_normal_known_points(self, basic_toroid_geometry, set_test_backend):
+        assert be.allclose(calculated_z, expected_z, rtol=1e-5, atol=1e-6)
+
+    def test_toroidal_normal_known_points(
+        self, basic_toroid_geometry, set_test_backend
+    ):
         """Test normal at specific points."""
         x = be.array([0.0, 10.0])
         y = be.array([10.0, 0.0])
@@ -921,47 +899,51 @@ class TestToroidalGeometry:
         assert be.allclose(nx, expected_nx, rtol=rtol, atol=atol)
         assert be.allclose(ny, expected_ny, rtol=rtol, atol=atol)
         assert be.allclose(nz, expected_nz, rtol=rtol, atol=atol)
-        
+
     def test_cylinder_x_sag(self, cylinder_x_geometry, set_test_backend):
         """Test sag for cylinder flat in X. Should only depend on y."""
         x = be.array([0.0, 10.0, 10.0])
         y = be.array([5.0, 5.0, 0.0])
-        
-        expected_z = be.array([-0.25062818, -0.25062818, 0.0]) 
+
+        expected_z = be.array([-0.25062818, -0.25062818, 0.0])
         calculated_z = cylinder_x_geometry.sag(x, y)
-        assert_allclose(calculated_z, expected_z, rtol=1e-5, atol=1e-6)    
-        
+        assert_allclose(calculated_z, expected_z, rtol=1e-5, atol=1e-6)
+
     def test_cylinder_y_sag(self, cylinder_y_geometry, set_test_backend):
         """Test sag for cylinder flat in Y. Should only depend on x."""
         x = be.array([5.0, 5.0, 0.0])
         y = be.array([0.0, 10.0, 10.0])
-        
-        expected_z = be.array([0.12507822, 0.12507822, 0.0]) 
+
+        expected_z = be.array([0.12507822, 0.12507822, 0.0])
         calculated_z = cylinder_y_geometry.sag(x, y)
         assert_allclose(calculated_z, expected_z, rtol=1e-5, atol=1e-6)
-    
-    def test_toroid_no_x_curvature_sag(self, toroid_no_x_curvature_coeffs, set_test_backend):
+
+    def test_toroid_no_x_curvature_sag(
+        self, toroid_no_x_curvature_coeffs, set_test_backend
+    ):
         """Test sag for toroid with R_rot = inf and YZ coeffs."""
-        
+
         geom = toroid_no_x_curvature_coeffs
-        x = be.array([0.0, 10.0, -5.0]) 
+        x = be.array([0.0, 10.0, -5.0])
         y = be.array([0.0, 1.0, 2.0])
 
         expected_z = be.array([0.0, 0.010007900001, 0.0400064001])
         calculated_z = geom.sag(x, y)
         assert_allclose(calculated_z, expected_z, rtol=1e-5, atol=1e-7)
 
-    def test_toroid_no_y_curvature_sag(self, toroid_no_y_curvature_coeffs, set_test_backend):
+    def test_toroid_no_y_curvature_sag(
+        self, toroid_no_y_curvature_coeffs, set_test_backend
+    ):
         """Test sag for toroid with R_yz = inf and YZ coeffs. Sag comes from X-rotation of YZ poly."""
-        
+
         geom = toroid_no_y_curvature_coeffs
         x = be.array([0.0, 10.0, 5.0])
-        y = be.array([0.0, 1.0, 2.0]) 
+        y = be.array([0.0, 1.0, 2.0])
 
         expected_z = be.array([0.0, 0.501264335, 0.12508686])
         calculated_z = geom.sag(x, y)
         assert_allclose(calculated_z, expected_z, rtol=1e-5, atol=1e-7)
-    
+
     def test_toroidal_sag_vs_zemax(self, basic_toroid_geometry, set_test_backend):
         """
         Compares sag values calculated by Optiland with
@@ -969,27 +951,28 @@ class TestToroidalGeometry:
         """
         geometry = basic_toroid_geometry
 
-        x_coords = be.array([0.0,  2.5, 0.0, -2.5,  5.0,  -5.0,  2.5, -2.5])
-        y_coords = be.array([0.0,  0.0, 2.5,  0.0,  2.5,  -2.5, -2.5,  2.5])
+        x_coords = be.array([0.0, 2.5, 0.0, -2.5, 5.0, -5.0, 2.5, -2.5])
+        y_coords = be.array([0.0, 0.0, 2.5, 0.0, 2.5, -2.5, -2.5, 2.5])
 
         # --- Zemax Sag Data ---
-        zemax_z_sag = be.array([
-            0.0,                      # (0, 0) - Vertex
-            3.125488433897521E-002,   # (2.5, 0) 
-            6.258204346657634E-002,   # (0, 2.5) 
-            3.125488433897521E-002,   # (-2.5, 0) 
-            1.877386899843393E-001,   # (5, 2.5) 
-            1.877386899843393E-001,   # (-5, -2.5) 
-            9.385650612446353E-002,   # (2.5, -2.5) 
-            9.385650612446353E-002    # (-2.5, 2.5) 
-        ]) 
+        zemax_z_sag = be.array(
+            [
+                0.0,  # (0, 0) - Vertex
+                3.125488433897521e-002,  # (2.5, 0)
+                6.258204346657634e-002,  # (0, 2.5)
+                3.125488433897521e-002,  # (-2.5, 0)
+                1.877386899843393e-001,  # (5, 2.5)
+                1.877386899843393e-001,  # (-5, -2.5)
+                9.385650612446353e-002,  # (2.5, -2.5)
+                9.385650612446353e-002,  # (-2.5, 2.5)
+            ]
+        )
 
         # Calculate sag using Optiland
         optiland_z_sag = geometry.sag(x_coords, y_coords)
 
-        assert be.allclose(optiland_z_sag,zemax_z_sag,rtol=1e-5,atol=1e-6)
+        assert be.allclose(optiland_z_sag, zemax_z_sag, rtol=1e-5, atol=1e-6)
 
-    
     def test_toroidal_ray_tracing_comparison(self, set_test_backend):
         """
         Traces rays through a single toroidal surface and compares output
@@ -1011,13 +994,13 @@ class TestToroidalGeometry:
         )
         lens.add_surface(index=2, thickness=10.0, material="air")
         lens.add_surface(index=3)
-        
+
         lens.set_aperture(aperture_type="EPD", value=10.0)
         lens.add_wavelength(value=0.550, is_primary=True)
         lens.set_field_type("angle")
         lens.add_field(y=0)
 
-        num_rays = 5 # Number of rays per fan
+        num_rays = 5  # Number of rays per fan
         wavelength = 0.550
         z_start = 0.0
 
@@ -1030,20 +1013,51 @@ class TestToroidalGeometry:
         M_in_yfan = be.zeros(num_rays)
         N_in_yfan = be.ones(num_rays)
         intensity_yfan = be.ones(num_rays)
-        rays_in_yfan = RealRays(x=x_in_yfan, y=y_in_yfan, z=z_in_yfan,
-                                L=L_in_yfan, M=M_in_yfan, N=N_in_yfan,
-                                wavelength=wavelength, intensity=intensity_yfan)
+        rays_in_yfan = RealRays(
+            x=x_in_yfan,
+            y=y_in_yfan,
+            z=z_in_yfan,
+            L=L_in_yfan,
+            M=M_in_yfan,
+            N=N_in_yfan,
+            wavelength=wavelength,
+            intensity=intensity_yfan,
+        )
 
         # Trace Y-Fan Rays
         rays_out_yfan = lens.surface_group.trace(rays_in_yfan)
         print("Y-Fan Rays:")
         print(rays_out_yfan.y)
-        zemax_x_out_yfan = be.array([0.0] * num_rays) 
-        zemax_y_out_yfan = be.array([-8.123193233401276E-001, -4.676255499616224E-001, 0, 4.676255499616224E-001, 8.123193233401276E-001]) 
-        zemax_z_out_yfan = be.array([15.0] * num_rays) 
-        zemax_L_out_yfan = be.array([0.0] * num_rays) 
-        zemax_M_out_yfan = be.array([3.251509839270260E-001, 1.537950377308984E-001, 0.0, -1.537950377308984E-001, -3.251509839270260E-001]) 
-        zemax_N_out_yfan = be.array([9.456621160072382E-001, 9.881027711576116E-001, 1.0, 9.881027711576116E-001, 9.456621160072382E-001]) 
+        zemax_x_out_yfan = be.array([0.0] * num_rays)
+        zemax_y_out_yfan = be.array(
+            [
+                -8.123193233401276e-001,
+                -4.676255499616224e-001,
+                0,
+                4.676255499616224e-001,
+                8.123193233401276e-001,
+            ]
+        )
+        zemax_z_out_yfan = be.array([15.0] * num_rays)
+        zemax_L_out_yfan = be.array([0.0] * num_rays)
+        zemax_M_out_yfan = be.array(
+            [
+                3.251509839270260e-001,
+                1.537950377308984e-001,
+                0.0,
+                -1.537950377308984e-001,
+                -3.251509839270260e-001,
+            ]
+        )
+        zemax_N_out_yfan = be.array(
+            [
+                9.456621160072382e-001,
+                9.881027711576116e-001,
+                1.0,
+                9.881027711576116e-001,
+                9.456621160072382e-001,
+            ]
+        )
 
         # Comparison Assertions for Y-Fan
         assert be.allclose(rays_out_yfan.x, zemax_x_out_yfan, rtol=1e-5, atol=1e-6)
@@ -1053,27 +1067,58 @@ class TestToroidalGeometry:
         assert be.allclose(rays_out_yfan.M, zemax_M_out_yfan, rtol=1e-5, atol=1e-6)
         assert be.allclose(rays_out_yfan.N, zemax_N_out_yfan, rtol=1e-5, atol=1e-6)
 
-         # --- Sagittal (X) Fan Test ---
-        x_coords = be.linspace(-5.0, 5.0, num_rays) 
+        # --- Sagittal (X) Fan Test ---
+        x_coords = be.linspace(-5.0, 5.0, num_rays)
         x_in_xfan = x_coords
         y_in_xfan = be.zeros(num_rays)
         z_in_xfan = be.array([z_start] * num_rays)
-        L_in_xfan = be.zeros(num_rays) 
-        M_in_xfan = be.zeros(num_rays) 
-        N_in_xfan = be.ones(num_rays)  
+        L_in_xfan = be.zeros(num_rays)
+        M_in_xfan = be.zeros(num_rays)
+        N_in_xfan = be.ones(num_rays)
         intensity_xfan = be.ones(num_rays)
-        rays_in_xfan = RealRays(x=x_in_xfan, y=y_in_xfan, z=z_in_xfan,
-                                L=L_in_xfan, M=M_in_xfan, N=N_in_xfan,
-                                wavelength=wavelength, intensity=intensity_xfan)
+        rays_in_xfan = RealRays(
+            x=x_in_xfan,
+            y=y_in_xfan,
+            z=z_in_xfan,
+            L=L_in_xfan,
+            M=M_in_xfan,
+            N=N_in_xfan,
+            wavelength=wavelength,
+            intensity=intensity_xfan,
+        )
 
         rays_out_xfan = lens.surface_group.trace(rays_in_xfan)
 
-        zemax_x_out_xfan = be.array([-4.668385225648558E+000, -2.333547899735358E+000, 0.0, 2.333547899735358E+000, 4.668385225648558E+000]) 
-        zemax_y_out_xfan = be.array([0.0] * num_rays)             
-        zemax_z_out_xfan = be.array([15.0] * num_rays)             
-        zemax_L_out_xfan = be.array([2.502086086422164E-002, 1.250260502601134E-002, 0.0, -1.250260502601134E-002, -2.502086086422164E-002]) 
-        zemax_M_out_xfan = be.array([0.0] * num_rays)             
-        zemax_N_out_xfan = be.array([9.996869292541608E-001, 9.999218393792406E-001, 1.0, 9.999218393792406E-001, 9.996869292541608E-001]) 
+        zemax_x_out_xfan = be.array(
+            [
+                -4.668385225648558e000,
+                -2.333547899735358e000,
+                0.0,
+                2.333547899735358e000,
+                4.668385225648558e000,
+            ]
+        )
+        zemax_y_out_xfan = be.array([0.0] * num_rays)
+        zemax_z_out_xfan = be.array([15.0] * num_rays)
+        zemax_L_out_xfan = be.array(
+            [
+                2.502086086422164e-002,
+                1.250260502601134e-002,
+                0.0,
+                -1.250260502601134e-002,
+                -2.502086086422164e-002,
+            ]
+        )
+        zemax_M_out_xfan = be.array([0.0] * num_rays)
+        zemax_N_out_xfan = be.array(
+            [
+                9.996869292541608e-001,
+                9.999218393792406e-001,
+                1.0,
+                9.999218393792406e-001,
+                9.996869292541608e-001,
+            ]
+        )
 
         assert be.allclose(rays_out_xfan.x, zemax_x_out_xfan, rtol=1e-5, atol=1e-6)
         assert be.allclose(rays_out_xfan.y, zemax_y_out_xfan, rtol=1e-5, atol=1e-6)
@@ -1085,13 +1130,13 @@ class TestToroidalGeometry:
     def test_toroidal_to_dict(self, basic_toroid_geometry, set_test_backend):
         """Test serialization to dictionary."""
         geom_dict = basic_toroid_geometry.to_dict()
-        assert geom_dict["type"] == "ToroidalGeometry" 
-        assert geom_dict["geometry_type"] == "Toroidal" 
+        assert geom_dict["type"] == "ToroidalGeometry"
+        assert geom_dict["geometry_type"] == "Toroidal"
         assert geom_dict["radius_rotation"] == 100.0
         assert geom_dict["radius_yz"] == 50.0
         assert geom_dict["conic_yz"] == -0.5
         assert geom_dict["coeffs_poly_y"] == pytest.approx([1e-5])
-        
+
         assert "radius" not in geom_dict
         assert "conic" not in geom_dict
         assert "coefficients" not in geom_dict
@@ -1099,9 +1144,9 @@ class TestToroidalGeometry:
     def test_toroidal_from_dict(self, basic_toroid_geometry, set_test_backend):
         """Test deserialization from dictionary."""
         geom_dict = basic_toroid_geometry.to_dict()
-        new_geometry = geometries.ToroidalGeometry.from_dict(geom_dict) 
+        new_geometry = geometries.ToroidalGeometry.from_dict(geom_dict)
         assert isinstance(new_geometry, geometries.ToroidalGeometry)
-        assert new_geometry.to_dict() == geom_dict 
+        assert new_geometry.to_dict() == geom_dict
 
     def test_toroidal_from_dict_invalid(self, set_test_backend):
         """Test deserialization with missing keys."""
@@ -1113,16 +1158,22 @@ class TestToroidalGeometry:
         }
         with pytest.raises(ValueError):
             geometries.ToroidalGeometry.from_dict(invalid_dict)
-            
-    def test_inf_radius_intersect_sphere_normal_incidence(self, cylinder_x_geometry, set_test_backend):
-            """Test _intersection_sphere for inf radius: ray normal to XY plane (z=0 plane)."""
-            
-            rays = RealRays(x=0.0, y=0.0, z=10.0, L=0.0, M=0.0, N=-1.0, intensity=1.0, wavelength=0.55)
-            ix, iy, iz = cylinder_x_geometry._intersection(rays)
-            
-            be.allclose(ix, be.array(0.0), rtol=1e-5, atol=1e-6)
-            be.allclose(iy, be.array(0.0), rtol=1e-5, atol=1e-6)
-            be.allclose(iz, be.array(0.0), rtol=1e-5, atol=1e-6) # Intersection coordinates should be (0, 0, 0)
+
+    def test_inf_radius_intersect_sphere_normal_incidence(
+        self, cylinder_x_geometry, set_test_backend
+    ):
+        """Test _intersection_sphere for inf radius: ray normal to XY plane (z=0 plane)."""
+
+        rays = RealRays(
+            x=0.0, y=0.0, z=10.0, L=0.0, M=0.0, N=-1.0, intensity=1.0, wavelength=0.55
+        )
+        ix, iy, iz = cylinder_x_geometry._intersection(rays)
+
+        be.allclose(ix, be.array(0.0), rtol=1e-5, atol=1e-6)
+        be.allclose(iy, be.array(0.0), rtol=1e-5, atol=1e-6)
+        be.allclose(
+            iz, be.array(0.0), rtol=1e-5, atol=1e-6
+        )  # Intersection coordinates should be (0, 0, 0)
 
 
 class TestBiconicGeometry:
@@ -1133,7 +1184,9 @@ class TestBiconicGeometry:
 
     def test_sag_vertex(self, set_test_backend):
         cs = CoordinateSystem()
-        geom = BiconicGeometry(cs, radius_x=10.0, radius_y=20.0, conic_x=0.5, conic_y=-0.5)
+        geom = BiconicGeometry(
+            cs, radius_x=10.0, radius_y=20.0, conic_x=0.5, conic_y=-0.5
+        )
         assert_allclose(geom.sag(0, 0), 0.0)
 
     def test_sag_finite_radii(self, set_test_backend):
@@ -1143,7 +1196,9 @@ class TestBiconicGeometry:
 
     def test_sag_rx_infinite(self, set_test_backend):
         cs = CoordinateSystem()
-        geom = BiconicGeometry(cs, radius_x=be.inf, radius_y=20.0) # Cylindrical along X
+        geom = BiconicGeometry(
+            cs, radius_x=be.inf, radius_y=20.0
+        )  # Cylindrical along X
         # Sag should only depend on y. zx(any) = 0
         # zy(1) for Ry=20 is 0.02501563183003138
         assert_allclose(geom.sag(x=10, y=1), 0.02501563183003138)
@@ -1151,7 +1206,9 @@ class TestBiconicGeometry:
 
     def test_sag_ry_infinite(self, set_test_backend):
         cs = CoordinateSystem()
-        geom = BiconicGeometry(cs, radius_x=10.0, radius_y=be.inf) # Cylindrical along Y
+        geom = BiconicGeometry(
+            cs, radius_x=10.0, radius_y=be.inf
+        )  # Cylindrical along Y
         # Sag should only depend on x. zy(any) = 0
         # zx(1) for Rx=10 is 0.05012562854249503
         assert_allclose(geom.sag(x=1, y=10), 0.05012562854249503)
@@ -1162,13 +1219,16 @@ class TestBiconicGeometry:
         geom = BiconicGeometry(cs, radius_x=be.inf, radius_y=be.inf)
         assert_allclose(geom.sag(x=10, y=20), 0.0)
         # Also test with conics, should still be zero
-        geom_conic = BiconicGeometry(cs, radius_x=be.inf, radius_y=be.inf, conic_x=0.5, conic_y=-1.0)
+        geom_conic = BiconicGeometry(
+            cs, radius_x=be.inf, radius_y=be.inf, conic_x=0.5, conic_y=-1.0
+        )
         assert_allclose(geom_conic.sag(x=10, y=20), 0.0)
 
-
     def test_sag_with_conics(self, set_test_backend):
         cs = CoordinateSystem()
-        geom = BiconicGeometry(cs, radius_x=10.0, radius_y=20.0, conic_x=-1.0, conic_y=0.5)
+        geom = BiconicGeometry(
+            cs, radius_x=10.0, radius_y=20.0, conic_x=-1.0, conic_y=0.5
+        )
         expected_sag = 0.05 + 0.02502345130203264
         assert_allclose(geom.sag(x=1, y=1), expected_sag)
 
@@ -1182,7 +1242,9 @@ class TestBiconicGeometry:
 
     def test_surface_normal_vertex(self, set_test_backend):
         cs = CoordinateSystem()
-        geom = BiconicGeometry(cs, radius_x=10.0, radius_y=20.0, conic_x=0.5, conic_y=-0.5)
+        geom = BiconicGeometry(
+            cs, radius_x=10.0, radius_y=20.0, conic_x=0.5, conic_y=-0.5
+        )
         nx, ny, nz = geom._surface_normal(x=0, y=0)
         assert_allclose(nx, 0.0)
         assert_allclose(ny, 0.0)
@@ -1199,41 +1261,50 @@ class TestBiconicGeometry:
 
     def test_surface_normal_cylindrical_rx_inf(self, set_test_backend):
         cs = CoordinateSystem()
-        geom = BiconicGeometry(cs, radius_x=be.inf, radius_y=10.0, conic_y=0.0) # Cylinder along X
-        nx, ny, nz = geom._surface_normal(x=5, y=1) # x shouldn't matter for dfdx=0
+        geom = BiconicGeometry(
+            cs, radius_x=be.inf, radius_y=10.0, conic_y=0.0
+        )  # Cylinder along X
+        nx, ny, nz = geom._surface_normal(x=5, y=1)  # x shouldn't matter for dfdx=0
         assert_allclose(nx, 0.0)
         assert_allclose(ny, 0.1)
         assert_allclose(nz, -0.99498743710662)
-        
+
     def test_surface_normal_array_input(self, set_test_backend):
         cs = CoordinateSystem()
-        geom = BiconicGeometry(cs, radius_x=10.0, radius_y=10.0, conic_x=0.0, conic_y=0.0) # Spherical R=10
+        geom = BiconicGeometry(
+            cs, radius_x=10.0, radius_y=10.0, conic_x=0.0, conic_y=0.0
+        )  # Spherical R=10
         x = be.array([0, 1])
         y = be.array([0, 1])
-        
+
         # Point (0,0): nx=0, ny=0, nz=-1
         # Point (1,1): nx=0.099503719, ny=0.099503719, nz=-0.990049373 (from spherical test)
         expected_nx = be.array([0.0, 0.09950371902099892])
         expected_ny = be.array([0.0, 0.09950371902099892])
         expected_nz = be.array([-1.0, -0.9900493732390136])
-        
+
         nx_calc, ny_calc, nz_calc = geom._surface_normal(x, y)
         assert_allclose(nx_calc, expected_nx)
         assert_allclose(ny_calc, expected_ny)
         assert_allclose(nz_calc, expected_nz)
 
-
     def test_distance_simple(self, set_test_backend):
         cs = CoordinateSystem()
         geom = BiconicGeometry(cs, radius_x=10.0, radius_y=20.0)
-        rays = RealRays(x=0.0, y=0.0, z=-5.0, L=0.0, M=0.0, N=1.0, wavelength=0.55, intensity=1.0)
+        rays = RealRays(
+            x=0.0, y=0.0, z=-5.0, L=0.0, M=0.0, N=1.0, wavelength=0.55, intensity=1.0
+        )
         # Sag at (0,0) is 0, so distance to surface at (0,0) should be 5.0
-        assert_allclose(geom.distance(rays), 5.0, atol=1e-9) # Newton-Raphson tolerance
+        assert_allclose(geom.distance(rays), 5.0, atol=1e-9)  # Newton-Raphson tolerance
 
     def test_distance_planar_biconic(self, set_test_backend):
         cs = CoordinateSystem()
-        geom = BiconicGeometry(cs, radius_x=be.inf, radius_y=be.inf, conic_x=0.0, conic_y=0.0) # Planar
-        rays = RealRays(x=1.0, y=1.0, z=-5.0, L=0.0, M=0.0, N=1.0, wavelength=0.55, intensity=1.0)
+        geom = BiconicGeometry(
+            cs, radius_x=be.inf, radius_y=be.inf, conic_x=0.0, conic_y=0.0
+        )  # Planar
+        rays = RealRays(
+            x=1.0, y=1.0, z=-5.0, L=0.0, M=0.0, N=1.0, wavelength=0.55, intensity=1.0
+        )
         # Distance to plane z=0 should be 5.0
         # The parent NewtonRaphsonGeometry's _intersection calls _intersection_plane if self.radius is inf.
         # In Biconic.__init__, self.radius is set to radius_x. So this path is tested.
@@ -1259,7 +1330,7 @@ class TestBiconicGeometry:
         assert geom_dict["conic_y"] == 0.2
         assert geom_dict["tol"] == 1e-9
         assert geom_dict["max_iter"] == 50
-        assert "radius" not in geom_dict 
+        assert "radius" not in geom_dict
         assert "conic" not in geom_dict
 
         reconstructed_geom = BiconicGeometry.from_dict(geom_dict)
@@ -1271,7 +1342,7 @@ class TestBiconicGeometry:
         assert reconstructed_geom.tol == original_geom.tol
         assert reconstructed_geom.max_iter == original_geom.max_iter
         assert reconstructed_geom.cs.to_dict() == original_geom.cs.to_dict()
-        
+
         # Check if the reconstructed dict is identical (it should be)
         assert reconstructed_geom.to_dict() == geom_dict
 
@@ -1284,23 +1355,29 @@ class TestBiconicGeometry:
             "radius_y": 20.0,
             # conic_x, conic_y, tol, max_iter will use defaults
         }
-        
+
         # Test missing radius_x
         invalid_dict_rx = minimal_valid_dict.copy()
         del invalid_dict_rx["radius_x"]
-        with pytest.raises(ValueError, match="Missing required BiconicGeometry keys: {'radius_x'}"):
+        with pytest.raises(
+            ValueError, match="Missing required BiconicGeometry keys: {'radius_x'}"
+        ):
             BiconicGeometry.from_dict(invalid_dict_rx)
 
         # Test missing radius_y
         invalid_dict_ry = minimal_valid_dict.copy()
         del invalid_dict_ry["radius_y"]
-        with pytest.raises(ValueError, match="Missing required BiconicGeometry keys: {'radius_y'}"):
+        with pytest.raises(
+            ValueError, match="Missing required BiconicGeometry keys: {'radius_y'}"
+        ):
             BiconicGeometry.from_dict(invalid_dict_ry)
-            
+
         # Test missing cs
         invalid_dict_cs = minimal_valid_dict.copy()
         del invalid_dict_cs["cs"]
-        with pytest.raises(ValueError, match="Missing required BiconicGeometry keys: {'cs'}"):
+        with pytest.raises(
+            ValueError, match="Missing required BiconicGeometry keys: {'cs'}"
+        ):
             BiconicGeometry.from_dict(invalid_dict_cs)
 
     def test_from_dict_default_conics_tol_max_iter(self, set_test_backend):
@@ -1316,5 +1393,9 @@ class TestBiconicGeometry:
         assert geom.Ry == 20.0
         assert geom.kx == 0.0  # Default
         assert geom.ky == 0.0  # Default
-        assert geom.tol == 1e-10  # Default from NewtonRaphsonGeometry via BiconicGeometry
-        assert geom.max_iter == 100 # Default from NewtonRaphsonGeometry via BiconicGeometry
+        assert (
+            geom.tol == 1e-10
+        )  # Default from NewtonRaphsonGeometry via BiconicGeometry
+        assert (
+            geom.max_iter == 100
+        )  # Default from NewtonRaphsonGeometry via BiconicGeometry
diff --git a/tests/test_jones.py b/tests/test_jones.py
index 40099e0d..7ee19555 100644
--- a/tests/test_jones.py
+++ b/tests/test_jones.py
@@ -207,36 +207,20 @@ def test_jones_linear_retarder(set_test_backend):
     jones_retarder = jones.JonesLinearRetarder(retardance=0.5, theta=0.0)
     jones_matrix = jones_retarder.calculate_matrix(rays)
     assert jones_matrix.shape == (1, 3, 3)
-    assert_allclose(be.real(jones_matrix[0, 0, 0]), 
-        0.9689124217106447
-    )
-    assert_allclose(be.imag(jones_matrix[0, 0, 0]), 
-        -0.2474039592545229
-    )
-    assert_allclose(be.real(jones_matrix[0, 1, 1]), 
-        0.9689124217106447
-    )
-    assert_allclose(be.imag(jones_matrix[0, 1, 1]), 
-        0.24740395925452294
-    )
+    assert_allclose(be.real(jones_matrix[0, 0, 0]), 0.9689124217106447)
+    assert_allclose(be.imag(jones_matrix[0, 0, 0]), -0.2474039592545229)
+    assert_allclose(be.real(jones_matrix[0, 1, 1]), 0.9689124217106447)
+    assert_allclose(be.imag(jones_matrix[0, 1, 1]), 0.24740395925452294)
     assert jones_matrix[0, 2, 2] == 1.0
 
     # Test with retardance = 0.5, theta = 0.5
     jones_retarder = jones.JonesLinearRetarder(retardance=0.5, theta=0.5)
     jones_matrix = jones_retarder.calculate_matrix(rays)
     assert jones_matrix.shape == (1, 3, 3)
-    assert_allclose(be.real(jones_matrix[0, 0, 0]), 
-        0.9689124217106448
-    )
-    assert_allclose(be.imag(jones_matrix[0, 0, 0]), 
-        -0.1336729296661260
-    )
-    assert_allclose(be.real(jones_matrix[0, 1, 1]), 
-        0.9689124217106448
-    )
-    assert_allclose(be.imag(jones_matrix[0, 1, 1]), 
-        0.13367292966612604
-    )
+    assert_allclose(be.real(jones_matrix[0, 0, 0]), 0.9689124217106448)
+    assert_allclose(be.imag(jones_matrix[0, 0, 0]), -0.1336729296661260)
+    assert_allclose(be.real(jones_matrix[0, 1, 1]), 0.9689124217106448)
+    assert_allclose(be.imag(jones_matrix[0, 1, 1]), 0.13367292966612604)
     assert jones_matrix[0, 2, 2] == 1.0
 
 
@@ -247,36 +231,20 @@ def test_jones_quarter_wave_retarder(set_test_backend):
     jones_retarder = jones.JonesQuarterWaveRetarder(theta=0.0)
     jones_matrix = jones_retarder.calculate_matrix(rays)
     assert jones_matrix.shape == (1, 3, 3)
-    assert_allclose(be.real(jones_matrix[0, 0, 0]), 
-        0.7071067811865476
-    )
-    assert_allclose(be.imag(jones_matrix[0, 0, 0]), 
-        -0.7071067811865476
-    )
-    assert_allclose(be.real(jones_matrix[0, 1, 1]), 
-        0.7071067811865476
-    )
-    assert_allclose(be.imag(jones_matrix[0, 1, 1]), 
-        0.7071067811865476
-    )
+    assert_allclose(be.real(jones_matrix[0, 0, 0]), 0.7071067811865476)
+    assert_allclose(be.imag(jones_matrix[0, 0, 0]), -0.7071067811865476)
+    assert_allclose(be.real(jones_matrix[0, 1, 1]), 0.7071067811865476)
+    assert_allclose(be.imag(jones_matrix[0, 1, 1]), 0.7071067811865476)
     assert jones_matrix[0, 2, 2] == 1.0
 
     # Test with theta = 0.5
     jones_retarder = jones.JonesQuarterWaveRetarder(theta=0.5)
     jones_matrix = jones_retarder.calculate_matrix(rays)
     assert jones_matrix.shape == (1, 3, 3)
-    assert_allclose(be.real(jones_matrix[0, 0, 0]), 
-        0.7071067811865476
-    )
-    assert_allclose(be.imag(jones_matrix[0, 0, 0]), 
-        -0.3820514243700897
-    )
-    assert_allclose(be.real(jones_matrix[0, 1, 1]), 
-        0.7071067811865476
-    )
-    assert_allclose(be.imag(jones_matrix[0, 1, 1]), 
-        0.38205142437008976
-    )
+    assert_allclose(be.real(jones_matrix[0, 0, 0]), 0.7071067811865476)
+    assert_allclose(be.imag(jones_matrix[0, 0, 0]), -0.3820514243700897)
+    assert_allclose(be.real(jones_matrix[0, 1, 1]), 0.7071067811865476)
+    assert_allclose(be.imag(jones_matrix[0, 1, 1]), 0.38205142437008976)
     assert jones_matrix[0, 2, 2] == 1.0
 
 
@@ -298,11 +266,7 @@ def test_jones_half_wave_retarder(set_test_backend):
     jones_matrix = jones_retarder.calculate_matrix(rays)
     assert jones_matrix.shape == (1, 3, 3)
     assert_allclose(be.real(jones_matrix[0, 0, 0]), 0.0)
-    assert_allclose(be.imag(jones_matrix[0, 0, 0]), 
-        -0.5403023058681398
-    )
+    assert_allclose(be.imag(jones_matrix[0, 0, 0]), -0.5403023058681398)
     assert_allclose(be.real(jones_matrix[0, 1, 1]), 0.0)
-    assert_allclose(be.imag(jones_matrix[0, 1, 1]), 
-        0.5403023058681398
-    )
+    assert_allclose(be.imag(jones_matrix[0, 1, 1]), 0.5403023058681398)
     assert jones_matrix[0, 2, 2] == 1.0
diff --git a/tests/test_monte_carlo.py b/tests/test_monte_carlo.py
index 366da9c8..4eabe91f 100644
--- a/tests/test_monte_carlo.py
+++ b/tests/test_monte_carlo.py
@@ -48,7 +48,7 @@ def test_run(monte_carlo):
         "0: f1",
         "1: f2",
         "Radius of Curvature, Surface 1",
-        "C0: Thickness, Surface 2"
+        "C0: Thickness, Surface 2",
     }
     assert set(monte_carlo._results.columns) == res
 
diff --git a/tests/test_mtf.py b/tests/test_mtf.py
index fb00a298..03b1a121 100644
--- a/tests/test_mtf.py
+++ b/tests/test_mtf.py
@@ -1,5 +1,6 @@
 import pytest
 import matplotlib
+
 matplotlib.use("Agg")  # ensure non-interactive backend for testing
 import matplotlib.pyplot as plt
 from unittest.mock import patch
diff --git a/tests/test_operand.py b/tests/test_operand.py
index f3701b85..f6cd9d7a 100644
--- a/tests/test_operand.py
+++ b/tests/test_operand.py
@@ -1,8 +1,11 @@
-import optiland.backend as be
 import pytest
 
+import optiland.backend as be
+from optiland.optic import Optic
 from optiland.optimization import operand
+from optiland.optimization.operand import RayOperand
 from optiland.samples.telescopes import HubbleTelescope
+
 from .utils import assert_allclose
 
 
@@ -312,6 +315,86 @@ class TestRayOperand:
         }
         assert_allclose(operand.RayOperand.OPD_difference(**data), 0.2211995620762635)
 
+    def create_test_optic(self):
+        lens = Optic(name="TMA")
+        lens.set_aperture(aperture_type="EPD", value=10)
+        lens.set_field_type(field_type="angle")
+        lens.add_field(y=0)
+        lens.add_field(y=+1.5)
+        lens.add_field(y=-1.5)
+        lens.add_wavelength(value=0.55, is_primary=True)
+
+        lens.add_surface(index=0, radius=be.inf, thickness=be.inf)
+        lens.add_surface(
+            index=1,
+            radius=-100,
+            thickness=-20,
+            conic=0,
+            material="mirror",
+            rx=be.radians(-15.0),
+            is_stop=True,
+        )
+        lens.add_surface(
+            index=2,
+            radius=-100,
+            thickness=+20,
+            conic=0,
+            material="mirror",
+            rx=be.radians(-10.0),
+            dy=-11.5,
+        )
+        lens.add_surface(
+            index=3,
+            radius=-100,
+            thickness=-19,
+            conic=0,
+            material="mirror",
+            rx=be.radians(-1.0),
+            dy=-15,
+        )
+        lens.add_surface(index=4, dy=-19.3)
+        return lens
+
+    def test_clearance(self, set_test_backend):
+        optic = self.create_test_optic()
+        wavelength = 0.55
+
+        dist1 = RayOperand.clearance(
+            optic=optic,
+            line_ray_surface_idx=1,
+            line_ray_field_coords=(0.0, 0.0),
+            line_ray_pupil_coords=(0.0, -1.0),
+            point_ray_surface_idx=3,
+            point_ray_field_coords=(0.0, 1.0),
+            point_ray_pupil_coords=(0.0, 0.0),
+            wavelength=wavelength,
+        )
+        assert_allclose(dist1, -7.412094834746042)
+
+        dist2 = RayOperand.clearance(
+            optic=optic,
+            line_ray_surface_idx=1,
+            line_ray_field_coords=(0.0, 0.0),
+            line_ray_pupil_coords=(0.0, 0.0),
+            point_ray_surface_idx=3,
+            point_ray_field_coords=(0.0, 0.0),
+            point_ray_pupil_coords=(0.0, 0.0),
+            wavelength=wavelength,
+        )
+        assert_allclose(dist2, -13.065596389231768)
+
+        dist3 = RayOperand.clearance(
+            optic=optic,
+            line_ray_surface_idx=1,
+            line_ray_field_coords=(0.0, -1.0),
+            line_ray_pupil_coords=(0.0, 1.0),
+            point_ray_surface_idx=3,
+            point_ray_field_coords=(0.0, 1.0),
+            point_ray_pupil_coords=(0.0, 0.0),
+            wavelength=wavelength,
+        )
+        assert_allclose(dist3, -15.730530102711754)
+
 
 class TestOperand:
     def test_get_value(self, set_test_backend, hubble):
diff --git a/tests/test_optic.py b/tests/test_optic.py
index 645b69c0..0a30da45 100644
--- a/tests/test_optic.py
+++ b/tests/test_optic.py
@@ -423,4 +423,6 @@ class TestOptic:
     def test_invalid_coordinate_system(self, set_test_backend):
         with pytest.raises(ValueError):
             self.optic.add_surface(index=0, radius=be.inf, z=-100)
-            self.optic.add_surface(index=1, radius=be.inf, z=0, dx=15)  # cannot use dx or dy with abs. z
+            self.optic.add_surface(
+                index=1, radius=be.inf, z=0, dx=15
+            )  # cannot use dx or dy with abs. z
diff --git a/tests/test_paraxial.py b/tests/test_paraxial.py
index 4712dde1..1016f7de 100644
--- a/tests/test_paraxial.py
+++ b/tests/test_paraxial.py
@@ -486,162 +486,97 @@ def test_paraxial_init(set_test_backend):
     assert paraxial.surfaces == optic.surface_group
 
 
-@pytest.mark.parametrize(
-    "optic_and_values",
-    get_optic_data(),
-    indirect=True
-)
+@pytest.mark.parametrize("optic_and_values", get_optic_data(), indirect=True)
 def test_calculate_f1(optic_and_values):
     optic_instance, values = optic_and_values
     assert_allclose(optic_instance.paraxial.f1(), values["f1"])
 
 
-@pytest.mark.parametrize(
-    "optic_and_values",
-    get_optic_data(),
-    indirect=True
-)
+@pytest.mark.parametrize("optic_and_values", get_optic_data(), indirect=True)
 def test_calculate_f2(optic_and_values):
     optic_instance, values = optic_and_values
     assert_allclose(optic_instance.paraxial.f2(), values["f2"])
 
 
-@pytest.mark.parametrize(
-    "optic_and_values",
-    get_optic_data(),
-    indirect=True
-)
+@pytest.mark.parametrize("optic_and_values", get_optic_data(), indirect=True)
 def test_calculate_F1(optic_and_values):
     optic_instance, values = optic_and_values
     assert_allclose(optic_instance.paraxial.F1(), values["F1"])
 
 
-@pytest.mark.parametrize(
-    "optic_and_values",
-    get_optic_data(),
-    indirect=True
-)
+@pytest.mark.parametrize("optic_and_values", get_optic_data(), indirect=True)
 def test_calculate_F2(optic_and_values):
     optic_instance, values = optic_and_values
     assert_allclose(optic_instance.paraxial.F2(), values["F2"])
 
 
-@pytest.mark.parametrize(
-    "optic_and_values",
-    get_optic_data(),
-    indirect=True
-)
+@pytest.mark.parametrize("optic_and_values", get_optic_data(), indirect=True)
 def test_calculate_P1(optic_and_values):
     optic_instance, values = optic_and_values
     assert_allclose(optic_instance.paraxial.P1(), values["P1"])
 
 
-@pytest.mark.parametrize(
-    "optic_and_values",
-    get_optic_data(),
-    indirect=True
-)
+@pytest.mark.parametrize("optic_and_values", get_optic_data(), indirect=True)
 def test_calculate_P2(optic_and_values):
     optic_instance, values = optic_and_values
     assert_allclose(optic_instance.paraxial.P2(), values["P2"])
 
 
-@pytest.mark.parametrize(
-    "optic_and_values",
-    get_optic_data(),
-    indirect=True
-)
+@pytest.mark.parametrize("optic_and_values", get_optic_data(), indirect=True)
 def test_calculate_P1anti(optic_and_values):
     optic_instance, values = optic_and_values
     assert_allclose(optic_instance.paraxial.P1anti(), values["P1anti"])
 
 
-
-@pytest.mark.parametrize(
-    "optic_and_values",
-    get_optic_data(),
-    indirect=True
-)
+@pytest.mark.parametrize("optic_and_values", get_optic_data(), indirect=True)
 def test_calculate_P2anti(optic_and_values):
     optic_instance, values = optic_and_values
     assert_allclose(optic_instance.paraxial.P2anti(), values["P2anti"])
 
 
-@pytest.mark.parametrize(
-    "optic_and_values",
-    get_optic_data(),
-    indirect=True
-)
+@pytest.mark.parametrize("optic_and_values", get_optic_data(), indirect=True)
 def test_calculate_N1(optic_and_values):
     optic_instance, values = optic_and_values
     assert_allclose(optic_instance.paraxial.N1(), values["N1"])
 
 
-@pytest.mark.parametrize(
-    "optic_and_values",
-    get_optic_data(),
-    indirect=True
-)
+@pytest.mark.parametrize("optic_and_values", get_optic_data(), indirect=True)
 def test_calculate_N2(optic_and_values):
     optic_instance, values = optic_and_values
     assert_allclose(optic_instance.paraxial.N2(), values["N2"])
 
 
-@pytest.mark.parametrize(
-    "optic_and_values",
-    get_optic_data(),
-    indirect=True
-)
+@pytest.mark.parametrize("optic_and_values", get_optic_data(), indirect=True)
 def test_calculate_EPL(optic_and_values):
     optic_instance, values = optic_and_values
     assert_allclose(optic_instance.paraxial.EPL(), values["EPL"])
 
 
-@pytest.mark.parametrize(
-    "optic_and_values",
-    get_optic_data(),
-    indirect=True
-)
+@pytest.mark.parametrize("optic_and_values", get_optic_data(), indirect=True)
 def test_calculate_EPD(optic_and_values):
     optic_instance, values = optic_and_values
     assert_allclose(optic_instance.paraxial.EPD(), values["EPD"])
 
 
-@pytest.mark.parametrize(
-    "optic_and_values",
-    get_optic_data(),
-    indirect=True
-)
+@pytest.mark.parametrize("optic_and_values", get_optic_data(), indirect=True)
 def test_calculate_XPL(optic_and_values):
     optic_instance, values = optic_and_values
     assert_allclose(optic_instance.paraxial.XPL(), values["XPL"])
 
 
-@pytest.mark.parametrize(
-    "optic_and_values",
-    get_optic_data(),
-    indirect=True
-)
+@pytest.mark.parametrize("optic_and_values", get_optic_data(), indirect=True)
 def test_calculate_XPD(optic_and_values):
     optic_instance, values = optic_and_values
     assert_allclose(optic_instance.paraxial.XPD(), values["XPD"])
 
 
-@pytest.mark.parametrize(
-    "optic_and_values",
-    get_optic_data(),
-    indirect=True
-)
+@pytest.mark.parametrize("optic_and_values", get_optic_data(), indirect=True)
 def test_calculate_FNO(optic_and_values):
     optic_instance, values = optic_and_values
     assert_allclose(optic_instance.paraxial.FNO(), values["FNO"])
 
 
-@pytest.mark.parametrize(
-    "optic_and_values",
-    get_optic_data(),
-    indirect=True
-)
+@pytest.mark.parametrize("optic_and_values", get_optic_data(), indirect=True)
 def test_calculate_invariant(optic_and_values):
     optic_instance, values = optic_and_values
     assert_allclose(optic_instance.paraxial.invariant(), values["invariant"])
@@ -679,7 +614,12 @@ def test_EPD_float_by_stop_size_finite(set_test_backend):
     lens.add_surface(index=0, radius=be.inf, thickness=be.inf)
     lens.add_surface(index=1, radius=22.01359, thickness=3.25896, material="SK16")
     lens.add_surface(index=2, radius=-435.76044, thickness=6.00755)
-    lens.add_surface(index=3, radius=-22.21328, thickness=0.99997, material=("F2", "schott"),)
+    lens.add_surface(
+        index=3,
+        radius=-22.21328,
+        thickness=0.99997,
+        material=("F2", "schott"),
+    )
     lens.add_surface(index=4, radius=20.29192, thickness=4.75041, is_stop=True)
     lens.add_surface(index=5, radius=79.68360, thickness=2.95208, material="SK16")
     lens.add_surface(index=6, radius=-18.39533, thickness=42.20778)
@@ -705,7 +645,12 @@ def test_EPD_float_by_stop_size_infinite(set_test_backend):
     lens.add_surface(index=0, radius=be.inf, thickness=10_000)
     lens.add_surface(index=1, radius=22.01359, thickness=3.25896, material="SK16")
     lens.add_surface(index=2, radius=-435.76044, thickness=6.00755)
-    lens.add_surface(index=3, radius=-22.21328, thickness=0.99997, material=("F2", "schott"),)
+    lens.add_surface(
+        index=3,
+        radius=-22.21328,
+        thickness=0.99997,
+        material=("F2", "schott"),
+    )
     lens.add_surface(index=4, radius=20.29192, thickness=4.75041, is_stop=True)
     lens.add_surface(index=5, radius=79.68360, thickness=2.95208, material="SK16")
     lens.add_surface(index=6, radius=-18.39533, thickness=42.20778)
diff --git a/tests/test_physical_apertures.py b/tests/test_physical_apertures.py
index 6f7d169f..72f2d9cb 100644
--- a/tests/test_physical_apertures.py
+++ b/tests/test_physical_apertures.py
@@ -415,11 +415,14 @@ class TestConfigureAperture:
     def test_none_input(self, set_test_backend):
         assert configure_aperture(None) is None
 
-    @pytest.mark.parametrize("scalar_input, expected_r_max", [
-        (2, 1.0),
-        (0.0, 0.0),
-        (3.5, 1.75),
-    ])
+    @pytest.mark.parametrize(
+        "scalar_input, expected_r_max",
+        [
+            (2, 1.0),
+            (0.0, 0.0),
+            (3.5, 1.75),
+        ],
+    )
     def test_scalar_input(self, set_test_backend, scalar_input, expected_r_max):
         result = configure_aperture(scalar_input)
         assert isinstance(result, RadialAperture)
@@ -429,20 +432,26 @@ class TestConfigureAperture:
         class DummyAperture(BaseAperture):
             def contains(self, x, y):
                 pass
+
             def extent(self):
                 pass
+
             def scale(self, scale_factor):
                 pass
+
         ap = DummyAperture()
         assert configure_aperture(ap) is ap
 
-    @pytest.mark.parametrize("invalid_input", [
-        "circle",
-        [1, 2],
-        {"r": 1},
-        object(),
-        set([1.0]),
-    ])
+    @pytest.mark.parametrize(
+        "invalid_input",
+        [
+            "circle",
+            [1, 2],
+            {"r": 1},
+            object(),
+            set([1.0]),
+        ],
+    )
     def test_invalid_input_raises_value_error(self, set_test_backend, invalid_input):
         with pytest.raises(ValueError, match="Invalid `aperture` provided"):
             configure_aperture(invalid_input)
@@ -451,12 +460,16 @@ class TestConfigureAperture:
         class CustomAperture(BaseAperture):
             def __init__(self):
                 self.special = True
+
             def contains(self, x, y):
                 pass
+
             def extent(self):
                 pass
+
             def scale(self, scale_factor):
                 pass
+
         ap = CustomAperture()
         result = configure_aperture(ap)
         assert result is ap
diff --git a/tests/test_psf.py b/tests/test_psf.py
index c00e7125..54eee416 100644
--- a/tests/test_psf.py
+++ b/tests/test_psf.py
@@ -25,6 +25,7 @@ def make_fftpsf(set_test_backend):
         if tweak_optic:
             tweak_optic(optic)
         return FFTPSF(optic, field, wavelength, num_rays, grid_size)
+
     return _factory
 
 
@@ -41,12 +42,15 @@ def test_strehl_ratio(make_fftpsf):
     assert 0 <= strehl_ratio <= 1
 
 
-@pytest.mark.parametrize("projection, log", [
-    ("2d", False),
-    ("3d", False),
-    ("2d", True),
-    ("3d", True),
-])
+@pytest.mark.parametrize(
+    "projection, log",
+    [
+        ("2d", False),
+        ("3d", False),
+        ("2d", True),
+        ("3d", True),
+    ],
+)
 @patch("matplotlib.pyplot.show")
 def test_view(mock_show, projection, log, make_fftpsf, set_test_backend):
     # Skip for torch since view isn't implemented there
diff --git a/tests/test_rays.py b/tests/test_rays.py
index bd352036..8de2f3c3 100644
--- a/tests/test_rays.py
+++ b/tests/test_rays.py
@@ -14,6 +14,7 @@ from optiland.samples.lithography import UVProjectionLens
 from optiland.samples.objectives import TessarLens
 from tests.utils import assert_allclose
 
+
 def test_translate(set_test_backend):
     rays = BaseRays()
     rays.x = 1.0
@@ -165,7 +166,7 @@ def test_rotate_x(set_test_backend):
     assert_allclose(rays.N[0], -1.0, rtol=0, atol=1e-10)
 
     rays.rotate_x(0.0)
-    
+
     assert_allclose(rays.x[0], 1.0, rtol=0, atol=1e-10)
     assert_allclose(rays.y[0], -2.0, rtol=0, atol=1e-10)
     assert_allclose(rays.z[0], -3.0, rtol=0, atol=1e-10)
@@ -588,7 +589,9 @@ class TestPolarizedRays:
         expected_jones_matrix = be.array(
             [[[1.0, 0.0, 0.0], [0.0, 0.99498744, -0.1], [0.0, 0.1, 0.99498744]]],
         )
-        assert_allclose(rays.p, expected_jones_matrix, atol=1e-8) # Reduced tolerance slightly for potential backend differences
+        assert_allclose(
+            rays.p, expected_jones_matrix, atol=1e-8
+        )  # Reduced tolerance slightly for potential backend differences
 
         # test case when jones = None
         rays.L0 = be.array([0.0])
@@ -641,6 +644,7 @@ class TestPolarizedRays:
         with pytest.raises(ValueError):
             rays._get_3d_electric_field(state)
 
+
 @pytest.mark.usefixtures("set_test_backend")
 class TestRayGenerator:
     def test_generate_rays(self):
diff --git a/tests/test_solves.py b/tests/test_solves.py
index 1121f76d..c7cccf17 100644
--- a/tests/test_solves.py
+++ b/tests/test_solves.py
@@ -70,14 +70,18 @@ class TestMarginalRayHeightSolve:
 class TestQuickfocusSolve:
     def test_quick_focus_solve_constructor(self, set_test_backend):
         optic = CookeTriplet()
-        optic.surface_group.surfaces[-1].geometry.cs.z = optic.surface_group.surfaces[-1].geometry.cs.z - be.array(10)
+        optic.surface_group.surfaces[-1].geometry.cs.z = optic.surface_group.surfaces[
+            -1
+        ].geometry.cs.z - be.array(10)
         solve = solves.QuickFocusSolve(optic)
 
         assert solve.optic == optic
 
     def test_quick_focus_solve_apply(self, set_test_backend):
         optic = CookeTriplet()
-        optic.surface_group.surfaces[-1].geometry.cs.z = optic.surface_group.surfaces[-1].geometry.cs.z - be.array(10)
+        optic.surface_group.surfaces[-1].geometry.cs.z = optic.surface_group.surfaces[
+            -1
+        ].geometry.cs.z - be.array(10)
         thickness = 42.21812063592369
         solve = solves.QuickFocusSolve(optic)
         solve.apply()
diff --git a/tests/test_tolerancing.py b/tests/test_tolerancing.py
index fd75b7c2..5922ead6 100644
--- a/tests/test_tolerancing.py
+++ b/tests/test_tolerancing.py
@@ -7,6 +7,7 @@ from optiland.tolerancing.core import Tolerancing
 from optiland.tolerancing.perturbation import Perturbation, ScalarSampler
 from tests.utils import assert_allclose
 
+
 @pytest.fixture
 def setup_tolerancing():
     optic = Edmund_49_847()
@@ -51,7 +52,7 @@ def test_add_operand_no_target(setup_tolerancing, set_test_backend):
     operand_type = "f2"
     input_data = {"optic": optic}
     weight = 2.0
-    target = optic.paraxial.f2() # Calculate target based on current state
+    target = optic.paraxial.f2()  # Calculate target based on current state
 
     # Ensure original optic data is not modified before adding operand
     original_f2 = optic.paraxial.f2()
@@ -79,17 +80,20 @@ def test_add_perturbation(setup_tolerancing, set_test_backend):
     surface_number = 1
 
     # Ensure original optic data is not modified before adding perturbation
-    original_radius = optic.surface_group.surfaces[surface_number].geometry.radius  
-    
-    tolerancing.add_perturbation(variable_type, perturbation, surface_number=surface_number)
+    original_radius = optic.surface_group.surfaces[surface_number].geometry.radius
+
+    tolerancing.add_perturbation(
+        variable_type, perturbation, surface_number=surface_number
+    )
     assert len(tolerancing.perturbations) == 1
     added_perturbation = tolerancing.perturbations[0]
     assert isinstance(added_perturbation, Perturbation)
     assert added_perturbation.optic == optic
     assert added_perturbation.type == variable_type
     # Check that adding the perturbation definition doesn't apply it yet
-    assert_allclose(optic.surface_group.surfaces[surface_number].geometry.radius, original_radius)
-
+    assert_allclose(
+        optic.surface_group.surfaces[surface_number].geometry.radius, original_radius
+    )
 
 
 def test_add_compensator(setup_tolerancing, set_test_backend):
diff --git a/tests/test_wavefront.py b/tests/test_wavefront.py
index af9df91c..901f6443 100644
--- a/tests/test_wavefront.py
+++ b/tests/test_wavefront.py
@@ -44,7 +44,9 @@ class TestWavefront:
         assert isinstance(w.data[((0.0, 0.0), 0.5876)].intensity, be.ndarray)
         assert isinstance(w.data[((0.0, 0.7), 0.5876)].opd, be.ndarray)
         assert isinstance(w.data[((0.0, 0.0), 0.5876)].pupil_x, be.ndarray)
-        assert be.size(w.data[((0.0, 1.0), 0.6563)].opd) == 469  # num points in the pupil
+        assert (
+            be.size(w.data[((0.0, 1.0), 0.6563)].opd) == 469
+        )  # num points in the pupil
 
     def test_trace_chief_ray(self, set_test_backend):
         optic = DoubleGauss()
@@ -171,6 +173,7 @@ class TestOPD:
         rms = opd.rms()
         assert_allclose(rms, 0.9709788038168692)
 
+
 class TestZernikeOPD:
     def test_zernike_opd_initialization(self, set_test_backend):
         optic = DoubleGauss()
diff --git a/tests/test_zernike.py b/tests/test_zernike.py
index f37f3588..20f1bdf2 100644
--- a/tests/test_zernike.py
+++ b/tests/test_zernike.py
@@ -720,7 +720,7 @@ class TestZernikeFit:
         plt.close()
 
     @patch("matplotlib.pyplot.show")
-    def test_view_residual_standard(self, mock_show,set_test_backend):
+    def test_view_residual_standard(self, mock_show, set_test_backend):
         zernike_fit_standard = zernike.ZernikeFit(
             self.x,
             self.y,
diff --git a/tests/utils.py b/tests/utils.py
index 27b1efaf..406e92c2 100644
--- a/tests/utils.py
+++ b/tests/utils.py
@@ -2,7 +2,7 @@ import numpy as np
 import optiland.backend as be
 
 
-def assert_allclose(a, b, rtol=1.e-5, atol=1.e-7):
+def assert_allclose(a, b, rtol=1.0e-5, atol=1.0e-7):
     """Assert that two arrays or tensors are element-wise equal within
     tolerance.
     """