From 7c550dceeb8257aaf7fa47bb215cd196ef2ec3ed Mon Sep 17 00:00:00 2001 From: Kyle Douglass Date: Mon, 6 Jul 2026 11:53:27 +0200 Subject: [PATCH 1/3] refactor: Enable per-path wavelength specs --- .../src/core/sequential_model/builder.rs | 139 ++++++++++++++---- .../src/core/sequential_model/mod.rs | 42 +++--- .../cherry-rs/src/examples/beam_splitter.rs | 3 +- .../cherry-rs/src/examples/wf_epi_emission.rs | 5 +- .../src/examples/wf_epi_excitation.rs | 2 +- .../src/examples/wf_epi_microscope.rs | 6 +- crates/cherry-rs/src/specs/paths.rs | 6 + crates/cherry-rs/src/views/paraxial.rs | 3 +- .../cherry-rs/src/views/ray_trace_3d/mod.rs | 97 +++++++++++- .../tests/beam_splitter_multipath.rs | 6 +- crates/cherry-rs/tests/wf_epi_microscope.rs | 55 ++++++- 11 files changed, 304 insertions(+), 60 deletions(-) diff --git a/crates/cherry-rs/src/core/sequential_model/builder.rs b/crates/cherry-rs/src/core/sequential_model/builder.rs index 7264c46c..bca6885e 100644 --- a/crates/cherry-rs/src/core/sequential_model/builder.rs +++ b/crates/cherry-rs/src/core/sequential_model/builder.rs @@ -60,12 +60,10 @@ impl SequentialModelBuilder { if self.paths.is_some() { let paths = self.paths.unwrap(); - let wavelengths = self.wavelengths.unwrap(); #[cfg(feature = "serde")] - let model = - SequentialModel::from_path_specs(paths, &wavelengths, self.registry.as_ref())?; + let model = SequentialModel::from_path_specs(paths, self.registry.as_ref())?; #[cfg(not(feature = "serde"))] - let model = SequentialModel::from_path_specs(paths, &wavelengths)?; + let model = SequentialModel::from_path_specs(paths)?; return Ok(BuildResult { model, gap_specs: vec![], @@ -166,10 +164,17 @@ impl SequentialModelBuilder { } if self.paths.is_some() { - if self.wavelengths.is_none() { - return Err(anyhow!("Wavelengths must be set")); - } else if self.wavelengths.as_ref().unwrap().is_empty() { - return Err(anyhow!("Wavelengths cannot be empty")); + if self.wavelengths.is_some() { + return Err(anyhow!( + "Cannot set both `paths` and `wavelengths` — set `wavelengths` on \ + each PathSpec instead" + )); + } + if self.stop_surface.is_some() { + return Err(anyhow!( + "Cannot set both `paths` and `stop_surface` — set `stop_surface` \ + on each PathSpec instead" + )); } return Ok(()); } @@ -605,6 +610,7 @@ mod tests { }], beam_splitter_arms: vec![], stop_surface: None, + wavelengths: vec![0.587], } } @@ -630,11 +636,9 @@ mod tests { }], beam_splitter_arms: vec![], stop_surface: None, + wavelengths: vec![0.587], }; - let result = SequentialModelBuilder::new() - .paths(vec![bad_path]) - .wavelengths(vec![0.587]) - .build(); + let result = SequentialModelBuilder::new().paths(vec![bad_path]).build(); assert!(result.is_err()); } @@ -652,11 +656,9 @@ mod tests { }], beam_splitter_arms: vec![], stop_surface: None, + wavelengths: vec![0.587], }; - let result = SequentialModelBuilder::new() - .paths(vec![bad_path]) - .wavelengths(vec![0.587]) - .build(); + let result = SequentialModelBuilder::new().paths(vec![bad_path]).build(); assert!(result.is_err()); } @@ -685,11 +687,9 @@ mod tests { ], beam_splitter_arms: vec![], stop_surface: None, + wavelengths: vec![0.587], }; - let result = SequentialModelBuilder::new() - .paths(vec![bad_path]) - .wavelengths(vec![0.587]) - .build(); + let result = SequentialModelBuilder::new().paths(vec![bad_path]).build(); assert!(result.is_err()); } @@ -712,10 +712,10 @@ mod tests { }], beam_splitter_arms: vec![], stop_surface: None, + wavelengths: vec![0.587], }; let result = SequentialModelBuilder::new() .paths(vec![path0, bad_path1]) - .wavelengths(vec![0.587]) .build(); assert!(result.is_err()); } @@ -755,11 +755,9 @@ mod tests { ], beam_splitter_arms: vec![], // missing arm declaration for the BS stop_surface: None, + wavelengths: vec![0.587], }; - let result = SequentialModelBuilder::new() - .paths(vec![path]) - .wavelengths(vec![0.587]) - .build(); + let result = SequentialModelBuilder::new().paths(vec![path]).build(); assert!(result.is_err()); } @@ -834,10 +832,10 @@ mod tests { ], beam_splitter_arms: vec![], stop_surface: None, + wavelengths: wls.to_vec(), }; let model_new = SequentialModelBuilder::new() .paths(vec![path]) - .wavelengths(wls.to_vec()) .build() .unwrap() .model; @@ -856,4 +854,93 @@ mod tests { epsilon = 1e-10 ); } + + // ── Per-path wavelengths ────────────────────────────────────────────── + + fn minimal_path_spec_with_wavelengths(wavelengths: Vec) -> PathSpec { + use crate::specs::paths::{PathSpec, PathSurfaceRef}; + PathSpec { + surface_refs: vec![ + PathSurfaceRef::New(SurfaceSpec::Object), + PathSurfaceRef::New(SurfaceSpec::Image { + rotation: Rotation3D::None, + decenter: Vec3::new(0.0, 0.0, 0.0), + rotation_offset: Rotation3D::None, + }), + ], + gaps: vec![GapSpec { + thickness: f64::INFINITY, + refractive_index: n!(1.0), + }], + beam_splitter_arms: vec![], + stop_surface: None, + wavelengths, + } + } + + #[test] + fn wavelengths_for_path_differ_across_paths() { + use crate::specs::paths::PathSurfaceRef; + + let path0 = minimal_path_spec_with_wavelengths(vec![0.488]); + // path1 shares path0's Object (store index 0) and adds its own Image, + // so it needs its own gap and a differing wavelength list/length. + let path1 = PathSpec { + surface_refs: vec![ + PathSurfaceRef::Shared(0), + PathSurfaceRef::New(SurfaceSpec::Image { + rotation: Rotation3D::None, + decenter: Vec3::new(0.0, 0.0, 0.0), + rotation_offset: Rotation3D::None, + }), + ], + gaps: vec![GapSpec { + thickness: 10.0, + refractive_index: n!(1.0), + }], + beam_splitter_arms: vec![], + stop_surface: None, + wavelengths: vec![0.500, 0.520, 0.540], + }; + + let model = SequentialModelBuilder::new() + .paths(vec![path0, path1]) + .build() + .unwrap() + .model; + + assert_eq!(model.wavelengths_for_path(0), &[0.488]); + assert_eq!(model.wavelengths_for_path(1), &[0.500, 0.520, 0.540]); + assert_eq!(model.wavelengths(), model.wavelengths_for_path(0)); + } + + #[test] + fn empty_wavelengths_on_one_path_is_rejected_even_if_others_are_nonempty() { + // Construct two independent PathSpecs directly rather than cloning — + // PathSpec derives no traits (not even Clone) as of this writing. + let good = minimal_path_spec_with_wavelengths(vec![0.5876]); + let bad = minimal_path_spec_with_wavelengths(vec![]); + let result = SequentialModelBuilder::new().paths(vec![good, bad]).build(); + assert!(result.is_err()); + } + + #[test] + fn build_fails_when_both_paths_and_wavelengths_are_set() { + let path = minimal_path_spec(); // now includes wavelengths: vec![0.587] + let result = SequentialModelBuilder::new() + .paths(vec![path]) + .wavelengths(vec![0.587]) + .build(); + assert!(result.is_err()); + } + + #[test] + fn build_fails_when_both_paths_and_stop_surface_are_set() { + let path = minimal_path_spec(); + let result = SequentialModelBuilder::new() + .paths(vec![path]) + .stop_surface(0) + .build(); + assert!(result.is_err()); + } } diff --git a/crates/cherry-rs/src/core/sequential_model/mod.rs b/crates/cherry-rs/src/core/sequential_model/mod.rs index 4a971042..13decc8a 100644 --- a/crates/cherry-rs/src/core/sequential_model/mod.rs +++ b/crates/cherry-rs/src/core/sequential_model/mod.rs @@ -69,6 +69,9 @@ struct OpticalPath { stop_surface: Option, /// Step-indexed cursor state, parallel to `surface_indices`. steps: Vec, + /// Wavelengths this path's `submodels` were built from, in the same + /// order (`submodels[i]` corresponds to `wavelengths[i]`). + wavelengths: Vec, } /// A gap between two surfaces in a sequential system. @@ -90,7 +93,6 @@ pub struct Gap { pub struct SequentialModel { store: SurfaceStore, paths: Vec, - wavelengths: Vec, } /// A submodel of a sequential optical system. @@ -403,11 +405,11 @@ impl SequentialModel { submodels, stop_surface, steps: cursor_placements, + wavelengths: wavelengths.to_vec(), }; Ok(Self { store, paths: vec![path], - wavelengths: wavelengths.to_vec(), }) } } @@ -447,11 +449,11 @@ impl SequentialModel { submodels, stop_surface, steps: cursor_placements, + wavelengths: wavelengths.to_vec(), }; Ok(Self { store, paths: vec![path], - wavelengths: wavelengths.to_vec(), }) } @@ -518,11 +520,11 @@ impl SequentialModel { submodels, stop_surface, steps: cursor_placements, + wavelengths: wavelengths.to_vec(), }; Ok(Self { store, paths: vec![path], - wavelengths: wavelengths.to_vec(), }) } @@ -533,13 +535,8 @@ impl SequentialModel { /// versions can supply the appropriate `surface_from_spec` variant. fn from_path_specs_with_builder( paths: Vec, - wavelengths: &[Float], mut build_surface: impl FnMut(&SurfaceSpec) -> Result>, ) -> Result { - if wavelengths.is_empty() { - return Err(anyhow!("At least one wavelength must be specified.")); - } - let mut store_surfaces: Vec> = Vec::new(); let mut store_placements: Vec = Vec::new(); let mut optical_paths: Vec = Vec::new(); @@ -845,8 +842,13 @@ impl SequentialModel { .cloned() .collect(); + if ps.wavelengths.is_empty() { + return Err(anyhow!( + "each PathSpec must have at least one wavelength; this path has none" + )); + } let mut submodels: Vec = Vec::new(); - for &wavelength in wavelengths.iter() { + for &wavelength in ps.wavelengths.iter() { let gaps = Self::gap_specs_to_gaps(&all_gap_specs, wavelength)?; submodels.push(SequentialSubModelBase::new(gaps)); } @@ -864,6 +866,7 @@ impl SequentialModel { submodels, stop_surface: ps.stop_surface, steps: path_steps, + wavelengths: ps.wavelengths.clone(), }); } @@ -874,7 +877,6 @@ impl SequentialModel { Ok(Self { store, paths: optical_paths, - wavelengths: wavelengths.to_vec(), }) } @@ -882,18 +884,15 @@ impl SequentialModel { #[cfg(feature = "serde")] pub(crate) fn from_path_specs( paths: Vec, - wavelengths: &[Float], registry: Option<&SurfaceRegistry>, ) -> Result { - Self::from_path_specs_with_builder(paths, wavelengths, |spec| { - surface_from_spec(spec, registry) - }) + Self::from_path_specs_with_builder(paths, |spec| surface_from_spec(spec, registry)) } /// Builds a multipath model from `PathSpec`s (non-serde variant). #[cfg(not(feature = "serde"))] - pub(crate) fn from_path_specs(paths: Vec, wavelengths: &[Float]) -> Result { - Self::from_path_specs_with_builder(paths, wavelengths, surface_from_spec) + pub(crate) fn from_path_specs(paths: Vec) -> Result { + Self::from_path_specs_with_builder(paths, surface_from_spec) } /// Number of optical paths in the model. @@ -1019,8 +1018,15 @@ impl SequentialModel { } /// Returns the wavelengths at which the system is modeled. + /// + /// Single-path shorthand; delegates to `paths[0]`. pub fn wavelengths(&self) -> &[Float] { - &self.wavelengths + self.wavelengths_for_path(0) + } + + /// Returns the wavelengths for path `path_id`, in submodel order. + pub fn wavelengths_for_path(&self, path_id: usize) -> &[Float] { + &self.paths[path_id].wavelengths } fn gap_specs_to_gaps(gap_specs: &[GapSpec], wavelength: Float) -> Result> { diff --git a/crates/cherry-rs/src/examples/beam_splitter.rs b/crates/cherry-rs/src/examples/beam_splitter.rs index c62821f3..5c44e820 100644 --- a/crates/cherry-rs/src/examples/beam_splitter.rs +++ b/crates/cherry-rs/src/examples/beam_splitter.rs @@ -58,6 +58,7 @@ pub fn two_path_model( ], beam_splitter_arms: vec![BeamSplitterPathKind::Transmitting], stop_surface: None, + wavelengths: wavelengths.to_vec(), }; // Path 1: reflected arm. @@ -77,11 +78,11 @@ pub fn two_path_model( ], beam_splitter_arms: vec![BeamSplitterPathKind::Reflecting], stop_surface: None, + wavelengths: wavelengths.to_vec(), }; SequentialModelBuilder::new() .paths(vec![path_t, path_r]) - .wavelengths(wavelengths.to_vec()) .build() .expect("beam splitter model builds") .model diff --git a/crates/cherry-rs/src/examples/wf_epi_emission.rs b/crates/cherry-rs/src/examples/wf_epi_emission.rs index a66a966f..f1866ed4 100644 --- a/crates/cherry-rs/src/examples/wf_epi_emission.rs +++ b/crates/cherry-rs/src/examples/wf_epi_emission.rs @@ -87,13 +87,12 @@ pub fn sequential_model( ], gaps: vec![gap_0, gap_1, gap_2, gap_3, gap_4], beam_splitter_arms: vec![BeamSplitterPathKind::Transmitting], - stop_surface: None, + stop_surface: Some(1), + wavelengths: wavelengths.to_vec(), }; SequentialModelBuilder::new() .paths(vec![path]) - .stop_surface(1) - .wavelengths(wavelengths.to_vec()) .build() .expect("wf_epi_emission model builds") .model diff --git a/crates/cherry-rs/src/examples/wf_epi_excitation.rs b/crates/cherry-rs/src/examples/wf_epi_excitation.rs index 0835db95..92195132 100644 --- a/crates/cherry-rs/src/examples/wf_epi_excitation.rs +++ b/crates/cherry-rs/src/examples/wf_epi_excitation.rs @@ -77,11 +77,11 @@ pub fn sequential_model( gaps: vec![gap_0, gap_1, gap_2, gap_3], beam_splitter_arms: vec![BeamSplitterPathKind::Reflecting], stop_surface: Some(3), + wavelengths: wavelengths.to_vec(), }; SequentialModelBuilder::new() .paths(vec![path]) - .wavelengths(wavelengths.to_vec()) .build() .expect("wf_epi_excitation model builds") .model diff --git a/crates/cherry-rs/src/examples/wf_epi_microscope.rs b/crates/cherry-rs/src/examples/wf_epi_microscope.rs index 32015904..9e6c6a9a 100644 --- a/crates/cherry-rs/src/examples/wf_epi_microscope.rs +++ b/crates/cherry-rs/src/examples/wf_epi_microscope.rs @@ -26,7 +26,8 @@ use crate::{ pub fn sequential_model( n_air: Rc, n_oil: Rc, - wavelengths: &[f64], + excitation_wavelengths: &[f64], + emission_wavelengths: &[f64], ) -> SequentialModel { // Excitation path gaps. let gap_0 = GapSpec { @@ -98,6 +99,7 @@ pub fn sequential_model( gaps: vec![gap_0, gap_1, gap_2, gap_3], beam_splitter_arms: vec![BeamSplitterPathKind::Reflecting], stop_surface: Some(3), + wavelengths: excitation_wavelengths.to_vec(), }; let path_emission = PathSpec { @@ -133,11 +135,11 @@ pub fn sequential_model( gaps: vec![gap_bs_mirror, gap_mirror_tube, gap_tube_camera], beam_splitter_arms: vec![BeamSplitterPathKind::Transmitting], stop_surface: Some(3), + wavelengths: emission_wavelengths.to_vec(), }; SequentialModelBuilder::new() .paths(vec![path_excitation, path_emission]) - .wavelengths(wavelengths.to_vec()) .build() .expect("wf_epi_microscope model builds") .model diff --git a/crates/cherry-rs/src/specs/paths.rs b/crates/cherry-rs/src/specs/paths.rs index 5b0fe236..03dff1e3 100644 --- a/crates/cherry-rs/src/specs/paths.rs +++ b/crates/cherry-rs/src/specs/paths.rs @@ -1,3 +1,4 @@ +use crate::core::Float; use crate::specs::{ gaps::GapSpec, surfaces::{BeamSplitterPathKind, SurfaceSpec}, @@ -26,6 +27,10 @@ pub struct PathSpec { /// User-specified aperture stop as a store index for this path, or `None` /// to fall back to heuristic aperture-stop selection. pub stop_surface: Option, + /// Wavelengths at which this path's submodels are built. Independent of + /// every other path's wavelength list — no cross-path length or value + /// constraint. + pub wavelengths: Vec, } /// One element of a [`PathSpec`]'s surface sequence. @@ -84,6 +89,7 @@ mod tests { }], beam_splitter_arms: vec![], stop_surface: None, + wavelengths: vec![0.5876], }; assert_eq!(ps.surface_refs.len(), 2); } diff --git a/crates/cherry-rs/src/views/paraxial.rs b/crates/cherry-rs/src/views/paraxial.rs index 1ed901f8..a4048eb1 100644 --- a/crates/cherry-rs/src/views/paraxial.rs +++ b/crates/cherry-rs/src/views/paraxial.rs @@ -1967,6 +1967,7 @@ mod test { gaps: vec![gap_inf(), gap_50(), gap_50()], beam_splitter_arms: vec![BeamSplitterPathKind::Transmitting], stop_surface: None, + wavelengths: vec![0.5876e-3], }; // Path 1: Object(Shared) → BS(Shared) → Iris2(SD=15) → Image_1 @@ -1981,11 +1982,11 @@ mod test { gaps: vec![gap_inf(), gap_50(), gap_50()], beam_splitter_arms: vec![BeamSplitterPathKind::Reflecting], stop_surface: None, + wavelengths: vec![0.5876e-3], }; let model = SequentialModelBuilder::new() .paths(vec![path_t, path_r]) - .wavelengths(vec![0.5876e-3]) .build() .unwrap() .model; 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 70f0ccef..e7cc9a44 100644 --- a/crates/cherry-rs/src/views/ray_trace_3d/mod.rs +++ b/crates/cherry-rs/src/views/ray_trace_3d/mod.rs @@ -176,9 +176,9 @@ pub fn ray_trace_3d_view( // Use Axis::U as the canonical submodel for all field/wavelength combinations. // 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 triples: Vec<(usize, usize, usize)> = (0..sequential_model.path_count()) .flat_map(|p| { + let n_wavelengths = sequential_model.wavelengths_for_path(p).len(); (0..field_specs.len()).flat_map(move |f| (0..n_wavelengths).map(move |w| (p, f, w))) }) .collect(); @@ -1448,6 +1448,101 @@ mod tests { ); } + #[test] + fn ray_trace_3d_view_handles_differing_wavelength_counts_per_path() { + use std::rc::Rc; + + use crate::RefractiveIndexSpec; + use crate::core::math::linalg::rotations::{EulerAngles, Rotation3D}; + use crate::core::sequential_model::builder::SequentialModelBuilder; + use crate::specs::gaps::{ConstantRefractiveIndex, GapSpec}; + use crate::specs::paths::{PathSpec, PathSurfaceRef}; + use crate::specs::surfaces::{BeamSplitterPathKind, SurfaceSpec}; + + let n_air: Rc = Rc::new(ConstantRefractiveIndex::new(1.0, 0.0)); + let bs_rotation = + Rotation3D::IntrinsicPassiveRUF(EulerAngles((-45_f64).to_radians(), 0.0, 0.0)); + let img = || SurfaceSpec::Image { + rotation: Rotation3D::None, + decenter: Vec3::new(0.0, 0.0, 0.0), + rotation_offset: Rotation3D::None, + }; + let gap_inf = || GapSpec { + thickness: f64::INFINITY, + refractive_index: n_air.clone(), + }; + + 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(img()), + ], + gaps: vec![ + gap_inf(), + GapSpec { + thickness: 10.0, + refractive_index: n_air.clone(), + }, + ], + beam_splitter_arms: vec![BeamSplitterPathKind::Transmitting], + stop_surface: None, + wavelengths: vec![0.5876], // 1 wavelength + }; + let path_r = PathSpec { + surface_refs: vec![ + PathSurfaceRef::Shared(0), + PathSurfaceRef::Shared(1), + PathSurfaceRef::New(img()), + ], + gaps: vec![ + gap_inf(), + GapSpec { + thickness: 10.0, + refractive_index: n_air.clone(), + }, + ], + beam_splitter_arms: vec![BeamSplitterPathKind::Reflecting], + stop_surface: None, + wavelengths: vec![0.4861, 0.5876, 0.6563], // 3 wavelengths + }; + + let model = SequentialModelBuilder::new() + .paths(vec![path_t, path_r]) + .build() + .expect("beam splitter model builds") + .model; + + 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, ¶xial_view, config) + .expect("ray trace should succeed with differing per-path wavelength counts"); + + assert_eq!( + results.iter().filter(|r| r.path_id() == 0).count(), + field_specs.len() + ); + assert_eq!( + results.iter().filter(|r| r.path_id() == 1).count(), + 3 * field_specs.len() + ); + } + #[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/tests/beam_splitter_multipath.rs b/crates/cherry-rs/tests/beam_splitter_multipath.rs index 44a04a22..1c5f5b3e 100644 --- a/crates/cherry-rs/tests/beam_splitter_multipath.rs +++ b/crates/cherry-rs/tests/beam_splitter_multipath.rs @@ -112,11 +112,9 @@ fn at13_repeated_store_index_is_permitted_via_pathspec() { ], beam_splitter_arms: vec![], stop_surface: None, + wavelengths: vec![0.587], }; - let result = SequentialModelBuilder::new() - .paths(vec![path]) - .wavelengths(vec![0.587]) - .build(); + let result = SequentialModelBuilder::new().paths(vec![path]).build(); assert!(result.is_ok()); let model = result.unwrap().model; assert_eq!(model.path_surface_indices(0), &[0, 1, 1, 2]); diff --git a/crates/cherry-rs/tests/wf_epi_microscope.rs b/crates/cherry-rs/tests/wf_epi_microscope.rs index 44bdabbd..552bb610 100644 --- a/crates/cherry-rs/tests/wf_epi_microscope.rs +++ b/crates/cherry-rs/tests/wf_epi_microscope.rs @@ -1,9 +1,10 @@ use approx::assert_abs_diff_eq; use cherry_rs::examples::wf_epi_microscope::sequential_model; -use cherry_rs::{FieldSpec, ParaxialView, SequentialModel, n}; +use cherry_rs::{ApertureSpec, FieldSpec, ParaxialView, SamplingConfig, SequentialModel, n}; -const WAVELENGTHS: [f64; 1] = [0.5876]; +const EXCITATION_WAVELENGTHS: [f64; 1] = [0.488]; +const EMISSION_WAVELENGTHS: [f64; 1] = [0.520]; const FIELD_SPECS: [FieldSpec; 1] = [FieldSpec::PointSource { x: 0.0, y: 1.5 }]; const EXC_EFFECTIVE_FOCAL_LENGTH: f64 = -1.8750; @@ -15,7 +16,12 @@ const EXC_LAGRANGE_INVARIANT: f64 = -0.1594; const EMI_IMAGE_LOCATION: f64 = 199.8800167976479; fn model() -> SequentialModel { - sequential_model(n!(1.0), n!(1.5), &WAVELENGTHS) + sequential_model( + n!(1.0), + n!(1.5), + &EXCITATION_WAVELENGTHS, + &EMISSION_WAVELENGTHS, + ) } #[test] @@ -135,3 +141,46 @@ fn at_emission_paraxial_image_location() { epsilon = 1e-4 ); } + +#[test] +fn at_excitation_and_emission_wavelengths_differ() { + let model = model(); + assert_eq!(model.wavelengths_for_path(0), &EXCITATION_WAVELENGTHS); + assert_eq!(model.wavelengths_for_path(1), &EMISSION_WAVELENGTHS); +} + +#[test] +fn at_ray_trace_succeeds_with_differing_wavelength_counts_per_path() { + use cherry_rs::ray_trace_3d_view; + + let model = sequential_model( + n!(1.0), + n!(1.5), + &[0.488], // excitation: 1 wavelength + &[0.500, 0.520, 0.540], // emission: 3 wavelengths + ); + let view = ParaxialView::new(&model, &FIELD_SPECS, false).unwrap(); + + // AT-7: subview counts per path match each path's own wavelength count. + let path0_subviews = view.iter().filter(|sv| sv.path_id() == 0).count(); + let path1_subviews = view.iter().filter(|sv| sv.path_id() == 1).count(); + assert_eq!(path1_subviews, 3 * path0_subviews); + + // AT-6: ray_trace_3d_view succeeds and produces the right result counts. + let aperture = ApertureSpec::EntrancePupil { semi_diameter: 1.5 }; + let config = SamplingConfig { + n_fan_rays: 5, + full_pupil_spacing: 0.1, + }; + let trace = ray_trace_3d_view(&aperture, &FIELD_SPECS, &model, &view, config) + .expect("ray trace should succeed with differing per-path wavelength counts"); + + assert_eq!( + trace.iter().filter(|r| r.path_id() == 0).count(), + FIELD_SPECS.len() + ); + assert_eq!( + trace.iter().filter(|r| r.path_id() == 1).count(), + 3 * FIELD_SPECS.len() + ); +} From 5933db076b728ac48794c103cbc97594a30d46ea Mon Sep 17 00:00:00 2001 From: Kyle Douglass Date: Tue, 7 Jul 2026 14:41:38 +0200 Subject: [PATCH 2/3] fix: Paraxial view resolve stop surfaces into step indexes, not store indexes --- crates/cherry-rs/src/views/paraxial.rs | 71 +++++++++++++++++++++++++- 1 file changed, 69 insertions(+), 2 deletions(-) diff --git a/crates/cherry-rs/src/views/paraxial.rs b/crates/cherry-rs/src/views/paraxial.rs index a4048eb1..b461556d 100644 --- a/crates/cherry-rs/src/views/paraxial.rs +++ b/crates/cherry-rs/src/views/paraxial.rs @@ -523,7 +523,7 @@ impl ParaxialSubView { )?; let aperture_stop = match data.stop_surface { - Some(i) => i, + Some(i) => resolve_stop_surface_step(i, surface_indices)?, None => Self::calc_aperture_stop( surfaces, placements, @@ -1432,6 +1432,38 @@ pub(crate) fn calc_marginal_ray( } } +/// Resolves a user-specified `stop_surface` (a store index, validated at +/// build time against the whole model's surface store — +/// `SequentialModel::validate_stop_surface`) to a step index within one +/// path's own traversal. +/// +/// Every paraxial computation that consumes an aperture stop +/// (`calc_marginal_ray`, the ray-height-based ratios in `calc_aperture_stop`) +/// indexes by step position within `surface_indices`, not by store position. +/// For a path whose surfaces are all `New` in store order, store index and step +/// index coincide by construction, so this only has visible effect once a path +/// contains `Shared`/`ObjectLinkedTo` steps ahead of the stop. +/// +/// Errors if `store_index` does not appear in `surface_indices` at all, or +/// appears more than once (ambiguous which occurrence is the stop). +fn resolve_stop_surface_step(store_index: usize, surface_indices: &[usize]) -> Result { + let mut occurrences = surface_indices + .iter() + .enumerate() + .filter(|&(_, &si)| si == store_index) + .map(|(step, _)| step); + match (occurrences.next(), occurrences.next()) { + (None, _) => Err(anyhow!( + "stop surface (store index {store_index}) is not visited by this path" + )), + (Some(_), Some(_)) => Err(anyhow!( + "stop surface (store index {store_index}) is visited more than once by this \ + path; ambiguous which occurrence is the aperture stop" + )), + (Some(step), None) => Ok(step), + } +} + /// Compute the paraxial marginal ray bundle for a given wavelength. /// /// Uses the first tangential direction `(0, 1, 0)`, valid for all rotationally @@ -1460,7 +1492,7 @@ pub(crate) fn marginal_ray_bundle( path_steps, )?; let stop = match model.stop_surface() { - Some(i) => i, + Some(i) => resolve_stop_surface_step(i, surface_indices)?, None => calc_aperture_stop( surfaces, placements, @@ -2016,6 +2048,41 @@ mod test { ); } + /// A path's user-specified `stop_surface` is a *store* index (validated + /// against the whole model's surface store at build time, + /// `SequentialModel::validate_stop_surface`), but every downstream + /// paraxial computation (`calc_marginal_ray`, the `None` branch's + /// `calc_aperture_stop`) indexes the aperture stop by *step* position + /// within that path's own traversal. For a path whose surfaces are all + /// `New` in store order (e.g. path 0 of most models, where store index + /// and step index coincide by construction), a mismatch between the two + /// numberspaces is invisible. `wf_epi_microscope`'s emission path uses + /// `Shared`/`ObjectLinkedTo` steps, so the two numberspaces genuinely + /// diverge there, catching it. + #[test] + fn stop_surface_step_resolves_from_store_index_on_a_reordered_path() { + use crate::examples::wf_epi_microscope::sequential_model; + + let model = sequential_model(n!(1.0), n!(1.5), &[0.488], &[0.520]); + + // path_emission's surface_indices are [5, 3, 2, 6, 7, 8] (synthesized + // Object=5, Shared objective=3, Shared beam splitter=2, mirror=6, + // tube lens=7, Image=8) — store index 3 (the objective, the value + // configured via `stop_surface`) sits at step 1, not step 3. + assert_eq!(model.path_surface_indices(1), &[5, 3, 2, 6, 7, 8]); + assert_eq!(model.stop_surface_for_path(1), Some(3)); + + let field_specs = vec![FieldSpec::PointSource { x: 0.0, y: 1.5 }]; + let view = ParaxialView::new(&model, &field_specs, false).unwrap(); + let sub = view.get_for_path(1, 0, 0).expect("path 1 subview"); + + assert_eq!( + *sub.aperture_stop(), + 1, + "the objective (store index 3) is path 1's step 1, not step 3" + ); + } + #[test] fn test_lagrange_invariant_is_conserved() { let (view, _) = setup(); From 80cd089b2ad30bcf1c29fabc4790adb15121875d Mon Sep 17 00:00:00 2001 From: Kyle Douglass Date: Tue, 7 Jul 2026 15:32:57 +0200 Subject: [PATCH 3/3] feat: Enable multipath solves --- .../src/core/sequential_model/builder.rs | 229 +++++++++- .../src/core/sequential_model/mod.rs | 45 +- .../src/core/sequential_model/solves/fno.rs | 392 ++++++++++++++++-- .../sequential_model/solves/marginal_ray.rs | 84 +++- .../src/core/sequential_model/solves/mod.rs | 35 +- .../src/examples/wf_epi_excitation.rs | 10 +- crates/cherry-rs/src/views/paraxial.rs | 26 +- 7 files changed, 749 insertions(+), 72 deletions(-) diff --git a/crates/cherry-rs/src/core/sequential_model/builder.rs b/crates/cherry-rs/src/core/sequential_model/builder.rs index bca6885e..3fe4df05 100644 --- a/crates/cherry-rs/src/core/sequential_model/builder.rs +++ b/crates/cherry-rs/src/core/sequential_model/builder.rs @@ -5,7 +5,10 @@ use anyhow::{Result, anyhow}; use crate::core::surfaces::SurfaceRegistry; use crate::specs::{gaps::GapSpec, paths::PathSpec, surfaces::SurfaceSpec}; -use super::{SequentialModel, solves::Solve}; +use super::{ + SequentialModel, + solves::{Solve, SolveKind}, +}; /// The output of a successful [`SequentialModelBuilder::build()`] call. /// Carries the model and the post-solve specs so callers can extract @@ -59,11 +62,30 @@ impl SequentialModelBuilder { self.validate()?; if self.paths.is_some() { - let paths = self.paths.unwrap(); + let mut paths = self.paths.unwrap(); + let solves = self.solves; #[cfg(feature = "serde")] - let model = SequentialModel::from_path_specs(paths, self.registry.as_ref())?; - #[cfg(not(feature = "serde"))] - let model = SequentialModel::from_path_specs(paths)?; + let registry = self.registry; + + let build = |paths: &[PathSpec]| -> Result { + #[cfg(feature = "serde")] + return SequentialModel::from_path_specs(paths, registry.as_ref()); + #[cfg(not(feature = "serde"))] + SequentialModel::from_path_specs(paths) + }; + + let mut model = build(&paths)?; + + if !solves.is_empty() { + validate_multipath_solves(&solves, &model)?; + let mut solves = solves; + solves.sort_by_key(|s| (s.parameter_kind(), s.surface_index(), s.path_id())); + for solve in &solves { + solve.apply_multipath(&model, &mut paths)?; + model = build(&paths)?; + } + } + return Ok(BuildResult { model, gap_specs: vec![], @@ -195,6 +217,47 @@ impl SequentialModelBuilder { } } +/// Validates every solve's `path_id()`/`surface_index()` against the +/// just-built (pre-solve) multipath `model`, and enforces FR-5's global +/// uniqueness rule for `Curvature`-kind solves. Runs once, before any solve +/// is applied, so a bad solve fails fast rather than mid-rebuild. +fn validate_multipath_solves(solves: &[Box], model: &SequentialModel) -> Result<()> { + let path_count = model.path_count(); + let n_surfaces = model.surfaces().len(); + for solve in solves { + let path_id = solve.path_id(); + if path_id >= path_count { + return Err(anyhow!( + "solve targets path_id {path_id} but the model has only {path_count} path(s)" + )); + } + if solve.parameter_kind() == SolveKind::Curvature { + let target = solve.surface_index(); + if target >= n_surfaces { + return Err(anyhow!( + "surface_index {target} does not exist in the model \ + (model has {n_surfaces} surface(s))" + )); + } + } + } + + let mut seen_curvature_targets: Vec = Vec::new(); + for solve in solves { + if solve.parameter_kind() == SolveKind::Curvature { + let target = solve.surface_index(); + if seen_curvature_targets.contains(&target) { + return Err(anyhow!( + "surface {target} has more than one Curvature-kind solve across the \ + model's paths; a shared surface may only be targeted by one" + )); + } + seen_curvature_targets.push(target); + } + } + Ok(()) +} + #[cfg(test)] mod tests { use super::*; @@ -943,4 +1006,160 @@ mod tests { .build(); assert!(result.is_err()); } + + // ── Per-path solves ─────────────────────────────────────────────────── + + #[test] + fn multipath_thickness_solve_on_one_path_leaves_other_unaffected() { + use crate::core::sequential_model::solves::MarginalRaySolve; + use crate::specs::fields::FieldSpec; + use crate::specs::paths::{PathSpec, PathSurfaceRef}; + use crate::views::paraxial::ParaxialView; + use approx::assert_abs_diff_eq; + + let img = || SurfaceSpec::Image { + rotation: Rotation3D::None, + decenter: Vec3::new(0.0, 0.0, 0.0), + rotation_offset: Rotation3D::None, + }; + let path0 = PathSpec { + surface_refs: vec![ + PathSurfaceRef::New(SurfaceSpec::Object), + PathSurfaceRef::New(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, + }), + PathSurfaceRef::New(img()), + ], + gaps: vec![ + GapSpec { + thickness: f64::INFINITY, + refractive_index: n!(1.0), + }, + GapSpec { + thickness: 1.0, + refractive_index: n!(1.5), + }, // untouched + ], + beam_splitter_arms: vec![], + stop_surface: None, + wavelengths: vec![0.5876], + }; + let path1 = PathSpec { + surface_refs: vec![ + PathSurfaceRef::Shared(0), // same Object + PathSurfaceRef::Shared(1), // same Sphere + PathSurfaceRef::New(img()), + ], + gaps: vec![ + GapSpec { + thickness: f64::INFINITY, + refractive_index: n!(1.0), + }, + GapSpec { + thickness: 1.0, + refractive_index: n!(1.5), + }, // will be solved + ], + beam_splitter_arms: vec![], + stop_surface: None, + wavelengths: vec![0.5876], + }; + + let model = SequentialModelBuilder::new() + .paths(vec![path0, path1]) + .solves(vec![Box::new( + MarginalRaySolve::new(1, 0.0, 0).with_path_id(1), + )]) + .build() + .expect("build should succeed") + .model; + + // Path 1's gap was solved to place its image at the paraxial focus. + let pv = ParaxialView::new( + &model, + &[FieldSpec::Angle { + chi: 0.0, + phi: 90.0, + }], + false, + ) + .unwrap(); + let sub1 = pv.get_for_path(1, 0, 0).unwrap(); + assert_abs_diff_eq!( + sub1.marginal_ray().rays_at_surface(2)[0].height, + 0.0, + epsilon = 1e-4 + ); + + // Path 0's own gap thickness is untouched (still 1.0, not solved). + // The shared Sphere lands at z = 0.0 (object-at-infinity placement + // rule), so path 0's image sits exactly one (unsolved) gap further. + assert_abs_diff_eq!(model.path_placement(0, 2).position.z(), 1.0, epsilon = 1e-9); + } + + #[test] + fn multipath_solve_with_out_of_range_path_id_is_rejected() { + use crate::core::sequential_model::solves::MarginalRaySolve; + + let path = minimal_path_spec(); // single path, path_count() == 1 + let result = SequentialModelBuilder::new() + .paths(vec![path]) + .solves(vec![Box::new( + MarginalRaySolve::new(0, 0.0, 0).with_path_id(1), // only path 0 exists + )]) + .build(); + assert!(result.is_err()); + } + + #[test] + fn multipath_curvature_solve_with_invalid_store_index_is_rejected_distinctly() { + use crate::core::sequential_model::solves::FNumberSolve; + + // Store index 99 doesn't exist anywhere in this model at all (only 2 + // surfaces exist: Object=0, Image=1) — distinct from the "exists but + // not visited by this path" case, covered in fno.rs's test module. + // `validate_multipath_solves` must reject this before any solve + // runs, with a message naming the invalid index rather than reusing + // `apply_multipath`'s "path does not visit surface" wording. + let path = minimal_path_spec(); + let result = SequentialModelBuilder::new() + .paths(vec![path]) + .solves(vec![Box::new(FNumberSolve::new(99, 4.0, 0))]) + .build(); + assert!(result.is_err()); + } + + #[test] + fn custom_solve_without_apply_multipath_errors_through_build() { + // A `Solve` implementation that only overrides `apply` (not + // `apply_multipath`) must hit the trait's default (error) body when + // run through the real multipath builder — not just when + // `apply_multipath` is called directly against a hand-built model. + struct SingleThicknessOnly; + impl Solve for SingleThicknessOnly { + fn apply( + &self, + _model: &SequentialModel, + _gap_specs: &mut Vec, + _surface_specs: &mut Vec, + ) -> Result<()> { + Ok(()) + } + fn surface_index(&self) -> usize { + 0 + } + } + + let path = minimal_path_spec(); + let result = SequentialModelBuilder::new() + .paths(vec![path]) + .solves(vec![Box::new(SingleThicknessOnly)]) + .build(); + assert!(result.is_err()); + } } diff --git a/crates/cherry-rs/src/core/sequential_model/mod.rs b/crates/cherry-rs/src/core/sequential_model/mod.rs index 13decc8a..3771a3f0 100644 --- a/crates/cherry-rs/src/core/sequential_model/mod.rs +++ b/crates/cherry-rs/src/core/sequential_model/mod.rs @@ -331,6 +331,43 @@ pub fn reversed_surface_id(num_surfaces: usize, surf_id: usize) -> usize { /// Only meaningful for `Reversed` paths; callers must gate on orientation /// before using the result (see /// [`SequentialModel::from_path_specs_with_builder`]). +/// Locates the `PathSpec`/step that introduced store index `store_index` via +/// `PathSurfaceRef::New`, using the identical counting rule +/// `from_path_specs_with_builder` uses to assign store indices: `New` and +/// `ObjectLinkedTo` steps each consume one index, in path/step order; +/// `Shared` steps do not. Returns `None` if `store_index` is out of range. +/// +/// An `ObjectLinkedTo` step also consumes a store index by this rule but +/// never owns a mutable `SurfaceSpec` (it always synthesizes a plain +/// `SurfaceSpec::Object` at build time) — callers needing a mutable target +/// must match the result against `PathSurfaceRef::New` and treat any other +/// variant as ineligible; no special-casing is done here. +/// +/// `paths`' topology (which steps are `New`/`Shared`/`ObjectLinkedTo`, and +/// their order) never changes across a builder's solve-rebuild iterations — +/// only the `SurfaceSpec` values inside `New` steps do — so a caller +/// applying multiple solves may compute this once and reuse it. +pub(crate) fn locate_surface_owner( + paths: &[PathSpec], + store_index: usize, +) -> Option<(usize, usize)> { + let mut counter = 0; + for (path_idx, ps) in paths.iter().enumerate() { + for (step_idx, sref) in ps.surface_refs.iter().enumerate() { + match sref { + PathSurfaceRef::New(_) | PathSurfaceRef::ObjectLinkedTo { .. } => { + if counter == store_index { + return Some((path_idx, step_idx)); + } + counter += 1; + } + PathSurfaceRef::Shared(_) => {} + } + } + } + None +} + fn count_leading_shared(refs: &[PathSurfaceRef]) -> usize { refs.iter() .skip(1) @@ -534,7 +571,7 @@ impl SequentialModel { /// from a [`SurfaceSpec`]. It is injected so that the serde and non-serde /// versions can supply the appropriate `surface_from_spec` variant. fn from_path_specs_with_builder( - paths: Vec, + paths: &[PathSpec], mut build_surface: impl FnMut(&SurfaceSpec) -> Result>, ) -> Result { let mut store_surfaces: Vec> = Vec::new(); @@ -652,7 +689,7 @@ impl SequentialModel { }; // Build the dense arm vec by consuming beam_splitter_arms in step order. - let mut bs_arms_iter = ps.beam_splitter_arms.into_iter(); + let mut bs_arms_iter = ps.beam_splitter_arms.iter().copied(); let mut dense_bs_arms: Vec> = Vec::with_capacity(n_refs); let mut cursor = match linked { @@ -883,7 +920,7 @@ impl SequentialModel { /// Builds a multipath model from `PathSpec`s (serde + registry variant). #[cfg(feature = "serde")] pub(crate) fn from_path_specs( - paths: Vec, + paths: &[PathSpec], registry: Option<&SurfaceRegistry>, ) -> Result { Self::from_path_specs_with_builder(paths, |spec| surface_from_spec(spec, registry)) @@ -891,7 +928,7 @@ impl SequentialModel { /// Builds a multipath model from `PathSpec`s (non-serde variant). #[cfg(not(feature = "serde"))] - pub(crate) fn from_path_specs(paths: Vec) -> Result { + pub(crate) fn from_path_specs(paths: &[PathSpec]) -> Result { Self::from_path_specs_with_builder(paths, surface_from_spec) } diff --git a/crates/cherry-rs/src/core/sequential_model/solves/fno.rs b/crates/cherry-rs/src/core/sequential_model/solves/fno.rs index 8597fd47..b99ebf5c 100644 --- a/crates/cherry-rs/src/core/sequential_model/solves/fno.rs +++ b/crates/cherry-rs/src/core/sequential_model/solves/fno.rs @@ -1,9 +1,16 @@ use anyhow::{Result, anyhow}; use crate::{ - core::{Float, sequential_model::SequentialSubModel}, - specs::{gaps::GapSpec, surfaces::SurfaceSpec}, - views::paraxial::marginal_ray_bundle, + core::{ + Float, + sequential_model::{SequentialSubModel, locate_surface_owner}, + }, + specs::{ + gaps::GapSpec, + paths::{PathSpec, PathSurfaceRef}, + surfaces::SurfaceSpec, + }, + views::paraxial::{marginal_ray_bundle, marginal_ray_bundle_for_path}, }; use super::super::SequentialModel; @@ -15,6 +22,7 @@ pub struct FNumberSolve { surface_index: usize, target_fno: Float, wavelength_id: usize, + path_id: usize, } impl FNumberSolve { @@ -23,11 +31,78 @@ impl FNumberSolve { surface_index, target_fno, wavelength_id, + path_id: 0, + } + } + + /// Targets path `path_id` in a multipath model. Defaults to 0; has no + /// effect on single-path models. + pub fn with_path_id(mut self, path_id: usize) -> Self { + self.path_id = path_id; + self + } + + /// Computes the radius of curvature that yields `target_fno`, given the + /// marginal ray's incoming angle/height and the refractive indices on + /// either side. Shared by `apply` and `apply_multipath`. + fn solved_roc( + prev_angle: Float, + height: Float, + n_0: Float, + n_1: Float, + target_fno: Float, + surface_index: usize, + ) -> Result { + let eps = Float::EPSILON * height.abs().max(1.0); + if height.abs() < eps { + return Err(anyhow!( + "marginal ray height at surface {surface_index} is effectively zero; \ + ROC is indeterminate" + )); + } + + let denom = n_1 / (2.0 * target_fno) + n_0 * prev_angle; + let roc = if denom.abs() < eps { + Float::INFINITY + } else { + (n_1 - n_0) * height / denom + }; + + if roc == 0.0 { + return Err(anyhow!( + "computed ROC is zero at surface {surface_index}; result is unphysical" + )); + } + Ok(roc) + } + + /// Writes `roc` into `spec` if it is a `Sphere` or `Conic`. Shared by + /// `apply` and `apply_multipath`. + fn write_roc(spec: &mut SurfaceSpec, roc: Float, surface_index: usize) -> Result<()> { + match spec { + SurfaceSpec::Sphere { + radius_of_curvature, + .. + } + | SurfaceSpec::Conic { + radius_of_curvature, + .. + } => { + *radius_of_curvature = roc; + Ok(()) + } + _ => Err(anyhow!( + "surface {surface_index} is not a Sphere or Conic; cannot set radius of curvature" + )), } } } impl Solve for FNumberSolve { + fn path_id(&self) -> usize { + self.path_id + } + fn apply( &self, model: &SequentialModel, @@ -59,15 +134,6 @@ impl Solve for FNumberSolve { let u = bundle.rays_at_surface(prev_idx)[0].angle; let y = bundle.rays_at_surface(self.surface_index)[0].height; - let eps = Float::EPSILON * y.abs().max(1.0); - if y.abs() < eps { - return Err(anyhow!( - "marginal ray height at surface {} is effectively zero; \ - ROC is indeterminate", - self.surface_index - )); - } - let submodel = model .submodel(self.wavelength_id) .ok_or_else(|| anyhow!("wavelength_id {} out of range", self.wavelength_id))?; @@ -83,37 +149,97 @@ impl Solve for FNumberSolve { .refractive_index .n(); - let denom = n_1 / (2.0 * self.target_fno) + n_0 * u; - let roc = if denom.abs() < eps { - Float::INFINITY - } else { - (n_1 - n_0) * y / denom - }; + let roc = Self::solved_roc(u, y, n_0, n_1, self.target_fno, self.surface_index)?; + Self::write_roc( + &mut surface_specs[self.surface_index], + roc, + self.surface_index, + ) + } - if roc == 0.0 { + fn apply_multipath(&self, model: &SequentialModel, paths: &mut [PathSpec]) -> Result<()> { + let path_id = self.path_id(); + if path_id >= paths.len() { return Err(anyhow!( - "computed ROC is zero at surface {}; result is unphysical", - self.surface_index + "solve targets path_id {path_id} but the model has only {} path(s)", + paths.len() )); } - match &mut surface_specs[self.surface_index] { - SurfaceSpec::Sphere { - radius_of_curvature, - .. + let store_index = self.surface_index; + let path_surface_indices = model.path_surface_indices(path_id); + let mut occurrences = path_surface_indices + .iter() + .enumerate() + .filter(|&(_, &si)| si == store_index) + .map(|(step, _)| step); + let eval_step = match (occurrences.next(), occurrences.next()) { + (None, _) => { + return Err(anyhow!( + "path {path_id} does not visit surface {store_index}" + )); } - | SurfaceSpec::Conic { - radius_of_curvature, - .. - } => { - *radius_of_curvature = roc; - Ok(()) + (Some(_), Some(_)) => { + return Err(anyhow!( + "surface {store_index} is visited more than once by path {path_id}; \ + ambiguous which occurrence to evaluate the solve against" + )); } - _ => Err(anyhow!( - "surface {} is not a Sphere or Conic; cannot set radius of curvature", - self.surface_index - )), + (Some(step), None) => step, + }; + let prev_step = eval_step.checked_sub(1).ok_or_else(|| { + anyhow!( + "surface {store_index} is the first step of path {path_id}; \ + cannot apply F/# solve at the object surface" + ) + })?; + + let n_wavelengths = model.wavelengths_for_path(path_id).len(); + if self.wavelength_id >= n_wavelengths { + return Err(anyhow!( + "wavelength_id {} is out of range (path {path_id} has {} wavelength(s))", + self.wavelength_id, + n_wavelengths + )); } + + let bundle = marginal_ray_bundle_for_path(model, path_id, self.wavelength_id)?; + let u = bundle.rays_at_surface(prev_step)[0].angle; + let y = bundle.rays_at_surface(eval_step)[0].height; + + let submodel = model + .submodels_for_path(path_id) + .get(self.wavelength_id) + .ok_or_else(|| anyhow!("wavelength_id {} out of range", self.wavelength_id))?; + let gaps = submodel.gaps(); + let n_0 = gaps + .get(prev_step) + .ok_or_else(|| anyhow!("no gap before step {prev_step} on path {path_id}"))? + .refractive_index + .n(); + let n_1 = gaps + .get(eval_step) + .ok_or_else(|| anyhow!("no gap after step {eval_step} on path {path_id}"))? + .refractive_index + .n(); + + let roc = Self::solved_roc(u, y, n_0, n_1, self.target_fno, store_index)?; + + let (owner_path, owner_step) = + locate_surface_owner(paths, store_index).ok_or_else(|| { + anyhow!("surface {store_index} is not introduced by any path (internal error)") + })?; + let spec = match &mut paths[owner_path].surface_refs[owner_step] { + PathSurfaceRef::New(spec) => spec, + _ => { + return Err(anyhow!( + "surface {store_index} (owned by path {owner_path}, step {owner_step}) is \ + not a New surface_ref; Shared/ObjectLinkedTo surfaces cannot be solved \ + directly" + )); + } + }; + Self::write_roc(spec, roc, store_index) } fn surface_index(&self) -> usize { @@ -321,4 +447,200 @@ mod tests { let bundle = marginal_ray_bundle(&model, 0).unwrap(); assert_abs_diff_eq!(bundle.rays_at_surface(3)[0].height, 0.0, epsilon = 1e-3); } + + // ── Per-path solves (multipath) ──────────────────────────────────────── + + /// Two paths sharing one refracting Sphere (store index 1). Path 0's own + /// Iris (SD=3) is more constraining than the Sphere (SD=10), so path 0's + /// aperture stop is its own Iris; path 1 has no Iris, so the shared + /// Sphere is path 1's aperture stop. + /// + /// Store: Object(0), Sphere(1), Iris(2), Image_0(3), Image_1(4). + fn shared_sphere_two_path_fixture() -> (PathSpec, PathSpec) { + use crate::specs::paths::PathSurfaceRef; + + let img = || SurfaceSpec::Image { + rotation: Rotation3D::None, + decenter: Vec3::new(0.0, 0.0, 0.0), + rotation_offset: Rotation3D::None, + }; + let path0 = PathSpec { + surface_refs: vec![ + PathSurfaceRef::New(SurfaceSpec::Object), + PathSurfaceRef::New(SurfaceSpec::Sphere { + semi_diameter: 10.0, + radius_of_curvature: 100.0, // placeholder; solved below + surf_kind: BoundaryKind::Refracting, + rotation: Rotation3D::None, + decenter: Vec3::new(0.0, 0.0, 0.0), + rotation_offset: Rotation3D::None, + }), + PathSurfaceRef::New(SurfaceSpec::Iris { + semi_diameter: 3.0, + rotation: Rotation3D::None, + decenter: Vec3::new(0.0, 0.0, 0.0), + rotation_offset: Rotation3D::None, + }), + PathSurfaceRef::New(img()), + ], + gaps: vec![ + GapSpec { + thickness: Float::INFINITY, + refractive_index: n!(1.0), + }, + GapSpec { + thickness: 50.0, + refractive_index: n!(1.5), + }, + GapSpec { + thickness: 50.0, + refractive_index: n!(1.5), + }, + ], + beam_splitter_arms: vec![], + stop_surface: None, + wavelengths: vec![0.5876], + }; + // Path 1's own gaps around the shared Sphere: glass before, air after + // (image space), so its image-space index matches the assumption + // `FNumberSolve::solved_roc` makes (exit angle == -1/(2*target_fno) + // directly, without an image-space index correction) — the same + // assumption the single-path convexplano fixture above satisfies by + // using air as its last gap's medium. + let path1 = PathSpec { + surface_refs: vec![ + PathSurfaceRef::Shared(0), + PathSurfaceRef::Shared(1), + PathSurfaceRef::New(img()), + ], + gaps: vec![ + GapSpec { + thickness: Float::INFINITY, + refractive_index: n!(1.5), + }, + GapSpec { + thickness: 50.0, + refractive_index: n!(1.0), + }, + ], + beam_splitter_arms: vec![], + stop_surface: None, + wavelengths: vec![0.5876], + }; + (path0, path1) + } + + #[test] + fn multipath_fno_solve_on_shared_surface_updates_both_paths() { + let (path0, path1) = shared_sphere_two_path_fixture(); + let target = 4.0; + + let model = SequentialModelBuilder::new() + .paths(vec![path0, path1]) + .solves(vec![Box::new( + FNumberSolve::new(1, target, 0).with_path_id(1), + )]) + .build() + .expect("build should succeed") + .model; + + let field_specs = vec![FieldSpec::Angle { + chi: 0.0, + phi: 90.0, + }]; + let pv = ParaxialView::new(&model, &field_specs, false).unwrap(); + + // Path 1 (the evaluation path): F/# matches the target directly. + let sub1 = pv.get_for_path(1, 0, 0).unwrap(); + assert_abs_diff_eq!(sub1.paraxial_fno().abs(), target, epsilon = 1e-3); + + // Path 0 sees the *same physical surface* — same store index, same + // curvature — but its own stop is its Iris, not the shared Sphere, + // so its own F/# generically differs from path 1's target (proving + // this was evaluated against path 1's own aperture, not path 0's). + let sub0 = pv.get_for_path(0, 0, 0).unwrap(); + assert!((sub0.paraxial_fno().abs() - target).abs() > 1e-2); + + // Both paths read the identical, updated ROC off the one physical surface. + let store_idx = model.path_surface_indices(1)[1]; + assert_eq!(model.path_surface_indices(0)[1], store_idx); + let roc = model.surfaces()[store_idx].roc(0.0); + assert!(roc.is_finite() && roc != 100.0); // changed from the placeholder + } + + #[test] + fn multipath_fno_solve_errors_if_path_does_not_visit_surface() { + let (path0, path1) = shared_sphere_two_path_fixture(); + // path0's own Iris (store index 2) is never visited by path 1. + let result = SequentialModelBuilder::new() + .paths(vec![path0, path1]) + .solves(vec![Box::new(FNumberSolve::new(2, 4.0, 0).with_path_id(1))]) + .build(); + assert!(result.is_err()); + } + + #[test] + fn multipath_fno_solve_errors_on_repeated_store_index_visit() { + use crate::specs::paths::PathSurfaceRef; + + // Path: [New(Obj), New(Sphere), Shared(1), New(Img)] — surface_indices = [0, 1, + // 1, 2]. + let path = PathSpec { + surface_refs: vec![ + PathSurfaceRef::New(SurfaceSpec::Object), + PathSurfaceRef::New(SurfaceSpec::Sphere { + semi_diameter: 12.7, + radius_of_curvature: 65.0, + surf_kind: BoundaryKind::Refracting, + rotation: Rotation3D::None, + decenter: Vec3::new(0.0, 0.0, 0.0), + rotation_offset: Rotation3D::None, + }), + PathSurfaceRef::Shared(1), + PathSurfaceRef::New(SurfaceSpec::Image { + rotation: Rotation3D::None, + decenter: Vec3::new(0.0, 0.0, 0.0), + rotation_offset: Rotation3D::None, + }), + ], + gaps: vec![ + GapSpec { + thickness: Float::INFINITY, + refractive_index: n!(1.0), + }, + GapSpec { + thickness: 10.0, + refractive_index: n!(1.5), + }, + GapSpec { + thickness: 10.0, + refractive_index: n!(1.5), + }, + ], + beam_splitter_arms: vec![], + stop_surface: None, + wavelengths: vec![0.5876], + }; + let result = SequentialModelBuilder::new() + .paths(vec![path]) + .solves(vec![Box::new(FNumberSolve::new(1, 4.0, 0))]) + .build(); + assert!(result.is_err()); + } + + #[test] + fn two_curvature_solves_on_the_same_store_index_are_rejected() { + // Two paths sharing one Sphere (store index 1); both declare an + // FNumberSolve on it — rejected regardless of differing path_id. + let (path0, path1) = shared_sphere_two_path_fixture(); + let result = SequentialModelBuilder::new() + .paths(vec![path0, path1]) + .solves(vec![ + Box::new(FNumberSolve::new(1, 4.0, 0)), // path_id defaults to 0 + Box::new(FNumberSolve::new(1, 8.0, 0).with_path_id(1)), /* different path_id, + * same target */ + ]) + .build(); + assert!(result.is_err()); + } } diff --git a/crates/cherry-rs/src/core/sequential_model/solves/marginal_ray.rs b/crates/cherry-rs/src/core/sequential_model/solves/marginal_ray.rs index 0043e8de..d7b49c97 100644 --- a/crates/cherry-rs/src/core/sequential_model/solves/marginal_ray.rs +++ b/crates/cherry-rs/src/core/sequential_model/solves/marginal_ray.rs @@ -2,8 +2,8 @@ use anyhow::{Result, anyhow}; use crate::{ core::Float, - specs::{gaps::GapSpec, surfaces::SurfaceSpec}, - views::paraxial::marginal_ray_bundle, + specs::{gaps::GapSpec, paths::PathSpec, surfaces::SurfaceSpec}, + views::paraxial::{ParaxialRayBundle, marginal_ray_bundle, marginal_ray_bundle_for_path}, }; use super::super::SequentialModel; @@ -19,6 +19,7 @@ pub struct MarginalRaySolve { gap_index: usize, target_height: Float, wavelength_id: usize, + path_id: usize, } impl MarginalRaySolve { @@ -27,8 +28,45 @@ impl MarginalRaySolve { gap_index, target_height, wavelength_id, + path_id: 0, } } + + /// Targets path `path_id` in a multipath model. Defaults to 0; has no + /// effect on single-path models (there is only ever path 0). + pub fn with_path_id(mut self, path_id: usize) -> Self { + self.path_id = path_id; + self + } + + /// Computes the gap thickness that places the marginal ray at + /// `target_height` at `gap_index` within `bundle`. Shared by `apply` + /// (single-path) and `apply_multipath`. + fn solved_thickness( + bundle: &ParaxialRayBundle, + gap_index: usize, + target_height: Float, + ) -> Result { + let ray = &bundle.rays_at_surface(gap_index)[0]; + let h = ray.height; + let u_prime = ray.angle; + + let eps = Float::EPSILON * h.abs().max(1.0); + if u_prime.abs() < eps { + return Err(anyhow!( + "marginal ray angle at surface {gap_index} is effectively zero; \ + thickness is indeterminate (collimated beam in gap space)" + )); + } + + let t = (target_height - h) / u_prime; + if t < 0.0 { + return Err(anyhow!( + "computed gap thickness {t} is negative for gap {gap_index}" + )); + } + Ok(t) + } } impl Solve for MarginalRaySolve { @@ -36,6 +74,10 @@ impl Solve for MarginalRaySolve { SolveKind::Thickness } + fn path_id(&self) -> usize { + self.path_id + } + fn apply( &self, model: &SequentialModel, @@ -58,28 +100,38 @@ impl Solve for MarginalRaySolve { } let bundle = marginal_ray_bundle(model, self.wavelength_id)?; - let ray = &bundle.rays_at_surface(self.gap_index)[0]; - let h = ray.height; - let u_prime = ray.angle; + let t = Self::solved_thickness(&bundle, self.gap_index, self.target_height)?; + gap_specs[self.gap_index].thickness = t; + Ok(()) + } - let eps = Float::EPSILON * h.abs().max(1.0); - if u_prime.abs() < eps { + fn apply_multipath(&self, model: &SequentialModel, paths: &mut [PathSpec]) -> Result<()> { + let path_id = self.path_id(); + if path_id >= paths.len() { return Err(anyhow!( - "marginal ray angle at surface {} is effectively zero; \ - thickness is indeterminate (collimated beam in gap space)", - self.gap_index + "solve targets path_id {path_id} but the model has only {} path(s)", + paths.len() )); } - - let t = (self.target_height - h) / u_prime; - if t < 0.0 { + if self.gap_index >= paths[path_id].gaps.len() { + return Err(anyhow!( + "gap_index {} is out of range (path {path_id} has {} gap(s))", + self.gap_index, + paths[path_id].gaps.len() + )); + } + let n_wavelengths = model.wavelengths_for_path(path_id).len(); + if self.wavelength_id >= n_wavelengths { return Err(anyhow!( - "computed gap thickness {t} is negative for gap {}", - self.gap_index + "wavelength_id {} is out of range (path {path_id} has {} wavelength(s))", + self.wavelength_id, + n_wavelengths )); } - gap_specs[self.gap_index].thickness = t; + let bundle = marginal_ray_bundle_for_path(model, path_id, self.wavelength_id)?; + let t = Self::solved_thickness(&bundle, self.gap_index, self.target_height)?; + paths[path_id].gaps[self.gap_index].thickness = t; Ok(()) } diff --git a/crates/cherry-rs/src/core/sequential_model/solves/mod.rs b/crates/cherry-rs/src/core/sequential_model/solves/mod.rs index 9f81fdb8..d47b3fd1 100644 --- a/crates/cherry-rs/src/core/sequential_model/solves/mod.rs +++ b/crates/cherry-rs/src/core/sequential_model/solves/mod.rs @@ -1,9 +1,9 @@ pub mod fno; pub mod marginal_ray; -use anyhow::Result; +use anyhow::{Result, anyhow}; -use crate::specs::{gaps::GapSpec, surfaces::SurfaceSpec}; +use crate::specs::{gaps::GapSpec, paths::PathSpec, surfaces::SurfaceSpec}; use super::SequentialModel; @@ -58,4 +58,35 @@ pub trait Solve { fn parameter_kind(&self) -> SolveKind { SolveKind::Curvature } + + /// Which path this solve applies to in a multipath model. Ignored by + /// single-path builds (`.gap_specs()`/`.surface_specs()`) — there is + /// only ever path 0. Defaults to 0 so every existing `Solve` + /// implementation keeps compiling and behaving identically without + /// change. + fn path_id(&self) -> usize { + 0 + } + + /// Applies this solve within a multipath build + /// (`SequentialModelBuilder::paths`). + /// + /// A multipath model has no flat `gap_specs`/`surface_specs` pair — + /// `paths` is the full, mutable set of `PathSpec`s instead. A + /// `Thickness`-kind solve mutates `paths[path_id()].gaps` directly + /// (gaps are never shared between paths, so no ownership resolution is + /// needed). A `Curvature`-kind solve's [`surface_index`] is a store + /// index; the implementation must resolve two independent things from + /// it: which path's own ray trace evaluates the constraint (`path_id()`), + /// and which path actually owns the mutable `SurfaceSpec` for that store + /// index (which may be a *different* path, when the surface is + /// `Shared`) — see `locate_surface_owner`. + /// + /// The default implementation errors; override to support multipath + /// models. + /// + /// [`surface_index`]: Solve::surface_index + fn apply_multipath(&self, _model: &SequentialModel, _paths: &mut [PathSpec]) -> Result<()> { + Err(anyhow!("this solve does not support multipath models")) + } } diff --git a/crates/cherry-rs/src/examples/wf_epi_excitation.rs b/crates/cherry-rs/src/examples/wf_epi_excitation.rs index 92195132..b9fdd3fb 100644 --- a/crates/cherry-rs/src/examples/wf_epi_excitation.rs +++ b/crates/cherry-rs/src/examples/wf_epi_excitation.rs @@ -16,7 +16,10 @@ use std::rc::Rc; use crate::{ BeamSplitterPathKind, EulerAngles, GapSpec, PathSpec, PathSurfaceRef, RefractiveIndexSpec, Rotation3D, SequentialModel, SurfaceSpec, Vec3, - core::{Float, sequential_model::builder::SequentialModelBuilder}, + core::{ + Float, + sequential_model::{builder::SequentialModelBuilder, solves::MarginalRaySolve}, + }, }; pub fn sequential_model( @@ -36,8 +39,10 @@ pub fn sequential_model( thickness: 50.0, refractive_index: n_air, }; + // Placeholder thickness; solved below to place the image plane at the + // paraxial focus (MarginalRaySolve targeting gap 3, target_height = 0.0). let gap_3 = GapSpec { - thickness: 5.0, + thickness: 1.0, refractive_index: n_oil, }; @@ -82,6 +87,7 @@ pub fn sequential_model( SequentialModelBuilder::new() .paths(vec![path]) + .solves(vec![Box::new(MarginalRaySolve::new(3, 0.0, 0))]) .build() .expect("wf_epi_excitation model builds") .model diff --git a/crates/cherry-rs/src/views/paraxial.rs b/crates/cherry-rs/src/views/paraxial.rs index b461556d..1c3c4e03 100644 --- a/crates/cherry-rs/src/views/paraxial.rs +++ b/crates/cherry-rs/src/views/paraxial.rs @@ -1464,22 +1464,24 @@ fn resolve_stop_surface_step(store_index: usize, surface_indices: &[usize]) -> R } } -/// Compute the paraxial marginal ray bundle for a given wavelength. +/// Compute the paraxial marginal ray bundle for a given path and wavelength. /// /// Uses the first tangential direction `(0, 1, 0)`, valid for all rotationally /// symmetric systems. -pub(crate) fn marginal_ray_bundle( +pub(crate) fn marginal_ray_bundle_for_path( model: &SequentialModel, + path_id: usize, wavelength_id: usize, ) -> Result { let submodel = model - .submodel(wavelength_id) - .ok_or_else(|| anyhow!("wavelength_id {wavelength_id} out of range"))?; + .submodels_for_path(path_id) + .get(wavelength_id) + .ok_or_else(|| anyhow!("wavelength_id {wavelength_id} out of range for path {path_id}"))?; 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 path_steps = model.path_steps(path_id); + let surface_indices = model.path_surface_indices(path_id); + let beam_splitter_arms = model.path_beam_splitter_arms(path_id); let v = Vec3::new(0.0, 1.0, 0.0); let per_surf_v = propagate_tangential_vec(v, surfaces, placements, surface_indices); @@ -1491,7 +1493,7 @@ pub(crate) fn marginal_ray_bundle( beam_splitter_arms, path_steps, )?; - let stop = match model.stop_surface() { + let stop = match model.stop_surface_for_path(path_id) { Some(i) => resolve_stop_surface_step(i, surface_indices)?, None => calc_aperture_stop( surfaces, @@ -1513,6 +1515,14 @@ pub(crate) fn marginal_ray_bundle( )) } +/// Single-path shorthand; delegates to path 0. +pub(crate) fn marginal_ray_bundle( + model: &SequentialModel, + wavelength_id: usize, +) -> Result { + marginal_ray_bundle_for_path(model, 0, wavelength_id) +} + fn argmin(ratios: &[Float]) -> usize { ratios .iter()