diff --git a/README.md b/README.md
index 759ddea2..d6c10756 100644
--- a/README.md
+++ b/README.md
@@ -66,7 +66,7 @@ Cherry might be for you if:
 
 ### Now
 - [ ] Systems with beam splitters
-- [ ] Paraxial surface types
+- [X] Thin lens surface types
 - [ ] Surface coatings
 - [ ] Encode designs into URLs and enable sharing via tiny URLs
 
diff --git a/crates/cherry-rs/src/core/math/vec3.rs b/crates/cherry-rs/src/core/math/vec3.rs
index df4c4454..c4d2e970 100644
--- a/crates/cherry-rs/src/core/math/vec3.rs
+++ b/crates/cherry-rs/src/core/math/vec3.rs
@@ -200,6 +200,11 @@ impl Vec3 {
         vecs
     }
 
+    /// Returns `true` if all three components are finite (not infinite or NaN).
+    pub fn is_finite(&self) -> bool {
+        self.e[0].is_finite() && self.e[1].is_finite() && self.e[2].is_finite()
+    }
+
     pub fn approx_eq(&self, rhs: &Self, tol: Float) -> bool {
         (self.e[0] - rhs.e[0]).abs() < tol
             && (self.e[1] - rhs.e[1]).abs() < tol
diff --git a/crates/cherry-rs/src/core/ray.rs b/crates/cherry-rs/src/core/ray.rs
index ac74a839..4267edef 100644
--- a/crates/cherry-rs/src/core/ray.rs
+++ b/crates/cherry-rs/src/core/ray.rs
@@ -1,7 +1,9 @@
 #[cfg(feature = "serde")]
 use serde::Serialize;
 
-use crate::core::{Float, PI, math::vec3::Vec3, sequential_model::placement::Placement};
+use crate::core::{
+    Float, PI, math::vec3::Vec3, sequential_model::surface_placement::SurfacePlacement,
+};
 
 /// A single ray to be traced through an optical system.
 ///
@@ -45,14 +47,14 @@ impl Ray {
 
     /// Transform a ray into the local coordinate system of a surface from the
     /// global system.
-    pub fn transform(&mut self, placement: &Placement) {
+    pub fn transform(&mut self, placement: &SurfacePlacement) {
         self.pos = placement.rotation_matrix * (self.pos - placement.position);
         self.dir = placement.rotation_matrix * self.dir;
     }
 
     /// Transform a ray from the local coordinate system of a surface into the
     /// global system.
-    pub fn i_transform(&mut self, placement: &Placement) {
+    pub fn i_transform(&mut self, placement: &SurfacePlacement) {
         self.pos = (placement.inv_rotation_matrix * self.pos) + placement.position;
         self.dir = placement.inv_rotation_matrix * self.dir;
     }
diff --git a/crates/cherry-rs/src/core/sequential_model/mod.rs b/crates/cherry-rs/src/core/sequential_model/mod.rs
index 1714d401..adebd583 100644
--- a/crates/cherry-rs/src/core/sequential_model/mod.rs
+++ b/crates/cherry-rs/src/core/sequential_model/mod.rs
@@ -1,22 +1,24 @@
 /// Data types for modeling sequential ray tracing systems.
 pub mod builder;
 pub(crate) mod cursor;
-pub mod placement;
 pub mod solves;
+pub mod surface_placement;
 
 use std::ops::Range;
 
 use anyhow::{Result, anyhow};
 
 use self::cursor::Cursor;
-use self::placement::Placement;
+use self::surface_placement::SurfacePlacement;
 #[cfg(feature = "serde")]
 use crate::core::surfaces::SurfaceRegistry;
 use crate::core::{
     Float,
     math::{linalg::mat3x3::Mat3x3, vec3::Vec3},
     refractive_index::RefractiveIndex,
-    surfaces::{BeamSplitter, Conic, Image, Iris, Object, Probe, Sphere, Surface, SurfaceKind},
+    surfaces::{
+        BeamSplitter, Conic, Image, Iris, Object, Probe, Sphere, Surface, SurfaceKind, ThinLens,
+    },
 };
 use crate::specs::surfaces::PlacementSpec;
 use crate::specs::{
@@ -25,37 +27,48 @@ use crate::specs::{
     surfaces::{BeamSplitterPathKind, BoundaryKind, SurfaceSpec},
 };
 
-/// Cursor forward direction at each surface.
-type AxisDirections = Vec<Vec3>;
-/// Nominal on-axis cursor position at each surface, before any decenter is
-/// applied.
-type CursorPositions = Vec<Vec3>;
-
-type SurfsPlacementsDirs = (
+type SurfaceStoreContents = (
     Vec<Box<dyn Surface>>,
-    Vec<Placement>,
-    AxisDirections,
-    CursorPositions,
+    Vec<SurfacePlacement>,
+    Vec<CursorPlacement>,
 );
 
 /// Owns all surface objects and their computed placements.
 #[derive(Debug)]
 struct SurfaceStore {
     surfaces: Vec<Box<dyn Surface>>,
-    placements: Vec<Placement>,
-    /// Cursor forward direction at each surface vertex (first-path walk).
-    axis_directions: Vec<Vec3>,
-    /// Nominal on-axis cursor position at each surface before any decenter
-    /// (first-path walk).
-    cursor_positions: Vec<Vec3>,
+    placements: Vec<SurfacePlacement>,
+}
+
+/// Per-step cursor state recorded along one optical path.
+///
+/// Captures the cursor orientation, position, and rotation matrix before any
+/// reflection at each surface, for both `New` and `Shared` steps. Carried by
+/// each iterator [`Step`] so that views have fully path-specific cursor data
+/// without separately querying the model.
+#[derive(Debug, Clone, Copy)]
+pub struct CursorPlacement {
+    /// Cursor forward direction (unit vector) as the beam approaches this
+    /// surface.
+    pub axis_direction: Vec3,
+    /// Nominal on-axis cursor position before any decenter.
+    pub cursor_position: Vec3,
+    /// Rotation from the global frame into the cursor frame at this step.
+    pub cursor_rotation_matrix: Mat3x3,
 }
 
 /// One optical path through the system.
 #[derive(Debug)]
 struct OpticalPath {
+    /// Ordered store indices visited by this path, one per step.
+    surface_indices: Vec<usize>,
+    /// Beam-splitter arm kind per step; parallel to `surface_indices`.
+    beam_splitter_arms: Vec<Option<BeamSplitterPathKind>>,
     submodels: Vec<SequentialSubModelBase>,
     /// User-specified aperture stop as a store index, or `None` for auto.
     stop_surface: Option<usize>,
+    /// Step-indexed cursor state, parallel to `surface_indices`.
+    steps: Vec<CursorPlacement>,
 }
 
 /// A gap between two surfaces in a sequential system.
@@ -160,19 +173,23 @@ pub trait SequentialSubModel {
     fn try_iter<'a>(
         &'a self,
         surfaces: &'a [Box<dyn Surface>],
-        placements: &'a [Placement],
+        placements: &'a [SurfacePlacement],
+        surface_indices: &'a [usize],
+        beam_splitter_arms: &'a [Option<BeamSplitterPathKind>],
+        path_steps: &'a [CursorPlacement],
     ) -> Result<SequentialSubModelIter<'a>>;
 
-    fn slice(&self, idx: Range<usize>) -> SequentialSubModelSlice<'_>;
+    fn slice<'a>(
+        &'a self,
+        idx: Range<usize>,
+        surface_indices: &'a [usize],
+        beam_splitter_arms: &'a [Option<BeamSplitterPathKind>],
+    ) -> SequentialSubModelSlice<'a>;
 }
 
 #[derive(Debug)]
 pub struct SequentialSubModelBase {
-    /// Ordered sequence of store indices visited by this path.
-    surface_indices: Vec<usize>,
     gaps: Vec<Gap>,
-    /// Dense arm kind per step; parallel to `surface_indices`.
-    beam_splitter_arms: Vec<Option<BeamSplitterPathKind>>,
 }
 
 /// A view of a single submodel in a sequential system.
@@ -183,6 +200,10 @@ pub struct SequentialSubModelSlice<'a> {
     surface_indices: &'a [usize],
     beam_splitter_arms: &'a [Option<BeamSplitterPathKind>],
     gaps: &'a [Gap],
+    /// First step index of this slice within the full path's step list.
+    /// Used by the iterator to correctly index into the caller-supplied
+    /// `path_steps` slice.
+    step_offset: usize,
 }
 
 /// An iterator over the surfaces and gaps in a submodel.
@@ -190,20 +211,24 @@ pub struct SequentialSubModelSlice<'a> {
 /// Most operations in sequential modeling involve use of this iterator.
 pub struct SequentialSubModelIter<'a> {
     surfaces: &'a [Box<dyn Surface>],
-    placements: &'a [Placement],
+    placements: &'a [SurfacePlacement],
     surface_indices: &'a [usize],
     beam_splitter_arms: &'a [Option<BeamSplitterPathKind>],
     gaps: &'a [Gap],
+    path_steps: &'a [CursorPlacement],
+    step_offset: usize,
     index: usize,
 }
 
 /// A reverse iterator over the surfaces and gaps in a submodel.
 pub struct SequentialSubModelReverseIter<'a> {
     surfaces: &'a [Box<dyn Surface>],
-    placements: &'a [Placement],
+    placements: &'a [SurfacePlacement],
     surface_indices: &'a [usize],
     beam_splitter_arms: &'a [Option<BeamSplitterPathKind>],
     gaps: &'a [Gap],
+    path_steps: &'a [CursorPlacement],
+    step_offset: usize,
     index: usize,
 }
 
@@ -215,10 +240,13 @@ pub struct Step<'a> {
     pub gap_before: &'a Gap,
     pub surface: &'a dyn Surface,
     pub gap_after: Option<&'a Gap>,
-    pub placement: &'a Placement,
+    pub surface_placement: &'a SurfacePlacement,
     /// The beam-splitter arm traversed at this step, or `None` for non-BS
     /// surfaces. Set from the path's `beam_splitter_arms` declaration.
     pub bs_arm: Option<BeamSplitterPathKind>,
+    /// Path-specific cursor state at this step (orientation, position, rotation
+    /// matrix of the cursor as it arrives at this surface).
+    pub cursor_placement: CursorPlacement,
 }
 
 /// Propagates a tangential direction unit vector through the mirror surfaces of
@@ -232,14 +260,16 @@ pub struct Step<'a> {
 pub(crate) fn propagate_tangential_vec(
     v_init: Vec3,
     surfaces: &[Box<dyn Surface>],
-    placements: &[Placement],
+    placements: &[SurfacePlacement],
+    surface_indices: &[usize],
 ) -> Vec<Vec3> {
     use crate::specs::surfaces::BoundaryKind;
     let mut v = v_init;
-    surfaces
+    surface_indices
         .iter()
-        .zip(placements.iter())
-        .map(|(surf, placement)| {
+        .map(|&idx| {
+            let surf = &surfaces[idx];
+            let placement = &placements[idx];
             let v_incident = v;
             if let BoundaryKind::Reflecting = surf.boundary_kind() {
                 // Normal in global frame derived from the *nominal* orientation
@@ -257,24 +287,34 @@ pub(crate) fn propagate_tangential_vec(
         .collect()
 }
 
-/// Returns the index of the first physical surface in the system.
+/// Returns the step index of the first physical surface visited by
+/// `surface_indices`.
 ///
-/// A physical surface is one that has a finite semi-diameter,
-/// i.e., a Conic or Iris. Object, Image, and Probe surfaces are excluded.
-pub(crate) fn first_physical_surface(surfaces: &[Box<dyn Surface>]) -> Option<usize> {
-    surfaces
+/// A physical surface has a finite semi-diameter (Conic or Iris). Object,
+/// Image, and Probe surfaces are excluded. Returns the position within
+/// `surface_indices` (i.e. a step index), not a store index.
+pub(crate) fn first_physical_step(
+    surface_indices: &[usize],
+    surfaces: &[Box<dyn Surface>],
+) -> Option<usize> {
+    surface_indices
         .iter()
-        .position(|surf| surf.mask().semi_diameter().is_finite())
+        .position(|&i| surfaces[i].mask().semi_diameter().is_finite())
 }
 
-/// Returns the index of the last physical surface in the system.
+/// Returns the step index of the last physical surface visited by
+/// `surface_indices`.
 ///
-/// A physical surface is one that limits the has a finite semi-diameter,
-/// i.e., a Conic or Iris. Object, Image, and Probe surfaces are excluded.
-pub fn last_physical_surface(surfaces: &[Box<dyn Surface>]) -> Option<usize> {
-    surfaces
+/// A physical surface has a finite semi-diameter (Conic or Iris). Object,
+/// Image, and Probe surfaces are excluded. Returns the position within
+/// `surface_indices` (i.e. a step index), not a store index.
+pub(crate) fn last_physical_step(
+    surface_indices: &[usize],
+    surfaces: &[Box<dyn Surface>],
+) -> Option<usize> {
+    surface_indices
         .iter()
-        .rposition(|surf| surf.mask().semi_diameter().is_finite())
+        .rposition(|&i| surfaces[i].mask().semi_diameter().is_finite())
 }
 
 /// Returns the id of a surface in a reversed system.
@@ -328,7 +368,7 @@ impl SequentialModel {
         #[cfg(not(feature = "serde"))]
         {
             Self::validate_specs(gap_specs, wavelengths)?;
-            let (surfaces, placements, axis_directions, cursor_positions) =
+            let (surfaces, placements, cursor_placements) =
                 Self::surf_specs_to_surfs(surface_specs, gap_specs)?;
             if let Some(i) = stop_surface {
                 Self::validate_stop_surface(&surfaces, i)?;
@@ -338,21 +378,18 @@ impl SequentialModel {
             let mut submodels: Vec<SequentialSubModelBase> = Vec::new();
             for &wavelength in wavelengths.iter() {
                 let gaps = Self::gap_specs_to_gaps(gap_specs, wavelength)?;
-                submodels.push(SequentialSubModelBase::new(
-                    surface_indices.clone(),
-                    gaps,
-                    bs_arms.clone(),
-                ));
+                submodels.push(SequentialSubModelBase::new(gaps));
             }
             let store = SurfaceStore {
                 surfaces,
                 placements,
-                axis_directions,
-                cursor_positions,
             };
             let path = OpticalPath {
+                surface_indices,
+                beam_splitter_arms: bs_arms,
                 submodels,
                 stop_surface,
+                steps: cursor_placements,
             };
             Ok(Self {
                 store,
@@ -375,7 +412,7 @@ impl SequentialModel {
         registry: Option<&SurfaceRegistry>,
     ) -> Result<Self> {
         Self::validate_specs(gap_specs, wavelengths)?;
-        let (surfaces, placements, axis_directions, cursor_positions) =
+        let (surfaces, placements, cursor_placements) =
             Self::surf_specs_to_surfs(surface_specs, gap_specs, registry)?;
         if let Some(i) = stop_surface {
             Self::validate_stop_surface(&surfaces, i)?;
@@ -385,21 +422,18 @@ impl SequentialModel {
         let mut submodels: Vec<SequentialSubModelBase> = Vec::new();
         for &wavelength in wavelengths.iter() {
             let gaps = Self::gap_specs_to_gaps(gap_specs, wavelength)?;
-            submodels.push(SequentialSubModelBase::new(
-                surface_indices.clone(),
-                gaps,
-                bs_arms.clone(),
-            ));
+            submodels.push(SequentialSubModelBase::new(gaps));
         }
         let store = SurfaceStore {
             surfaces,
             placements,
-            axis_directions,
-            cursor_positions,
         };
         let path = OpticalPath {
+            surface_indices,
+            beam_splitter_arms: bs_arms,
             submodels,
             stop_surface,
+            steps: cursor_placements,
         };
         Ok(Self {
             store,
@@ -446,7 +480,7 @@ impl SequentialModel {
         }
         Self::validate_specs(gap_specs, wavelengths)?;
 
-        let (placements, axis_directions, cursor_positions) =
+        let (placements, cursor_placements) =
             Self::build_placements_and_directions(&surfaces, placement_specs, gap_specs);
 
         if let Some(i) = stop_surface {
@@ -458,22 +492,19 @@ impl SequentialModel {
         let mut submodels: Vec<SequentialSubModelBase> = Vec::new();
         for &wavelength in wavelengths.iter() {
             let gaps = Self::gap_specs_to_gaps(gap_specs, wavelength)?;
-            submodels.push(SequentialSubModelBase::new(
-                surface_indices.clone(),
-                gaps,
-                bs_arms.clone(),
-            ));
+            submodels.push(SequentialSubModelBase::new(gaps));
         }
 
         let store = SurfaceStore {
             surfaces,
             placements,
-            axis_directions,
-            cursor_positions,
         };
         let path = OpticalPath {
+            surface_indices,
+            beam_splitter_arms: bs_arms,
             submodels,
             stop_surface,
+            steps: cursor_placements,
         };
         Ok(Self {
             store,
@@ -498,9 +529,7 @@ impl SequentialModel {
         }
 
         let mut store_surfaces: Vec<Box<dyn Surface>> = Vec::new();
-        let mut store_placements: Vec<Placement> = Vec::new();
-        let mut store_axis_directions: Vec<Vec3> = Vec::new();
-        let mut store_cursor_positions: Vec<Vec3> = Vec::new();
+        let mut store_placements: Vec<SurfacePlacement> = Vec::new();
         let mut optical_paths: Vec<OpticalPath> = Vec::new();
 
         for ps in paths {
@@ -524,6 +553,7 @@ impl SequentialModel {
 
             let mut cursor = Cursor::new(-ps.gaps[0].thickness);
             let mut surface_indices: Vec<usize> = Vec::new();
+            let mut path_steps: Vec<CursorPlacement> = Vec::new();
 
             for (step, sref) in ps.surface_refs.iter().enumerate() {
                 let is_first = step == 0;
@@ -553,6 +583,13 @@ impl SequentialModel {
                 };
                 dense_bs_arms.push(bs_arm);
 
+                // Record cursor state before any reflection at this surface.
+                path_steps.push(CursorPlacement {
+                    axis_direction: cursor.forward(),
+                    cursor_position: cursor.pos(),
+                    cursor_rotation_matrix: cursor.rotation_matrix(),
+                });
+
                 match sref {
                     PathSurfaceRef::New(spec) => {
                         let surface = build_surface(spec)?;
@@ -571,12 +608,9 @@ impl SequentialModel {
                             ));
                         }
 
-                        store_axis_directions.push(cursor.forward());
-                        store_cursor_positions.push(cursor.pos());
-
                         let nominal_rot = spec.rotation().rotation_matrix();
                         let actual_rot = spec.rotation_offset().rotation_matrix() * nominal_rot;
-                        let placement = Placement::from_decenter_and_rotation(
+                        let placement = SurfacePlacement::from_decenter_and_rotation(
                             spec.decenter(),
                             actual_rot,
                             nominal_rot,
@@ -654,15 +688,14 @@ impl SequentialModel {
             let mut submodels: Vec<SequentialSubModelBase> = Vec::new();
             for &wavelength in wavelengths.iter() {
                 let gaps = Self::gap_specs_to_gaps(&ps.gaps, wavelength)?;
-                submodels.push(SequentialSubModelBase::new(
-                    surface_indices.clone(),
-                    gaps,
-                    dense_bs_arms.clone(),
-                ));
+                submodels.push(SequentialSubModelBase::new(gaps));
             }
             optical_paths.push(OpticalPath {
+                surface_indices,
+                beam_splitter_arms: dense_bs_arms,
                 submodels,
                 stop_surface,
+                steps: path_steps,
             });
         }
 
@@ -673,8 +706,6 @@ impl SequentialModel {
         let store = SurfaceStore {
             surfaces: store_surfaces,
             placements: store_placements,
-            axis_directions: store_axis_directions,
-            cursor_positions: store_cursor_positions,
         };
         Ok(Self {
             store,
@@ -711,16 +742,20 @@ impl SequentialModel {
         self.paths.len()
     }
 
-    /// Ordered store indices visited by path `path_id` (identical for every
-    /// wavelength submodel within that path).
+    /// Ordered store indices visited by path `path_id`.
     pub fn path_surface_indices(&self, path_id: usize) -> &[usize] {
-        &self.paths[path_id].submodels[0].surface_indices
+        &self.paths[path_id].surface_indices
+    }
+
+    /// Beam-splitter arm kind per step for path `path_id`.
+    pub fn path_beam_splitter_arms(&self, path_id: usize) -> &[Option<BeamSplitterPathKind>] {
+        &self.paths[path_id].beam_splitter_arms
     }
 
     /// Placement of the surface at `step` in path `path_id`, looked up from
     /// the store via that path's `surface_indices`.
-    pub fn path_placement(&self, path_id: usize, step: usize) -> &Placement {
-        let store_idx = self.paths[path_id].submodels[0].surface_indices[step];
+    pub fn path_placement(&self, path_id: usize, step: usize) -> &SurfacePlacement {
+        let store_idx = self.paths[path_id].surface_indices[step];
         &self.store.placements[store_idx]
     }
 
@@ -759,6 +794,20 @@ impl SequentialModel {
         self.paths[path_id].stop_surface
     }
 
+    /// Returns the per-step cursor data for path `path_id`.
+    ///
+    /// Each entry corresponds to one surface in that path's traversal order,
+    /// and records the cursor state as the beam *approaches* that surface
+    /// (before any reflection).
+    pub fn path_steps(&self, path_id: usize) -> &[CursorPlacement] {
+        &self.paths[path_id].steps
+    }
+
+    /// Returns all wavelength submodels for path `path_id`.
+    pub fn submodels_for_path(&self, path_id: usize) -> &[SequentialSubModelBase] {
+        &self.paths[path_id].submodels
+    }
+
     /// Returns the largest semi-diameter of any surface in the system.
     ///
     /// This ignores surfaces without any size, such as object, probe, and image
@@ -786,7 +835,7 @@ impl SequentialModel {
     ///
     /// The i-th placement corresponds to the i-th surface returned by
     /// [`surfaces()`](Self::surfaces).
-    pub fn placements(&self) -> &[Placement] {
+    pub fn placements(&self) -> &[SurfacePlacement] {
         &self.store.placements
     }
 
@@ -814,15 +863,6 @@ impl SequentialModel {
         &self.wavelengths
     }
 
-    /// Returns the optical axis directions at each surface vertex.
-    pub fn axis_directions(&self) -> &[Vec3] {
-        &self.store.axis_directions
-    }
-
-    pub fn cursor_positions(&self) -> &[Vec3] {
-        &self.store.cursor_positions
-    }
-
     fn gap_specs_to_gaps(gap_specs: &[GapSpec], wavelength: Float) -> Result<Vec<Gap>> {
         let mut gaps = Vec::new();
         for gap_spec in gap_specs.iter() {
@@ -835,15 +875,24 @@ impl SequentialModel {
     /// Returns true if the system is rotationally symmetric about the optical
     /// axis.
     ///
-    /// A system is rotationally symmetric if no physical surface has a tilt
-    /// relative to the optical axis, i.e., the surface-tilt rotation equals
-    /// the cursor rotation at every physical surface.
-    pub fn is_rotationally_symmetric(placements: &[Placement]) -> bool {
-        !placements.iter().any(|p| {
-            // R_surf = surface_tilt × cursor = global_to_local · cursor_to_global
-            let r_surf = p.rotation_matrix * p.cursor_rotation_matrix.transpose();
-            !r_surf.approx_eq(&Mat3x3::identity(), 1e-10)
-        })
+    /// A system is rotationally symmetric if no surface has a tilt relative to
+    /// the cursor approaching it on any path, i.e., the surface-tilt rotation
+    /// equals the cursor rotation at every step across all paths.
+    pub fn is_rotationally_symmetric(&self) -> bool {
+        let placements = &self.store.placements;
+        for path_id in 0..self.path_count() {
+            let steps = self.path_steps(path_id);
+            let indices = self.path_surface_indices(path_id);
+            for (step, &idx) in steps.iter().zip(indices.iter()) {
+                let p = &placements[idx];
+                // R_surf = surface_tilt × cursor = global_to_local · cursor_to_global
+                let r_surf = p.rotation_matrix * step.cursor_rotation_matrix.transpose();
+                if !r_surf.approx_eq(&Mat3x3::identity(), 1e-10) {
+                    return false;
+                }
+            }
+        }
+        true
     }
 
     /// Walks the cursor through the system, building placements and axis
@@ -855,10 +904,9 @@ impl SequentialModel {
         surfaces: &[Box<dyn Surface>],
         surface_placements: &[PlacementSpec],
         gap_specs: &[GapSpec],
-    ) -> (Vec<Placement>, Vec<Vec3>, Vec<Vec3>) {
+    ) -> (Vec<SurfacePlacement>, Vec<CursorPlacement>) {
         let mut placements = Vec::new();
-        let mut axis_directions = Vec::new();
-        let mut cursor_positions = Vec::new();
+        let mut cursor_placements = Vec::new();
         let mut cursor = Cursor::new(-gap_specs[0].thickness);
 
         // Surfaces 0 to N-2 (each paired with a gap that follows it).
@@ -867,12 +915,15 @@ impl SequentialModel {
             .zip(surface_placements.iter())
             .zip(gap_specs.iter())
         {
-            axis_directions.push(cursor.forward());
-            cursor_positions.push(cursor.pos());
+            cursor_placements.push(CursorPlacement {
+                axis_direction: cursor.forward(),
+                cursor_position: cursor.pos(),
+                cursor_rotation_matrix: cursor.rotation_matrix(),
+            });
 
             let nominal_rot = sp.rotation.rotation_matrix();
             let actual_rot = sp.rotation_offset.rotation_matrix() * nominal_rot;
-            let placement = Placement::from_decenter_and_rotation(
+            let placement = SurfacePlacement::from_decenter_and_rotation(
                 sp.decenter,
                 actual_rot,
                 nominal_rot,
@@ -893,19 +944,22 @@ impl SequentialModel {
         }
 
         // Last surface - no gap after it.
-        axis_directions.push(cursor.forward());
-        cursor_positions.push(cursor.pos());
+        cursor_placements.push(CursorPlacement {
+            axis_direction: cursor.forward(),
+            cursor_position: cursor.pos(),
+            cursor_rotation_matrix: cursor.rotation_matrix(),
+        });
         let sp = surface_placements.last().expect("at least one surface");
         let nominal_rot = sp.rotation.rotation_matrix();
         let actual_rot = sp.rotation_offset.rotation_matrix() * nominal_rot;
-        placements.push(Placement::from_decenter_and_rotation(
+        placements.push(SurfacePlacement::from_decenter_and_rotation(
             sp.decenter,
             actual_rot,
             nominal_rot,
             &cursor,
         ));
 
-        (placements, axis_directions, cursor_positions)
+        (placements, cursor_placements)
     }
 
     #[cfg(feature = "serde")]
@@ -913,7 +967,7 @@ impl SequentialModel {
         surf_specs: &[SurfaceSpec],
         gap_specs: &[GapSpec],
         registry: Option<&SurfaceRegistry>,
-    ) -> Result<SurfsPlacementsDirs> {
+    ) -> Result<SurfaceStoreContents> {
         let surfaces: Vec<Box<dyn Surface>> = surf_specs
             .iter()
             .map(|s| surface_from_spec(s, registry))
@@ -926,16 +980,16 @@ impl SequentialModel {
                 rotation_offset: spec.rotation_offset(),
             })
             .collect();
-        let (placements, axis_directions, cursor_positions) =
+        let (placements, cursor_placements) =
             Self::build_placements_and_directions(&surfaces, &surface_placements, gap_specs);
-        Ok((surfaces, placements, axis_directions, cursor_positions))
+        Ok((surfaces, placements, cursor_placements))
     }
 
     #[cfg(not(feature = "serde"))]
     fn surf_specs_to_surfs(
         surf_specs: &[SurfaceSpec],
         gap_specs: &[GapSpec],
-    ) -> Result<SurfsPlacementsDirs> {
+    ) -> Result<SurfaceStoreContents> {
         let surfaces: Vec<Box<dyn Surface>> = surf_specs
             .iter()
             .map(surface_from_spec)
@@ -948,9 +1002,9 @@ impl SequentialModel {
                 rotation_offset: spec.rotation_offset(),
             })
             .collect();
-        let (placements, axis_directions, cursor_positions) =
+        let (placements, cursor_placements) =
             Self::build_placements_and_directions(&surfaces, &surface_placements, gap_specs);
-        Ok((surfaces, placements, axis_directions, cursor_positions))
+        Ok((surfaces, placements, cursor_placements))
     }
 
     fn validate_gaps(gaps: &[GapSpec]) -> Result<()> {
@@ -976,16 +1030,8 @@ impl SequentialModel {
 }
 
 impl SequentialSubModelBase {
-    pub(crate) fn new(
-        surface_indices: Vec<usize>,
-        gaps: Vec<Gap>,
-        beam_splitter_arms: Vec<Option<BeamSplitterPathKind>>,
-    ) -> Self {
-        Self {
-            surface_indices,
-            gaps,
-            beam_splitter_arms,
-        }
+    pub(crate) fn new(gaps: Vec<Gap>) -> Self {
+        Self { gaps }
     }
 }
 
@@ -1005,23 +1051,34 @@ impl SequentialSubModel for SequentialSubModelBase {
     fn try_iter<'a>(
         &'a self,
         surfaces: &'a [Box<dyn Surface>],
-        placements: &'a [Placement],
+        placements: &'a [SurfacePlacement],
+        surface_indices: &'a [usize],
+        beam_splitter_arms: &'a [Option<BeamSplitterPathKind>],
+        path_steps: &'a [CursorPlacement],
     ) -> Result<SequentialSubModelIter<'a>> {
         SequentialSubModelIter::new(
             surfaces,
             placements,
-            &self.surface_indices,
-            &self.beam_splitter_arms,
+            surface_indices,
+            beam_splitter_arms,
             &self.gaps,
+            path_steps,
+            0,
         )
     }
 
-    fn slice(&self, idx: Range<usize>) -> SequentialSubModelSlice<'_> {
+    fn slice<'a>(
+        &'a self,
+        idx: Range<usize>,
+        surface_indices: &'a [usize],
+        beam_splitter_arms: &'a [Option<BeamSplitterPathKind>],
+    ) -> SequentialSubModelSlice<'a> {
         let si_range = idx.start..=idx.end;
         SequentialSubModelSlice {
-            surface_indices: &self.surface_indices[si_range.clone()],
-            beam_splitter_arms: &self.beam_splitter_arms[si_range],
-            gaps: &self.gaps[idx],
+            surface_indices: &surface_indices[si_range.clone()],
+            beam_splitter_arms: &beam_splitter_arms[si_range],
+            gaps: &self.gaps[idx.clone()],
+            step_offset: idx.start,
         }
     }
 }
@@ -1042,7 +1099,10 @@ impl SequentialSubModel for SequentialSubModelSlice<'_> {
     fn try_iter<'b>(
         &'b self,
         surfaces: &'b [Box<dyn Surface>],
-        placements: &'b [Placement],
+        placements: &'b [SurfacePlacement],
+        _surface_indices: &'b [usize],
+        _beam_splitter_arms: &'b [Option<BeamSplitterPathKind>],
+        path_steps: &'b [CursorPlacement],
     ) -> Result<SequentialSubModelIter<'b>> {
         SequentialSubModelIter::new(
             surfaces,
@@ -1050,15 +1110,23 @@ impl SequentialSubModel for SequentialSubModelSlice<'_> {
             self.surface_indices,
             self.beam_splitter_arms,
             self.gaps,
+            path_steps,
+            self.step_offset,
         )
     }
 
-    fn slice(&self, idx: Range<usize>) -> SequentialSubModelSlice<'_> {
+    fn slice<'b>(
+        &'b self,
+        idx: Range<usize>,
+        _surface_indices: &'b [usize],
+        _beam_splitter_arms: &'b [Option<BeamSplitterPathKind>],
+    ) -> SequentialSubModelSlice<'b> {
         let si_range = idx.start..=idx.end;
         SequentialSubModelSlice {
             surface_indices: &self.surface_indices[si_range.clone()],
             beam_splitter_arms: &self.beam_splitter_arms[si_range],
-            gaps: &self.gaps[idx],
+            gaps: &self.gaps[idx.clone()],
+            step_offset: self.step_offset + idx.start,
         }
     }
 }
@@ -1066,10 +1134,12 @@ impl SequentialSubModel for SequentialSubModelSlice<'_> {
 impl<'a> SequentialSubModelIter<'a> {
     fn new(
         surfaces: &'a [Box<dyn Surface>],
-        placements: &'a [Placement],
+        placements: &'a [SurfacePlacement],
         surface_indices: &'a [usize],
         beam_splitter_arms: &'a [Option<BeamSplitterPathKind>],
         gaps: &'a [Gap],
+        path_steps: &'a [CursorPlacement],
+        step_offset: usize,
     ) -> Result<Self> {
         if surface_indices.len() != gaps.len() + 1 {
             return Err(anyhow!(
@@ -1083,6 +1153,8 @@ impl<'a> SequentialSubModelIter<'a> {
             surface_indices,
             beam_splitter_arms,
             gaps,
+            path_steps,
+            step_offset,
             index: 0,
         })
     }
@@ -1094,6 +1166,8 @@ impl<'a> SequentialSubModelIter<'a> {
             self.surface_indices,
             self.beam_splitter_arms,
             self.gaps,
+            self.path_steps,
+            self.step_offset,
         )
     }
 }
@@ -1108,22 +1182,25 @@ impl<'a> Iterator for SequentialSubModelIter<'a> {
         let path_pos = self.index + 1;
         let store_idx = self.surface_indices[path_pos];
         let bs_arm = self.beam_splitter_arms[path_pos];
+        let cursor_placement = self.path_steps[self.step_offset + path_pos];
         let result = if self.index == self.gaps.len() - 1 {
             // We are at the image space gap
             Step {
                 gap_before: &self.gaps[self.index],
                 surface: self.surfaces[store_idx].as_ref(),
                 gap_after: None,
-                placement: &self.placements[store_idx],
+                surface_placement: &self.placements[store_idx],
                 bs_arm,
+                cursor_placement,
             }
         } else {
             Step {
                 gap_before: &self.gaps[self.index],
                 surface: self.surfaces[store_idx].as_ref(),
                 gap_after: Some(&self.gaps[self.index + 1]),
-                placement: &self.placements[store_idx],
+                surface_placement: &self.placements[store_idx],
                 bs_arm,
+                cursor_placement,
             }
         };
         self.index += 1;
@@ -1140,10 +1217,12 @@ impl ExactSizeIterator for SequentialSubModelIter<'_> {
 impl<'a> SequentialSubModelReverseIter<'a> {
     fn new(
         surfaces: &'a [Box<dyn Surface>],
-        placements: &'a [Placement],
+        placements: &'a [SurfacePlacement],
         surface_indices: &'a [usize],
         beam_splitter_arms: &'a [Option<BeamSplitterPathKind>],
         gaps: &'a [Gap],
+        path_steps: &'a [CursorPlacement],
+        step_offset: usize,
     ) -> Result<Self> {
         // Note that this requirement is different than the forward iterator.
         if surface_indices.len() != gaps.len() + 1 {
@@ -1158,6 +1237,8 @@ impl<'a> SequentialSubModelReverseIter<'a> {
             surface_indices,
             beam_splitter_arms,
             gaps,
+            path_steps,
+            step_offset,
             // We will never iterate from the image space surface in reverse.
             index: 1,
         })
@@ -1174,13 +1255,15 @@ impl<'a> Iterator for SequentialSubModelReverseIter<'a> {
         if self.index < n {
             let store_idx = self.surface_indices[forward_index];
             let bs_arm = self.beam_splitter_arms[forward_index];
+            let cursor_placement = self.path_steps[self.step_offset + forward_index];
             // We are somewhere in the middle of the system or at the object space gap.
             let result = Some(Step {
                 gap_before: &self.gaps[forward_index],
                 surface: self.surfaces[store_idx].as_ref(),
                 gap_after: Some(&self.gaps[forward_index - 1]),
-                placement: &self.placements[store_idx],
+                surface_placement: &self.placements[store_idx],
                 bs_arm,
+                cursor_placement,
             });
             self.index += 1;
             result
@@ -1219,6 +1302,11 @@ pub(crate) fn surface_from_spec(
             *radius_of_curvature,
             *surf_kind,
         ))),
+        SurfaceSpec::ThinLens {
+            semi_diameter,
+            focal_length,
+            ..
+        } => Ok(Box::new(ThinLens::new(*semi_diameter, *focal_length))),
         SurfaceSpec::Custom {
             type_id, params, ..
         } => registry
@@ -1265,6 +1353,11 @@ pub(crate) fn surface_from_spec(spec: &SurfaceSpec) -> Result<Box<dyn Surface>>
             *radius_of_curvature,
             *surf_kind,
         ))),
+        SurfaceSpec::ThinLens {
+            semi_diameter,
+            focal_length,
+            ..
+        } => Ok(Box::new(ThinLens::new(*semi_diameter, *focal_length))),
         SurfaceSpec::Image { .. } => Ok(Box::new(Image::new())),
         SurfaceSpec::Object => Ok(Box::new(Object::new())),
         SurfaceSpec::Probe { .. } => Ok(Box::new(Probe::new())),
@@ -1283,36 +1376,37 @@ mod tests {
         specs::surfaces::BoundaryKind,
     };
 
-    // Helper: build a Placement for a surface with the given rotation, in an
-    // identity cursor frame (cursor aligned with global axes, origin at (0,0,0)).
-    fn placement_with_rotation(rotation: Rotation3D) -> Placement {
+    // Helper: build a SurfacePlacement for a surface with the given rotation, in
+    // an identity cursor frame (cursor aligned with global axes, origin at
+    // (0,0,0)).
+    fn placement_with_rotation(rotation: Rotation3D) -> (SurfacePlacement, Mat3x3) {
         let cursor_rotation_matrix = Mat3x3::identity();
         let rotation_matrix = rotation.rotation_matrix() * cursor_rotation_matrix;
-        Placement::new(
+        let sp = SurfacePlacement::new(
             Vec3::new(0.0, 0.0, 0.0),
             0.0,
             rotation_matrix,
             rotation_matrix,
-            cursor_rotation_matrix,
-        )
+        );
+        (sp, cursor_rotation_matrix)
     }
 
     #[test]
     fn projected_sd_untilted_surface() {
         let r = 10.0;
-        let placement = placement_with_rotation(Rotation3D::None);
+        let (placement, crm) = placement_with_rotation(Rotation3D::None);
         let tol = 1e-12;
         let v_u = Vec3::new(0.0, 1.0, 0.0);
         let v_r = Vec3::new(1.0, 0.0, 0.0);
         assert!(
-            (placement.projected_semi_diameter(r, v_u) - r).abs() < tol,
+            (placement.projected_semi_diameter(crm, r, v_u) - r).abs() < tol,
             "U axis: expected {r}, got {}",
-            placement.projected_semi_diameter(r, v_u)
+            placement.projected_semi_diameter(crm, r, v_u)
         );
         assert!(
-            (placement.projected_semi_diameter(r, v_r) - r).abs() < tol,
+            (placement.projected_semi_diameter(crm, r, v_r) - r).abs() < tol,
             "R axis: expected {r}, got {}",
-            placement.projected_semi_diameter(r, v_r)
+            placement.projected_semi_diameter(crm, r, v_r)
         );
     }
 
@@ -1321,22 +1415,22 @@ mod tests {
         // 45° rotation about cursor-R; foreshortens only the U axis.
         let r = 10.0;
         let theta = 45.0_f64.to_radians();
-        let placement = placement_with_rotation(Rotation3D::IntrinsicPassiveRUF(EulerAngles(
-            theta, 0.0, 0.0,
-        )));
+        let (placement, crm) = placement_with_rotation(Rotation3D::IntrinsicPassiveRUF(
+            EulerAngles(theta, 0.0, 0.0),
+        ));
         let tol = 1e-10;
         let v_u = Vec3::new(0.0, 1.0, 0.0);
         let v_r = Vec3::new(1.0, 0.0, 0.0);
         assert!(
-            (placement.projected_semi_diameter(r, v_u) - r * theta.cos()).abs() < tol,
+            (placement.projected_semi_diameter(crm, r, v_u) - r * theta.cos()).abs() < tol,
             "U axis: expected {}, got {}",
             r * theta.cos(),
-            placement.projected_semi_diameter(r, v_u)
+            placement.projected_semi_diameter(crm, r, v_u)
         );
         assert!(
-            (placement.projected_semi_diameter(r, v_r) - r).abs() < tol,
+            (placement.projected_semi_diameter(crm, r, v_r) - r).abs() < tol,
             "R axis: expected {r}, got {}",
-            placement.projected_semi_diameter(r, v_r)
+            placement.projected_semi_diameter(crm, r, v_r)
         );
     }
 
@@ -1345,21 +1439,21 @@ mod tests {
         // 30° rotation about cursor-U; foreshortens only the R axis.
         let r = 10.0;
         let psi = 30.0_f64.to_radians();
-        let placement =
+        let (placement, crm) =
             placement_with_rotation(Rotation3D::IntrinsicPassiveRUF(EulerAngles(0.0, psi, 0.0)));
         let tol = 1e-10;
         let v_u = Vec3::new(0.0, 1.0, 0.0);
         let v_r = Vec3::new(1.0, 0.0, 0.0);
         assert!(
-            (placement.projected_semi_diameter(r, v_r) - r * psi.cos()).abs() < tol,
+            (placement.projected_semi_diameter(crm, r, v_r) - r * psi.cos()).abs() < tol,
             "R axis: expected {}, got {}",
             r * psi.cos(),
-            placement.projected_semi_diameter(r, v_r)
+            placement.projected_semi_diameter(crm, r, v_r)
         );
         assert!(
-            (placement.projected_semi_diameter(r, v_u) - r).abs() < tol,
+            (placement.projected_semi_diameter(crm, r, v_u) - r).abs() < tol,
             "U axis: expected {r}, got {}",
-            placement.projected_semi_diameter(r, v_u)
+            placement.projected_semi_diameter(crm, r, v_u)
         );
     }
 
@@ -1373,6 +1467,7 @@ mod tests {
         let model = mirrors_figure_z::sequential_model(air, &wavelengths);
         let surfaces = model.surfaces();
         let placements = model.placements();
+        let path_steps = model.path_steps(0);
         let r = 12.7_f64;
         let expected_u = r * (30.0_f64.to_radians()).cos();
         let tol = 1e-10;
@@ -1383,15 +1478,16 @@ mod tests {
         for &mirror_idx in &[1usize, 2usize] {
             let sd = surfaces[mirror_idx].mask().semi_diameter();
             let placement = &placements[mirror_idx];
+            let crm = path_steps[mirror_idx].cursor_rotation_matrix;
             assert!(
-                (placement.projected_semi_diameter(sd, v_u) - expected_u).abs() < tol,
+                (placement.projected_semi_diameter(crm, sd, v_u) - expected_u).abs() < tol,
                 "Mirror {mirror_idx} U: expected {expected_u}, got {}",
-                placement.projected_semi_diameter(sd, v_u)
+                placement.projected_semi_diameter(crm, sd, v_u)
             );
             assert!(
-                (placement.projected_semi_diameter(sd, v_r) - r).abs() < tol,
+                (placement.projected_semi_diameter(crm, sd, v_r) - r).abs() < tol,
                 "Mirror {mirror_idx} R: expected {r}, got {}",
-                placement.projected_semi_diameter(sd, v_r)
+                placement.projected_semi_diameter(crm, sd, v_r)
             );
         }
     }
@@ -1412,7 +1508,13 @@ mod tests {
         let model = mirrors_figure_z::sequential_model(n!(1.0), &[0.5876]);
         let v_init = Vec3::new(0.0, 1.0, 0.0); // phi = 90°
 
-        let vecs = propagate_tangential_vec(v_init, model.surfaces(), model.placements());
+        let surface_indices: Vec<usize> = (0..model.surfaces().len()).collect();
+        let vecs = propagate_tangential_vec(
+            v_init,
+            model.surfaces(),
+            model.placements(),
+            &surface_indices,
+        );
 
         let sqrt3_over_2 = (3.0_f64 / 4.0_f64).sqrt();
 
@@ -1429,27 +1531,54 @@ mod tests {
 
     #[test]
     fn is_rotationally_symmetric() {
-        // A system with identity rotations is rotationally symmetric.
-        let id = Mat3x3::identity();
-        let placements = vec![
-            Placement::new(Vec3::new(0.0, 0.0, 0.0), 0.0, id, id, id),
-            Placement::new(Vec3::new(0.0, 0.0, 0.0), 0.0, id, id, id),
+        // A simple on-axis system (Object → Sphere → Image) is symmetric.
+        use crate::{GapSpec, SurfaceSpec, specs::surfaces::BoundaryKind};
+        let air = n!(1.0);
+        let glass = n!(1.5);
+        let gaps = vec![
+            GapSpec {
+                thickness: f64::INFINITY,
+                refractive_index: air.clone(),
+            },
+            GapSpec {
+                thickness: 5.0,
+                refractive_index: glass,
+            },
+            GapSpec {
+                thickness: 50.0,
+                refractive_index: air,
+            },
         ];
-        assert!(SequentialModel::is_rotationally_symmetric(&placements));
+        let surfaces = vec![
+            SurfaceSpec::Object,
+            SurfaceSpec::Sphere {
+                semi_diameter: 12.5,
+                radius_of_curvature: 25.8,
+                surf_kind: BoundaryKind::Refracting,
+                rotation: Rotation3D::None,
+                decenter: Vec3::new(0.0, 0.0, 0.0),
+                rotation_offset: Rotation3D::None,
+            },
+            SurfaceSpec::Image {
+                rotation: Rotation3D::None,
+                decenter: Vec3::new(0.0, 0.0, 0.0),
+                rotation_offset: Rotation3D::None,
+            },
+        ];
+        let simple =
+            SequentialModel::from_surface_specs(&gaps, &surfaces, &[0.5876], None).unwrap();
+        assert!(simple.is_rotationally_symmetric());
 
         // A system with tilted surfaces is not rotationally symmetric.
         use crate::examples::mirrors_figure_z;
-        let air = n!(1.0);
-        let wavelengths = [0.5876];
-        let figure_z = mirrors_figure_z::sequential_model(air, &wavelengths);
-        assert!(!SequentialModel::is_rotationally_symmetric(
-            figure_z.placements()
-        ));
+        let air2 = n!(1.0);
+        let figure_z = mirrors_figure_z::sequential_model(air2, &[0.5876]);
+        assert!(!figure_z.is_rotationally_symmetric());
     }
 
     #[test]
-    fn test_first_physical_surface() {
-        // Object(0), Probe(1), Sphere(2), Sphere(3), Image(4) — first physical is index
+    fn test_first_physical_step() {
+        // Object(0), Probe(1), Sphere(2), Sphere(3), Image(4) — first physical step is
         // 2.
         let surfaces: Vec<Box<dyn Surface>> = vec![
             Box::new(Object::new()),
@@ -1458,14 +1587,15 @@ mod tests {
             Box::new(Sphere::new(1.0, 1.0, BoundaryKind::Refracting)),
             Box::new(Image::new()),
         ];
+        let surface_indices: Vec<usize> = (0..surfaces.len()).collect();
 
-        let result = first_physical_surface(&surfaces);
+        let result = first_physical_step(&surface_indices, &surfaces);
         assert_eq!(result, Some(2));
     }
 
     #[test]
-    fn test_last_physical_surface() {
-        // Object(0), Sphere(1), Sphere(2), Probe(3), Image(4) — last physical is index
+    fn test_last_physical_step() {
+        // Object(0), Sphere(1), Sphere(2), Probe(3), Image(4) — last physical step is
         // 2.
         let surfaces: Vec<Box<dyn Surface>> = vec![
             Box::new(Object::new()),
@@ -1474,8 +1604,9 @@ mod tests {
             Box::new(Probe::new()),
             Box::new(Image::new()),
         ];
+        let surface_indices: Vec<usize> = (0..surfaces.len()).collect();
 
-        let result = last_physical_surface(&surfaces);
+        let result = last_physical_step(&surface_indices, &surfaces);
         assert_eq!(result, Some(2));
     }
 
@@ -1494,19 +1625,19 @@ mod tests {
         let id = Mat3x3::identity();
 
         // z-coordinate infinite
-        let p = Placement::new(Vec3::new(0.0, 0.0, Float::INFINITY), 0.0, id, id, id);
+        let p = SurfacePlacement::new(Vec3::new(0.0, 0.0, Float::INFINITY), 0.0, id, id);
         assert!(p.is_infinite());
 
         // y-coordinate infinite
-        let p = Placement::new(Vec3::new(0.0, Float::INFINITY, 0.0), 0.0, id, id, id);
+        let p = SurfacePlacement::new(Vec3::new(0.0, Float::INFINITY, 0.0), 0.0, id, id);
         assert!(p.is_infinite());
 
         // x-coordinate infinite
-        let p = Placement::new(Vec3::new(Float::INFINITY, 0.0, 0.0), 0.0, id, id, id);
+        let p = SurfacePlacement::new(Vec3::new(Float::INFINITY, 0.0, 0.0), 0.0, id, id);
         assert!(p.is_infinite());
 
         // finite
-        let p = Placement::new(Vec3::new(0.0, 0.0, 0.0), 0.0, id, id, id);
+        let p = SurfacePlacement::new(Vec3::new(0.0, 0.0, 0.0), 0.0, id, id);
         assert!(!p.is_infinite());
     }
 
@@ -1531,15 +1662,14 @@ mod tests {
 
     /// For a straight system, axis_direction should equal (0, 0, 1) everywhere.
     #[test]
-    fn placement_axis_direction_straight_system() {
+    fn cursor_placement_axis_direction_straight_system() {
         use crate::examples::convexplano_lens;
         use approx::assert_abs_diff_eq;
         let model = convexplano_lens::sequential_model(n!(1.0), n!(1.515), &[0.5876]);
-        for placement in model.placements() {
-            let axis = placement.axis_direction();
-            assert_abs_diff_eq!(axis.x(), 0.0, epsilon = 1e-12);
-            assert_abs_diff_eq!(axis.y(), 0.0, epsilon = 1e-12);
-            assert_abs_diff_eq!(axis.z(), 1.0, epsilon = 1e-12);
+        for cp in model.path_steps(0) {
+            assert_abs_diff_eq!(cp.axis_direction.x(), 0.0, epsilon = 1e-12);
+            assert_abs_diff_eq!(cp.axis_direction.y(), 0.0, epsilon = 1e-12);
+            assert_abs_diff_eq!(cp.axis_direction.z(), 1.0, epsilon = 1e-12);
         }
     }
 
@@ -1733,11 +1863,18 @@ mod tests {
         ];
         let id = Mat3x3::identity();
         let placements = vec![
-            Placement::new(Vec3::new(0., 0., 0.), 0., id, id, id),
-            Placement::new(Vec3::new(0., 0., 5.), 5., id, id, id),
-            Placement::new(Vec3::new(0., 0., 10.), 10., id, id, id),
-            Placement::new(Vec3::new(0., 0., 15.), 15., id, id, id),
+            SurfacePlacement::new(Vec3::new(0., 0., 0.), 0., id, id),
+            SurfacePlacement::new(Vec3::new(0., 0., 5.), 5., id, id),
+            SurfacePlacement::new(Vec3::new(0., 0., 10.), 10., id, id),
+            SurfacePlacement::new(Vec3::new(0., 0., 15.), 15., id, id),
         ];
+        let path_steps: Vec<CursorPlacement> = (0..4)
+            .map(|_| CursorPlacement {
+                axis_direction: Vec3::new(0.0, 0.0, 1.0),
+                cursor_position: Vec3::new(0.0, 0.0, 0.0),
+                cursor_rotation_matrix: id,
+            })
+            .collect();
         let gaps = vec![
             Gap {
                 thickness: 5.0,
@@ -1752,8 +1889,18 @@ mod tests {
                 refractive_index: RefractiveIndex::try_from_spec(n!(1.0).as_ref(), 0.5876).unwrap(),
             },
         ];
-        let submodel = SequentialSubModelBase::new(vec![0, 2, 1, 3], gaps, vec![None; 4]);
-        let mut iter = submodel.try_iter(&surfaces, &placements).unwrap();
+        let surface_indices = vec![0, 2, 1, 3];
+        let bs_arms = vec![None; 4];
+        let submodel = SequentialSubModelBase::new(gaps);
+        let mut iter = submodel
+            .try_iter(
+                &surfaces,
+                &placements,
+                &surface_indices,
+                &bs_arms,
+                &path_steps,
+            )
+            .unwrap();
 
         let step0 = iter.next().unwrap();
         // step0 surface should be Sphere at index 2 (sd = 20.0)
@@ -1775,10 +1922,17 @@ mod tests {
         ];
         let id = Mat3x3::identity();
         let placements = vec![
-            Placement::new(Vec3::new(0., 0., 0.), 0., id, id, id),
-            Placement::new(Vec3::new(0., 0., 5.), 5., id, id, id),
-            Placement::new(Vec3::new(0., 0., 10.), 10., id, id, id),
+            SurfacePlacement::new(Vec3::new(0., 0., 0.), 0., id, id),
+            SurfacePlacement::new(Vec3::new(0., 0., 5.), 5., id, id),
+            SurfacePlacement::new(Vec3::new(0., 0., 10.), 10., id, id),
         ];
+        let path_steps: Vec<CursorPlacement> = (0..4)
+            .map(|_| CursorPlacement {
+                axis_direction: Vec3::new(0.0, 0.0, 1.0),
+                cursor_position: Vec3::new(0.0, 0.0, 0.0),
+                cursor_rotation_matrix: id,
+            })
+            .collect();
         let gaps = vec![
             Gap {
                 thickness: 5.0,
@@ -1793,9 +1947,20 @@ mod tests {
                 refractive_index: RefractiveIndex::try_from_spec(n!(1.0).as_ref(), 0.5876).unwrap(),
             },
         ];
-        let submodel = SequentialSubModelBase::new(vec![0, 1, 1, 2], gaps, vec![None; 4]);
+        let surface_indices = vec![0, 1, 1, 2];
+        let bs_arms = vec![None; 4];
+        let submodel = SequentialSubModelBase::new(gaps);
         // Must not panic or error — iteration visits index 1 twice.
-        let count = submodel.try_iter(&surfaces, &placements).unwrap().count();
+        let count = submodel
+            .try_iter(
+                &surfaces,
+                &placements,
+                &surface_indices,
+                &bs_arms,
+                &path_steps,
+            )
+            .unwrap()
+            .count();
         assert_eq!(count, 3); // 3 gaps → 3 steps
     }
 
@@ -1813,11 +1978,18 @@ mod tests {
         ];
         let id = Mat3x3::identity();
         let placements = vec![
-            Placement::new(Vec3::new(0., 0., 0.), 0., id, id, id),
-            Placement::new(Vec3::new(0., 0., 10.), 10., id, id, id),
-            Placement::new(Vec3::new(0., 0., 20.), 20., id, id, id),
-            Placement::new(Vec3::new(0., 0., 30.), 30., id, id, id),
+            SurfacePlacement::new(Vec3::new(0., 0., 0.), 0., id, id),
+            SurfacePlacement::new(Vec3::new(0., 0., 10.), 10., id, id),
+            SurfacePlacement::new(Vec3::new(0., 0., 20.), 20., id, id),
+            SurfacePlacement::new(Vec3::new(0., 0., 30.), 30., id, id),
         ];
+        let path_steps: Vec<CursorPlacement> = (0..4)
+            .map(|_| CursorPlacement {
+                axis_direction: Vec3::new(0.0, 0.0, 1.0),
+                cursor_position: Vec3::new(0.0, 0.0, 0.0),
+                cursor_rotation_matrix: id,
+            })
+            .collect();
         let gaps = vec![
             Gap {
                 thickness: 10.0,
@@ -1832,9 +2004,18 @@ mod tests {
                 refractive_index: RefractiveIndex::try_from_spec(n!(1.0).as_ref(), 0.5876).unwrap(),
             },
         ];
+        let surface_indices = vec![0, 1, 2, 3];
         let bs_arms = vec![None, None, Some(BeamSplitterPathKind::Transmitting), None];
-        let submodel = SequentialSubModelBase::new(vec![0, 1, 2, 3], gaps, bs_arms);
-        let mut iter = submodel.try_iter(&surfaces, &placements).unwrap();
+        let submodel = SequentialSubModelBase::new(gaps);
+        let mut iter = submodel
+            .try_iter(
+                &surfaces,
+                &placements,
+                &surface_indices,
+                &bs_arms,
+                &path_steps,
+            )
+            .unwrap();
 
         let step0 = iter.next().unwrap(); // Sphere
         assert_eq!(step0.bs_arm, None);
@@ -1953,4 +2134,78 @@ mod tests {
             img_with.z()
         );
     }
+
+    // AT: path 0 cursor positions lie along +Z; path 1 cursor positions
+    // diverge to +Y after the beam splitter.
+    #[test]
+    fn path_steps_cursor_positions_two_path_model() {
+        use crate::examples::beam_splitter;
+        use crate::specs::gaps::ConstantRefractiveIndex;
+        use std::rc::Rc;
+
+        let n_air: Rc<dyn crate::RefractiveIndexSpec> =
+            Rc::new(ConstantRefractiveIndex::new(1.0, 0.0));
+        let model = beam_splitter::two_path_model(n_air, &[0.5876], 10.0, 10.0);
+
+        let tol = 1e-10;
+
+        // Path 0: Object(z=−inf), BS(z=0), Image_T(z=10)
+        let steps0 = model.path_steps(0);
+        assert_eq!(steps0.len(), 3);
+        // Object cursor is at −∞ along Z.
+        assert!(steps0[0].cursor_position.z().is_infinite());
+        // BS cursor at z=0.
+        assert!((steps0[1].cursor_position.z()).abs() < tol);
+        // Image_T cursor at z=10.
+        assert!((steps0[2].cursor_position.z() - 10.0).abs() < tol);
+        // All path-0 cursor positions have zero x and y.
+        for s in steps0.iter().filter(|s| s.cursor_position.z().is_finite()) {
+            assert!(s.cursor_position.x().abs() < tol);
+            assert!(s.cursor_position.y().abs() < tol);
+        }
+
+        // Path 1: Object(shared), BS(shared, reflected to −Y), Image_R
+        // The BS uses a −45° rotation about R, so the reflected arm travels in −Y.
+        let steps1 = model.path_steps(1);
+        assert_eq!(steps1.len(), 3);
+        // Image_R cursor position diverges to −Y: x=0, y=−10, z=0.
+        assert!(
+            steps1[2].cursor_position.x().abs() < tol,
+            "x={}",
+            steps1[2].cursor_position.x()
+        );
+        assert!(
+            (steps1[2].cursor_position.y() + 10.0).abs() < tol,
+            "y={}",
+            steps1[2].cursor_position.y()
+        );
+        assert!(
+            steps1[2].cursor_position.z().abs() < tol,
+            "z={}",
+            steps1[2].cursor_position.z()
+        );
+    }
+
+    // AT: path 1's axis direction at the Image_R step points along −Y after a
+    // −45° BS rotation about R.
+    #[test]
+    fn path_steps_axis_direction_reflected_arm() {
+        use crate::examples::beam_splitter;
+        use crate::specs::gaps::ConstantRefractiveIndex;
+        use std::rc::Rc;
+
+        let n_air: Rc<dyn crate::RefractiveIndexSpec> =
+            Rc::new(ConstantRefractiveIndex::new(1.0, 0.0));
+        let model = beam_splitter::two_path_model(n_air, &[0.5876], 10.0, 10.0);
+
+        let tol = 1e-10;
+        let steps1 = model.path_steps(1);
+
+        // The BS example uses a −45° rotation about R, so the reflected arm
+        // travels in the −Y direction.
+        let d = steps1[2].axis_direction;
+        assert!(d.x().abs() < tol, "x={}", d.x());
+        assert!((d.y() + 1.0).abs() < tol, "y={}", d.y());
+        assert!(d.z().abs() < tol, "z={}", d.z());
+    }
 }
diff --git a/crates/cherry-rs/src/core/sequential_model/placement.rs b/crates/cherry-rs/src/core/sequential_model/surface_placement.rs
similarity index 78%
rename from crates/cherry-rs/src/core/sequential_model/placement.rs
rename to crates/cherry-rs/src/core/sequential_model/surface_placement.rs
index 867a9038..675641f2 100644
--- a/crates/cherry-rs/src/core/sequential_model/placement.rs
+++ b/crates/cherry-rs/src/core/sequential_model/surface_placement.rs
@@ -1,11 +1,17 @@
 /// Placement of a surface in a sequential optical system.
 ///
-/// A [`Placement`] describes *where* a surface sits in 3D space and how the
-/// optical axis (cursor) is oriented when it arrives at that surface. It is
-/// intentionally separate from surface geometry ([`Surface`]) so that
-/// coordinate-system operations and intrinsic-geometry operations do not mix.
+/// A [`SurfacePlacement`] describes *where* a surface sits in 3D space and
+/// how its rotation relates to the global frame. It is intentionally separate
+/// from surface geometry ([`Surface`]) so that coordinate-system operations
+/// and intrinsic-geometry operations do not mix.
+///
+/// The cursor-frame state (axis direction, cursor position, cursor rotation
+/// matrix) is path-specific and is carried in [`CursorPlacement`] on each
+/// iterator [`Step`] rather than stored here.
 ///
 /// [`Surface`]: crate::core::surfaces::Surface
+/// [`CursorPlacement`]: crate::core::sequential_model::CursorPlacement
+/// [`Step`]: crate::core::sequential_model::Step
 use crate::core::{
     Float,
     math::{linalg::mat3x3::Mat3x3, vec3::Vec3},
@@ -15,7 +21,7 @@ use super::cursor::Cursor;
 
 /// Position and orientation of a surface in the global coordinate system.
 #[derive(Debug, Clone)]
-pub struct Placement {
+pub struct SurfacePlacement {
     /// Vertex position in the global coordinate system.
     pub position: Vec3,
 
@@ -39,18 +45,12 @@ pub struct Placement {
     /// `rotation_offset` is `None`. Used by `propagate_tangential_vec` so
     /// that the paraxial tangential axis follows the nominal system.
     pub nominal_inv_rotation_matrix: Mat3x3,
-
-    /// Rotation from the global frame into the optical-axis (cursor) frame
-    /// only, without any surface tilt applied.
-    ///
-    /// Needed for aperture-projection calculations (`projected_semi_diameter`)
-    /// and for `is_rotationally_symmetric`.
-    pub cursor_rotation_matrix: Mat3x3,
 }
 
-impl Placement {
-    /// Build a [`Placement`] from a pre-composed surface rotation matrix, the
-    /// nominal rotation matrix, a decenter, and the current cursor state.
+impl SurfacePlacement {
+    /// Build a [`SurfacePlacement`] from a pre-composed surface rotation
+    /// matrix, the nominal rotation matrix, a decenter, and the current cursor
+    /// state.
     ///
     /// `actual_rotation_matrix` is
     /// `rotation_offset.rotation_matrix() * rotation.rotation_matrix()`.
@@ -69,22 +69,15 @@ impl Placement {
         let nom_rot_matrix = nominal_rotation_matrix * cursor_rotation_matrix;
         let offset_global = cursor_rotation_matrix.transpose() * decenter;
         let position = cursor.pos() + offset_global;
-        Self::new(
-            position,
-            cursor.track(),
-            rotation_matrix,
-            nom_rot_matrix,
-            cursor_rotation_matrix,
-        )
+        Self::new(position, cursor.track(), rotation_matrix, nom_rot_matrix)
     }
 
-    /// Create a new [`Placement`] from its constituent parts.
+    /// Create a new [`SurfacePlacement`] from its constituent parts.
     pub fn new(
         position: Vec3,
         track: Float,
         rotation_matrix: Mat3x3,
         nominal_rotation_matrix: Mat3x3,
-        cursor_rotation_matrix: Mat3x3,
     ) -> Self {
         let inv_rotation_matrix = rotation_matrix.transpose();
         let nominal_inv_rotation_matrix = nominal_rotation_matrix.transpose();
@@ -94,7 +87,6 @@ impl Placement {
             rotation_matrix,
             inv_rotation_matrix,
             nominal_inv_rotation_matrix,
-            cursor_rotation_matrix,
         }
     }
 
@@ -113,40 +105,32 @@ impl Placement {
             || self.position.z().is_infinite()
     }
 
-    /// Returns the unit vector pointing along the optical axis (cursor forward
-    /// direction) at this surface, expressed in the global frame.
-    pub fn axis_direction(&self) -> Vec3 {
-        // The third row of cursor_rotation_matrix is the forward direction in
-        // global coords when the matrix is global-to-cursor.
-        // Transposing maps it back to global, giving the forward vector.
-        self.cursor_rotation_matrix.transpose() * Vec3::new(0.0, 0.0, 1.0)
-    }
-
     /// Returns the semi-diameter as seen by a paraxial ray travelling along
     /// the cursor axis in the tangential plane defined by `v`.
     ///
-    /// `r` is the surface's clear-aperture semi-diameter (from
-    /// [`Surface::semi_diameter`]). `v` is a unit vector in the global
-    /// frame that lies in the transverse plane and defines the meridional
-    /// plane of interest.
-    ///
-    /// For a tilted surface the effective limit on cursor height is
-    /// `r · |n_F| / sqrt(n_φ² + n_F²)` where `(n_R, n_U, n_F)` are the
-    /// surface-normal components in the cursor frame and
-    /// `n_φ = n_R · v_x + n_U · v_y` is the component along `v`.
+    /// `cursor_rotation_matrix` is the path-specific cursor orientation at
+    /// this step (from [`CursorPlacement::cursor_rotation_matrix`]).
+    /// `r` is the surface's clear-aperture semi-diameter. `v` is a unit
+    /// vector in the global frame that lies in the transverse plane and
+    /// defines the meridional plane of interest.
     ///
     /// Returns [`Float::INFINITY`] when `r` is infinite (non-aperture
     /// surfaces).
     ///
-    /// [`Surface::semi_diameter`]: crate::core::surfaces::Surface::semi_diameter
-    pub fn projected_semi_diameter(&self, r: Float, v: Vec3) -> Float {
+    /// [`CursorPlacement::cursor_rotation_matrix`]: crate::core::sequential_model::CursorPlacement::cursor_rotation_matrix
+    pub fn projected_semi_diameter(
+        &self,
+        cursor_rotation_matrix: Mat3x3,
+        r: Float,
+        v: Vec3,
+    ) -> Float {
         if r.is_infinite() {
             return Float::INFINITY;
         }
 
         // R_surf = cursor_to_local = global_to_local · cursor_to_global
         //        = rotation_matrix · cursor_rotation_matrix.transpose()
-        let r_surf = self.rotation_matrix * self.cursor_rotation_matrix.transpose();
+        let r_surf = self.rotation_matrix * cursor_rotation_matrix.transpose();
 
         // Third row of R_surf is the surface normal expressed in cursor frame:
         // (n_R, n_U, n_F)
@@ -185,7 +169,8 @@ mod tests {
     fn at1_zero_displacement_identity_cursor() {
         let id = Mat3x3::identity();
         let cursor = identity_cursor();
-        let p = Placement::from_decenter_and_rotation(Vec3::new(0.0, 0.0, 0.0), id, id, &cursor);
+        let p =
+            SurfacePlacement::from_decenter_and_rotation(Vec3::new(0.0, 0.0, 0.0), id, id, &cursor);
         assert!(
             p.position.approx_eq(&Vec3::new(0.0, 0.0, 0.0), 1e-15),
             "position: {:?}",
@@ -193,7 +178,6 @@ mod tests {
         );
         assert!(p.rotation_matrix.approx_eq(&id, 1e-15));
         assert!(p.nominal_inv_rotation_matrix.approx_eq(&id, 1e-15));
-        assert!(p.cursor_rotation_matrix.approx_eq(&id, 1e-15));
         assert!((p.track - 0.0).abs() < 1e-15);
     }
 
@@ -203,7 +187,8 @@ mod tests {
         let id = Mat3x3::identity();
         let dx = 2.5;
         let cursor = cursor_at(10.0);
-        let p = Placement::from_decenter_and_rotation(Vec3::new(dx, 0.0, 0.0), id, id, &cursor);
+        let p =
+            SurfacePlacement::from_decenter_and_rotation(Vec3::new(dx, 0.0, 0.0), id, id, &cursor);
         let tol = 1e-15;
         assert!((p.position.x() - dx).abs() < tol, "x: {}", p.position.x());
         assert!(p.position.y().abs() < tol, "y: {}", p.position.y());
@@ -216,7 +201,8 @@ mod tests {
         let id = Mat3x3::identity();
         let dy = -1.3;
         let cursor = cursor_at(5.0);
-        let p = Placement::from_decenter_and_rotation(Vec3::new(0.0, dy, 0.0), id, id, &cursor);
+        let p =
+            SurfacePlacement::from_decenter_and_rotation(Vec3::new(0.0, dy, 0.0), id, id, &cursor);
         let tol = 1e-15;
         assert!(p.position.x().abs() < tol, "x: {}", p.position.x());
         assert!((p.position.y() - dy).abs() < tol, "y: {}", p.position.y());
@@ -230,7 +216,8 @@ mod tests {
         let dz = 3.0;
         let nominal_z = 7.0;
         let cursor = cursor_at(nominal_z);
-        let p = Placement::from_decenter_and_rotation(Vec3::new(0.0, 0.0, dz), id, id, &cursor);
+        let p =
+            SurfacePlacement::from_decenter_and_rotation(Vec3::new(0.0, 0.0, dz), id, id, &cursor);
         let tol = 1e-15;
         assert!(p.position.x().abs() < tol, "x: {}", p.position.x());
         assert!(p.position.y().abs() < tol, "y: {}", p.position.y());
@@ -244,7 +231,8 @@ mod tests {
     // AT-5: F-axis decenter after a 90° fold shifts along the post-fold axis,
     // not along global Z. After a 90° fold about the R-axis the cursor
     // forward direction is (0, -1, 0), so a cursor-F decenter of dz moves the
-    // vertex by (0, -dz, 0) in global coordinates.
+    // vertex by (0, -dz, 0) in global coords relative to the mirror without
+    // decenter.
     #[test]
     fn at5_axial_decenter_after_fold() {
         // Build a 3-surface system: Object — flat mirror (45° theta fold) — Image
diff --git a/crates/cherry-rs/src/core/surfaces/kinds/mod.rs b/crates/cherry-rs/src/core/surfaces/kinds/mod.rs
index 4afa12f9..42bdf98a 100644
--- a/crates/cherry-rs/src/core/surfaces/kinds/mod.rs
+++ b/crates/cherry-rs/src/core/surfaces/kinds/mod.rs
@@ -5,6 +5,7 @@ pub mod iris;
 pub mod object;
 pub mod probe;
 pub mod sphere;
+pub mod thin_lens;
 
 pub use beam_splitter::BeamSplitter;
 pub use conic::Conic;
@@ -13,3 +14,4 @@ pub use iris::Iris;
 pub use object::Object;
 pub use probe::Probe;
 pub use sphere::Sphere;
+pub use thin_lens::ThinLens;
diff --git a/crates/cherry-rs/src/core/surfaces/kinds/thin_lens.rs b/crates/cherry-rs/src/core/surfaces/kinds/thin_lens.rs
new file mode 100644
index 00000000..1ac06e3a
--- /dev/null
+++ b/crates/cherry-rs/src/core/surfaces/kinds/thin_lens.rs
@@ -0,0 +1,127 @@
+/// A thin lens surface type for 3D ray tracing.
+use crate::core::{
+    Float,
+    math::vec3::Vec3,
+    surfaces::{BoundaryKind, Mask, Surface, SurfaceKind, solvers::flat_surface},
+};
+
+/// A thin lens for 3D ray tracing.
+///
+/// The focal length of the thin lens is specified as that when the refractive
+/// index of both sides of the lens is in air (1.0).
+#[derive(Debug)]
+pub struct ThinLens {
+    pub focal_length: Float,
+    mask: Mask,
+}
+
+impl ThinLens {
+    pub fn new(semi_diameter: Float, focal_length: Float) -> Self {
+        Self {
+            focal_length,
+            mask: Mask::Circular { semi_diameter },
+        }
+    }
+}
+
+impl Surface for ThinLens {
+    fn boundary_kind(&self) -> BoundaryKind {
+        BoundaryKind::Refracting
+    }
+
+    fn interact(&self, ray: &mut crate::Ray, n_0: Float, n_1: Float, _norm: crate::Vec3) {
+        let power = self.power(0.0, n_0, n_1);
+        let ux_0 = ray.dir().l() / ray.dir().n();
+        let uy_0 = ray.dir().m() / ray.dir().n();
+
+        // Paraxial ray transer equations for a thin lens
+        let ux_1 = (n_0 * ux_0 - power * ray.pos().x()) / n_1;
+        let uy_1 = (n_0 * uy_0 - power * ray.pos().y()) / n_1;
+
+        // l_1 = ux_1 * dir_cos_n_1, m_1 = uy_1 * dir_cos_n_1, and
+        // l_1^2 + m_1^2 + dir_cos_n_1^2 = 1, so:
+        // dir_cos_n_1^2 * (ux_1^2 + uy_1^2 + 1) = 1
+        let dir_cos_n_1 = (1.0 + ux_1 * ux_1 + uy_1 * uy_1).powf(-0.5);
+        let dir_cos_l_1 = ux_1 * dir_cos_n_1;
+        let dir_cos_m_1 = uy_1 * dir_cos_n_1;
+
+        ray.set_dir(Vec3::new(dir_cos_l_1, dir_cos_m_1, dir_cos_n_1).normalize());
+    }
+
+    fn intersect(
+        &self,
+        ray: &crate::core::ray::Ray,
+        max_iter: usize,
+    ) -> anyhow::Result<(Vec3, Vec3)> {
+        flat_surface(ray, self, max_iter)
+    }
+
+    fn mask(&self) -> &Mask {
+        &self.mask
+    }
+
+    fn power(&self, _azimuth_rad: Float, _n_0: Float, _n_1: Float) -> Float {
+        1.0 / self.focal_length
+    }
+
+    fn norm(&self, _pos: Vec3) -> Vec3 {
+        Vec3::new(0.0, 0.0, 1.0)
+    }
+
+    fn sag(&self, _pos: crate::Vec3) -> Float {
+        0.0
+    }
+
+    fn surface_kind(&self) -> SurfaceKind {
+        SurfaceKind::ThinLens
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::core::ray::Ray;
+    use approx::assert_abs_diff_eq;
+
+    /// Tests a ray incident at 10 degrees at a height of 10 mm on a thin lens
+    /// with a focal length of 100 mm.
+    #[test]
+    fn test_ray_interact() {
+        let focal_length = 100.0;
+        let semi_diameter = 25.0;
+
+        let lens = ThinLens::new(semi_diameter, focal_length);
+        let mut ray = Ray::new(
+            Vec3::new(0.0, 10.0, 0.0),
+            Vec3::new(
+                0.0,
+                (80.0_f64).to_radians().cos(),
+                (10.0_f64).to_radians().cos(),
+            )
+            .normalize(),
+        );
+
+        let norm = lens.norm(ray.pos());
+        assert_abs_diff_eq!(norm.x(), 0.0, epsilon = 1e-6);
+        assert_abs_diff_eq!(norm.y(), 0.0, epsilon = 1e-6);
+        assert_abs_diff_eq!(norm.z(), 1.0, epsilon = 1e-6);
+
+        lens.interact(&mut ray, 1.0, 1.0, norm);
+        let dir = ray.dir();
+        assert_abs_diff_eq!(dir.x(), 0.0, epsilon = 1e-6);
+        assert_abs_diff_eq!(dir.y(), 0.07610561, epsilon = 1e-6);
+        assert_abs_diff_eq!(dir.z(), 0.99709976, epsilon = 1e-6);
+    }
+
+    /// A thin lens's power is `1/focal_length` regardless of the surrounding
+    /// media — unlike a curved refracting surface, it isn't derived from
+    /// `roc()` and `n_0`/`n_1`.
+    #[test]
+    fn test_power_is_media_independent() {
+        let lens = ThinLens::new(25.0, 100.0);
+        let expected = 1.0 / 100.0;
+        assert_abs_diff_eq!(lens.power(0.0, 1.0, 1.0), expected, epsilon = 1e-12);
+        assert_abs_diff_eq!(lens.power(0.0, 1.0, 1.5), expected, epsilon = 1e-12);
+        assert_abs_diff_eq!(lens.power(1.23, 1.0, 1.0), expected, epsilon = 1e-12);
+    }
+}
diff --git a/crates/cherry-rs/src/core/surfaces/mod.rs b/crates/cherry-rs/src/core/surfaces/mod.rs
index ed6d9bd1..c07f84cf 100644
--- a/crates/cherry-rs/src/core/surfaces/mod.rs
+++ b/crates/cherry-rs/src/core/surfaces/mod.rs
@@ -10,7 +10,7 @@ pub mod solvers;
 #[cfg(feature = "serde")]
 pub mod surface_registry;
 
-pub use kinds::{BeamSplitter, Conic, Image, Iris, Object, Probe, Sphere};
+pub use kinds::{BeamSplitter, Conic, Image, Iris, Object, Probe, Sphere, ThinLens};
 #[cfg(feature = "serde")]
 pub use surface_registry::{SurfaceConstructor, SurfaceRegistry};
 
@@ -31,6 +31,7 @@ pub enum SurfaceKind {
     Object,
     Probe,
     Sphere,
+    ThinLens,
     Custom,
 }
 
@@ -83,6 +84,19 @@ pub trait Surface: std::fmt::Debug + Send + Sync {
         Float::INFINITY
     }
 
+    /// Returns the surface's contribution to paraxial refracting power.
+    ///
+    /// The default implementation is the standard single-interface power
+    /// formula, `(n_1 - n_0) / roc(azimuth_rad)`, used by [`ParaxialView`] to
+    /// build ray transfer matrices. Surfaces whose power is not determined by
+    /// curvature and surrounding media (e.g. an idealized thin lens) override
+    /// this directly instead of trying to encode their power into `roc()`.
+    ///
+    /// [`ParaxialView`]: crate::views::paraxial::ParaxialView
+    fn power(&self, azimuth_rad: Float, n_0: Float, n_1: Float) -> Float {
+        (n_1 - n_0) / self.roc(azimuth_rad)
+    }
+
     /// Returns the surface sag at a given position in local coordinates.
     fn sag(&self, pos: Vec3) -> Float;
 
diff --git a/crates/cherry-rs/src/examples/mod.rs b/crates/cherry-rs/src/examples/mod.rs
index b88ee675..63db1fc8 100644
--- a/crates/cherry-rs/src/examples/mod.rs
+++ b/crates/cherry-rs/src/examples/mod.rs
@@ -7,3 +7,4 @@ pub mod f_theta_scan_lens;
 pub mod galvo_mirror;
 pub mod mirrors_figure_z;
 pub mod petzval_lens;
+pub mod thin_lens_singlet;
diff --git a/crates/cherry-rs/src/examples/thin_lens_singlet.rs b/crates/cherry-rs/src/examples/thin_lens_singlet.rs
new file mode 100644
index 00000000..09bec35d
--- /dev/null
+++ b/crates/cherry-rs/src/examples/thin_lens_singlet.rs
@@ -0,0 +1,36 @@
+//! A single idealized thin lens with f = 100 mm, object at infinity.
+use std::rc::Rc;
+
+use crate::{GapSpec, RefractiveIndexSpec, Rotation3D, SequentialModel, SurfaceSpec, Vec3};
+
+pub fn sequential_model(
+    n_air: Rc<dyn RefractiveIndexSpec>,
+    wavelengths: &[f64],
+) -> SequentialModel {
+    let gap_0 = GapSpec {
+        thickness: f64::INFINITY,
+        refractive_index: n_air.clone(),
+    };
+    let gap_1 = GapSpec {
+        thickness: 100.0,
+        refractive_index: n_air,
+    };
+    let gaps = vec![gap_0, gap_1];
+
+    let surf_0 = SurfaceSpec::Object;
+    let surf_1 = SurfaceSpec::ThinLens {
+        semi_diameter: 12.5,
+        focal_length: 100.0,
+        rotation: Rotation3D::None,
+        decenter: Vec3::new(0.0, 0.0, 0.0),
+        rotation_offset: Rotation3D::None,
+    };
+    let surf_2 = SurfaceSpec::Image {
+        rotation: Rotation3D::None,
+        decenter: Vec3::new(0.0, 0.0, 0.0),
+        rotation_offset: Rotation3D::None,
+    };
+    let surfaces = vec![surf_0, surf_1, surf_2];
+
+    SequentialModel::from_surface_specs(&gaps, &surfaces, wavelengths, None).unwrap()
+}
diff --git a/crates/cherry-rs/src/gui/app.rs b/crates/cherry-rs/src/gui/app.rs
index c4d372b0..0e2ab284 100644
--- a/crates/cherry-rs/src/gui/app.rs
+++ b/crates/cherry-rs/src/gui/app.rs
@@ -474,6 +474,10 @@ impl eframe::App for CherryApp {
 
                 ui.menu_button("Examples", |ui| {
                     ui.label("Simple");
+                    if ui.button("Thin Lens").clicked() {
+                        self.load_specs(examples::thin_lens());
+                        ui.close();
+                    }
                     if ui.button("Convexplano Lens").clicked() {
                         self.load_specs(SystemSpecs::default());
                         ui.close();
diff --git a/crates/cherry-rs/src/gui/compute.rs b/crates/cherry-rs/src/gui/compute.rs
index 926d09d3..e6b6f4d8 100644
--- a/crates/cherry-rs/src/gui/compute.rs
+++ b/crates/cherry-rs/src/gui/compute.rs
@@ -272,6 +272,7 @@ fn build_surface_descs(seq: &SequentialModel) -> Vec<SurfaceDesc> {
                 SurfaceKind::Probe => "Probe",
                 SurfaceKind::Iris => "Iris",
                 SurfaceKind::Sphere => "Sphere",
+                SurfaceKind::ThinLens => "Thin Lens",
                 SurfaceKind::Custom => "Custom",
             };
             SurfaceDesc {
@@ -279,7 +280,6 @@ fn build_surface_descs(seq: &SequentialModel) -> Vec<SurfaceDesc> {
                 label: format!("{name} [{i}]"),
                 pos: p.position,
                 rot_mat: p.rotation_matrix,
-                cursor_rot_mat: p.cursor_rotation_matrix,
             }
         })
         .collect()
diff --git a/crates/cherry-rs/src/gui/convert.rs b/crates/cherry-rs/src/gui/convert.rs
index 53b41bc3..e42bc9e1 100644
--- a/crates/cherry-rs/src/gui/convert.rs
+++ b/crates/cherry-rs/src/gui/convert.rs
@@ -143,6 +143,19 @@ fn convert_specs_inner(
                     rotation_offset: Rotation3D::None,
                 }
             }
+            SurfaceVariant::ThinLens => {
+                let semi_diameter = parse_float(&row.semi_diameter)
+                    .with_context(|| format!("surface {i}: semi-diameter"))?;
+                let focal_length = parse_float(&row.focal_length)
+                    .with_context(|| format!("surface {i}: focal length"))?;
+                SurfaceSpec::ThinLens {
+                    semi_diameter,
+                    focal_length,
+                    rotation: Rotation3D::None,
+                    decenter: Vec3::new(0.0, 0.0, 0.0),
+                    rotation_offset: Rotation3D::None,
+                }
+            }
             SurfaceVariant::Iris => {
                 let semi_diameter = parse_float(&row.semi_diameter)
                     .with_context(|| format!("surface {i}: semi-diameter"))?;
@@ -313,6 +326,7 @@ fn apply_group_transforms(
                     Component::Element { surf_idxs } => all_surfs.extend(surf_idxs),
                     Component::Iris { stop_idx } => all_surfs.push(*stop_idx),
                     Component::Mirror { surf_idx } => all_surfs.push(*surf_idx),
+                    Component::ThinLens { surf_idx } => all_surfs.push(*surf_idx),
                     Component::UnpairedSurface { surf_idx } => all_surfs.push(*surf_idx),
                 }
             }
@@ -327,9 +341,10 @@ fn apply_group_transforms(
         // The first surface in the group is the pivot / coordinate-frame origin.
         let s1 = *all_surfs.first().unwrap();
         let p = placements[s1].position; // pivot vertex, global frame
-        let c_s1 = placements[s1].cursor_rotation_matrix; // passive global→cursor at s1
+        let c_s1 = placements[s1].rotation_matrix; // passive global→surface local at s1
 
-        // Convert group rotation (degrees, cursor frame at s1) to passive matrix.
+        // Convert group rotation (degrees, surface-local frame at s1) to passive
+        // matrix.
         let [theta_deg, psi_deg, phi_deg] = group.rotation;
         let r_cursor_passive = Rotation3D::IntrinsicPassiveRUF(EulerAngles(
             theta_deg.to_radians(),
@@ -339,11 +354,12 @@ fn apply_group_transforms(
         .rotation_matrix();
 
         // R_group: passive rotation in global frame.
-        // R_group = C_{s1}^T · R_cursor_passive · C_{s1}
+        // R_group = C_{s1}^T · R_local_passive · C_{s1}
         let c_s1_t = c_s1.transpose();
         let r_group = c_s1_t * r_cursor_passive * c_s1;
 
-        // d_global: group decenter converted from cursor frame of s1 to global frame.
+        // d_global: group decenter converted from surface-local frame of s1 to global
+        // frame.
         let [dr, du, df] = group.decenter;
         let d_user = Vec3::new(dr, du, df);
         let d_global = c_s1_t * d_user;
@@ -353,16 +369,16 @@ fn apply_group_transforms(
                 continue;
             }
             let v_i = placements[i].position; // nominal vertex, global frame
-            let c_i = placements[i].cursor_rotation_matrix;
+            let c_i = placements[i].rotation_matrix;
 
             // Rotate about pivot (active = r_group^T), then translate.
             let rotated = r_group.transpose() * (v_i - p);
             let v_i_prime = p + rotated + d_global;
 
-            // Per-surface decenter in cursor frame i.
+            // Per-surface decenter in surface-local frame i.
             let decenter_i = c_i * (v_i_prime - v_i);
 
-            // Per-surface rotation_offset in cursor frame i (passive).
+            // Per-surface rotation_offset in surface-local frame i (passive).
             let rot_off_mat = c_i * r_group * c_i.transpose();
             let rotation_offset_i = mat3x3_to_rotation3d(rot_off_mat);
 
@@ -378,6 +394,7 @@ fn component_first_idx(c: &Component) -> usize {
         Component::Element { surf_idxs } => *surf_idxs.first().unwrap_or(&usize::MAX),
         Component::Iris { stop_idx } => *stop_idx,
         Component::Mirror { surf_idx } => *surf_idx,
+        Component::ThinLens { surf_idx } => *surf_idx,
         Component::UnpairedSurface { surf_idx } => *surf_idx,
     }
 }
@@ -419,6 +436,11 @@ fn set_surface_displacement(
             rotation_offset: ro,
             ..
         }
+        | SurfaceSpec::ThinLens {
+            decenter: d,
+            rotation_offset: ro,
+            ..
+        }
         | SurfaceSpec::Iris {
             decenter: d,
             rotation_offset: ro,
diff --git a/crates/cherry-rs/src/gui/examples.rs b/crates/cherry-rs/src/gui/examples.rs
index ea5ed850..581caf0a 100644
--- a/crates/cherry-rs/src/gui/examples.rs
+++ b/crates/cherry-rs/src/gui/examples.rs
@@ -1,4 +1,43 @@
-use super::model::{BoundaryVariant, FieldMode, FieldRow, SurfaceRow, SurfaceVariant, SystemSpecs};
+use super::model::{
+    BoundaryVariant, FieldMode, FieldRow, SolveSpec, SurfaceRow, SurfaceVariant, SystemSpecs,
+};
+
+/// Simple thin lens: f = 100 mm, object at infinity.
+///
+/// Carries an M (marginal ray height) solve on the lens-to-image gap with a
+/// target height of 0, demonstrating how to keep the image plane at the
+/// lens's focus instead of hand-entering the back focal distance.
+pub fn thin_lens() -> SystemSpecs {
+    SystemSpecs {
+        surfaces: vec![
+            SurfaceRow::new_object("Infinity"),
+            SurfaceRow::new_thin_lens("12.5", "100.0", "50.0", "1.0"),
+            SurfaceRow::new_image(),
+        ],
+        fields: vec![FieldRow {
+            chi: "0.0".into(),
+            phi: "90.0".into(),
+            x: "0.0".into(),
+        }],
+        aperture_semi_diameter: "12.5".into(),
+        wavelengths: vec!["0.5876".into()],
+        field_mode: FieldMode::Angle,
+        use_materials: false,
+        selected_materials: Vec::new(),
+        cross_section_n_rays: 3,
+        full_pupil_spacing: "0.1".into(),
+        n_fan_rays: 65,
+        background_n: "1.0".into(),
+        background_material_key: None,
+        stop_surface: None,
+        solves: vec![SolveSpec::MarginalRayHeight {
+            gap_index: 1,
+            target_height: 0.0,
+            wavelength_id: 0,
+        }],
+        lens_groups: Vec::new(),
+    }
+}
 
 /// Figure-Z two-mirror system: two flat mirrors at 30° tilt, separated by 100
 /// mm, returning the beam parallel to the z-axis.
@@ -14,6 +53,7 @@ pub fn mirrors_figure_z() -> SystemSpecs {
                 semi_diameter: "12.7".into(),
                 radius_of_curvature: "Infinity".into(),
                 conic_constant: String::new(),
+                focal_length: String::new(),
                 theta: "30".into(),
                 psi: "0".into(),
                 material_key: None,
@@ -26,6 +66,7 @@ pub fn mirrors_figure_z() -> SystemSpecs {
                 semi_diameter: "12.7".into(),
                 radius_of_curvature: "Infinity".into(),
                 conic_constant: String::new(),
+                focal_length: String::new(),
                 theta: "30".into(),
                 psi: "0".into(),
                 material_key: None,
@@ -148,6 +189,7 @@ pub fn convexplano_lens_with_materials() -> SystemSpecs {
                 semi_diameter: "12.5".into(),
                 radius_of_curvature: "Infinity".into(),
                 conic_constant: "0.0".into(),
+                focal_length: String::new(),
                 theta: "0".into(),
                 psi: "0".into(),
                 material_key: Some("other:air:Ciddor".into()),
@@ -160,6 +202,7 @@ pub fn convexplano_lens_with_materials() -> SystemSpecs {
                 semi_diameter: "12.5".into(),
                 radius_of_curvature: "25.8".into(),
                 conic_constant: String::new(),
+                focal_length: String::new(),
                 theta: "0".into(),
                 psi: "0".into(),
                 material_key: Some("popular_glass:BK7:SCHOTT".into()),
@@ -172,6 +215,7 @@ pub fn convexplano_lens_with_materials() -> SystemSpecs {
                 semi_diameter: "12.5".into(),
                 radius_of_curvature: "Infinity".into(),
                 conic_constant: String::new(),
+                focal_length: String::new(),
                 theta: "0".into(),
                 psi: "0".into(),
                 material_key: Some("other:air:Ciddor".into()),
@@ -221,6 +265,7 @@ pub fn f_theta_scan_lens() -> SystemSpecs {
                 semi_diameter: "12.5".into(),
                 radius_of_curvature: "Infinity".into(),
                 conic_constant: "0.0".into(),
+                focal_length: String::new(),
                 theta: "0".into(),
                 psi: "0".into(),
                 material_key: Some("other:air:Ciddor".into()),
@@ -233,6 +278,7 @@ pub fn f_theta_scan_lens() -> SystemSpecs {
                 semi_diameter: "0.5".into(),
                 radius_of_curvature: "Infinity".into(),
                 conic_constant: "0.0".into(),
+                focal_length: String::new(),
                 theta: "0".into(),
                 psi: "0".into(),
                 material_key: Some("other:air:Ciddor".into()),
@@ -245,6 +291,7 @@ pub fn f_theta_scan_lens() -> SystemSpecs {
                 semi_diameter: "2".into(),
                 radius_of_curvature: "-2.2136".into(),
                 conic_constant: String::new(),
+                focal_length: String::new(),
                 theta: "0".into(),
                 psi: "0".into(),
                 material_key: Some("specs:SCHOTT-optical:N-SF57".into()),
@@ -257,6 +304,7 @@ pub fn f_theta_scan_lens() -> SystemSpecs {
                 semi_diameter: "2".into(),
                 radius_of_curvature: "-2.6575".into(),
                 conic_constant: String::new(),
+                focal_length: String::new(),
                 theta: "0".into(),
                 psi: "0".into(),
                 material_key: Some("other:air:Ciddor".into()),
@@ -269,6 +317,7 @@ pub fn f_theta_scan_lens() -> SystemSpecs {
                 semi_diameter: "2".into(),
                 radius_of_curvature: "-5.5022".into(),
                 conic_constant: String::new(),
+                focal_length: String::new(),
                 theta: "0".into(),
                 psi: "0".into(),
                 material_key: Some("specs:SCHOTT-optical:N-SF57".into()),
@@ -281,6 +330,7 @@ pub fn f_theta_scan_lens() -> SystemSpecs {
                 semi_diameter: "2".into(),
                 radius_of_curvature: "-3.8129".into(),
                 conic_constant: String::new(),
+                focal_length: String::new(),
                 theta: "0".into(),
                 psi: "0".into(),
                 material_key: Some("other:air:Ciddor".into()),
@@ -293,6 +343,7 @@ pub fn f_theta_scan_lens() -> SystemSpecs {
                 semi_diameter: "3".into(),
                 radius_of_curvature: "7.9951".into(),
                 conic_constant: String::new(),
+                focal_length: String::new(),
                 theta: "0".into(),
                 psi: "0".into(),
                 material_key: Some("specs:SCHOTT-optical:N-SF57".into()),
@@ -305,6 +356,7 @@ pub fn f_theta_scan_lens() -> SystemSpecs {
                 semi_diameter: "3".into(),
                 radius_of_curvature: "8.3651".into(),
                 conic_constant: String::new(),
+                focal_length: String::new(),
                 theta: "0".into(),
                 psi: "0".into(),
                 material_key: Some("other:air:Ciddor".into()),
@@ -351,6 +403,7 @@ pub fn galvo_scan_lens_negrean_mansvelder() -> SystemSpecs {
                 semi_diameter: String::new(),
                 radius_of_curvature: String::new(),
                 conic_constant: String::new(),
+                focal_length: String::new(),
                 theta: "0".into(),
                 psi: "0".into(),
                 material_key: Some("other:air:Ciddor".into()),
@@ -363,6 +416,7 @@ pub fn galvo_scan_lens_negrean_mansvelder() -> SystemSpecs {
                 semi_diameter: "2".into(),
                 radius_of_curvature: "Infinity".into(),
                 conic_constant: String::new(),
+                focal_length: String::new(),
                 theta: "-45".into(),
                 psi: "0".into(),
                 material_key: Some("other:air:Ciddor".into()),
@@ -375,6 +429,7 @@ pub fn galvo_scan_lens_negrean_mansvelder() -> SystemSpecs {
                 semi_diameter: "9".into(),
                 radius_of_curvature: "21.423".into(),
                 conic_constant: String::new(),
+                focal_length: String::new(),
                 theta: "0".into(),
                 psi: "0".into(),
                 material_key: Some("specs:SCHOTT-optical:N-KZFS5".into()),
@@ -387,6 +442,7 @@ pub fn galvo_scan_lens_negrean_mansvelder() -> SystemSpecs {
                 semi_diameter: "8".into(),
                 radius_of_curvature: "13.471".into(),
                 conic_constant: String::new(),
+                focal_length: String::new(),
                 theta: "0".into(),
                 psi: "0".into(),
                 material_key: Some("other:air:Ciddor".into()),
@@ -399,6 +455,7 @@ pub fn galvo_scan_lens_negrean_mansvelder() -> SystemSpecs {
                 semi_diameter: "15".into(),
                 radius_of_curvature: "88.222".into(),
                 conic_constant: String::new(),
+                focal_length: String::new(),
                 theta: "0".into(),
                 psi: "0".into(),
                 material_key: Some("specs:SCHOTT-optical:N-PK51".into()),
@@ -411,6 +468,7 @@ pub fn galvo_scan_lens_negrean_mansvelder() -> SystemSpecs {
                 semi_diameter: "15".into(),
                 radius_of_curvature: "-23.392".into(),
                 conic_constant: String::new(),
+                focal_length: String::new(),
                 theta: "0".into(),
                 psi: "0".into(),
                 material_key: Some("other:air:Ciddor".into()),
@@ -423,6 +481,7 @@ pub fn galvo_scan_lens_negrean_mansvelder() -> SystemSpecs {
                 semi_diameter: "15".into(),
                 radius_of_curvature: "211.304".into(),
                 conic_constant: String::new(),
+                focal_length: String::new(),
                 theta: "0".into(),
                 psi: "0".into(),
                 material_key: Some("specs:OHARA-optical:S-FPM2".into()),
@@ -435,6 +494,7 @@ pub fn galvo_scan_lens_negrean_mansvelder() -> SystemSpecs {
                 semi_diameter: "15".into(),
                 radius_of_curvature: "-20.385".into(),
                 conic_constant: String::new(),
+                focal_length: String::new(),
                 theta: "0".into(),
                 psi: "0".into(),
                 material_key: Some("specs:SCHOTT-optical:N-KZFS11".into()),
@@ -447,6 +507,7 @@ pub fn galvo_scan_lens_negrean_mansvelder() -> SystemSpecs {
                 semi_diameter: "15".into(),
                 radius_of_curvature: "Infinity".into(),
                 conic_constant: String::new(),
+                focal_length: String::new(),
                 theta: "0".into(),
                 psi: "0".into(),
                 material_key: Some("other:air:Ciddor".into()),
@@ -498,6 +559,7 @@ pub fn concave_mirror() -> SystemSpecs {
                 semi_diameter: "12.5".into(),
                 radius_of_curvature: "-200.0".into(),
                 conic_constant: String::new(),
+                focal_length: String::new(),
                 theta: "0".into(),
                 psi: "0".into(),
                 material_key: None,
@@ -527,7 +589,59 @@ pub fn concave_mirror() -> SystemSpecs {
         background_n: "1.0".into(),
         background_material_key: None,
         stop_surface: None,
-        solves: Vec::new(),
+        solves: vec![SolveSpec::MarginalRayHeight {
+            gap_index: 1,
+            target_height: 0.0,
+            wavelength_id: 0,
+        }],
         lens_groups: Vec::new(),
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::{SequentialModelBuilder, gui::convert};
+
+    fn parse(specs: &SystemSpecs) -> convert::ParsedSpecs {
+        #[cfg(not(feature = "ri-info"))]
+        return convert::convert_specs(specs).expect("convert");
+        #[cfg(feature = "ri-info")]
+        return convert::convert_specs(specs, &Default::default()).expect("convert");
+    }
+
+    /// The thin lens preset must parse into a valid, buildable model — the
+    /// bug this guards against: a stubbed-out field (e.g. an empty
+    /// `focal_length`) that compiles but fails at conversion time.
+    #[test]
+    fn thin_lens_example_converts_to_valid_model() {
+        let specs = thin_lens();
+        let parsed = parse(&specs);
+        SequentialModelBuilder::new()
+            .gap_specs(parsed.gaps)
+            .surface_specs(parsed.surfaces)
+            .wavelengths(parsed.wavelengths)
+            .build()
+            .expect("model");
+    }
+
+    /// The M solve on the lens-to-image gap must resolve to the lens's back
+    /// focal distance (== focal length, for a thin lens in air with an
+    /// object at infinity), demonstrating that the image plane tracks the
+    /// lens's focus rather than a hand-entered distance.
+    #[test]
+    fn thin_lens_example_solve_places_image_at_focus() {
+        let specs = thin_lens();
+        let parsed = parse(&specs);
+        let build_result = SequentialModelBuilder::new()
+            .gap_specs(parsed.gaps)
+            .surface_specs(parsed.surfaces)
+            .wavelengths(parsed.wavelengths)
+            .solves(parsed.solves)
+            .build()
+            .expect("model with solve applied");
+
+        let solved_thickness = build_result.gap_specs[1].thickness;
+        approx::assert_abs_diff_eq!(solved_thickness, 100.0, epsilon = 1e-6);
+    }
+}
diff --git a/crates/cherry-rs/src/gui/model.rs b/crates/cherry-rs/src/gui/model.rs
index 8cac3d82..93f4806b 100644
--- a/crates/cherry-rs/src/gui/model.rs
+++ b/crates/cherry-rs/src/gui/model.rs
@@ -142,6 +142,7 @@ pub enum SurfaceVariant {
     Object,
     Sphere,
     Conic,
+    ThinLens,
     Iris,
     Probe,
     Image,
@@ -153,6 +154,7 @@ impl SurfaceVariant {
     pub const SELECTABLE: &[SurfaceVariant] = &[
         SurfaceVariant::Sphere,
         SurfaceVariant::Conic,
+        SurfaceVariant::ThinLens,
         SurfaceVariant::Iris,
         SurfaceVariant::Probe,
     ];
@@ -164,6 +166,7 @@ impl std::fmt::Display for SurfaceVariant {
             SurfaceVariant::Object => write!(f, "Object"),
             SurfaceVariant::Sphere => write!(f, "Sphere"),
             SurfaceVariant::Conic => write!(f, "Conic"),
+            SurfaceVariant::ThinLens => write!(f, "Thin Lens"),
             SurfaceVariant::Iris => write!(f, "Iris"),
             SurfaceVariant::Probe => write!(f, "Probe"),
             SurfaceVariant::Image => write!(f, "Image"),
@@ -202,6 +205,9 @@ pub struct SurfaceRow {
     pub semi_diameter: String,
     pub radius_of_curvature: String,
     pub conic_constant: String,
+    /// Focal length. Only meaningful for `SurfaceVariant::ThinLens`.
+    #[serde(default)]
+    pub focal_length: String,
     /// Tilt in UF plane (about cursor-R axis), degrees. Only meaningful for
     /// reflecting Conic surfaces.
     #[serde(default = "default_zero")]
@@ -226,6 +232,7 @@ impl SurfaceRow {
             semi_diameter: String::new(),
             radius_of_curvature: String::new(),
             conic_constant: String::new(),
+            focal_length: String::new(),
             theta: "0".into(),
             psi: "0".into(),
             material_key: None,
@@ -247,6 +254,7 @@ impl SurfaceRow {
             semi_diameter: semi_diameter.into(),
             radius_of_curvature: radius_of_curvature.into(),
             conic_constant: conic_constant.into(),
+            focal_length: String::new(),
             theta: "0".into(),
             psi: "0".into(),
             material_key: None,
@@ -267,6 +275,28 @@ impl SurfaceRow {
             semi_diameter: semi_diameter.into(),
             radius_of_curvature: radius_of_curvature.into(),
             conic_constant: String::new(),
+            focal_length: String::new(),
+            theta: "0".into(),
+            psi: "0".into(),
+            material_key: None,
+        }
+    }
+
+    pub fn new_thin_lens(
+        semi_diameter: &str,
+        focal_length: &str,
+        thickness: &str,
+        refractive_index: &str,
+    ) -> Self {
+        Self {
+            variant: SurfaceVariant::ThinLens,
+            boundary_variant: BoundaryVariant::Refracting,
+            refractive_index: refractive_index.into(),
+            thickness: thickness.into(),
+            semi_diameter: semi_diameter.into(),
+            radius_of_curvature: String::new(),
+            conic_constant: String::new(),
+            focal_length: focal_length.into(),
             theta: "0".into(),
             psi: "0".into(),
             material_key: None,
@@ -282,6 +312,7 @@ impl SurfaceRow {
             semi_diameter: semi_diameter.into(),
             radius_of_curvature: String::new(),
             conic_constant: String::new(),
+            focal_length: String::new(),
             theta: "0".into(),
             psi: "0".into(),
             material_key: None,
@@ -297,6 +328,7 @@ impl SurfaceRow {
             semi_diameter: String::new(),
             radius_of_curvature: String::new(),
             conic_constant: String::new(),
+            focal_length: String::new(),
             theta: "0".into(),
             psi: "0".into(),
             material_key: None,
diff --git a/crates/cherry-rs/src/gui/panels/lens_overlay.rs b/crates/cherry-rs/src/gui/panels/lens_overlay.rs
index 9175a8e1..92c354d5 100644
--- a/crates/cherry-rs/src/gui/panels/lens_overlay.rs
+++ b/crates/cherry-rs/src/gui/panels/lens_overlay.rs
@@ -26,7 +26,9 @@ fn component_first_idx(c: &Component) -> usize {
             .first()
             .expect("Element must have at least one surface"),
         Component::Iris { stop_idx } => *stop_idx,
-        Component::Mirror { surf_idx } | Component::UnpairedSurface { surf_idx } => *surf_idx,
+        Component::Mirror { surf_idx }
+        | Component::ThinLens { surf_idx }
+        | Component::UnpairedSurface { surf_idx } => *surf_idx,
     }
 }
 
@@ -43,6 +45,7 @@ fn default_group_name(c: &Component) -> String {
         }
         Component::Iris { stop_idx } => format!("Iris ({stop_idx})"),
         Component::Mirror { surf_idx } => format!("Mirror ({surf_idx})"),
+        Component::ThinLens { surf_idx } => format!("Thin Lens ({surf_idx})"),
         Component::UnpairedSurface { surf_idx } => format!("Surface ({surf_idx})"),
     }
 }
diff --git a/crates/cherry-rs/src/gui/panels/surfaces.rs b/crates/cherry-rs/src/gui/panels/surfaces.rs
index db3a1502..678df01f 100644
--- a/crates/cherry-rs/src/gui/panels/surfaces.rs
+++ b/crates/cherry-rs/src/gui/panels/surfaces.rs
@@ -33,6 +33,10 @@ pub fn surfaces_panel(
         .surfaces
         .iter()
         .any(|s| s.variant == SurfaceVariant::Conic);
+    let has_thin_lens = specs
+        .surfaces
+        .iter()
+        .any(|s| s.variant == SurfaceVariant::ThinLens);
 
     egui::ScrollArea::horizontal().show(ui, |ui| {
         let ctx = ui.ctx().clone();
@@ -53,6 +57,12 @@ pub fn surfaces_panel(
             table
         };
 
+        let table = if has_thin_lens {
+            table.column(Column::initial(90.0).resizable(true)) // Focal Length
+        } else {
+            table
+        };
+
         let table = table
             .column(Column::initial(80.0).resizable(true)) // Thickness
             .column(Column::initial(n_col_width).resizable(true)); // n / Material
@@ -81,6 +91,9 @@ pub fn surfaces_panel(
                 if has_conic {
                     header.col(|ui| header_cell(ui, None, "Conic"));
                 }
+                if has_thin_lens {
+                    header.col(|ui| header_cell(ui, None, "Focal Length"));
+                }
                 header.col(|ui| header_cell(ui, None, "Thickness"));
                 header.col(|ui| header_cell(ui, None, "n"));
                 if has_reflecting {
@@ -115,6 +128,7 @@ pub fn surfaces_panel(
                         let is_image = surf.variant == SurfaceVariant::Image;
                         let is_conic = surf.variant == SurfaceVariant::Conic;
                         let is_sphere = surf.variant == SurfaceVariant::Sphere;
+                        let is_thin_lens = surf.variant == SurfaceVariant::ThinLens;
                         let is_curved = is_conic || is_sphere;
                         let is_locked = is_object || is_image;
 
@@ -263,6 +277,30 @@ pub fn surfaces_panel(
                             });
                         }
 
+                        // Focal Length (only when the system has Thin Lens surfaces)
+                        if has_thin_lens {
+                            row.col(|ui| {
+                                if is_thin_lens {
+                                    // Normalize empty string to a sensible nonzero
+                                    // default so the stored value is always valid
+                                    // (unlike RoC/conic constant, 0 is not a legal
+                                    // focal length).
+                                    if surf.focal_length.is_empty() {
+                                        surf.focal_length = "100".into();
+                                        changed = true;
+                                    }
+                                    changed |= drag_inf(
+                                        ui,
+                                        &mut surf.focal_length,
+                                        row_idx,
+                                        "fl",
+                                        f64::NEG_INFINITY..=f64::INFINITY,
+                                        1.0,
+                                    );
+                                }
+                            });
+                        }
+
                         // Thickness
                         row.col(|ui| {
                             if !is_image {
@@ -638,6 +676,80 @@ mod tests {
         harness.get_by_label("Nominal Rotation");
     }
 
+    /// The "Focal Length" column is absent without a thin lens.
+    #[test]
+    fn focal_length_column_absent_without_thin_lens() {
+        let mut specs = minimal_specs();
+        let mut harness = Harness::builder()
+            .with_size(egui::vec2(2000.0, 600.0))
+            .build_ui(|ui| {
+                default_panel(ui, &mut specs);
+            });
+        harness.run();
+        assert!(
+            harness.query_all_by_label("Focal Length").next().is_none(),
+            "Focal Length column should not appear without a thin lens"
+        );
+    }
+
+    /// The "Focal Length" column appears, and is editable, when a thin lens
+    /// surface exists.
+    #[test]
+    fn focal_length_column_present_with_thin_lens() {
+        let mut specs = SystemSpecs {
+            surfaces: vec![
+                SurfaceRow::new_object("Infinity"),
+                SurfaceRow::new_thin_lens("12.5", "100.0", "100.0", "1.0"),
+                SurfaceRow::new_image(),
+            ],
+            ..Default::default()
+        };
+        let mut harness = Harness::builder()
+            .with_size(egui::vec2(2000.0, 600.0))
+            .build_ui(|ui| {
+                default_panel(ui, &mut specs);
+            });
+        harness.run();
+        harness.get_by_label("Focal Length");
+    }
+
+    /// Regression: a thin lens row with an empty `focal_length` (e.g. just
+    /// switched from another variant via the dropdown, which doesn't reset
+    /// fields) must self-heal to a nonzero default on render, the same way
+    /// the Conic Constant column self-heals to "0". An empty string is
+    /// rejected downstream by `convert.rs::parse_float`, unlike "0" which
+    /// would be silently accepted but physically invalid for a focal length.
+    #[test]
+    fn focal_length_empty_string_self_heals_to_nonzero_default() {
+        let mut specs = SystemSpecs {
+            surfaces: vec![
+                SurfaceRow::new_object("Infinity"),
+                SurfaceRow::new_thin_lens("12.5", "", "100.0", "1.0"),
+                SurfaceRow::new_image(),
+            ],
+            ..Default::default()
+        };
+        {
+            let mut harness = Harness::builder()
+                .with_size(egui::vec2(2000.0, 600.0))
+                .build_ui(|ui| {
+                    default_panel(ui, &mut specs);
+                });
+            harness.run();
+        }
+
+        let focal_length = &specs.surfaces[1].focal_length;
+        assert!(
+            !focal_length.is_empty(),
+            "focal_length should self-heal to a nonzero default, not stay empty"
+        );
+        assert_ne!(
+            focal_length.as_str(),
+            "0",
+            "focal_length must not default to 0 (infinite power, divide-by-zero downstream)"
+        );
+    }
+
     fn lens_specs() -> SystemSpecs {
         SystemSpecs {
             surfaces: vec![
diff --git a/crates/cherry-rs/src/gui/result_package.rs b/crates/cherry-rs/src/gui/result_package.rs
index 471dd4c2..4109dbda 100644
--- a/crates/cherry-rs/src/gui/result_package.rs
+++ b/crates/cherry-rs/src/gui/result_package.rs
@@ -24,9 +24,6 @@ pub struct SurfaceDesc {
     pub pos: Vec3,
     /// Rotation matrix from global into the surface's local coordinate system.
     pub rot_mat: Mat3x3,
-    /// Cursor rotation matrix C_i at this surface (passive, global→cursor
-    /// frame).
-    pub cursor_rot_mat: Mat3x3,
 }
 
 /// Lightweight description of a field point for display.
diff --git a/crates/cherry-rs/src/gui/windows/cross_section.rs b/crates/cherry-rs/src/gui/windows/cross_section.rs
index f5c79c53..7cc60189 100644
--- a/crates/cherry-rs/src/gui/windows/cross_section.rs
+++ b/crates/cherry-rs/src/gui/windows/cross_section.rs
@@ -200,6 +200,13 @@ impl CrossSectionWindow {
             draw_element(&painter, elem, &w2s, ui.visuals());
         }
 
+        // Draw optical axis before rays so rays appear on top.
+        if self.annotations.show_axis {
+            for path in &geom.axis_paths {
+                draw_axis(&painter, path, &w2s);
+            }
+        }
+
         // Draw rays.
         for (wl_idx, paths) in geom.ray_paths.iter().enumerate() {
             let color = wavelengths
@@ -210,10 +217,7 @@ impl CrossSectionWindow {
             draw_rays(&painter, paths, &w2s, color);
         }
 
-        // Draw annotations.
-        if self.annotations.show_axis {
-            draw_axis(&painter, &geom.axis_path, &w2s);
-        }
+        // Draw remaining annotations.
         if self.annotations.show_scalebar {
             draw_scalebar(&painter, rect, &geom.bounding_box);
         }
@@ -229,9 +233,11 @@ impl CrossSectionWindow {
             d.remove::<usize>(hover_id);
             v
         });
+        // surface_frames is path × step; for the hover annotation use path 0
+        // (single-path systems: step == store idx; multipath GUI TBD).
         if self.annotations.show_ruf_axes
             && let Some(idx) = hover_idx
-            && let Some(Some(frame)) = geom.surface_frames.get(idx)
+            && let Some(Some(frame)) = geom.surface_frames.first().and_then(|p| p.get(idx))
         {
             draw_ruf_axes(&painter, frame, &w2s, self.cutting_plane);
         }
@@ -396,6 +402,93 @@ fn draw_element(
         DrawElement::FlatPlane { p1, p2, kind } => {
             draw_flat_plane(painter, *p1, *p2, *kind, w2s);
         }
+        DrawElement::ThinLens {
+            center_z,
+            center_t,
+            fwd_z,
+            fwd_t,
+            half_gap,
+            converging,
+        } => {
+            draw_thin_lens(
+                painter,
+                *center_z as f32,
+                *center_t as f32,
+                *fwd_z as f32,
+                *fwd_t as f32,
+                *half_gap as f32,
+                *converging,
+                w2s,
+            );
+        }
+    }
+}
+
+/// Color used for the thin-lens glyph, distinct from the rest of the
+/// cross-section palette (blue = glass lens, orange = mirror/unpaired
+/// profile, green = image, yellow = probe, gray = object).
+const THIN_LENS_COLOR: egui::Color32 = egui::Color32::from_rgb(150, 80, 200);
+
+#[allow(clippy::too_many_arguments)]
+fn draw_thin_lens(
+    painter: &egui::Painter,
+    center_z: f32,
+    center_t: f32,
+    fwd_z: f32,
+    fwd_t: f32,
+    half_gap: f32,
+    converging: bool,
+    w2s: &WorldToScreen,
+) {
+    let stroke = egui::Stroke::new(1.5, THIN_LENS_COLOR);
+    // Perpendicular to (fwd_z, fwd_t) is (-fwd_t, fwd_z) — direction along the
+    // lens surface in the 2-D (z, transverse) plot, same convention as
+    // `draw_stop`.
+    let perp_z = -fwd_t;
+    let perp_t = fwd_z;
+    let center = w2s.map(center_z, center_t);
+    let top = w2s.map(center_z + perp_z * half_gap, center_t + perp_t * half_gap);
+    let bot = w2s.map(center_z - perp_z * half_gap, center_t - perp_t * half_gap);
+    painter.line_segment([top, bot], stroke);
+
+    // Compute outward directions in screen space (not world space) so the
+    // arrowheads are correct regardless of axis scaling/flips in `w2s`.
+    draw_arrowhead(
+        painter,
+        top,
+        (top - center).normalized(),
+        converging,
+        stroke,
+    );
+    draw_arrowhead(
+        painter,
+        bot,
+        (bot - center).normalized(),
+        converging,
+        stroke,
+    );
+}
+
+/// Draws a small "V" arrowhead at `tip`. `outward` is the unit direction
+/// pointing away from the lens center along the shaft. When `converging` is
+/// true the arrowhead opens outward (the standard symbol for a positive,
+/// converging focal length); when false it opens inward (negative,
+/// diverging).
+fn draw_arrowhead(
+    painter: &egui::Painter,
+    tip: egui::Pos2,
+    outward: egui::Vec2,
+    converging: bool,
+    stroke: egui::Stroke,
+) {
+    const LEN: f32 = 8.0;
+    const SPREAD: f32 = 0.45; // radians from the shaft direction
+    let dir = if converging { outward } else { -outward };
+    for sign in [-1.0_f32, 1.0] {
+        let angle = sign * SPREAD;
+        let (sin, cos) = angle.sin_cos();
+        let rotated = egui::vec2(dir.x * cos - dir.y * sin, dir.x * sin + dir.y * cos);
+        painter.line_segment([tip, tip - rotated * LEN], stroke);
     }
 }
 
@@ -852,7 +945,9 @@ fn render_svg(
         r#"<rect width="{w}" height="{h}" fill="none" stroke="{border}" stroke-width="1"/>"#
     ));
 
-    svg_axis(&mut s, &geom.axis_path, &w2s, scalebar_color);
+    for path in &geom.axis_paths {
+        svg_axis(&mut s, path, &w2s, scalebar_color);
+    }
 
     for elem in &geom.elements {
         match elem {
@@ -890,6 +985,25 @@ fn render_svg(
             DrawElement::FlatPlane { p1, p2, kind } => {
                 svg_flat_plane(&mut s, *p1, *p2, *kind, &w2s);
             }
+            DrawElement::ThinLens {
+                center_z,
+                center_t,
+                fwd_z,
+                fwd_t,
+                half_gap,
+                converging,
+            } => {
+                svg_thin_lens(
+                    &mut s,
+                    *center_z,
+                    *center_t,
+                    *fwd_z,
+                    *fwd_t,
+                    *half_gap,
+                    *converging,
+                    &w2s,
+                );
+            }
         }
     }
 
@@ -1028,6 +1142,62 @@ fn svg_stop(
     ));
 }
 
+/// Hex color for the thin-lens glyph, matching `THIN_LENS_COLOR`.
+const THIN_LENS_SVG_COLOR: &str = "#9650c8";
+
+#[allow(clippy::too_many_arguments)]
+fn svg_thin_lens(
+    s: &mut String,
+    center_z: f64,
+    center_t: f64,
+    fwd_z: f64,
+    fwd_t: f64,
+    half_gap: f64,
+    converging: bool,
+    w2s: &WorldToSvg,
+) {
+    let perp_z = -fwd_t;
+    let perp_t = fwd_z;
+    let center = w2s.map(center_z, center_t);
+    let top = w2s.map(center_z + perp_z * half_gap, center_t + perp_t * half_gap);
+    let bot = w2s.map(center_z - perp_z * half_gap, center_t - perp_t * half_gap);
+    s.push_str(&format!(
+        r#"<line x1="{:.2}" y1="{:.2}" x2="{:.2}" y2="{:.2}" stroke="{THIN_LENS_SVG_COLOR}" stroke-width="1.5"/>"#,
+        top.0, top.1, bot.0, bot.1
+    ));
+    svg_arrowhead(s, top, outward_dir(top, center), converging);
+    svg_arrowhead(s, bot, outward_dir(bot, center), converging);
+}
+
+fn outward_dir(tip: (f64, f64), center: (f64, f64)) -> (f64, f64) {
+    let (dx, dy) = (tip.0 - center.0, tip.1 - center.1);
+    let len = (dx * dx + dy * dy).sqrt().max(f64::EPSILON);
+    (dx / len, dy / len)
+}
+
+/// Draws a small "V" arrowhead at `tip`, mirroring the painter-side
+/// `draw_arrowhead`. `outward` points away from the lens center; the
+/// arrowhead opens outward when `converging`, inward otherwise.
+fn svg_arrowhead(s: &mut String, tip: (f64, f64), outward: (f64, f64), converging: bool) {
+    const LEN: f64 = 8.0;
+    const SPREAD: f64 = 0.45;
+    let dir = if converging {
+        outward
+    } else {
+        (-outward.0, -outward.1)
+    };
+    for sign in [-1.0_f64, 1.0] {
+        let angle = sign * SPREAD;
+        let (sin, cos) = angle.sin_cos();
+        let rotated = (dir.0 * cos - dir.1 * sin, dir.0 * sin + dir.1 * cos);
+        let end = (tip.0 - rotated.0 * LEN, tip.1 - rotated.1 * LEN);
+        s.push_str(&format!(
+            r#"<line x1="{:.2}" y1="{:.2}" x2="{:.2}" y2="{:.2}" stroke="{THIN_LENS_SVG_COLOR}" stroke-width="1.5"/>"#,
+            tip.0, tip.1, end.0, end.1
+        ));
+    }
+}
+
 fn svg_flat_plane(
     s: &mut String,
     p1: [f64; 2],
@@ -1196,7 +1366,7 @@ mod tests {
                 },
                 elements: Vec::new(),
                 ray_paths: Vec::new(),
-                axis_path: Vec::new(),
+                axis_paths: Vec::new(),
                 surface_frames: Vec::new(),
             },
             xz: PlaneGeometry {
@@ -1206,7 +1376,7 @@ mod tests {
                 },
                 elements: Vec::new(),
                 ray_paths: Vec::new(),
-                axis_path: Vec::new(),
+                axis_paths: Vec::new(),
                 surface_frames: Vec::new(),
             },
         };
@@ -1232,4 +1402,63 @@ mod tests {
         harness.step();
         harness.get_by_label_contains("The optical axis leaves both coordinate planes");
     }
+
+    /// Rendering a `DrawElement::ThinLens` (painter path) must not panic, for
+    /// both converging and diverging glyphs.
+    #[test]
+    fn thin_lens_element_renders_without_panicking() {
+        use crate::views::cross_section::{Bounds2D, CrossSectionView, PlaneGeometry};
+
+        for converging in [true, false] {
+            let mut window = CrossSectionWindow::default();
+            let plane = PlaneGeometry {
+                bounding_box: Bounds2D {
+                    z: (-10.0, 10.0),
+                    transverse: (-15.0, 15.0),
+                },
+                elements: vec![DrawElement::ThinLens {
+                    center_z: 0.0,
+                    center_t: 0.0,
+                    fwd_z: 1.0,
+                    fwd_t: 0.0,
+                    half_gap: 12.5,
+                    converging,
+                }],
+                ray_paths: Vec::new(),
+                axis_paths: Vec::new(),
+                surface_frames: Vec::new(),
+            };
+            let cs = CrossSectionView {
+                wavelengths: vec![0.5876],
+                yz_valid: true,
+                xz_valid: true,
+                yz: plane,
+                xz: PlaneGeometry {
+                    bounding_box: Bounds2D {
+                        z: (-1.0, 1.0),
+                        transverse: (-1.0, 1.0),
+                    },
+                    elements: Vec::new(),
+                    ray_paths: Vec::new(),
+                    axis_paths: Vec::new(),
+                    surface_frames: Vec::new(),
+                },
+            };
+            let result = ResultPackage {
+                id: 1,
+                wavelengths: vec![0.5876],
+                surfaces: Vec::new(),
+                fields: Vec::new(),
+                field_specs: Vec::new(),
+                paraxial: None,
+                ray_trace: None,
+                cross_section: Some(cs),
+                error: None,
+                solved_values: Default::default(),
+                components: Vec::new(),
+            };
+            let mut harness = Harness::new(|ctx| show_window(&mut window, Some(&result), ctx));
+            harness.step();
+        }
+    }
 }
diff --git a/crates/cherry-rs/src/gui/windows/ray_fan.rs b/crates/cherry-rs/src/gui/windows/ray_fan.rs
index bbd5812e..4ff10a72 100644
--- a/crates/cherry-rs/src/gui/windows/ray_fan.rs
+++ b/crates/cherry-rs/src/gui/windows/ray_fan.rs
@@ -542,6 +542,7 @@ mod tests {
                     SurfaceKind::Probe => "Probe",
                     SurfaceKind::Iris => "Iris",
                     SurfaceKind::Sphere => "Sphere",
+                    SurfaceKind::ThinLens => "Thin Lens",
                     SurfaceKind::Custom => "Custom",
                 };
                 SurfaceDesc {
@@ -549,7 +550,6 @@ mod tests {
                     label: format!("{name} [{i}]"),
                     pos: p.position,
                     rot_mat: p.rotation_matrix,
-                    cursor_rot_mat: p.cursor_rotation_matrix,
                 }
             })
             .collect();
diff --git a/crates/cherry-rs/src/gui/windows/spot_diagram.rs b/crates/cherry-rs/src/gui/windows/spot_diagram.rs
index 8ce92c6f..b752f6e3 100644
--- a/crates/cherry-rs/src/gui/windows/spot_diagram.rs
+++ b/crates/cherry-rs/src/gui/windows/spot_diagram.rs
@@ -463,19 +463,19 @@ mod tests {
         let result = ResultPackage {
             id: 1,
             wavelengths: seq.wavelengths().to_vec(),
-            surfaces: seq
-                .surfaces()
-                .iter()
-                .zip(seq.placements().iter())
-                .enumerate()
-                .map(|(i, (_s, p))| crate::gui::result_package::SurfaceDesc {
-                    index: i,
-                    label: format!("S{i}"),
-                    pos: p.position,
-                    rot_mat: p.rotation_matrix,
-                    cursor_rot_mat: p.cursor_rotation_matrix,
-                })
-                .collect(),
+            surfaces: {
+                seq.surfaces()
+                    .iter()
+                    .zip(seq.placements().iter())
+                    .enumerate()
+                    .map(|(i, (_s, p))| crate::gui::result_package::SurfaceDesc {
+                        index: i,
+                        label: format!("S{i}"),
+                        pos: p.position,
+                        rot_mat: p.rotation_matrix,
+                    })
+                    .collect()
+            },
             fields: vec![
                 FieldDesc {
                     label: "0.000\u{00b0}".into(),
diff --git a/crates/cherry-rs/src/lib.rs b/crates/cherry-rs/src/lib.rs
index 2ca9e523..886a05f4 100644
--- a/crates/cherry-rs/src/lib.rs
+++ b/crates/cherry-rs/src/lib.rs
@@ -147,12 +147,14 @@ pub use core::{
     math::vec3::Vec3,
     ray::Ray,
     sequential_model::{
-        SequentialModel, SequentialSubModel, Step,
+        CursorPlacement, SequentialModel, SequentialSubModel, Step,
         builder::{BuildResult, SequentialModelBuilder},
-        placement::Placement,
         solves::{FNumberSolve, MarginalRaySolve, Solve, SolveKind},
+        surface_placement::SurfacePlacement,
+    },
+    surfaces::{
+        BeamSplitter, Conic, Image, Iris, Object, Probe, Sphere, Surface, SurfaceKind, ThinLens,
     },
-    surfaces::{BeamSplitter, Conic, Image, Iris, Object, Probe, Sphere, Surface, SurfaceKind},
 };
 pub use specs::{
     aperture::ApertureSpec,
diff --git a/crates/cherry-rs/src/specs/surfaces.rs b/crates/cherry-rs/src/specs/surfaces.rs
index 54bf221c..6261624d 100644
--- a/crates/cherry-rs/src/specs/surfaces.rs
+++ b/crates/cherry-rs/src/specs/surfaces.rs
@@ -121,6 +121,17 @@ pub enum SurfaceSpec {
         #[cfg_attr(feature = "serde", serde(default = "default_rotation3d_none"))]
         rotation_offset: Rotation3D,
     },
+    /// An idealized thin lens with a focal length independent of the
+    /// surrounding media.
+    ThinLens {
+        semi_diameter: Float,
+        focal_length: Float,
+        rotation: Rotation3D,
+        #[cfg_attr(feature = "serde", serde(default = "default_zero_vec3"))]
+        decenter: Vec3,
+        #[cfg_attr(feature = "serde", serde(default = "default_rotation3d_none"))]
+        rotation_offset: Rotation3D,
+    },
     /// A user-defined surface type registered with a [`SurfaceRegistry`].
     ///
     /// `type_id` must match a key registered via
@@ -191,6 +202,7 @@ impl SurfaceSpec {
         match self {
             SurfaceSpec::Conic { rotation, .. }
             | SurfaceSpec::Sphere { rotation, .. }
+            | SurfaceSpec::ThinLens { rotation, .. }
             | SurfaceSpec::Image { rotation, .. }
             | SurfaceSpec::Probe { rotation, .. }
             | SurfaceSpec::Iris { rotation, .. }
@@ -210,6 +222,9 @@ impl SurfaceSpec {
             | SurfaceSpec::Sphere {
                 rotation_offset, ..
             }
+            | SurfaceSpec::ThinLens {
+                rotation_offset, ..
+            }
             | SurfaceSpec::Image {
                 rotation_offset, ..
             }
@@ -234,6 +249,7 @@ impl SurfaceSpec {
         match self {
             SurfaceSpec::Conic { decenter, .. }
             | SurfaceSpec::Sphere { decenter, .. }
+            | SurfaceSpec::ThinLens { decenter, .. }
             | SurfaceSpec::Image { decenter, .. }
             | SurfaceSpec::Probe { decenter, .. }
             | SurfaceSpec::Iris { decenter, .. }
@@ -320,6 +336,43 @@ mod tests {
         }
     }
 
+    // ThinLens: serialize then deserialize a ThinLens spec with non-zero
+    // decenter and rotation_offset; assert the round-tripped values match,
+    // and that focal_length survives the round trip.
+    #[cfg(feature = "serde")]
+    #[test]
+    fn thin_lens_serde_round_trip_preserves_focal_length_and_decenter() {
+        use crate::core::math::linalg::rotations::EulerAngles;
+
+        let phi = 0.2_f64;
+        let spec = SurfaceSpec::ThinLens {
+            semi_diameter: 25.0,
+            focal_length: 100.0,
+            rotation: Rotation3D::None,
+            decenter: Vec3::new(0.1, 0.2, 0.0),
+            rotation_offset: Rotation3D::IntrinsicPassiveRUF(EulerAngles(phi, 0.0, 0.0)),
+        };
+
+        let json = serde_json::to_string(&spec).expect("serialize");
+        let back: SurfaceSpec = serde_json::from_str(&json).expect("deserialize");
+
+        match back {
+            SurfaceSpec::ThinLens {
+                focal_length,
+                semi_diameter,
+                ..
+            } => {
+                assert!((focal_length - 100.0).abs() < 1e-15);
+                assert!((semi_diameter - 25.0).abs() < 1e-15);
+            }
+            other => panic!("unexpected spec variant: {:?}", other),
+        }
+
+        let d = back.decenter();
+        assert!((d.x() - 0.1).abs() < 1e-15, "decenter x: {}", d.x());
+        assert!((d.y() - 0.2).abs() < 1e-15, "decenter y: {}", d.y());
+    }
+
     // AT-10: deserializing a JSON string that omits decenter and rotation_offset
     // applies the correct defaults (zero vector and None).
     #[cfg(feature = "serde")]
diff --git a/crates/cherry-rs/src/views/components/mod.rs b/crates/cherry-rs/src/views/components/mod.rs
index a5d0b6ea..57848291 100644
--- a/crates/cherry-rs/src/views/components/mod.rs
+++ b/crates/cherry-rs/src/views/components/mod.rs
@@ -35,6 +35,11 @@ pub enum Component {
     Mirror {
         surf_idx: usize,
     },
+    /// An idealized thin lens, detected by `surface_kind()` rather than by
+    /// gap-pairing (it normally sits between two equal-index gaps).
+    ThinLens {
+        surf_idx: usize,
+    },
     UnpairedSurface {
         surf_idx: usize,
     },
@@ -89,6 +94,9 @@ pub fn components_view(
         } else if kind == SurfaceKind::Iris {
             non_elements.push(Component::Iris { stop_idx: i });
             claimed.insert(i);
+        } else if kind == SurfaceKind::ThinLens {
+            non_elements.push(Component::ThinLens { surf_idx: i });
+            claimed.insert(i);
         }
     }
 
@@ -203,6 +211,7 @@ pub fn components_view(
         Component::Element { surf_idxs } => *surf_idxs.first().unwrap_or(&usize::MAX),
         Component::Iris { stop_idx } => *stop_idx,
         Component::Mirror { surf_idx } => *surf_idx,
+        Component::ThinLens { surf_idx } => *surf_idx,
         Component::UnpairedSurface { surf_idx } => *surf_idx,
     });
     Ok(result)
@@ -523,6 +532,50 @@ mod tests {
         }));
     }
 
+    pub fn thin_lens_singlet() -> SequentialModel {
+        // A thin lens borders equal-index (air) gaps on both sides, so it
+        // can't be detected by the gap-pairing logic used for glass
+        // elements — it must be classified by surface_kind() instead.
+        let air = n!(1.0);
+
+        let surf_0 = SurfaceSpec::Object;
+        let gap_0 = GapSpec {
+            thickness: Float::INFINITY,
+            refractive_index: air.clone(),
+        };
+        let surf_1 = SurfaceSpec::ThinLens {
+            semi_diameter: 12.5,
+            focal_length: 100.0,
+            rotation: Rotation3D::None,
+            decenter: Vec3::new(0.0, 0.0, 0.0),
+            rotation_offset: Rotation3D::None,
+        };
+        let gap_1 = GapSpec {
+            thickness: 100.0,
+            refractive_index: air,
+        };
+        let surf_2 = SurfaceSpec::Image {
+            rotation: Rotation3D::None,
+            decenter: Vec3::new(0.0, 0.0, 0.0),
+            rotation_offset: Rotation3D::None,
+        };
+
+        let surfaces = vec![surf_0, surf_1, surf_2];
+        let gaps = vec![gap_0, gap_1];
+        let wavelengths = vec![0.5876];
+
+        SequentialModel::from_surface_specs(&gaps, &surfaces, &wavelengths, None).unwrap()
+    }
+
+    #[test]
+    fn test_thin_lens_is_standalone_component() {
+        let sequential_model = thin_lens_singlet();
+        let components = components_view(&sequential_model, n!(1.0)).unwrap();
+
+        assert_eq!(components.len(), 1);
+        assert!(components.contains(&Component::ThinLens { surf_idx: 1 }));
+    }
+
     #[test]
     fn test_silly_single_surface_and_stop() {
         // Sphere1 borders a glass gap on its right side even though the iris
diff --git a/crates/cherry-rs/src/views/cross_section.rs b/crates/cherry-rs/src/views/cross_section.rs
index 28b6b441..e5da619b 100644
--- a/crates/cherry-rs/src/views/cross_section.rs
+++ b/crates/cherry-rs/src/views/cross_section.rs
@@ -2,7 +2,10 @@
 
 use crate::{
     SequentialModel, SurfaceKind,
-    core::{Float, math::vec3::Vec3, sequential_model::placement::Placement, surfaces::Surface},
+    core::{
+        Float, math::vec3::Vec3, sequential_model::surface_placement::SurfacePlacement,
+        surfaces::Surface,
+    },
     views::{components::Component, ray_trace_3d::RayBundle},
 };
 
@@ -16,6 +19,19 @@ pub enum GlobalAxis {
 const N_PTS: usize = 64;
 const EPS: f64 = 1e-6;
 
+/// Project a 3D cursor position onto the cross-section plane.
+/// Returns `None` for positions at infinity (e.g. object at infinity).
+fn to_plot_point(p: Vec3, axis: GlobalAxis) -> Option<[f64; 2]> {
+    if !p.is_finite() {
+        return None;
+    }
+    let t = match axis {
+        GlobalAxis::Y => p.y(),
+        GlobalAxis::X => p.x(),
+    };
+    Some([p.z(), t])
+}
+
 /// The complete 2D cross-section view of a sequential optical system.
 pub struct CrossSectionView {
     pub wavelengths: Vec<f64>,
@@ -54,14 +70,14 @@ pub struct PlaneGeometry {
     pub elements: Vec<DrawElement>,
     /// ray_paths[wavelength_idx][path_idx] = Vec<[z, transverse]>
     pub ray_paths: Vec<Vec<Vec<[f64; 2]>>>,
-    /// On-axis positions [z, transverse] from the first finite surface to the
-    /// image surface. Starts at the object surface for finite-conjugate
-    /// systems, or at the first lens surface for infinite-conjugate systems
-    /// (where the object placement is infinite).
-    pub axis_path: Vec<[f64; 2]>,
-    /// Per-surface local RUF frame in 2D plot coordinates, indexed by surface
-    /// index. `None` for surfaces with an infinite vertex position.
-    pub surface_frames: Vec<Option<SurfaceFrame2D>>,
+    /// On-axis positions [z, transverse] per optical path.
+    /// Each inner `Vec` lists [z, transverse] points from the first finite
+    /// surface to the image surface for that path.
+    pub axis_paths: Vec<Vec<[f64; 2]>>,
+    /// Per-path, per-step local RUF frame in 2D plot coordinates.
+    /// `surface_frames[path_id][step_id]` is `None` for steps whose surface
+    /// has an infinite vertex position (e.g. object/image at infinity).
+    pub surface_frames: Vec<Vec<Option<SurfaceFrame2D>>>,
 }
 
 /// Axis-aligned bounding box in the (z, transverse) 2D coordinate system.
@@ -103,6 +119,20 @@ pub enum DrawElement {
         p2: [f64; 2],
         kind: FlatPlaneKind,
     },
+    ThinLens {
+        center_z: f64,
+        /// Transverse position of the surface vertex (non-zero when the
+        /// lens has a decenter along the cross-section's transverse axis).
+        center_t: f64,
+        /// Forward (optical-axis) direction at this surface in (z, t) plot
+        /// space.
+        fwd_z: f64,
+        fwd_t: f64,
+        half_gap: f64,
+        /// True for a positive (converging) focal length — drawn with
+        /// outward-pointing arrowheads; false (diverging) draws inward.
+        converging: bool,
+    },
 }
 
 #[derive(Debug, Clone, Copy, PartialEq, Eq)]
@@ -126,7 +156,6 @@ pub fn cross_section_view(
     components: &[Component],
 ) -> CrossSectionView {
     let wavelengths = model.wavelengths().to_vec();
-    let axis_dirs = model.axis_directions();
     let placements = model.placements();
 
     // A plane is valid when (a) the optical axis lies in it, (b) every surface
@@ -135,16 +164,24 @@ pub fn cross_section_view(
     // the plane). The surface normal in global coords is
     // `rotation_matrix.transpose() * local_z` because `rotation_matrix` maps
     // global→local.
-    let yz_valid = axis_dirs.iter().all(|d| d.x().abs() < EPS)
-        && placements.iter().all(|p| {
-            let n = p.rotation_matrix.transpose() * Vec3::new(0.0, 0.0, 1.0);
-            p.position.x().abs() < EPS && n.x().abs() < EPS
-        });
-    let xz_valid = axis_dirs.iter().all(|d| d.y().abs() < EPS)
-        && placements.iter().all(|p| {
-            let n = p.rotation_matrix.transpose() * Vec3::new(0.0, 0.0, 1.0);
-            p.position.y().abs() < EPS && n.y().abs() < EPS
-        });
+    let yz_valid = (0..model.path_count()).all(|pid| {
+        model
+            .path_steps(pid)
+            .iter()
+            .all(|s| s.axis_direction.x().abs() < EPS)
+    }) && placements.iter().all(|p| {
+        let n = p.rotation_matrix.transpose() * Vec3::new(0.0, 0.0, 1.0);
+        p.position.x().abs() < EPS && n.x().abs() < EPS
+    });
+    let xz_valid = (0..model.path_count()).all(|pid| {
+        model
+            .path_steps(pid)
+            .iter()
+            .all(|s| s.axis_direction.y().abs() < EPS)
+    }) && placements.iter().all(|p| {
+        let n = p.rotation_matrix.transpose() * Vec3::new(0.0, 0.0, 1.0);
+        p.position.y().abs() < EPS && n.y().abs() < EPS
+    });
 
     let yz = build_plane_geometry(model, cross_section_rays, GlobalAxis::Y, components);
     let xz = build_plane_geometry(model, cross_section_rays, GlobalAxis::X, components);
@@ -221,6 +258,32 @@ fn build_plane_geometry(
                     extent: largest_sd * 1.5,
                 });
             }
+            Component::ThinLens { surf_idx } => {
+                let placement = &placements[*surf_idx];
+                let center_z = placement.position.z();
+                let center_t = match axis {
+                    GlobalAxis::Y => placement.position.y(),
+                    GlobalAxis::X => placement.position.x(),
+                };
+                let fwd = placement.inv_rotation_matrix * Vec3::new(0.0, 0.0, 1.0);
+                let fwd_z = fwd.z();
+                let fwd_t = match axis {
+                    GlobalAxis::Y => fwd.y(),
+                    GlobalAxis::X => fwd.x(),
+                };
+                let sd = surfaces[*surf_idx].mask().semi_diameter();
+                // n_0 = n_1 = 1.0 is safe here: ThinLens::power() ignores
+                // both arguments, so this just recovers 1/focal_length's sign.
+                let converging = surfaces[*surf_idx].power(0.0, 1.0, 1.0).is_sign_positive();
+                elements.push(DrawElement::ThinLens {
+                    center_z,
+                    center_t,
+                    fwd_z,
+                    fwd_t,
+                    half_gap: sd,
+                    converging,
+                });
+            }
             Component::Mirror { surf_idx } => {
                 let pts = sample_surface(
                     surfaces[*surf_idx].as_ref(),
@@ -332,65 +395,68 @@ fn build_plane_geometry(
         }
     }
 
-    let axis_path: Vec<[f64; 2]> = model
-        .cursor_positions()
-        .iter()
-        .filter(|p| p.x().is_finite() && p.y().is_finite() && p.z().is_finite())
-        .map(|p| {
-            let t = match axis {
-                GlobalAxis::Y => p.y(),
-                GlobalAxis::X => p.x(),
-            };
-            [p.z(), t]
+    let axis_paths: Vec<Vec<[f64; 2]>> = (0..model.path_count())
+        .map(|path_id| {
+            model
+                .path_steps(path_id)
+                .iter()
+                .filter_map(|s| to_plot_point(s.cursor_position, axis))
+                .collect()
         })
         .collect();
 
     let bounding_box = compute_bounds(&elements, &ray_paths);
 
-    // Use cursor_positions (axis position before any decenter) rather than
+    // Use cursor_position (axis position before any decenter) rather than
     // placement.position (physical vertex) so that group tilts and decenters
     // don't displace the annotation away from the optical axis.
-    let cursor_positions = model.cursor_positions();
-    let surface_frames: Vec<Option<SurfaceFrame2D>> = placements
-        .iter()
-        .zip(cursor_positions.iter())
-        .map(|(p, cursor_pos)| {
-            if p.is_infinite() {
-                return None;
-            }
-            let crm_t = p.cursor_rotation_matrix.transpose();
-            let f = crm_t * Vec3::new(0.0, 0.0, 1.0);
-            let r = crm_t * Vec3::new(1.0, 0.0, 0.0);
-            let u = crm_t * Vec3::new(0.0, 1.0, 0.0);
-            let (vertex_z, vertex_t, f_z, f_t, t_z, t_t, oop_out_of_screen) = match axis {
-                GlobalAxis::Y => (
-                    cursor_pos.z(),
-                    cursor_pos.y(),
-                    f.z(),
-                    f.y(),
-                    u.z(),
-                    u.y(),
-                    r.x() < 0.0,
-                ),
-                GlobalAxis::X => (
-                    cursor_pos.z(),
-                    cursor_pos.x(),
-                    f.z(),
-                    f.x(),
-                    r.z(),
-                    r.x(),
-                    u.y() > 0.0,
-                ),
-            };
-            Some(SurfaceFrame2D {
-                vertex_z,
-                vertex_t,
-                f_z,
-                f_t,
-                t_z,
-                t_t,
-                oop_out_of_screen,
-            })
+    // Frames are per-path and per-step; cursor data is path-specific.
+    let surface_frames: Vec<Vec<Option<SurfaceFrame2D>>> = (0..model.path_count())
+        .map(|path_id| {
+            model
+                .path_steps(path_id)
+                .iter()
+                .zip(model.path_surface_indices(path_id).iter())
+                .map(|(step, &idx)| {
+                    if placements[idx].is_infinite() {
+                        return None;
+                    }
+                    let crm_t = step.cursor_rotation_matrix.transpose();
+                    let cursor_pos = step.cursor_position;
+                    let f = crm_t * Vec3::new(0.0, 0.0, 1.0);
+                    let r = crm_t * Vec3::new(1.0, 0.0, 0.0);
+                    let u = crm_t * Vec3::new(0.0, 1.0, 0.0);
+                    let (vertex_z, vertex_t, f_z, f_t, t_z, t_t, oop_out_of_screen) = match axis {
+                        GlobalAxis::Y => (
+                            cursor_pos.z(),
+                            cursor_pos.y(),
+                            f.z(),
+                            f.y(),
+                            u.z(),
+                            u.y(),
+                            r.x() < 0.0,
+                        ),
+                        GlobalAxis::X => (
+                            cursor_pos.z(),
+                            cursor_pos.x(),
+                            f.z(),
+                            f.x(),
+                            r.z(),
+                            r.x(),
+                            u.y() > 0.0,
+                        ),
+                    };
+                    Some(SurfaceFrame2D {
+                        vertex_z,
+                        vertex_t,
+                        f_z,
+                        f_t,
+                        t_z,
+                        t_t,
+                        oop_out_of_screen,
+                    })
+                })
+                .collect()
         })
         .collect();
 
@@ -398,7 +464,7 @@ fn build_plane_geometry(
         bounding_box,
         elements,
         ray_paths,
-        axis_path,
+        axis_paths,
         surface_frames,
     }
 }
@@ -410,7 +476,7 @@ fn build_plane_geometry(
 /// (z, x) pairs.
 fn sample_surface(
     surf: &dyn Surface,
-    placement: &Placement,
+    placement: &SurfacePlacement,
     axis: GlobalAxis,
     n_pts: usize,
 ) -> Vec<[f64; 2]> {
@@ -514,6 +580,33 @@ fn compute_bounds(elements: &[DrawElement], ray_paths: &[Vec<Vec<[f64; 2]>>]) ->
                 update(p1[0], p1[1], &mut z_min, &mut z_max, &mut t_min, &mut t_max);
                 update(p2[0], p2[1], &mut z_min, &mut z_max, &mut t_min, &mut t_max);
             }
+            DrawElement::ThinLens {
+                center_z,
+                center_t,
+                fwd_z,
+                fwd_t,
+                half_gap,
+                ..
+            } => {
+                let perp_z = -fwd_t;
+                let perp_t = fwd_z;
+                update(
+                    center_z + perp_z * half_gap,
+                    center_t + perp_t * half_gap,
+                    &mut z_min,
+                    &mut z_max,
+                    &mut t_min,
+                    &mut t_max,
+                );
+                update(
+                    center_z - perp_z * half_gap,
+                    center_t - perp_t * half_gap,
+                    &mut z_min,
+                    &mut z_max,
+                    &mut t_min,
+                    &mut t_max,
+                );
+            }
         }
     }
 
@@ -848,6 +941,74 @@ mod tests {
         );
     }
 
+    fn thin_lens_model(focal_length: Float) -> SequentialModel {
+        let air = n!(1.0);
+        let gaps = vec![
+            GapSpec {
+                thickness: Float::INFINITY,
+                refractive_index: air.clone(),
+            },
+            GapSpec {
+                thickness: 100.0,
+                refractive_index: air,
+            },
+        ];
+        let surfs = vec![
+            SurfaceSpec::Object,
+            SurfaceSpec::ThinLens {
+                semi_diameter: 12.5,
+                focal_length,
+                rotation: Rotation3D::None,
+                decenter: Vec3::new(0.0, 0.0, 0.0),
+                rotation_offset: Rotation3D::None,
+            },
+            SurfaceSpec::Image {
+                rotation: Rotation3D::None,
+                decenter: Vec3::new(0.0, 0.0, 0.0),
+                rotation_offset: Rotation3D::None,
+            },
+        ];
+        SequentialModel::from_surface_specs(&gaps, &surfs, &[0.5876], None)
+            .expect("build thin lens model")
+    }
+
+    fn find_thin_lens(cs: &CrossSectionView) -> Option<(f64, bool)> {
+        cs.yz.elements.iter().find_map(|e| match e {
+            DrawElement::ThinLens {
+                half_gap,
+                converging,
+                ..
+            } => Some((*half_gap, *converging)),
+            _ => None,
+        })
+    }
+
+    #[test]
+    fn thin_lens_draws_as_standalone_element_with_correct_semi_diameter() {
+        let model = thin_lens_model(100.0);
+        let components = components_view(&model, n!(1.0)).unwrap();
+        let cs = cross_section_view(&model, None, &components);
+
+        let (half_gap, converging) =
+            find_thin_lens(&cs).expect("ThinLens DrawElement not found in YZ plane");
+        assert!((half_gap - 12.5).abs() < 1e-9, "half_gap: {half_gap}");
+        assert!(converging, "positive focal length should be converging");
+    }
+
+    #[test]
+    fn thin_lens_diverging_for_negative_focal_length() {
+        let model = thin_lens_model(-100.0);
+        let components = components_view(&model, n!(1.0)).unwrap();
+        let cs = cross_section_view(&model, None, &components);
+
+        let (_, converging) =
+            find_thin_lens(&cs).expect("ThinLens DrawElement not found in YZ plane");
+        assert!(
+            !converging,
+            "negative focal length should be diverging (converging == false)"
+        );
+    }
+
     #[test]
     fn axis_path_starts_at_first_lens_for_infinite_object() {
         // convexplano_lens uses INFINITY for the first gap, so the object
@@ -859,21 +1020,24 @@ mod tests {
         let components = components_view(&model, air).unwrap();
         let cs = cross_section_view(&model, None, &components);
 
-        // Object is infinite → axis_path must not contain an infinite coordinate.
+        // Object is infinite → axis_paths[0] must not contain an infinite coordinate.
         assert!(
-            cs.yz
-                .axis_path
+            cs.yz.axis_paths[0]
                 .iter()
                 .all(|&[z, t]| z.is_finite() && t.is_finite()),
-            "axis_path must not contain infinite coordinates"
+            "axis_paths[0] must not contain infinite coordinates"
         );
         // For a straight system every on-axis transverse coordinate is zero.
-        for &[_z, t] in &cs.yz.axis_path {
+        for &[_z, t] in &cs.yz.axis_paths[0] {
             assert!(t.abs() < EPS, "expected on-axis transverse ≈ 0, got {t}");
         }
         // The path should include the two lens surfaces plus the image surface
         // (3 finite surfaces: front lens, back lens, image).
-        assert_eq!(cs.yz.axis_path.len(), 3, "expected 3 points in axis_path");
+        assert_eq!(
+            cs.yz.axis_paths[0].len(),
+            3,
+            "expected 3 points in axis_paths[0]"
+        );
     }
 
     #[test]
@@ -886,18 +1050,17 @@ mod tests {
 
         // Object(finite) + Sphere + Image = 3 surfaces.
         assert_eq!(
-            cs.yz.axis_path.len(),
+            cs.yz.axis_paths[0].len(),
             3,
             "expected 3 points: object, sphere, image; got {}",
-            cs.yz.axis_path.len()
+            cs.yz.axis_paths[0].len()
         );
         // All coordinates must be finite.
         assert!(
-            cs.yz
-                .axis_path
+            cs.yz.axis_paths[0]
                 .iter()
                 .all(|&[z, t]| z.is_finite() && t.is_finite()),
-            "axis_path must not contain infinite coordinates"
+            "axis_paths[0] must not contain infinite coordinates"
         );
     }
 
@@ -908,8 +1071,13 @@ mod tests {
         let cs = cross_section_view(&model, None, &components);
         assert_eq!(
             cs.yz.surface_frames.len(),
-            model.surfaces().len(),
-            "surface_frames must have one entry per surface"
+            model.path_count(),
+            "surface_frames outer length must equal path count"
+        );
+        assert_eq!(
+            cs.yz.surface_frames[0].len(),
+            model.path_steps(0).len(),
+            "surface_frames[0] inner length must equal step count for path 0"
         );
     }
 
@@ -924,7 +1092,7 @@ mod tests {
         let components = components_view(&model, air).unwrap();
         let cs = cross_section_view(&model, None, &components);
         assert!(
-            cs.yz.surface_frames[0].is_none(),
+            cs.yz.surface_frames[0][0].is_none(),
             "object at infinity must produce a None surface frame"
         );
     }
@@ -940,7 +1108,7 @@ mod tests {
 
         // Surface 1 is the sphere. Check directions; don't assert specific z since
         // Cherry places the first refracting surface at z=0.
-        let frame = cs.yz.surface_frames[1]
+        let frame = cs.yz.surface_frames[0][1]
             .as_ref()
             .expect("sphere surface frame must be Some");
         assert!(
@@ -1037,7 +1205,7 @@ mod tests {
 
         // Iris is at index 3, after the 45° fold — its cursor F should point in
         // the transverse direction (f_z ≈ 0, |f_t| ≈ 1).
-        let frame = cs.yz.surface_frames[3]
+        let frame = cs.yz.surface_frames[0][3]
             .as_ref()
             .expect("iris surface frame must be Some");
         assert!(
@@ -1071,10 +1239,10 @@ mod tests {
         let components = components_view(&model, n!(1.0)).unwrap();
         let cs = cross_section_view(&model, None, &components);
 
-        for &[_z, t] in &cs.yz.axis_path {
+        for &[_z, t] in &cs.yz.axis_paths[0] {
             assert!(
                 t.abs() < EPS,
-                "axis_path must stay on-axis despite decenter, got t={t}"
+                "axis_paths[0] must stay on-axis despite decenter, got t={t}"
             );
         }
     }
@@ -1086,7 +1254,7 @@ mod tests {
         let model = straight_sphere_model(Vec3::new(0.0, 2.0, 0.0), Rotation3D::None);
         let components = components_view(&model, n!(1.0)).unwrap();
         let cs = cross_section_view(&model, None, &components);
-        let frame = cs.yz.surface_frames[1]
+        let frame = cs.yz.surface_frames[0][1]
             .as_ref()
             .expect("sphere surface frame must be Some");
         assert!(
@@ -1158,10 +1326,10 @@ mod tests {
 
         // After the fold, at least one point in axis_path must have non-zero
         // transverse.
-        let has_nonzero_t = cs.yz.axis_path.iter().any(|&[_z, t]| t.abs() > 0.1);
+        let has_nonzero_t = cs.yz.axis_paths[0].iter().any(|&[_z, t]| t.abs() > 0.1);
         assert!(
             has_nonzero_t,
-            "expected non-zero transverse in axis_path after 45° fold"
+            "expected non-zero transverse in axis_paths[0] after 45° fold"
         );
     }
 
@@ -1247,4 +1415,23 @@ mod tests {
             "inner surface profile must be non-empty"
         );
     }
+
+    #[test]
+    fn multipath_produces_two_axis_paths() {
+        use crate::examples::beam_splitter;
+        use crate::specs::gaps::ConstantRefractiveIndex;
+        use std::rc::Rc;
+
+        let n_air: Rc<dyn crate::RefractiveIndexSpec> =
+            Rc::new(ConstantRefractiveIndex::new(1.0, 0.0));
+        let model = beam_splitter::two_path_model(n_air.clone(), &[0.5876], 10.0, 10.0);
+        let components = components_view(&model, n_air).unwrap();
+        let cs = cross_section_view(&model, None, &components);
+
+        assert_eq!(
+            cs.yz.axis_paths.len(),
+            2,
+            "expected one axis_path per optical path"
+        );
+    }
 }
diff --git a/crates/cherry-rs/src/views/paraxial.rs b/crates/cherry-rs/src/views/paraxial.rs
index 3f8f76dc..ba8fb918 100644
--- a/crates/cherry-rs/src/views/paraxial.rs
+++ b/crates/cherry-rs/src/views/paraxial.rs
@@ -17,13 +17,16 @@ use crate::{
         Float,
         math::{linalg::mat2x2::Mat2x2, vec3::Vec3},
         sequential_model::{
-            SequentialModel, SequentialSubModel, Step, first_physical_surface,
-            last_physical_surface, placement::Placement, propagate_tangential_vec,
-            reversed_surface_id,
+            CursorPlacement, SequentialModel, SequentialSubModel, Step, first_physical_step,
+            last_physical_step, propagate_tangential_vec, reversed_surface_id,
+            surface_placement::SurfacePlacement,
         },
         surfaces::Surface,
     },
-    specs::{fields::unique_tangential_vecs, surfaces::BoundaryKind},
+    specs::{
+        fields::unique_tangential_vecs,
+        surfaces::{BeamSplitterPathKind, BoundaryKind},
+    },
 };
 
 const DEFAULT_THICKNESS: Float = 0.0;
@@ -118,11 +121,12 @@ pub struct ParaxialViewDescription {
 
 /// A paraxial subview of an optical system.
 ///
-/// A paraxial subview is identified by a wavelength index and a tangential
-/// direction index. It is not created by the user, but rather by instantiating
-/// a new ParaxialView struct.
+/// A paraxial subview is identified by a path index, a wavelength index, and a
+/// tangential direction index. It is not created by the user, but rather by
+/// instantiating a new ParaxialView struct.
 #[derive(Debug)]
 pub struct ParaxialSubView {
+    path_id: usize,
     wavelength_id: usize,
     tangential_vec_id: usize,
     is_obj_space_telecentric: bool,
@@ -148,6 +152,7 @@ pub struct ParaxialSubView {
 #[derive(Debug)]
 #[cfg_attr(feature = "serde", derive(Serialize))]
 pub struct ParaxialSubViewDescription {
+    path_id: usize,
     wavelength_id: usize,
     tangential_vec_id: usize,
     aperture_stop: usize,
@@ -277,27 +282,45 @@ impl ParaxialView {
     ) -> Result<Self> {
         let surfaces = sequential_model.surfaces();
         let placements = sequential_model.placements();
-        let tangential_vecs: Vec<TangentialVector> =
-            if SequentialModel::is_rotationally_symmetric(placements) {
-                vec![Vec3::new(0.0, 1.0, 0.0)]
-            } else {
-                unique_tangential_vecs(field_specs)
-            };
+        let tangential_vecs: Vec<TangentialVector> = if sequential_model.is_rotationally_symmetric()
+        {
+            vec![Vec3::new(0.0, 1.0, 0.0)]
+        } else {
+            unique_tangential_vecs(field_specs)
+        };
 
-        let stop_surface = sequential_model.stop_surface();
         let mut subviews = Vec::new();
-        for (wav_idx, submodel) in sequential_model.submodels().iter().enumerate() {
-            for (v_idx, &v) in tangential_vecs.iter().enumerate() {
-                let data = SubModelData {
-                    sequential_sub_model: submodel as &dyn SequentialSubModel,
-                    surfaces,
-                    placements,
-                    field_specs,
-                    stop_surface,
-                };
-                let subview =
-                    ParaxialSubView::new(wav_idx, v_idx, &data, v, is_obj_space_telecentric)?;
-                subviews.push(subview);
+        for path_id in 0..sequential_model.path_count() {
+            let stop_surface = sequential_model.stop_surface_for_path(path_id);
+            let path_steps = sequential_model.path_steps(path_id);
+            let surface_indices = sequential_model.path_surface_indices(path_id);
+            let beam_splitter_arms = sequential_model.path_beam_splitter_arms(path_id);
+            for (wav_idx, submodel) in sequential_model
+                .submodels_for_path(path_id)
+                .iter()
+                .enumerate()
+            {
+                for (v_idx, &v) in tangential_vecs.iter().enumerate() {
+                    let data = SubModelData {
+                        sequential_sub_model: submodel as &dyn SequentialSubModel,
+                        surfaces,
+                        placements,
+                        surface_indices,
+                        beam_splitter_arms,
+                        path_steps,
+                        field_specs,
+                        stop_surface,
+                    };
+                    let subview = ParaxialSubView::new(
+                        path_id,
+                        wav_idx,
+                        v_idx,
+                        &data,
+                        v,
+                        is_obj_space_telecentric,
+                    )?;
+                    subviews.push(subview);
+                }
             }
         }
 
@@ -318,11 +341,24 @@ impl ParaxialView {
         }
     }
 
-    /// Returns the subview for the given wavelength and tangential-direction
-    /// indices, or `None` if no such subview exists.
+    /// Returns the subview for path 0, the given wavelength, and
+    /// tangential-direction indices, or `None` if no such subview exists.
     pub fn get(&self, wavelength_id: usize, tangential_vec_id: usize) -> Option<&ParaxialSubView> {
+        self.get_for_path(0, wavelength_id, tangential_vec_id)
+    }
+
+    /// Returns the subview for the given path, wavelength, and
+    /// tangential-direction indices, or `None` if no such subview exists.
+    pub fn get_for_path(
+        &self,
+        path_id: usize,
+        wavelength_id: usize,
+        tangential_vec_id: usize,
+    ) -> Option<&ParaxialSubView> {
         self.subviews.iter().find(|sv| {
-            sv.wavelength_id == wavelength_id && sv.tangential_vec_id == tangential_vec_id
+            sv.path_id == path_id
+                && sv.wavelength_id == wavelength_id
+                && sv.tangential_vec_id == tangential_vec_id
         })
     }
 
@@ -423,7 +459,10 @@ impl ParaxialView {
 struct SubModelData<'a> {
     sequential_sub_model: &'a dyn SequentialSubModel,
     surfaces: &'a [Box<dyn Surface>],
-    placements: &'a [Placement],
+    placements: &'a [SurfacePlacement],
+    surface_indices: &'a [usize],
+    beam_splitter_arms: &'a [Option<BeamSplitterPathKind>],
+    path_steps: &'a [CursorPlacement],
     field_specs: &'a [FieldSpec],
     stop_surface: Option<usize>,
 }
@@ -435,6 +474,7 @@ impl ParaxialSubView {
     /// plane (e.g. `(0,1,0)` for phi=90°). It is propagated through mirror
     /// surfaces internally to compute per-surface foreshortening.
     fn new(
+        path_id: usize,
         wavelength_id: usize,
         tangential_vec_id: usize,
         data: &SubModelData<'_>,
@@ -444,28 +484,64 @@ impl ParaxialSubView {
         let sequential_sub_model = data.sequential_sub_model;
         let surfaces = data.surfaces;
         let placements = data.placements;
+        let surface_indices = data.surface_indices;
+        let beam_splitter_arms = data.beam_splitter_arms;
+        let path_steps = data.path_steps;
         let field_specs = data.field_specs;
-        // Propagate v through mirror surfaces to get per-surface tangential vectors.
-        let per_surf_v: Vec<TangentialVector> = propagate_tangential_vec(v, surfaces, placements);
+        // Propagate v through this path's mirror surfaces to get per-step tangential
+        // vectors.
+        let per_surf_v: Vec<TangentialVector> =
+            propagate_tangential_vec(v, surfaces, placements, surface_indices);
 
-        let pseudo_marginal_ray =
-            Self::calc_pseudo_marginal_ray(sequential_sub_model, surfaces, placements)?;
-        let parallel_ray = Self::calc_parallel_ray(sequential_sub_model, surfaces, placements)?;
-        let reverse_parallel_ray =
-            Self::calc_reverse_parallel_ray(sequential_sub_model, surfaces, placements)?;
+        let pseudo_marginal_ray = Self::calc_pseudo_marginal_ray(
+            sequential_sub_model,
+            surfaces,
+            placements,
+            surface_indices,
+            beam_splitter_arms,
+            path_steps,
+        )?;
+        let parallel_ray = Self::calc_parallel_ray(
+            sequential_sub_model,
+            surfaces,
+            placements,
+            surface_indices,
+            beam_splitter_arms,
+            path_steps,
+        )?;
+        let reverse_parallel_ray = Self::calc_reverse_parallel_ray(
+            sequential_sub_model,
+            surfaces,
+            placements,
+            surface_indices,
+            beam_splitter_arms,
+            path_steps,
+        )?;
 
         let aperture_stop = match data.stop_surface {
             Some(i) => i,
-            None => {
-                Self::calc_aperture_stop(surfaces, placements, &pseudo_marginal_ray, &per_surf_v)
-            }
+            None => Self::calc_aperture_stop(
+                surfaces,
+                placements,
+                surface_indices,
+                path_steps,
+                &pseudo_marginal_ray,
+                &per_surf_v,
+            ),
         };
-        let back_focal_distance = Self::calc_back_focal_distance(surfaces, &parallel_ray)?;
-        let front_focal_distance =
-            Self::calc_front_focal_distance(surfaces, &reverse_parallel_ray)?;
+        let back_focal_distance =
+            Self::calc_back_focal_distance(surfaces, surface_indices, &parallel_ray)?;
+        let front_focal_distance = Self::calc_front_focal_distance(
+            sequential_sub_model,
+            surfaces,
+            surface_indices,
+            &reverse_parallel_ray,
+        )?;
         let marginal_ray = Self::calc_marginal_ray(
             surfaces,
             placements,
+            surface_indices,
+            path_steps,
             &pseudo_marginal_ray,
             &aperture_stop,
             &per_surf_v,
@@ -474,6 +550,9 @@ impl ParaxialSubView {
             sequential_sub_model,
             surfaces,
             placements,
+            surface_indices,
+            beam_splitter_arms,
+            path_steps,
             is_obj_space_telecentric,
             &aperture_stop,
             &per_surf_v,
@@ -483,6 +562,9 @@ impl ParaxialSubView {
             sequential_sub_model,
             surfaces,
             placements,
+            surface_indices,
+            beam_splitter_arms,
+            path_steps,
             &aperture_stop,
             &marginal_ray,
         )?;
@@ -497,14 +579,22 @@ impl ParaxialSubView {
             surfaces,
             sequential_sub_model,
             placements,
+            surface_indices,
+            beam_splitter_arms,
+            path_steps,
             v,
             field_specs,
             &entrance_pupil,
         )?;
-        let paraxial_image_plane =
-            Self::calc_paraxial_image_plane(surfaces, placements, &marginal_ray, &chief_ray)?;
+        let paraxial_image_plane = Self::calc_paraxial_image_plane(
+            surfaces,
+            placements,
+            surface_indices,
+            &marginal_ray,
+            &chief_ray,
+        )?;
 
-        let last_phys_id = last_physical_surface(surfaces)
+        let last_phys_id = last_physical_step(surface_indices, surfaces)
             .ok_or_else(|| anyhow!("There are no physical surfaces"))?;
         let n_image = sequential_sub_model
             .gaps()
@@ -517,6 +607,7 @@ impl ParaxialSubView {
         let image_space_fno = effective_focal_length / (2.0 * entrance_pupil.semi_diameter);
 
         Ok(Self {
+            path_id,
             wavelength_id,
             tangential_vec_id,
             is_obj_space_telecentric,
@@ -539,6 +630,7 @@ impl ParaxialSubView {
 
     fn describe(&self) -> ParaxialSubViewDescription {
         ParaxialSubViewDescription {
+            path_id: self.path_id,
             wavelength_id: self.wavelength_id,
             tangential_vec_id: self.tangential_vec_id,
             aperture_stop: self.aperture_stop,
@@ -557,6 +649,10 @@ impl ParaxialSubView {
         }
     }
 
+    pub fn path_id(&self) -> usize {
+        self.path_id
+    }
+
     pub fn wavelength_id(&self) -> usize {
         self.wavelength_id
     }
@@ -623,21 +719,30 @@ impl ParaxialSubView {
 
     fn calc_aperture_stop(
         surfaces: &[Box<dyn Surface>],
-        placements: &[Placement],
+        placements: &[SurfacePlacement],
+        surface_indices: &[usize],
+        path_steps: &[CursorPlacement],
         pseudo_marginal_ray: &ParaxialRayBundle,
         per_surf_v: &[TangentialVector],
     ) -> usize {
-        calc_aperture_stop(surfaces, placements, pseudo_marginal_ray, per_surf_v)
+        calc_aperture_stop(
+            surfaces,
+            placements,
+            surface_indices,
+            path_steps,
+            pseudo_marginal_ray,
+            per_surf_v,
+        )
     }
 
     fn calc_back_focal_distance(
         surfaces: &[Box<dyn Surface>],
+        surface_indices: &[usize],
         parallel_ray: &ParaxialRayBundle,
     ) -> Result<Float> {
-        let last_physical_surface_index =
-            last_physical_surface(surfaces).ok_or(anyhow!("There are no physical surfaces"))?;
-        let intercepts =
-            axis_intercepts(parallel_ray.rays_at_surface(last_physical_surface_index))?;
+        let last_physical_step_index = last_physical_step(surface_indices, surfaces)
+            .ok_or(anyhow!("There are no physical surfaces"))?;
+        let intercepts = axis_intercepts(parallel_ray.rays_at_surface(last_physical_step_index))?;
 
         let bfd = intercepts[0];
 
@@ -672,10 +777,14 @@ impl ParaxialSubView {
     /// Only field specs whose phi angle matches `v` are used. This ensures each
     /// submodel's chief ray is computed from the fields that lie in its
     /// meridional plane.
+    #[allow(clippy::too_many_arguments)]
     fn calc_chief_ray(
         surfaces: &[Box<dyn Surface>],
         sequential_sub_model: &dyn SequentialSubModel,
-        placements: &[Placement],
+        placements: &[SurfacePlacement],
+        surface_indices: &[usize],
+        beam_splitter_arms: &[Option<BeamSplitterPathKind>],
+        path_steps: &[CursorPlacement],
         v: TangentialVector,
         field_specs: &[FieldSpec],
         entrance_pupil: &Pupil,
@@ -716,6 +825,9 @@ impl ParaxialSubView {
             sequential_sub_model,
             surfaces,
             placements,
+            surface_indices,
+            beam_splitter_arms,
+            path_steps,
             false,
         )
     }
@@ -737,10 +849,14 @@ impl ParaxialSubView {
         efl.abs()
     }
 
+    #[allow(clippy::too_many_arguments)]
     fn calc_entrance_pupil(
         sequential_sub_model: &dyn SequentialSubModel,
         surfaces: &[Box<dyn Surface>],
-        placements: &[Placement],
+        placements: &[SurfacePlacement],
+        surface_indices: &[usize],
+        beam_splitter_arms: &[Option<BeamSplitterPathKind>],
+        path_steps: &[CursorPlacement],
         is_obj_space_telecentric: bool,
         aperture_stop: &usize,
         per_surf_v: &[TangentialVector],
@@ -754,12 +870,17 @@ impl ParaxialSubView {
             });
         }
 
-        // In case the aperture stop is the first surface.
+        // In case the aperture stop is the first surface (step 1).
         if *aperture_stop == 1usize {
+            let store_idx = surface_indices[1];
+            let crm = path_steps[1].cursor_rotation_matrix;
             return Ok(Pupil {
                 location: 0.0,
-                semi_diameter: placements[1]
-                    .projected_semi_diameter(surfaces[1].mask().semi_diameter(), per_surf_v[1]),
+                semi_diameter: placements[store_idx].projected_semi_diameter(
+                    crm,
+                    surfaces[store_idx].mask().semi_diameter(),
+                    per_surf_v[1],
+                ),
             });
         }
 
@@ -771,9 +892,12 @@ impl ParaxialSubView {
         }];
         let results = Self::trace(
             ray,
-            &sequential_sub_model.slice(0..*aperture_stop),
+            &sequential_sub_model.slice(0..*aperture_stop, surface_indices, beam_splitter_arms),
             surfaces,
             placements,
+            surface_indices,
+            beam_splitter_arms,
+            path_steps,
             true,
         )?;
         let location = axis_intercepts(results.last_surface().unwrap())?[0];
@@ -798,19 +922,24 @@ impl ParaxialSubView {
         })
     }
 
+    #[allow(clippy::too_many_arguments)]
     fn calc_exit_pupil(
         sequential_sub_model: &dyn SequentialSubModel,
         surfaces: &[Box<dyn Surface>],
-        placements: &[Placement],
+        placements: &[SurfacePlacement],
+        surface_indices: &[usize],
+        beam_splitter_arms: &[Option<BeamSplitterPathKind>],
+        path_steps: &[CursorPlacement],
         aperture_stop: &usize,
         marginal_ray: &ParaxialRayBundle,
     ) -> Result<Pupil> {
-        let last_physical_surface_id =
-            last_physical_surface(surfaces).ok_or(anyhow!("There are no physical surfaces"))?;
-        if last_physical_surface_id == *aperture_stop {
+        let last_physical_step_id = last_physical_step(surface_indices, surfaces)
+            .ok_or(anyhow!("There are no physical surfaces"))?;
+        if last_physical_step_id == *aperture_stop {
+            let store_idx = surface_indices[last_physical_step_id];
             return Ok(Pupil {
                 location: 0.0,
-                semi_diameter: surfaces[last_physical_surface_id].mask().semi_diameter(),
+                semi_diameter: surfaces[store_idx].mask().semi_diameter(),
             });
         }
 
@@ -822,20 +951,26 @@ impl ParaxialSubView {
 
         let results = Self::trace(
             ray,
-            &sequential_sub_model.slice(*aperture_stop..sequential_sub_model.len()),
+            &sequential_sub_model.slice(
+                *aperture_stop..sequential_sub_model.len(),
+                surface_indices,
+                beam_splitter_arms,
+            ),
             surfaces,
             placements,
+            surface_indices,
+            beam_splitter_arms,
+            path_steps,
             false,
         )?;
 
-        // Distance is relative to the last physical surface
-        let sliced_last_physical_surface_id = last_physical_surface_id - aperture_stop;
-        let distance =
-            axis_intercepts(results.rays_at_surface(sliced_last_physical_surface_id))?[0];
+        // Distance is relative to the last physical surface (step-indexed in the slice)
+        let sliced_last_physical_step_id = last_physical_step_id - aperture_stop;
+        let distance = axis_intercepts(results.rays_at_surface(sliced_last_physical_step_id))?[0];
 
         // Propagate the marginal ray to the exit pupil location and find its height
         let semi_diameter = propagate(
-            marginal_ray.rays_at_surface(last_physical_surface_id),
+            marginal_ray.rays_at_surface(last_physical_step_id),
             distance,
         )[0]
         .height;
@@ -847,12 +982,14 @@ impl ParaxialSubView {
     }
 
     fn calc_front_focal_distance(
+        sequential_sub_model: &dyn SequentialSubModel,
         surfaces: &[Box<dyn Surface>],
+        surface_indices: &[usize],
         reverse_parallel_ray: &ParaxialRayBundle,
     ) -> Result<Float> {
-        let first_physical_surface_index =
-            first_physical_surface(surfaces).ok_or(anyhow!("There are no physical surfaces"))?;
-        let index = reversed_surface_id(surfaces.len(), first_physical_surface_index);
+        let first_physical_step_index = first_physical_step(surface_indices, surfaces)
+            .ok_or(anyhow!("There are no physical surfaces"))?;
+        let index = reversed_surface_id(sequential_sub_model.len() + 1, first_physical_step_index);
         let intercepts = axis_intercepts(reverse_parallel_ray.rays_at_surface(index))?;
 
         let ffd = intercepts[0];
@@ -880,7 +1017,9 @@ impl ParaxialSubView {
 
     fn calc_marginal_ray(
         surfaces: &[Box<dyn Surface>],
-        placements: &[Placement],
+        placements: &[SurfacePlacement],
+        surface_indices: &[usize],
+        path_steps: &[CursorPlacement],
         pseudo_marginal_ray: &ParaxialRayBundle,
         aperture_stop: &usize,
         per_surf_v: &[TangentialVector],
@@ -888,6 +1027,8 @@ impl ParaxialSubView {
         calc_marginal_ray(
             surfaces,
             placements,
+            surface_indices,
+            path_steps,
             pseudo_marginal_ray,
             aperture_stop,
             per_surf_v,
@@ -898,37 +1039,51 @@ impl ParaxialSubView {
     fn calc_parallel_ray(
         sequential_sub_model: &dyn SequentialSubModel,
         surfaces: &[Box<dyn Surface>],
-        placements: &[Placement],
+        placements: &[SurfacePlacement],
+        surface_indices: &[usize],
+        beam_splitter_arms: &[Option<BeamSplitterPathKind>],
+        path_steps: &[CursorPlacement],
     ) -> Result<ParaxialRayBundle> {
         let ray = vec![ParaxialRay {
             height: 1.0,
             angle: 0.0,
         }];
 
-        Self::trace(ray, sequential_sub_model, surfaces, placements, false)
+        Self::trace(
+            ray,
+            sequential_sub_model,
+            surfaces,
+            placements,
+            surface_indices,
+            beam_splitter_arms,
+            path_steps,
+            false,
+        )
     }
 
     /// Compute the paraxial image plane.
     fn calc_paraxial_image_plane(
         surfaces: &[Box<dyn Surface>],
-        placements: &[Placement],
+        placements: &[SurfacePlacement],
+        surface_indices: &[usize],
         marginal_ray: &ParaxialRayBundle,
         chief_ray: &ParaxialRayBundle,
     ) -> Result<ImagePlane> {
-        let last_physical_surface_id =
-            last_physical_surface(surfaces).ok_or(anyhow!("There are no physical surfaces"))?;
+        let last_physical_step_id = last_physical_step(surface_indices, surfaces)
+            .ok_or(anyhow!("There are no physical surfaces"))?;
+        let store_idx = surface_indices[last_physical_step_id];
 
-        let d_axis = axis_intercepts(marginal_ray.rays_at_surface(last_physical_surface_id))?[0];
+        let d_axis = axis_intercepts(marginal_ray.rays_at_surface(last_physical_step_id))?[0];
         let location = if d_axis.is_infinite() {
             // Ensure positive infinity is returned for infinite image planes
             Float::INFINITY
         } else {
-            placements[last_physical_surface_id].track + d_axis
+            placements[store_idx].track + d_axis
         };
 
         // Propagate the chief ray from the last physical surface to the image plane to
         // determine its semi-diameter.
-        let propagated = propagate(chief_ray.rays_at_surface(last_physical_surface_id), d_axis);
+        let propagated = propagate(chief_ray.rays_at_surface(last_physical_step_id), d_axis);
         let semi_diameter = propagated[0].height.abs();
 
         Ok(ImagePlane {
@@ -941,30 +1096,55 @@ impl ParaxialSubView {
     fn calc_pseudo_marginal_ray(
         sequential_sub_model: &dyn SequentialSubModel,
         surfaces: &[Box<dyn Surface>],
-        placements: &[Placement],
+        placements: &[SurfacePlacement],
+        surface_indices: &[usize],
+        beam_splitter_arms: &[Option<BeamSplitterPathKind>],
+        path_steps: &[CursorPlacement],
     ) -> Result<ParaxialRayBundle> {
-        calc_pseudo_marginal_ray(sequential_sub_model, surfaces, placements)
+        calc_pseudo_marginal_ray(
+            sequential_sub_model,
+            surfaces,
+            placements,
+            surface_indices,
+            beam_splitter_arms,
+            path_steps,
+        )
     }
 
     /// Compute the reverse parallel ray.
     fn calc_reverse_parallel_ray(
         sequential_sub_model: &dyn SequentialSubModel,
         surfaces: &[Box<dyn Surface>],
-        placements: &[Placement],
+        placements: &[SurfacePlacement],
+        surface_indices: &[usize],
+        beam_splitter_arms: &[Option<BeamSplitterPathKind>],
+        path_steps: &[CursorPlacement],
     ) -> Result<ParaxialRayBundle> {
         let ray = vec![ParaxialRay {
             height: 1.0,
             angle: 0.0,
         }];
 
-        Self::trace(ray, sequential_sub_model, surfaces, placements, true)
+        Self::trace(
+            ray,
+            sequential_sub_model,
+            surfaces,
+            placements,
+            surface_indices,
+            beam_splitter_arms,
+            path_steps,
+            true,
+        )
     }
 
     /// Compute the ray transfer matrix for each gap/surface pair.
     fn rtms(
         sequential_sub_model: &dyn SequentialSubModel,
         surfaces: &[Box<dyn Surface>],
-        placements: &[Placement],
+        placements: &[SurfacePlacement],
+        surface_indices: &[usize],
+        beam_splitter_arms: &[Option<BeamSplitterPathKind>],
+        path_steps: &[CursorPlacement],
         reverse: bool,
     ) -> Result<Vec<RayTransferMatrix>> {
         let mut txs: Vec<RayTransferMatrix> = Vec::new();
@@ -972,11 +1152,23 @@ impl ParaxialSubView {
         let mut reverse_iter;
         let steps: &mut dyn Iterator<Item = Step> = if reverse {
             reverse_iter = sequential_sub_model
-                .try_iter(surfaces, placements)?
+                .try_iter(
+                    surfaces,
+                    placements,
+                    surface_indices,
+                    beam_splitter_arms,
+                    path_steps,
+                )?
                 .try_reverse()?;
             &mut reverse_iter
         } else {
-            forward_iter = sequential_sub_model.try_iter(surfaces, placements)?;
+            forward_iter = sequential_sub_model.try_iter(
+                surfaces,
+                placements,
+                surface_indices,
+                beam_splitter_arms,
+                path_steps,
+            )?;
             &mut forward_iter
         };
         for Step {
@@ -996,8 +1188,6 @@ impl ParaxialSubView {
                 gap_0.thickness
             };
 
-            let roc = surface.roc(0.0);
-
             let n_0 = gap_0.refractive_index.n();
             let n_1 = if let Some(gap_1) = gap_1 {
                 gap_1.refractive_index.n()
@@ -1005,21 +1195,33 @@ impl ParaxialSubView {
                 gap_0.refractive_index.n()
             };
 
-            let rtm = surface_to_rtm(surface, t, roc, n_0, n_1);
+            let rtm = surface_to_rtm(surface, t, n_0, n_1);
             txs.push(rtm);
         }
 
         Ok(txs)
     }
 
+    #[allow(clippy::too_many_arguments)]
     fn trace(
         initial_rays: Vec<ParaxialRay>,
         sequential_sub_model: &dyn SequentialSubModel,
         surfaces: &[Box<dyn Surface>],
-        placements: &[Placement],
+        placements: &[SurfacePlacement],
+        surface_indices: &[usize],
+        beam_splitter_arms: &[Option<BeamSplitterPathKind>],
+        path_steps: &[CursorPlacement],
         reverse: bool,
     ) -> Result<ParaxialRayBundle> {
-        let txs = Self::rtms(sequential_sub_model, surfaces, placements, reverse)?;
+        let txs = Self::rtms(
+            sequential_sub_model,
+            surfaces,
+            placements,
+            surface_indices,
+            beam_splitter_arms,
+            path_steps,
+            reverse,
+        )?;
         let num_surfaces = txs.len() + 1;
         let num_rays = initial_rays.len();
         let mut flat: Vec<ParaxialRay> = Vec::with_capacity(num_surfaces * num_rays);
@@ -1047,21 +1249,16 @@ impl ParaxialSubView {
 
 /// Compute the ray transfer matrix for propagation to and interaction with a
 /// surface.
-fn surface_to_rtm(
-    surface: &dyn Surface,
-    t: Float,
-    roc: Float,
-    n_0: Float,
-    n_1: Float,
-) -> RayTransferMatrix {
+fn surface_to_rtm(surface: &dyn Surface, t: Float, n_0: Float, n_1: Float) -> RayTransferMatrix {
     match surface.boundary_kind() {
-        BoundaryKind::Refracting => Mat2x2::new(
-            1.0,
-            t,
-            (n_0 - n_1) / n_1 / roc,
-            t * (n_0 - n_1) / n_1 / roc + n_0 / n_1,
-        ),
-        BoundaryKind::Reflecting => Mat2x2::new(1.0, t, 2.0 / roc, 2.0 * t / roc + 1.0),
+        BoundaryKind::Refracting => {
+            let phi = surface.power(0.0, n_0, n_1);
+            Mat2x2::new(1.0, t, -phi / n_1, -phi * t / n_1 + n_0 / n_1)
+        }
+        BoundaryKind::Reflecting => {
+            let roc = surface.roc(0.0);
+            Mat2x2::new(1.0, t, 2.0 / roc, 2.0 * t / roc + 1.0)
+        }
         BoundaryKind::NoOp => Mat2x2::new(1.0, t, 0.0, 1.0),
     }
 }
@@ -1070,7 +1267,10 @@ fn surface_to_rtm(
 pub(crate) fn calc_pseudo_marginal_ray(
     sequential_sub_model: &dyn SequentialSubModel,
     surfaces: &[Box<dyn Surface>],
-    placements: &[Placement],
+    placements: &[SurfacePlacement],
+    surface_indices: &[usize],
+    beam_splitter_arms: &[Option<BeamSplitterPathKind>],
+    path_steps: &[CursorPlacement],
 ) -> Result<ParaxialRayBundle> {
     let ray = if sequential_sub_model.is_obj_at_inf() {
         vec![ParaxialRay {
@@ -1083,44 +1283,70 @@ pub(crate) fn calc_pseudo_marginal_ray(
             angle: 1.0,
         }]
     };
-    ParaxialSubView::trace(ray, sequential_sub_model, surfaces, placements, false)
+    ParaxialSubView::trace(
+        ray,
+        sequential_sub_model,
+        surfaces,
+        placements,
+        surface_indices,
+        beam_splitter_arms,
+        path_steps,
+        false,
+    )
 }
 
-/// Compute the aperture stop surface index using the minimum aperture-ratio
+/// Compute the aperture stop step index using the minimum aperture-ratio
 /// heuristic.
+///
+/// Returns a step index (position within the path's traversal), not a store
+/// index. `per_surf_v` must be step-indexed, as returned by
+/// `propagate_tangential_vec` with the path's `surface_indices`.
 pub(crate) fn calc_aperture_stop(
     surfaces: &[Box<dyn Surface>],
-    placements: &[Placement],
+    placements: &[SurfacePlacement],
+    surface_indices: &[usize],
+    path_steps: &[CursorPlacement],
     pseudo_marginal_ray: &ParaxialRayBundle,
     per_surf_v: &[TangentialVector],
 ) -> usize {
-    let ratios: Vec<Float> = surfaces
+    let ratios: Vec<Float> = surface_indices
         .iter()
-        .zip(placements.iter())
+        .zip(path_steps.iter())
         .zip(pseudo_marginal_ray.iter_surfaces())
         .zip(per_surf_v.iter())
-        .map(|(((s, p), rays), &v)| {
-            (p.projected_semi_diameter(s.mask().semi_diameter(), v) / rays[0].height).abs()
+        .map(|(((idx, step), rays), &v)| {
+            let s = &surfaces[*idx];
+            let p = &placements[*idx];
+            let crm = step.cursor_rotation_matrix;
+            (p.projected_semi_diameter(crm, s.mask().semi_diameter(), v) / rays[0].height).abs()
         })
         .collect();
     argmin(&ratios[1..ratios.len() - 1]) + 1
 }
 
 /// Scale the pseudo-marginal ray to match the aperture stop semi-diameter.
+///
+/// `aperture_stop` is a step index. `per_surf_v` must be step-indexed, as
+/// returned by `propagate_tangential_vec` with the path's `surface_indices`.
 pub(crate) fn calc_marginal_ray(
     surfaces: &[Box<dyn Surface>],
-    placements: &[Placement],
+    placements: &[SurfacePlacement],
+    surface_indices: &[usize],
+    path_steps: &[CursorPlacement],
     pseudo_marginal_ray: &ParaxialRayBundle,
     aperture_stop: &usize,
     per_surf_v: &[TangentialVector],
 ) -> ParaxialRayBundle {
-    let ratios: Vec<Float> = surfaces
+    let ratios: Vec<Float> = surface_indices
         .iter()
-        .zip(placements.iter())
+        .zip(path_steps.iter())
         .zip(pseudo_marginal_ray.iter_surfaces())
         .zip(per_surf_v.iter())
-        .map(|(((s, p), rays), &v)| {
-            p.projected_semi_diameter(s.mask().semi_diameter(), v) / rays[0].height
+        .map(|(((idx, step), rays), &v)| {
+            let s = &surfaces[*idx];
+            let p = &placements[*idx];
+            let crm = step.cursor_rotation_matrix;
+            p.projected_semi_diameter(crm, s.mask().semi_diameter(), v) / rays[0].height
         })
         .collect();
     let scale_factor = ratios[*aperture_stop];
@@ -1153,17 +1379,36 @@ pub(crate) fn marginal_ray_bundle(
         .ok_or_else(|| anyhow!("wavelength_id {wavelength_id} out of range"))?;
     let surfaces = model.surfaces();
     let placements = model.placements();
+    let path_steps = model.path_steps(0);
+    let surface_indices = model.path_surface_indices(0);
+    let beam_splitter_arms = model.path_beam_splitter_arms(0);
 
     let v = Vec3::new(0.0, 1.0, 0.0);
-    let per_surf_v = propagate_tangential_vec(v, surfaces, placements);
-    let pseudo = calc_pseudo_marginal_ray(submodel, surfaces, placements)?;
+    let per_surf_v = propagate_tangential_vec(v, surfaces, placements, surface_indices);
+    let pseudo = calc_pseudo_marginal_ray(
+        submodel,
+        surfaces,
+        placements,
+        surface_indices,
+        beam_splitter_arms,
+        path_steps,
+    )?;
     let stop = match model.stop_surface() {
         Some(i) => i,
-        None => calc_aperture_stop(surfaces, placements, &pseudo, &per_surf_v),
+        None => calc_aperture_stop(
+            surfaces,
+            placements,
+            surface_indices,
+            path_steps,
+            &pseudo,
+            &per_surf_v,
+        ),
     };
     Ok(calc_marginal_ray(
         surfaces,
         placements,
+        surface_indices,
+        path_steps,
         &pseudo,
         &stop,
         &per_surf_v,
@@ -1296,11 +1541,15 @@ mod test {
             sequential_sub_model: seq_sub_model as &dyn SequentialSubModel,
             surfaces: sequential_model.surfaces(),
             placements: sequential_model.placements(),
+            surface_indices: sequential_model.path_surface_indices(0),
+            beam_splitter_arms: sequential_model.path_beam_splitter_arms(0),
+            path_steps: sequential_model.path_steps(0),
             field_specs: &field_specs,
             stop_surface: None,
         };
         (
             ParaxialSubView::new(
+                0, // path_id
                 0,
                 0,
                 &data,
@@ -1370,6 +1619,9 @@ mod test {
             seq_sub_model,
             sequential_model.surfaces(),
             sequential_model.placements(),
+            sequential_model.path_surface_indices(0),
+            sequential_model.path_beam_splitter_arms(0),
+            sequential_model.path_steps(0),
         )
         .unwrap();
 
@@ -1400,6 +1652,9 @@ mod test {
             seq_sub_model,
             sequential_model.surfaces(),
             sequential_model.placements(),
+            sequential_model.path_surface_indices(0),
+            sequential_model.path_beam_splitter_arms(0),
+            sequential_model.path_steps(0),
         )
         .unwrap();
 
@@ -1484,11 +1739,14 @@ mod test {
             sequential_sub_model: seq_sub_model as &dyn SequentialSubModel,
             surfaces: sequential_model.surfaces(),
             placements: sequential_model.placements(),
+            surface_indices: sequential_model.path_surface_indices(0),
+            beam_splitter_arms: sequential_model.path_beam_splitter_arms(0),
+            path_steps: sequential_model.path_steps(0),
             field_specs: &field_specs,
             stop_surface: None,
         };
 
-        let view = ParaxialSubView::new(0, 0, &data, Vec3::new(0.0, 1.0, 0.0), false).unwrap();
+        let view = ParaxialSubView::new(0, 0, 0, &data, Vec3::new(0.0, 1.0, 0.0), false).unwrap();
 
         assert_eq!(*view.aperture_stop(), 2);
     }
@@ -1567,4 +1825,125 @@ mod test {
         let expected = *sub.effective_focal_length() / (2.0 * sub.entrance_pupil().semi_diameter);
         assert_abs_diff_eq!(sub.image_space_fno(), expected, epsilon = 1e-6);
     }
+
+    #[test]
+    fn paraxial_view_two_path_model_has_subviews_for_both_paths() {
+        use std::rc::Rc;
+
+        use crate::examples::beam_splitter::two_path_model;
+        use crate::specs::gaps::ConstantRefractiveIndex;
+
+        let n_air = Rc::new(ConstantRefractiveIndex::new(1.0, 0.0));
+        let model = two_path_model(n_air, &[0.5876e-3], 10.0, 10.0);
+        let field = vec![FieldSpec::Angle {
+            chi: 0.0,
+            phi: 90.0,
+        }];
+        let pv = ParaxialView::new(&model, &field, false).unwrap();
+
+        assert!(pv.get_for_path(0, 0, 0).is_some(), "path 0 subview missing");
+        assert!(pv.get_for_path(1, 0, 0).is_some(), "path 1 subview missing");
+    }
+
+    /// The aperture stop for each path must be computed from that path's own
+    /// surfaces, not the store-indexed surface list. This test uses two paths
+    /// with different irises in each arm: path 0 (Iris SD=5) and path 1
+    /// (Iris SD=15), with a shared BS (SD=10). The BS is the most constraining
+    /// surface on path 1, so its aperture stop must be at step 1 (BS), not
+    /// step 2.
+    #[test]
+    fn paraxial_aperture_stop_uses_path_specific_surfaces() {
+        use std::rc::Rc;
+
+        use crate::{
+            BeamSplitterPathKind, EulerAngles, GapSpec, PathSpec, PathSurfaceRef, Rotation3D,
+            SurfaceSpec, Vec3, core::sequential_model::builder::SequentialModelBuilder,
+            specs::gaps::ConstantRefractiveIndex,
+        };
+
+        let n_air = Rc::new(ConstantRefractiveIndex::new(1.0, 0.0));
+        let bs_rotation =
+            Rotation3D::IntrinsicPassiveRUF(EulerAngles((-45_f64 as Float).to_radians(), 0.0, 0.0));
+        let gap_inf = || GapSpec {
+            thickness: Float::INFINITY,
+            refractive_index: n_air.clone(),
+        };
+        let gap_50 = || GapSpec {
+            thickness: 50.0,
+            refractive_index: n_air.clone(),
+        };
+        let img = || SurfaceSpec::Image {
+            rotation: Rotation3D::None,
+            decenter: Vec3::new(0.0, 0.0, 0.0),
+            rotation_offset: Rotation3D::None,
+        };
+        let iris = |sd: Float| SurfaceSpec::Iris {
+            semi_diameter: sd,
+            rotation: Rotation3D::None,
+            decenter: Vec3::new(0.0, 0.0, 0.0),
+            rotation_offset: Rotation3D::None,
+        };
+
+        // Path 0: Object → BS(SD=10) → Iris1(SD=5) → Image_0
+        // Store: Object(0), BS(1), Iris1(2), Image_0(3)
+        let path_t = PathSpec {
+            surface_refs: vec![
+                PathSurfaceRef::New(SurfaceSpec::Object),
+                PathSurfaceRef::New(SurfaceSpec::BeamSplitter {
+                    semi_diameter: 10.0,
+                    rotation: bs_rotation,
+                    decenter: Vec3::new(0.0, 0.0, 0.0),
+                    rotation_offset: Rotation3D::None,
+                }),
+                PathSurfaceRef::New(iris(5.0)),
+                PathSurfaceRef::New(img()),
+            ],
+            gaps: vec![gap_inf(), gap_50(), gap_50()],
+            beam_splitter_arms: vec![BeamSplitterPathKind::Transmitting],
+        };
+
+        // Path 1: Object(Shared) → BS(Shared) → Iris2(SD=15) → Image_1
+        // Store: ..., Iris2(4), Image_1(5)
+        let path_r = PathSpec {
+            surface_refs: vec![
+                PathSurfaceRef::Shared(0),
+                PathSurfaceRef::Shared(1),
+                PathSurfaceRef::New(iris(15.0)),
+                PathSurfaceRef::New(img()),
+            ],
+            gaps: vec![gap_inf(), gap_50(), gap_50()],
+            beam_splitter_arms: vec![BeamSplitterPathKind::Reflecting],
+        };
+
+        let model = SequentialModelBuilder::new()
+            .paths(vec![path_t, path_r])
+            .wavelengths(vec![0.5876e-3])
+            .build()
+            .unwrap()
+            .model;
+
+        let field = vec![FieldSpec::Angle {
+            chi: 0.0,
+            phi: 90.0,
+        }];
+        let pv = ParaxialView::new(&model, &field, false).unwrap();
+
+        // Path 0: BS(SD=10) at step 1, Iris1(SD=5) at step 2. Iris1 is the stop.
+        let path0 = pv.get_for_path(0, 0, 0).expect("path 0 subview");
+        assert_eq!(
+            *path0.aperture_stop(),
+            2,
+            "path 0 stop should be Iris1 (step 2)"
+        );
+
+        // Path 1: BS(SD=10) at step 1, Iris2(SD=15) at step 2. BS is the stop.
+        // Bug: if calc_aperture_stop uses store-indexed surfaces, it sees Iris1
+        // (SD=5) at step 2 position and incorrectly identifies step 2 as the stop.
+        let path1 = pv.get_for_path(1, 0, 0).expect("path 1 subview");
+        assert_eq!(
+            *path1.aperture_stop(),
+            1,
+            "path 1 stop should be BS (step 1)"
+        );
+    }
 }
diff --git a/crates/cherry-rs/src/views/ray_trace_3d/mod.rs b/crates/cherry-rs/src/views/ray_trace_3d/mod.rs
index 7adc28d9..57aed369 100644
--- a/crates/cherry-rs/src/views/ray_trace_3d/mod.rs
+++ b/crates/cherry-rs/src/views/ray_trace_3d/mod.rs
@@ -13,12 +13,16 @@ use crate::{
         Float, PI,
         math::vec3::Vec3,
         ray::Ray,
-        sequential_model::{SequentialModel, SequentialSubModel, placement::Placement},
+        sequential_model::{
+            CursorPlacement, SequentialModel, SequentialSubModel,
+            surface_placement::SurfacePlacement,
+        },
         surfaces::Surface,
     },
     specs::{
         aperture::ApertureSpec,
         fields::{FieldSpec, PupilSampling},
+        surfaces::BeamSplitterPathKind,
     },
 };
 
@@ -62,11 +66,13 @@ pub struct TraceResultsCollection {
 /// The results of a 3D ray trace.
 ///
 /// This represents the results of a 3D ray trace for a single set of values of
-/// 1. wavelength ID and
-/// 2. field ID.
+/// 1. path ID,
+/// 2. wavelength ID, and
+/// 3. field ID.
 #[derive(Debug)]
 #[cfg_attr(feature = "serde", derive(Serialize))]
 pub struct TraceResults {
+    path_id: usize,
     wavelength_id: usize,
     field_id: usize,
 
@@ -133,6 +139,9 @@ pub fn trace_ray_bundle(
                     sequential_submodel,
                     sequential_model.surfaces(),
                     sequential_model.placements(),
+                    sequential_model.path_surface_indices(0),
+                    sequential_model.path_beam_splitter_arms(0),
+                    sequential_model.path_steps(0),
                     aperture_spec,
                     &field_specs[field_id],
                     paraxial_subview,
@@ -168,92 +177,116 @@ pub fn ray_trace_3d_view(
     // For rotationally symmetric systems only Axis::U exists; for non-symmetric
     // systems this is the fallback that is always present.
     let n_wavelengths = sequential_model.wavelengths().len();
-    let pairs: Vec<(usize, usize)> = (0..field_specs.len())
-        .flat_map(|f| (0..n_wavelengths).map(move |w| (f, w)))
+    let triples: Vec<(usize, usize, usize)> = (0..sequential_model.path_count())
+        .flat_map(|p| {
+            (0..field_specs.len()).flat_map(move |f| (0..n_wavelengths).map(move |w| (p, f, w)))
+        })
         .collect();
 
-    let results: Vec<TraceResults> = pairs
+    let results: Vec<TraceResults> = triples
         .into_par_iter()
-        .map(|(field_id, wavelength_id)| -> Result<TraceResults> {
-            tracing::trace!(
-                "Tracing rays for field_id={}, wavelength_id={}",
-                field_id,
-                wavelength_id,
-            );
+        .map(
+            |(path_id, field_id, wavelength_id)| -> Result<TraceResults> {
+                tracing::trace!(
+                    "Tracing rays for path_id={}, field_id={}, wavelength_id={}",
+                    path_id,
+                    field_id,
+                    wavelength_id,
+                );
 
-            let sequential_submodel = sequential_model
-                .submodel(wavelength_id)
-                .ok_or_else(|| anyhow!("Submodel not found"))?;
-            let tangential_vec_id =
-                paraxial_view.tangential_vec_id_for_phi(field_specs[field_id].tangential_fan_phi());
-            let paraxial_subview = paraxial_view
-                .get(wavelength_id, tangential_vec_id)
-                .ok_or_else(|| anyhow!("Submodel not found"))?;
-
-            let field_spec = &field_specs[field_id];
-            let surfaces = sequential_model.surfaces();
-            let placements = sequential_model.placements();
-
-            let chief_ray = ray_trace_submodel(
-                sequential_submodel,
-                surfaces,
-                placements,
-                aperture_spec,
-                field_spec,
-                paraxial_subview,
-                PupilSampling::ChiefRay,
-            )?;
-            let full_pupil = ray_trace_submodel(
-                sequential_submodel,
-                surfaces,
-                placements,
-                aperture_spec,
-                field_spec,
-                paraxial_subview,
-                PupilSampling::SquareGrid {
-                    spacing: config.full_pupil_spacing,
-                },
-            )?;
-            let tangential_fan = ray_trace_submodel(
-                sequential_submodel,
-                surfaces,
-                placements,
-                aperture_spec,
-                field_spec,
-                paraxial_subview,
-                PupilSampling::TangentialRayFan {
-                    n: config.n_fan_rays,
-                },
-            )?;
-            let sagittal_fan = ray_trace_submodel(
-                sequential_submodel,
-                surfaces,
-                placements,
-                aperture_spec,
-                field_spec,
-                paraxial_subview,
-                PupilSampling::SagittalRayFan {
-                    n: config.n_fan_rays,
-                },
-            )?;
-
-            trace!(
-                field_id,
-                wavelength_id,
-                "Finished tracing all bundles for field {}, wavelength {}",
-                field_id,
-                wavelength_id
-            );
+                let sequential_submodel = sequential_model
+                    .submodels_for_path(path_id)
+                    .get(wavelength_id)
+                    .ok_or_else(|| anyhow!("Submodel not found"))?;
+                let tangential_vec_id = paraxial_view
+                    .tangential_vec_id_for_phi(field_specs[field_id].tangential_fan_phi());
+                let paraxial_subview = paraxial_view
+                    .get_for_path(path_id, wavelength_id, tangential_vec_id)
+                    .ok_or_else(|| anyhow!("Paraxial subview not found"))?;
 
-            Ok(TraceResults {
-                wavelength_id,
-                field_id,
-                chief_ray,
-                full_pupil,
-                tangential_fan,
-                sagittal_fan,
-            })
-        })
+                let field_spec = &field_specs[field_id];
+                let surfaces = sequential_model.surfaces();
+                let placements = sequential_model.placements();
+                let surface_indices = sequential_model.path_surface_indices(path_id);
+                let beam_splitter_arms = sequential_model.path_beam_splitter_arms(path_id);
+                let path_steps = sequential_model.path_steps(path_id);
+
+                let chief_ray = ray_trace_submodel(
+                    sequential_submodel,
+                    surfaces,
+                    placements,
+                    surface_indices,
+                    beam_splitter_arms,
+                    path_steps,
+                    aperture_spec,
+                    field_spec,
+                    paraxial_subview,
+                    PupilSampling::ChiefRay,
+                )?;
+                let full_pupil = ray_trace_submodel(
+                    sequential_submodel,
+                    surfaces,
+                    placements,
+                    surface_indices,
+                    beam_splitter_arms,
+                    path_steps,
+                    aperture_spec,
+                    field_spec,
+                    paraxial_subview,
+                    PupilSampling::SquareGrid {
+                        spacing: config.full_pupil_spacing,
+                    },
+                )?;
+                let tangential_fan = ray_trace_submodel(
+                    sequential_submodel,
+                    surfaces,
+                    placements,
+                    surface_indices,
+                    beam_splitter_arms,
+                    path_steps,
+                    aperture_spec,
+                    field_spec,
+                    paraxial_subview,
+                    PupilSampling::TangentialRayFan {
+                        n: config.n_fan_rays,
+                    },
+                )?;
+                let sagittal_fan = ray_trace_submodel(
+                    sequential_submodel,
+                    surfaces,
+                    placements,
+                    surface_indices,
+                    beam_splitter_arms,
+                    path_steps,
+                    aperture_spec,
+                    field_spec,
+                    paraxial_subview,
+                    PupilSampling::SagittalRayFan {
+                        n: config.n_fan_rays,
+                    },
+                )?;
+
+                trace!(
+                    path_id,
+                    field_id,
+                    wavelength_id,
+                    "Finished tracing all bundles for path {}, field {}, wavelength {}",
+                    path_id,
+                    field_id,
+                    wavelength_id
+                );
+
+                Ok(TraceResults {
+                    path_id,
+                    wavelength_id,
+                    field_id,
+                    chief_ray,
+                    full_pupil,
+                    tangential_fan,
+                    sagittal_fan,
+                })
+            },
+        )
         .collect::<Result<Vec<_>>>()?;
 
     Ok(TraceResultsCollection::new(results))
@@ -264,11 +297,21 @@ impl TraceResultsCollection {
         Self { results }
     }
 
-    /// Get results for a specific field and wavelength.
+    /// Get results for path 0, a specific field, and wavelength.
     pub fn get(&self, field_id: usize, wavelength_id: usize) -> Option<&TraceResults> {
-        self.results
-            .iter()
-            .find(|r| r.field_id == field_id && r.wavelength_id == wavelength_id)
+        self.get_for_path(0, field_id, wavelength_id)
+    }
+
+    /// Get results for a specific path, field, and wavelength.
+    pub fn get_for_path(
+        &self,
+        path_id: usize,
+        field_id: usize,
+        wavelength_id: usize,
+    ) -> Option<&TraceResults> {
+        self.results.iter().find(|r| {
+            r.path_id == path_id && r.field_id == field_id && r.wavelength_id == wavelength_id
+        })
     }
 
     /// Get all results for a given wavelength.
@@ -309,6 +352,11 @@ impl TraceResultsCollection {
 }
 
 impl TraceResults {
+    /// Returns the path ID of this result.
+    pub fn path_id(&self) -> usize {
+        self.path_id
+    }
+
     // Returns the field ID of the ray bundle.
     pub fn field_id(&self) -> usize {
         self.field_id
@@ -346,10 +394,14 @@ impl TraceResults {
     }
 }
 
+#[allow(clippy::too_many_arguments)]
 fn ray_trace_submodel(
     sequential_submodel: &impl SequentialSubModel,
     surfaces: &[Box<dyn Surface>],
-    placements: &[Placement],
+    placements: &[SurfacePlacement],
+    surface_indices: &[usize],
+    beam_splitter_arms: &[Option<BeamSplitterPathKind>],
+    path_steps: &[CursorPlacement],
     aperture_spec: &ApertureSpec,
     field_spec: &FieldSpec,
     paraxial_subview: &ParaxialSubView,
@@ -363,7 +415,13 @@ fn ray_trace_submodel(
         pupil_sampling,
     )?;
 
-    let mut sequential_sub_model_iter = sequential_submodel.try_iter(surfaces, placements)?;
+    let mut sequential_sub_model_iter = sequential_submodel.try_iter(
+        surfaces,
+        placements,
+        surface_indices,
+        beam_splitter_arms,
+        path_steps,
+    )?;
     Ok(trace(&mut sequential_sub_model_iter, rays))
 }
 
@@ -378,7 +436,7 @@ fn ray_trace_submodel(
 /// * `field_spec` - The field specification.
 /// * `sampling` - The pupil sampling method.
 fn rays(
-    placements: &[Placement],
+    placements: &[SurfacePlacement],
     aperture_spec: &ApertureSpec,
     paraxial_subview: &ParaxialSubView,
     field_spec: &FieldSpec,
@@ -479,7 +537,7 @@ fn rays(
 /// * `phi` - The azimuthal angle of the ray in the x-y plane, radians.
 /// * `chi` - The zenith angle of the ray w.r.t. the z-axis, radians.
 fn chief_ray_from_angle(
-    placements: &[Placement],
+    placements: &[SurfacePlacement],
     aperture_spec: &ApertureSpec,
     paraxial_subview: &ParaxialSubView,
     phi: Float,
@@ -525,7 +583,7 @@ fn chief_ray_from_pos(
 /// * `chi` - The zenith angle of the ray w.r.t. the z-axis, radians.
 #[allow(clippy::too_many_arguments)]
 fn parallel_ray_fan(
-    placements: &[Placement],
+    placements: &[SurfacePlacement],
     aperture_spec: &ApertureSpec,
     paraxial_subview: &ParaxialSubView,
     num_rays: usize,
@@ -586,7 +644,7 @@ fn parallel_ray_fan(
 ///   (marginal rays).
 /// * `chi` - The zenith angle of the ray bundle w.r.t. the z-axis in radians.
 fn parallel_ray_bundle_on_sq_grid(
-    placements: &[Placement],
+    placements: &[SurfacePlacement],
     aperture_spec: &ApertureSpec,
     paraxial_subview: &ParaxialSubView,
     spacing: Float,
@@ -770,7 +828,7 @@ fn axial_launch_point(obj_z: Float, sur_z: Float, enp_z: Float) -> Float {
 /// * `phi` - The azimuthal angle of the ray fan in the x-y plane, radians.
 /// * `chi` - The zenith angle of the ray w.r.t. the z-axis, radians.
 fn parallel_ray_bundle_origin(
-    placements: &[Placement],
+    placements: &[SurfacePlacement],
     aperture_spec: &ApertureSpec,
     paraxial_subview: &ParaxialSubView,
     phi: Float,
@@ -1322,6 +1380,40 @@ mod tests {
         approx::assert_abs_diff_eq!(dy / dx, 1.0, epsilon = 1e-6);
     }
 
+    #[test]
+    fn ray_trace_two_path_model_produces_results_for_both_paths() {
+        use std::rc::Rc;
+
+        use crate::examples::beam_splitter::two_path_model;
+        use crate::specs::gaps::ConstantRefractiveIndex;
+
+        let n_air = Rc::new(ConstantRefractiveIndex::new(1.0, 0.0));
+        let model = two_path_model(n_air, &[0.5876e-3], 10.0, 10.0);
+        let field_specs = vec![FieldSpec::Angle {
+            chi: 0.0,
+            phi: 90.0,
+        }];
+        let aperture_spec = ApertureSpec::EntrancePupil { semi_diameter: 5.0 };
+        let paraxial_view = ParaxialView::new(&model, &field_specs, false).unwrap();
+        let config = SamplingConfig {
+            n_fan_rays: 3,
+            full_pupil_spacing: 0.5,
+        };
+
+        let results =
+            ray_trace_3d_view(&aperture_spec, &field_specs, &model, &paraxial_view, config)
+                .unwrap();
+
+        assert!(
+            results.get_for_path(0, 0, 0).is_some(),
+            "path 0 results missing"
+        );
+        assert!(
+            results.get_for_path(1, 0, 0).is_some(),
+            "path 1 results missing"
+        );
+    }
+
     #[test]
     fn chief_ray_reached_image_on_axis() {
         // On-axis chief ray should always reach the image surface in a standard
diff --git a/crates/cherry-rs/src/views/ray_trace_3d/trace.rs b/crates/cherry-rs/src/views/ray_trace_3d/trace.rs
index 92796235..872079d9 100644
--- a/crates/cherry-rs/src/views/ray_trace_3d/trace.rs
+++ b/crates/cherry-rs/src/views/ray_trace_3d/trace.rs
@@ -52,7 +52,7 @@ pub fn trace(sequential_submodel: &mut SequentialSubModelIter, mut rays: Vec<Ray
             }
 
             // Transform into coordinate system of the surface
-            ray.transform(step.placement);
+            ray.transform(step.surface_placement);
 
             // Find the ray intersection with the surface.
             // Errors if the intersection point does not converge.
@@ -91,7 +91,7 @@ pub fn trace(sequential_submodel: &mut SequentialSubModelIter, mut rays: Vec<Ray
             step.surface.interact(ray, n_0, n_1, norm);
 
             // Transform back to the global coordinate system
-            ray.i_transform(step.placement);
+            ray.i_transform(step.surface_placement);
 
             rays_at_surface[ray_id] = ray.clone();
         }
diff --git a/crates/cherry-rs/tests/thin_lens_singlet.rs b/crates/cherry-rs/tests/thin_lens_singlet.rs
new file mode 100644
index 00000000..a0847a26
--- /dev/null
+++ b/crates/cherry-rs/tests/thin_lens_singlet.rs
@@ -0,0 +1,39 @@
+use cherry_rs::{FieldSpec, ParaxialView, examples::thin_lens_singlet, n};
+
+const WAVELENGTHS: [f64; 1] = [0.5876];
+const FIELD_SPECS: [FieldSpec; 1] = [FieldSpec::Angle { chi: 0.0, phi: 0.0 }];
+
+// An idealized thin lens has no thickness, so its principal planes coincide
+// with the lens itself — back focal distance equals the effective focal
+// length exactly, for an object at infinity.
+
+#[test]
+fn thin_lens_efl_equals_focal_length() {
+    let model = thin_lens_singlet::sequential_model(n!(1.0), &WAVELENGTHS);
+    let view = ParaxialView::new(&model, &FIELD_SPECS, false).expect("paraxial view");
+    for sub_view in view.iter() {
+        approx::assert_abs_diff_eq!(*sub_view.effective_focal_length(), 100.0, epsilon = 1e-9);
+    }
+}
+
+#[test]
+fn thin_lens_bfd_equals_focal_length() {
+    let model = thin_lens_singlet::sequential_model(n!(1.0), &WAVELENGTHS);
+    let view = ParaxialView::new(&model, &FIELD_SPECS, false).expect("paraxial view");
+    for sub_view in view.iter() {
+        approx::assert_abs_diff_eq!(*sub_view.back_focal_distance(), 100.0, epsilon = 1e-9);
+    }
+}
+
+#[test]
+fn thin_lens_marginal_ray_crosses_axis_at_focal_plane() {
+    let model = thin_lens_singlet::sequential_model(n!(1.0), &WAVELENGTHS);
+    let view = ParaxialView::new(&model, &FIELD_SPECS, false).expect("paraxial view");
+    for sub_view in view.iter() {
+        let marginal = sub_view.marginal_ray();
+        // Surface 0 = Object (at infinity), 1 = ThinLens, 2 = Image (placed at
+        // the back focal distance by the example's gap thickness of 100 mm).
+        let h_image = marginal.rays_at_surface(2)[0].height;
+        approx::assert_abs_diff_eq!(h_image, 0.0, epsilon = 1e-9);
+    }
+}