From d57633157ebc6398c3836f0006adc34049b5d2f5 Mon Sep 17 00:00:00 2001 From: Enrico Sita Date: Wed, 21 Oct 2026 10:08:33 +0200 Subject: [PATCH] M2: calculation engine (grids, tilt, vertical illuminance, statistics) Luminaires with heading and tilt (IES tilt factors applied), direct illuminance on any surface, area and full-width route grids summed over luminaires, min/mean/max and uniformity statistics that are unverifiable when any point is; photometry plausibility warnings. Isotropic and Lambertian cases within 1e-6; 100 luminaires x 10,000 points in ~55 ms. Co-Authored-By: Claude Opus 5.5 --- ROADMAP.md | 5 +- apps/luxpath-cli/src/main.rs | 4 + apps/luxpath-cli/tests/photometry.rs | 14 ++ crates/calc/src/grid.rs | 210 ++++++++++++++++++++++++ crates/calc/src/lib.rs | 236 ++++++++++++++++----------- crates/calc/src/stats.rs | 65 ++++++++ crates/calc/src/tests.rs | 176 ++++++++++++++++++++ crates/calc/tests/performance.rs | 43 +++++ crates/geom/src/lib.rs | 115 +++++++++++++ crates/photometry/src/ldt.rs | 6 +- crates/photometry/src/lib.rs | 7 +- crates/photometry/src/model.rs | 77 +++++++++ docs/gaps.md | 7 +- docs/roadmap.md | 4 +- docs/status.md | 15 +- parity/features.toml | 22 ++- xtask/src/main.rs | 5 +- 17 files changed, 892 insertions(+), 119 deletions(-) create mode 100644 crates/calc/src/grid.rs create mode 100644 crates/calc/src/stats.rs create mode 100644 crates/calc/src/tests.rs create mode 100644 crates/calc/tests/performance.rs diff --git a/ROADMAP.md b/ROADMAP.md index 4620888..17f67b8 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -2,7 +2,7 @@ **Stage: pre-alpha** · next: alpha (M1–M4 + M6 validation) · ~150–250 agent-assisted hours away -> **Last updated:** 2026-10-21 · **Change:** M1 (photometry readers) done +> **Last updated:** 2026-10-21 · **Change:** M2 (engine) done luxpath's one-page summary. Details live in: @@ -24,7 +24,7 @@ luxpath's one-page summary. Details live in: ## Next up -1. **M2:** grids, sums over luminaires, tilt, vertical illuminance, photometry plausibility warnings (gap 13). +1. **M3:** the scene format and `luxpath check`; obstruction by walls (gap 16) before M4. 2. **Real photometric files** from the owner's practice, kept in `plan/` (gap 12). 3. **Open question:** buy EN 1838:2024 (around €100 from a national standards body) to verify the requirement values. @@ -32,6 +32,7 @@ luxpath's one-page summary. Details live in: | Date | Change | |---|---| +| 2026-10-21 | M2: luminaires with heading and tilt (IES tilt factors applied), illuminance on any surface, area and full-width route grids summed over luminaires, statistics that are unverifiable when any point is; isotropic and Lambertian cases within 1e-6; 100 luminaires × 10,000 points in about 55 ms (release, in `cargo xtask ci`); photometry plausibility warnings. | | 2026-10-21 | M1: IES LM-63 (1986–2019) and EULUMDAT readers, symmetry expansion, tilt data, emergency flux, located errors, property tests, a 10-minute fuzz run (22.8 million inputs) without a panic; `luxpath photometry`. Existing crates evaluated first: none usable (docs/architecture.md). | | 2026-10-16 | M7 (early): `luxpath-ifc` and `luxpath inspect-ifc`; 9 public models inspected, findings in docs/research.md §7. | | 2026-10-16 | M0: workspace, layering and parity gates, C-γ intensity table, direct illuminance, demo CLI, docs. | diff --git a/apps/luxpath-cli/src/main.rs b/apps/luxpath-cli/src/main.rs index be05467..ecd2b8d 100644 --- a/apps/luxpath-cli/src/main.rs +++ b/apps/luxpath-cli/src/main.rs @@ -150,6 +150,9 @@ fn describe(p: &Photometry) -> String { if let Some(tilt) = p.tilt() { line("tilt", format!("{} angles, lamp geometry {}", tilt.angles_deg().len(), tilt.lamp_geometry)); } + for warning in p.warnings() { + line("warning", warning.to_string()); + } out } @@ -171,6 +174,7 @@ fn point(args: &[String]) -> Result { let luminaire = Luminaire { position: Vec3::new(0.0, 0.0, height), rotation_deg: 0.0, + tilt_deg: 0.0, distribution, flux_lm: flux, maintenance_factor: 1.0, diff --git a/apps/luxpath-cli/tests/photometry.rs b/apps/luxpath-cli/tests/photometry.rs index 1b362da..8d163f1 100644 --- a/apps/luxpath-cli/tests/photometry.rs +++ b/apps/luxpath-cli/tests/photometry.rs @@ -24,9 +24,23 @@ fn describes_ies_and_eulumdat_files_alike() { assert!(out.contains("1000.0 lm (lamp flux)"), "{out}"); assert!(out.contains("rotational"), "{out}"); assert!(out.contains("300.0 cd/klm at C0° γ0°"), "{out}"); + assert!(!out.contains("warning"), "{out}"); } } +#[test] +fn implausible_files_are_described_with_a_warning() { + let dir = std::env::temp_dir().join(format!("luxpath-cli-warn-{}", std::process::id())); + std::fs::create_dir_all(&dir).unwrap(); + let path = dir.join("bright.ies"); + // 1000 cd all round from a 1000 lm lamp: 4π × 1000 lm, 1256.6 % of the lamp flux. + std::fs::write(&path, "IESNA:LM-63-2002\nTILT=NONE\n1 1000 1 2 1 1 2 0 0 0\n1 1 0\n0 180\n0\n1000 1000\n").unwrap(); + let (code, out, _) = luxpath(&["photometry", path.to_str().unwrap()]); + assert_eq!(code, 0); + assert!(out.contains("warning") && out.contains("1256.6 %"), "{out}"); + std::fs::remove_dir_all(&dir).ok(); +} + #[test] fn a_broken_file_is_an_error_with_its_line() { let dir = std::env::temp_dir().join(format!("luxpath-cli-{}", std::process::id())); diff --git a/crates/calc/src/grid.rs b/crates/calc/src/grid.rs new file mode 100644 index 0000000..447e9d5 --- /dev/null +++ b/crates/calc/src/grid.rs @@ -0,0 +1,210 @@ +//! Calculation grids: points where illuminance is computed, all on surfaces with one normal. + +use luxpath_geom::{Vec2, Vec3, distance_to_boundary, point_in_polygon}; + +/// Grids larger than this are refused (a 100 m × 100 m floor at 0.1 m is a million points). +pub const MAX_GRID_POINTS: usize = 1_000_000; + +#[derive(Debug, Clone, PartialEq)] +pub struct Grid { + points: Vec, + normal: Vec3, +} + +#[derive(Debug, Clone, PartialEq)] +pub enum GridError { + /// A spacing, width, border or height is not a finite number in range. + BadDimension(&'static str), + /// The polygon has fewer than three vertices or a non-finite one; the route fewer than two. + BadOutline, + /// More than [`MAX_GRID_POINTS`]. + TooManyPoints, + /// No point is left (the area is smaller than its border, or the spacing): nothing to check. + Empty, +} + +impl std::fmt::Display for GridError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + Self::BadDimension(what) => write!(f, "the grid's {what} is not a usable number"), + Self::BadOutline => write!(f, "the outline needs at least three finite corners (two for a route)"), + Self::TooManyPoints => write!(f, "the grid would have more than {MAX_GRID_POINTS} points"), + Self::Empty => write!(f, "the grid has no points: the area is smaller than its border or spacing"), + } + } +} + +impl std::error::Error for GridError {} + +fn positive(value: f64, what: &'static str) -> Result { + if value.is_finite() && value > 0.0 { Ok(value) } else { Err(GridError::BadDimension(what)) } +} + +impl Grid { + /// Points given one by one, on surfaces facing `normal` (vertical illuminance at a fire + /// extinguisher: the point on the equipment, the normal towards the room). + pub fn new(points: Vec, normal: Vec3) -> Result { + if points.len() > MAX_GRID_POINTS { + return Err(GridError::TooManyPoints); + } + if points.is_empty() { + return Err(GridError::Empty); + } + if !points.iter().all(|p| p.is_finite()) || normal.normalized().is_none() { + return Err(GridError::BadDimension("point or normal")); + } + Ok(Self { points, normal }) + } + + /// A horizontal grid facing up over the `polygon` at `height_m`: the centres of `spacing_m` + /// cells across its bounding box that lie inside it and at least `border_m` from its edges. + pub fn area(polygon: &[Vec2], height_m: f64, spacing_m: f64, border_m: f64) -> Result { + let spacing = positive(spacing_m, "spacing")?; + if !height_m.is_finite() { + return Err(GridError::BadDimension("height")); + } + if !border_m.is_finite() || border_m < 0.0 { + return Err(GridError::BadDimension("border")); + } + if polygon.len() < 3 || !polygon.iter().all(|p| p.is_finite()) { + return Err(GridError::BadOutline); + } + let (min, max) = polygon.iter().fold( + (Vec2::new(f64::INFINITY, f64::INFINITY), Vec2::new(f64::NEG_INFINITY, f64::NEG_INFINITY)), + |(lo, hi), p| (Vec2::new(lo.x.min(p.x), lo.y.min(p.y)), Vec2::new(hi.x.max(p.x), hi.y.max(p.y))), + ); + let (nx, ny) = (((max.x - min.x) / spacing).ceil(), ((max.y - min.y) / spacing).ceil()); + if !(nx * ny).is_finite() || nx * ny > MAX_GRID_POINTS as f64 { + return Err(GridError::TooManyPoints); + } + // Finite and at most MAX_GRID_POINTS: the casts are exact. + let (nx, ny) = (nx as usize, ny as usize); + let mut points = Vec::new(); + for j in 0..ny { + for i in 0..nx { + let p = Vec2::new(min.x + (i as f64 + 0.5) * spacing, min.y + (j as f64 + 0.5) * spacing); + if point_in_polygon(p, polygon) && distance_to_boundary(p, polygon) >= border_m { + points.push(p.at(height_m)); + } + } + } + Self::new(points, Vec3::UP) + } + + /// A horizontal grid facing up along a route: its centre line `centre` at `height_m`, sampled + /// every `spacing_m` along each segment (both ends included), with `across` points evenly over + /// the full `width_m` (edges included; one point is the centre line). + pub fn route( + centre: &[Vec2], + width_m: f64, + height_m: f64, + spacing_m: f64, + across: usize, + ) -> Result { + let spacing = positive(spacing_m, "spacing")?; + let width = positive(width_m, "width")?; + if !height_m.is_finite() { + return Err(GridError::BadDimension("height")); + } + if across == 0 || across > MAX_GRID_POINTS { + return Err(GridError::BadDimension("number of points across")); + } + if centre.len() < 2 || !centre.iter().all(|p| p.is_finite()) { + return Err(GridError::BadOutline); + } + let offsets: Vec = if across == 1 { + vec![0.0] + } else { + (0..across).map(|k| -width / 2.0 + width * k as f64 / (across - 1) as f64).collect() + }; + let mut points = Vec::new(); + for (index, pair) in centre.windows(2).enumerate() { + let (Some(a), Some(b)) = (pair.first(), pair.get(1)) else { continue }; + let along = *b - *a; + let length = along.length(); + if length == 0.0 { + continue; + } + let (ux, uy) = (along.x / length, along.y / length); + let steps = (length / spacing).ceil(); + if !steps.is_finite() || steps * across as f64 + points.len() as f64 > MAX_GRID_POINTS as f64 { + return Err(GridError::TooManyPoints); + } + // Finite and bounded by MAX_GRID_POINTS: the cast is exact. + let steps = steps as usize; + // Each segment's start is the previous segment's end: include it only once. + let first = if index == 0 { 0 } else { 1 }; + for s in first..=steps { + let t = length * s as f64 / steps as f64; + for offset in &offsets { + // The perpendicular to the left of the direction of travel is (−uy, ux). + let p = Vec2::new(a.x + ux * t - uy * offset, a.y + uy * t + ux * offset); + points.push(p.at(height_m)); + } + } + } + Self::new(points, Vec3::UP) + } + + pub fn points(&self) -> &[Vec3] { + &self.points + } + + pub fn normal(&self) -> Vec3 { + self.normal + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn square(side: f64) -> Vec { + vec![Vec2::new(0.0, 0.0), Vec2::new(side, 0.0), Vec2::new(side, side), Vec2::new(0.0, side)] + } + + #[test] + fn an_area_grid_takes_cell_centres_inside_less_the_border() { + // A 4 m square at 1 m spacing: 16 centres at 0.5, 1.5, 2.5, 3.5. + let g = Grid::area(&square(4.0), 0.0, 1.0, 0.0).unwrap(); + assert_eq!(g.points().len(), 16); + assert_eq!(g.points()[0], Vec3::new(0.5, 0.5, 0.0)); + // A 0.5 m border keeps every centre (they are exactly 0.5 m in); 0.6 m keeps the inner 4. + assert_eq!(Grid::area(&square(4.0), 0.0, 1.0, 0.5).unwrap().points().len(), 16); + assert_eq!(Grid::area(&square(4.0), 0.0, 1.0, 0.6).unwrap().points().len(), 4); + assert_eq!(Grid::area(&square(1.0), 0.0, 1.0, 0.6), Err(GridError::Empty)); + } + + #[test] + fn a_route_grid_covers_the_full_width_and_each_corner_once() { + // An L-shaped route, 2 m then 2 m, 1 m wide, every 1 m, three points across. + let route = [Vec2::new(0.0, 0.0), Vec2::new(2.0, 0.0), Vec2::new(2.0, 2.0)]; + let g = Grid::route(&route, 1.0, 0.0, 1.0, 3).unwrap(); + // Stations: 0, 1, 2 on the first leg, then 1, 2 on the second: five, times three across. + assert_eq!(g.points().len(), 15); + // The first station spans y = −0.5 to 0.5 (left of travel along +x is +y). + assert_eq!(&g.points()[..3], &[Vec3::new(0.0, -0.5, 0.0), Vec3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.5, 0.0)]); + // On the second leg (along +y), left is −x. + assert_eq!(g.points()[14], Vec3::new(1.5, 2.0, 0.0)); + assert_eq!(Grid::route(&route, 1.0, 0.0, 1.0, 1).unwrap().points().len(), 5); + } + + #[test] + fn bad_dimensions_and_huge_grids_are_refused() { + assert_eq!(Grid::area(&square(4.0), 0.0, 0.0, 0.0), Err(GridError::BadDimension("spacing"))); + assert_eq!(Grid::area(&square(4.0), 0.0, f64::NAN, 0.0), Err(GridError::BadDimension("spacing"))); + assert_eq!(Grid::area(&square(4.0), 0.0, 1.0, -1.0), Err(GridError::BadDimension("border"))); + assert_eq!(Grid::area(&square(4.0)[..2], 0.0, 1.0, 0.0), Err(GridError::BadOutline)); + assert_eq!(Grid::area(&square(1e6), 0.0, 0.01, 0.0), Err(GridError::TooManyPoints)); + assert_eq!(Grid::route(&[Vec2::default()], 1.0, 0.0, 1.0, 3), Err(GridError::BadOutline)); + assert_eq!( + Grid::route(&[Vec2::default(), Vec2::new(1e9, 0.0)], 1.0, 0.0, 1e-3, 3), + Err(GridError::TooManyPoints) + ); + assert_eq!( + Grid::route(&[Vec2::default(), Vec2::new(1.0, 0.0)], 1.0, 0.0, 1.0, 0), + Err(GridError::BadDimension("number of points across")) + ); + assert_eq!(Grid::new(vec![Vec3::default()], Vec3::default()), Err(GridError::BadDimension("point or normal"))); + } +} diff --git a/crates/calc/src/lib.rs b/crates/calc/src/lib.rs index 9457744..aff95b9 100644 --- a/crates/calc/src/lib.rs +++ b/crates/calc/src/lib.rs @@ -1,136 +1,186 @@ //! Illuminance calculations. //! -//! M0 has the core of the engine: **direct** horizontal illuminance at a point from one -//! luminaire, by the inverse-square cosine law: +//! The engine computes **direct** illuminance on a surface with normal `n` at a point, from each +//! luminaire, by the inverse-square cosine law, and sums it over luminaires: //! //! ```text //! E = I(C, γ) · cos θ / d² //! ``` //! //! where `I` is the luminaire's intensity towards the point (cd), `d` the distance (m) and `θ` the -//! angle between the light's direction and the surface normal. Reflected light is not included: -//! emergency-lighting tools leave it out by default (see `docs/research.md`), and inter-reflection -//! is a later milestone (M9). +//! angle between the arriving light and `n` (light from behind the surface counts 0). Reflected +//! light is not included: emergency-lighting tools leave it out by default (see `docs/research.md`), +//! and inter-reflection is a later milestone (M9). //! -//! Next (M2): luminaire tilt, grids of points, sums over many luminaires, vertical planes. +//! - [`Luminaire`]: a placed, oriented (heading and tilt) distribution with its flux; +//! [`Luminaire::from_photometry`] applies the file's tilt factors. +//! - [`Grid`]: points with one surface normal: over an area (a polygon, less a border), along a +//! route (its full width), or given one by one (vertical illuminance at equipment). +//! - [`illuminance`] at every grid point, summed over luminaires; [`Stats`] (min, mean, max, the +//! uniformity ratios). +//! +//! An undefined result (a point at a luminaire's centre, a non-finite input) is `None` for that +//! point, and statistics over it are `None`: a check reports "cannot verify", never a value. //! //! Layer L3: depends on `luxpath-geom` and `luxpath-photometry` (and later `luxpath-scene`). +mod grid; +mod stats; + +pub use grid::{Grid, GridError, MAX_GRID_POINTS}; +pub use stats::Stats; + use luxpath_geom::{Vec3, wrap_degrees}; -use luxpath_photometry::IntensityTable; +use luxpath_photometry::{IntensityTable, Photometry}; /// One luminaire as the calculation sees it. /// -/// M0 supports luminaires pointing straight down (the common case for emergency luminaires), -/// turned about the vertical axis by `rotation_deg` (C0 then points along that heading from +x). +/// Orientation: the luminaire's photometric frame (C0 along its x axis, C90 along y, γ0 straight +/// down its −z axis) is first **tilted** by `tilt_deg` about its own y axis, so that the nadir swings +/// towards C0 (90° aims the light horizontally, as for a wall-mounted luminaire whose photometry +/// was measured pointing down), then turned by `rotation_deg` about the vertical (C0 then points +/// along that heading, counter-clockwise from +x). #[derive(Debug, Clone)] pub struct Luminaire { pub position: Vec3, pub rotation_deg: f64, + pub tilt_deg: f64, pub distribution: IntensityTable, - /// Luminous flux in emergency mode, in lumens (often far below the mains-mode flux). + /// Luminous flux in the mode being checked, in lumens (in emergency mode often far below the + /// mains-mode flux). pub flux_lm: f64, /// Maintenance factor (0–1] applied to the intensity. pub maintenance_factor: f64, } -/// Direct illuminance (lux) on a horizontal surface facing up at `point`. +/// Why a luminaire cannot be placed. +#[derive(Debug, Clone, PartialEq)] +pub enum CalcError { + /// The file's tilt table does not cover this tilt, so the light output there is unknown. + TiltOutsideTable { tilt_deg: f64 }, + /// A factor scaled the table out of range. + Photometry(luxpath_photometry::PhotometryError), +} + +impl std::fmt::Display for CalcError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + Self::TiltOutsideTable { tilt_deg } => { + write!(f, "the photometric file's tilt table does not cover a tilt of {tilt_deg}°") + } + Self::Photometry(e) => write!(f, "{e}"), + } + } +} + +impl std::error::Error for CalcError {} + +impl Luminaire { + /// A luminaire from a file's photometry, with the file's tilt factor (IES `TILT=INCLUDE`) for + /// `tilt_deg` applied to the distribution. A tilt the table does not cover is an error, not a + /// guess. + pub fn from_photometry( + photometry: &Photometry, + position: Vec3, + rotation_deg: f64, + tilt_deg: f64, + flux_lm: f64, + maintenance_factor: f64, + ) -> Result { + let distribution = match photometry.tilt() { + None => photometry.table().clone(), + Some(tilt) => { + let factor = tilt.factor(tilt_deg).ok_or(CalcError::TiltOutsideTable { tilt_deg })?; + photometry.table().scaled(factor).map_err(CalcError::Photometry)? + } + }; + Ok(Self { position, rotation_deg, tilt_deg, distribution, flux_lm, maintenance_factor }) + } + + /// The precomputed numbers a calculation needs, or `None` if an input is not finite. + fn frame(&self) -> Option { + let ok = self.position.is_finite() + && self.rotation_deg.is_finite() + && self.tilt_deg.is_finite() + && self.flux_lm.is_finite() + && self.maintenance_factor.is_finite(); + let (sin_r, cos_r) = self.rotation_deg.to_radians().sin_cos(); + let (sin_t, cos_t) = self.tilt_deg.to_radians().sin_cos(); + ok.then_some(Frame { sin_r, cos_r, sin_t, cos_t, scale: self.flux_lm / 1000.0 * self.maintenance_factor }) + } +} + +struct Frame { + sin_r: f64, + cos_r: f64, + sin_t: f64, + cos_t: f64, + /// cd per cd/klm: flux / 1000 × maintenance factor. + scale: f64, +} + +impl Frame { + /// A world direction in the luminaire's own frame: undo the heading (about z), then the tilt + /// (about the luminaire's y axis). + fn to_local(&self, d: Vec3) -> Vec3 { + // Rz(−rotation) + let x = d.x * self.cos_r + d.y * self.sin_r; + let y = -d.x * self.sin_r + d.y * self.cos_r; + // Tilting by t maps the nadir (0, 0, −1) to (sin t, 0, −cos t); undo it with Ry(t): + // x' = x cos t + z sin t, z' = −x sin t + z cos t. + Vec3::new(x * self.cos_t + d.z * self.sin_t, y, -x * self.sin_t + d.z * self.cos_t) + } +} + +/// Direct illuminance (lux) at `point` on a surface facing `normal` (need not be unit length). /// -/// Returns `None` when the result is undefined: the point coincides with the luminaire, or an -/// input is not finite. Points above the luminaire's light (θ ≥ 90°) receive 0 lux. +/// `None` when the result is undefined: the point coincides with the luminaire, the normal is zero, +/// or an input is not finite. Light arriving from behind the surface gives 0 lux. +pub fn direct_illuminance(luminaire: &Luminaire, point: Vec3, normal: Vec3) -> Option { + let normal = normal.normalized()?; + direct(luminaire, &luminaire.frame()?, point, normal) +} + +/// Direct illuminance (lux) on a horizontal surface facing up at `point`. pub fn direct_horizontal_illuminance(luminaire: &Luminaire, point: Vec3) -> Option { - if !point.is_finite() - || !luminaire.position.is_finite() - || !luminaire.flux_lm.is_finite() - || !luminaire.maintenance_factor.is_finite() - { + direct_illuminance(luminaire, point, Vec3::UP) +} + +fn direct(luminaire: &Luminaire, frame: &Frame, point: Vec3, unit_normal: Vec3) -> Option { + if !point.is_finite() { return None; } let to_point = point - luminaire.position; let d2 = to_point.length_squared(); let dir = to_point.normalized()?; - // Light reaches an upward-facing surface only when it travels downwards. - let cos_theta = -dir.z; + // The light travels along `dir`; it lights the surface when it arrives against the normal. + let cos_theta = -dir.dot(unit_normal); if cos_theta <= 0.0 { return Some(0.0); } - // γ is measured from the nadir, which for a downward-pointing luminaire is the same angle as θ. - let gamma = cos_theta.clamp(-1.0, 1.0).acos().to_degrees(); - let c = if dir.x == 0.0 && dir.y == 0.0 { - 0.0 - } else { - wrap_degrees(dir.y.atan2(dir.x).to_degrees() - luminaire.rotation_deg) - }; - let candela = - luminaire.distribution.intensity(c, gamma) * luminaire.flux_lm / 1000.0 * luminaire.maintenance_factor; + let local = frame.to_local(dir); + let gamma = (-local.z).clamp(-1.0, 1.0).acos().to_degrees(); + let c = if local.x == 0.0 && local.y == 0.0 { 0.0 } else { wrap_degrees(local.y.atan2(local.x).to_degrees()) }; + let candela = luminaire.distribution.intensity(c, gamma) * frame.scale; let e = candela * cos_theta / d2; e.is_finite().then_some(e) } -#[cfg(test)] -mod tests { - use super::*; - - fn isotropic_at(height: f64) -> Luminaire { - Luminaire { - position: Vec3::new(0.0, 0.0, height), - rotation_deg: 0.0, - distribution: IntensityTable::isotropic(1000.0).unwrap(), - flux_lm: 1000.0, - maintenance_factor: 1.0, +/// Illuminance (lux) at every point of `grid`, summed over `luminaires`; `None` at a point where +/// any luminaire's contribution is undefined. +pub fn illuminance(luminaires: &[Luminaire], grid: &Grid) -> Vec> { + let mut sums: Vec> = vec![Some(0.0); grid.points().len()]; + let Some(normal) = grid.normal().normalized() else { return vec![None; grid.points().len()] }; + for luminaire in luminaires { + let Some(frame) = luminaire.frame() else { return vec![None; grid.points().len()] }; + for (sum, point) in sums.iter_mut().zip(grid.points()) { + if let Some(total) = sum { + *sum = direct(luminaire, &frame, *point, normal).map(|e| *total + e); + } } } - - #[test] - fn directly_below_is_intensity_over_distance_squared() { - // 1000 cd at 2 m: E = 1000 / 4 = 250 lx. - let e = direct_horizontal_illuminance(&isotropic_at(2.0), Vec3::new(0.0, 0.0, 0.0)).unwrap(); - assert!((e - 250.0).abs() < 1e-9); - } - - #[test] - fn off_axis_points_follow_the_cosine_law() { - // Point 2 m across, 2 m down: d² = 8, cos θ = 1/√2 → E = 1000 · 0.7071 / 8. - let e = direct_horizontal_illuminance(&isotropic_at(2.0), Vec3::new(2.0, 0.0, 0.0)).unwrap(); - assert!((e - 1000.0 * std::f64::consts::FRAC_1_SQRT_2 / 8.0).abs() < 1e-9); - } - - #[test] - fn flux_and_maintenance_factor_scale_linearly() { - let mut lum = isotropic_at(2.0); - lum.flux_lm = 200.0; - lum.maintenance_factor = 0.8; - let e = direct_horizontal_illuminance(&lum, Vec3::default()).unwrap(); - assert!((e - 250.0 * 0.2 * 0.8).abs() < 1e-9); - } - - #[test] - fn rotation_selects_the_c_plane() { - // C0 = 200 cd/klm, C180 = 100 cd/klm, flat in γ. - let table = IntensityTable::new(vec![0.0, 180.0], vec![0.0, 180.0], vec![200.0, 200.0, 100.0, 100.0]).unwrap(); - let mut lum = isotropic_at(1.0); - lum.distribution = table; - let east = Vec3::new(1.0, 0.0, 0.0); - let unrotated = direct_horizontal_illuminance(&lum, east).unwrap(); - lum.rotation_deg = 180.0; - let rotated = direct_horizontal_illuminance(&lum, east).unwrap(); - assert!((unrotated - 2.0 * rotated).abs() < 1e-9); - } - - #[test] - fn points_level_with_or_above_the_luminaire_get_nothing() { - let lum = isotropic_at(2.0); - assert_eq!(direct_horizontal_illuminance(&lum, Vec3::new(5.0, 0.0, 2.0)), Some(0.0)); - assert_eq!(direct_horizontal_illuminance(&lum, Vec3::new(0.0, 0.0, 3.0)), Some(0.0)); - } - - #[test] - fn degenerate_inputs_are_none_not_a_panic() { - let lum = isotropic_at(2.0); - assert_eq!(direct_horizontal_illuminance(&lum, Vec3::new(0.0, 0.0, 2.0)), None); - assert_eq!(direct_horizontal_illuminance(&lum, Vec3::new(f64::NAN, 0.0, 0.0)), None); - let mut bad = isotropic_at(2.0); - bad.flux_lm = f64::INFINITY; - assert_eq!(direct_horizontal_illuminance(&bad, Vec3::default()), None); - } + sums } + +#[cfg(test)] +mod tests; diff --git a/crates/calc/src/stats.rs b/crates/calc/src/stats.rs new file mode 100644 index 0000000..f9ce9e6 --- /dev/null +++ b/crates/calc/src/stats.rs @@ -0,0 +1,65 @@ +//! Statistics over a grid's illuminances. + +/// Minimum, mean and maximum illuminance (lux) over a grid's points. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Stats { + pub min_lx: f64, + pub mean_lx: f64, + pub max_lx: f64, + pub points: usize, +} + +impl Stats { + /// The statistics of `values`; `None` if there are none or any is undefined (a check over them + /// cannot be verified). + pub fn of(values: &[Option]) -> Option { + let mut stats = Self { min_lx: f64::INFINITY, mean_lx: 0.0, max_lx: f64::NEG_INFINITY, points: 0 }; + let mut sum = 0.0; + for v in values { + let v = (*v).filter(|v| v.is_finite())?; + stats.min_lx = stats.min_lx.min(v); + stats.max_lx = stats.max_lx.max(v); + sum += v; + stats.points += 1; + } + if stats.points == 0 { + return None; + } + stats.mean_lx = sum / stats.points as f64; + Some(stats) + } + + /// Max over min (EN 1838 limits it, e.g. 40:1); `None` when the minimum is 0 (unbounded). + pub fn max_to_min(&self) -> Option { + (self.min_lx > 0.0).then(|| self.max_lx / self.min_lx) + } + + /// Min over mean (U₀); `None` when the mean is 0. + pub fn min_to_mean(&self) -> Option { + (self.mean_lx > 0.0).then(|| self.min_lx / self.mean_lx) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn statistics_by_hand() { + // 1, 2, 3, 6: min 1, mean 3, max 6; max/min 6, min/mean 1/3. + let s = Stats::of(&[Some(1.0), Some(2.0), Some(3.0), Some(6.0)]).unwrap(); + assert_eq!((s.min_lx, s.mean_lx, s.max_lx, s.points), (1.0, 3.0, 6.0, 4)); + assert_eq!(s.max_to_min(), Some(6.0)); + assert_eq!(s.min_to_mean(), Some(1.0 / 3.0)); + } + + #[test] + fn undefined_or_empty_values_give_no_statistics() { + assert_eq!(Stats::of(&[]), None); + assert_eq!(Stats::of(&[Some(1.0), None]), None); + assert_eq!(Stats::of(&[Some(f64::NAN)]), None); + let dark = Stats::of(&[Some(0.0), Some(2.0)]).unwrap(); + assert_eq!(dark.max_to_min(), None); + assert_eq!(Stats::of(&[Some(0.0)]).unwrap().min_to_mean(), None); + } +} diff --git a/crates/calc/src/tests.rs b/crates/calc/src/tests.rs new file mode 100644 index 0000000..f305c90 --- /dev/null +++ b/crates/calc/src/tests.rs @@ -0,0 +1,176 @@ +use super::*; +use luxpath_geom::Vec2; + +fn isotropic_at(height: f64) -> Luminaire { + Luminaire { + position: Vec3::new(0.0, 0.0, height), + rotation_deg: 0.0, + tilt_deg: 0.0, + distribution: IntensityTable::isotropic(1000.0).unwrap(), + flux_lm: 1000.0, + maintenance_factor: 1.0, + } +} + +/// I = 1000 cd/klm × cos γ below the horizon, none above; γ every 0.1° so that linear +/// interpolation stays within 4e-7 of the cosine (error ≤ Δγ²/8 = (0.1° in rad)²/8). +fn lambertian_at(height: f64) -> Luminaire { + let gammas: Vec = (0..=1800).map(|i| f64::from(i) / 10.0).collect(); + let values = gammas.iter().map(|g| if *g <= 90.0 { 1000.0 * g.to_radians().cos() } else { 0.0 }).collect(); + Luminaire { distribution: IntensityTable::new(vec![0.0], gammas, values).unwrap(), ..isotropic_at(height) } +} + +fn close(a: f64, b: f64) -> bool { + (a - b).abs() <= 1e-6 * b.abs().max(1e-12) +} + +#[test] +fn directly_below_is_intensity_over_distance_squared() { + // 1000 cd at 2 m: E = 1000 / 4 = 250 lx. + let e = direct_horizontal_illuminance(&isotropic_at(2.0), Vec3::new(0.0, 0.0, 0.0)).unwrap(); + assert!((e - 250.0).abs() < 1e-9); +} + +#[test] +fn off_axis_points_follow_the_cosine_law() { + // Point 2 m across, 2 m down: d² = 8, cos θ = 1/√2 → E = 1000 · 0.7071 / 8. + let e = direct_horizontal_illuminance(&isotropic_at(2.0), Vec3::new(2.0, 0.0, 0.0)).unwrap(); + assert!((e - 1000.0 * std::f64::consts::FRAC_1_SQRT_2 / 8.0).abs() < 1e-9); +} + +#[test] +fn a_lambertian_source_matches_its_closed_form() { + // I = I₀ cos γ, and on the floor cos γ = cos θ = h/d: E = I₀ h² / d⁴. + let lum = lambertian_at(2.0); + for r in [0.0, 0.5, 1.5, 3.0, 8.0] { + let d2: f64 = 4.0 + r * r; + let expected = 1000.0 * 4.0 / (d2 * d2); + let e = direct_horizontal_illuminance(&lum, Vec3::new(r, 0.0, 0.0)).unwrap(); + assert!(close(e, expected), "r {r}: {e} against {expected}"); + } +} + +#[test] +fn vertical_illuminance_follows_the_cosine_law_on_the_wall() { + // Isotropic 1000 cd at 2 m; a wall surface at (3, 0, 0) facing back towards the luminaire + // (normal −x): cos θ = 3/d, so E = 1000 · 3 / d³ with d = √13. + let lum = isotropic_at(2.0); + let d = 13f64.sqrt(); + let e = direct_illuminance(&lum, Vec3::new(3.0, 0.0, 0.0), Vec3::new(-1.0, 0.0, 0.0)).unwrap(); + assert!(close(e, 3000.0 / (d * d * d))); + // Facing away, it gets nothing. + assert_eq!(direct_illuminance(&lum, Vec3::new(3.0, 0.0, 0.0), Vec3::new(1.0, 0.0, 0.0)), Some(0.0)); + assert_eq!(direct_illuminance(&lum, Vec3::new(3.0, 0.0, 0.0), Vec3::default()), None); +} + +/// 1000 cd/klm over the lower hemisphere of its own frame, none above. +fn downlight() -> IntensityTable { + IntensityTable::new(vec![0.0], vec![0.0, 90.0], vec![1000.0, 1000.0]).unwrap() +} + +#[test] +fn a_tilted_luminaire_aims_its_light_towards_c0() { + // Tilted 90°: its nadir points along +x. A wall 2 m ahead facing back gets I/d² = 1000/4; + // behind it, nothing (its upper hemisphere is dark). + let lum = Luminaire { tilt_deg: 90.0, distribution: downlight(), ..isotropic_at(2.0) }; + let ahead = direct_illuminance(&lum, Vec3::new(2.0, 0.0, 2.0), Vec3::new(-1.0, 0.0, 0.0)).unwrap(); + assert!(close(ahead, 250.0)); + let behind = direct_illuminance(&lum, Vec3::new(-2.0, 0.0, 2.0), Vec3::new(1.0, 0.0, 0.0)).unwrap(); + assert_eq!(behind, 0.0); + // Turned to a heading of 90°, it aims along +y instead. + let turned = Luminaire { rotation_deg: 90.0, ..lum }; + let north = direct_illuminance(&turned, Vec3::new(0.0, 2.0, 2.0), Vec3::new(0.0, -1.0, 0.0)).unwrap(); + assert!(close(north, 250.0)); + // Behind it (−y) is its dark upper hemisphere. + assert_eq!(direct_illuminance(&turned, Vec3::new(0.0, -2.0, 2.0), Vec3::new(0.0, 1.0, 0.0)), Some(0.0)); +} + +#[test] +fn flux_and_maintenance_factor_scale_linearly() { + let mut lum = isotropic_at(2.0); + lum.flux_lm = 200.0; + lum.maintenance_factor = 0.8; + let e = direct_horizontal_illuminance(&lum, Vec3::default()).unwrap(); + assert!((e - 250.0 * 0.2 * 0.8).abs() < 1e-9); +} + +#[test] +fn rotation_selects_the_c_plane() { + // C0 = 200 cd/klm, C180 = 100 cd/klm, flat in γ. + let table = IntensityTable::new(vec![0.0, 180.0], vec![0.0, 180.0], vec![200.0, 200.0, 100.0, 100.0]).unwrap(); + let mut lum = isotropic_at(1.0); + lum.distribution = table; + let east = Vec3::new(1.0, 0.0, 0.0); + let unrotated = direct_horizontal_illuminance(&lum, east).unwrap(); + lum.rotation_deg = 180.0; + let rotated = direct_horizontal_illuminance(&lum, east).unwrap(); + assert!((unrotated - 2.0 * rotated).abs() < 1e-9); +} + +#[test] +fn points_level_with_or_above_the_luminaire_get_nothing() { + let lum = isotropic_at(2.0); + assert_eq!(direct_horizontal_illuminance(&lum, Vec3::new(5.0, 0.0, 2.0)), Some(0.0)); + assert_eq!(direct_horizontal_illuminance(&lum, Vec3::new(0.0, 0.0, 3.0)), Some(0.0)); +} + +#[test] +fn degenerate_inputs_are_none_not_a_panic() { + let lum = isotropic_at(2.0); + assert_eq!(direct_horizontal_illuminance(&lum, Vec3::new(0.0, 0.0, 2.0)), None); + assert_eq!(direct_horizontal_illuminance(&lum, Vec3::new(f64::NAN, 0.0, 0.0)), None); + let mut bad = isotropic_at(2.0); + bad.flux_lm = f64::INFINITY; + assert_eq!(direct_horizontal_illuminance(&bad, Vec3::default()), None); + bad = Luminaire { tilt_deg: f64::NAN, ..isotropic_at(2.0) }; + assert_eq!(direct_horizontal_illuminance(&bad, Vec3::default()), None); +} + +#[test] +fn grids_sum_every_luminaire() { + // Two isotropic 1000 cd sources at 2 m, at x = ±2; the point between them gets 2 × 1000 · (2/√8) / 8. + let a = Luminaire { position: Vec3::new(-2.0, 0.0, 2.0), ..isotropic_at(2.0) }; + let b = Luminaire { position: Vec3::new(2.0, 0.0, 2.0), ..isotropic_at(2.0) }; + let grid = Grid::new(vec![Vec3::new(0.0, 0.0, 0.0), Vec3::new(-2.0, 0.0, 0.0)], Vec3::UP).unwrap(); + let e = illuminance(&[a.clone(), b.clone()], &grid); + let one = 1000.0 * std::f64::consts::FRAC_1_SQRT_2 / 8.0; + assert!(close(e[0].unwrap(), 2.0 * one)); + // Under a: 250 lx from a, plus b from 4 m across: 1000 · (2/√20) / 20. + assert!(close(e[1].unwrap(), 250.0 + 1000.0 * 2.0 / 20f64.sqrt() / 20.0)); + // No luminaires: 0 lx everywhere (dark, not unknown). + assert_eq!(illuminance(&[], &grid), vec![Some(0.0), Some(0.0)]); +} + +#[test] +fn an_undefined_point_makes_the_statistics_unverifiable() { + let lum = isotropic_at(2.0); + let grid = Grid::new(vec![Vec3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 2.0)], Vec3::UP).unwrap(); + let e = illuminance(&[lum], &grid); + assert_eq!(e[1], None); + assert_eq!(Stats::of(&e), None); +} + +#[test] +fn an_area_under_one_luminaire_matches_the_closed_form_point_by_point() { + let lum = lambertian_at(2.5); + let room = [Vec2::new(-3.0, -2.0), Vec2::new(3.0, -2.0), Vec2::new(3.0, 2.0), Vec2::new(-3.0, 2.0)]; + let grid = Grid::area(&room, 0.0, 0.25, 0.5).unwrap(); + let e = illuminance(&[lum], &grid); + for (p, e) in grid.points().iter().zip(&e) { + let d2 = p.x * p.x + p.y * p.y + 2.5 * 2.5; + assert!(close(e.unwrap(), 1000.0 * 2.5 * 2.5 / (d2 * d2))); + } + let stats = Stats::of(&e).unwrap(); + assert!(stats.max_to_min().unwrap() > 1.0); +} + +#[test] +fn tilt_factors_from_the_file_scale_the_light_and_unknown_tilts_are_errors() { + let ies = "IESNA:LM-63-2002\nTILT=INCLUDE\n1\n3\n0 45 90\n1 0.9 0.8\n1 1000 1 2 1 1 2 0 0 0\n1 1 0\n0 180\n0\n1000 1000\n"; + let p = luxpath_photometry::read_ies(ies.as_bytes()).unwrap(); + let at = |tilt: f64| Luminaire::from_photometry(&p, Vec3::new(0.0, 0.0, 2.0), 0.0, tilt, 1000.0, 1.0); + let level = direct_horizontal_illuminance(&at(0.0).unwrap(), Vec3::default()).unwrap(); + let tilted = direct_horizontal_illuminance(&at(45.0).unwrap(), Vec3::default()).unwrap(); + assert!(close(tilted, 0.9 * level)); + assert_eq!(at(100.0).unwrap_err(), CalcError::TiltOutsideTable { tilt_deg: 100.0 }); +} diff --git a/crates/calc/tests/performance.rs b/crates/calc/tests/performance.rs new file mode 100644 index 0000000..c2b06bb --- /dev/null +++ b/crates/calc/tests/performance.rs @@ -0,0 +1,43 @@ +//! M2's performance budget: 100 luminaires × 10,000 points in under 100 ms, with a realistic table +//! (C every 5°, γ every 1°). Timing means something only in release builds, so this test is +//! ignored by default and `cargo xtask ci` runs it with `--release -- --ignored`. + +#![allow(clippy::unwrap_used)] + +use std::time::{Duration, Instant}; + +use luxpath_calc::{Grid, Luminaire, Stats, illuminance}; +use luxpath_geom::{Vec2, Vec3}; +use luxpath_photometry::IntensityTable; + +#[test] +#[ignore = "timing: run in release (cargo xtask ci does)"] +fn a_hundred_luminaires_over_ten_thousand_points_within_budget() { + let c: Vec = (0..72).map(|i| f64::from(i) * 5.0).collect(); + let g: Vec = (0..=180).map(f64::from).collect(); + let values = c.iter().flat_map(|c| g.iter().map(move |g| 100.0 + c / 10.0 + (180.0 - g))).collect(); + let table = IntensityTable::new(c, g, values).unwrap(); + let luminaires: Vec = (0..100) + .map(|i| Luminaire { + position: Vec3::new(f64::from(i % 10) * 5.0, f64::from(i / 10) * 5.0, 3.0), + rotation_deg: f64::from(i) * 7.0, + tilt_deg: f64::from(i % 3) * 10.0, + distribution: table.clone(), + flux_lm: 200.0, + maintenance_factor: 0.8, + }) + .collect(); + let floor = [Vec2::new(0.0, 0.0), Vec2::new(50.0, 0.0), Vec2::new(50.0, 50.0), Vec2::new(0.0, 50.0)]; + let grid = Grid::area(&floor, 0.0, 0.5, 0.0).unwrap(); + assert_eq!(grid.points().len(), 10_000); + + let mut best = Duration::MAX; + for _ in 0..5 { + let start = Instant::now(); + let e = illuminance(&luminaires, &grid); + best = best.min(start.elapsed()); + assert!(Stats::of(&e).is_some()); + } + eprintln!("100 luminaires × 10,000 points: {best:?} (best of 5)"); + assert!(best < Duration::from_millis(100), "{best:?} is over the 100 ms budget"); +} diff --git a/crates/geom/src/lib.rs b/crates/geom/src/lib.rs index de795c0..7717721 100644 --- a/crates/geom/src/lib.rs +++ b/crates/geom/src/lib.rs @@ -51,6 +51,22 @@ impl Vec3 { } } +impl std::ops::Add for Vec3 { + type Output = Vec3; + + fn add(self, other: Vec3) -> Vec3 { + Vec3::new(self.x + other.x, self.y + other.y, self.z + other.z) + } +} + +impl std::ops::Neg for Vec3 { + type Output = Vec3; + + fn neg(self) -> Vec3 { + Vec3::new(-self.x, -self.y, -self.z) + } +} + impl std::ops::Sub for Vec3 { type Output = Vec3; @@ -59,6 +75,76 @@ impl std::ops::Sub for Vec3 { } } +/// A point or vector in plan (the floor plane), in metres. +#[derive(Debug, Clone, Copy, PartialEq, Default)] +pub struct Vec2 { + pub x: f64, + pub y: f64, +} + +impl Vec2 { + pub const fn new(x: f64, y: f64) -> Self { + Self { x, y } + } + + pub fn length(self) -> f64 { + self.x.hypot(self.y) + } + + pub fn is_finite(self) -> bool { + self.x.is_finite() && self.y.is_finite() + } + + /// This point at height `z`. + pub fn at(self, z: f64) -> Vec3 { + Vec3::new(self.x, self.y, z) + } +} + +impl std::ops::Sub for Vec2 { + type Output = Vec2; + + fn sub(self, other: Vec2) -> Vec2 { + Vec2::new(self.x - other.x, self.y - other.y) + } +} + +/// The edges of a closed polygon: each vertex with the next, the last with the first. +fn edges(polygon: &[Vec2]) -> impl Iterator + '_ { + polygon.iter().zip(polygon.iter().cycle().skip(1)).map(|(a, b)| (*a, *b)) +} + +/// Whether `p` lies inside the closed `polygon` (even-odd rule; points exactly on an edge may go +/// either way). Fewer than three vertices enclose nothing. +pub fn point_in_polygon(p: Vec2, polygon: &[Vec2]) -> bool { + if polygon.len() < 3 { + return false; + } + let mut inside = false; + for (a, b) in edges(polygon) { + if (a.y > p.y) != (b.y > p.y) { + let x = a.x + (p.y - a.y) / (b.y - a.y) * (b.x - a.x); + if p.x < x { + inside = !inside; + } + } + } + inside +} + +/// The distance from `p` to the segment from `a` to `b`. +pub fn distance_to_segment(p: Vec2, a: Vec2, b: Vec2) -> f64 { + let ab = b - a; + let len2 = ab.x * ab.x + ab.y * ab.y; + let t = if len2 > 0.0 { (((p.x - a.x) * ab.x + (p.y - a.y) * ab.y) / len2).clamp(0.0, 1.0) } else { 0.0 }; + (p - Vec2::new(a.x + ab.x * t, a.y + ab.y * t)).length() +} + +/// The distance from `p` to the nearest edge of the closed `polygon` (infinite for no edges). +pub fn distance_to_boundary(p: Vec2, polygon: &[Vec2]) -> f64 { + edges(polygon).map(|(a, b)| distance_to_segment(p, a, b)).fold(f64::INFINITY, f64::min) +} + /// Wraps an angle in degrees into `[0, 360)`. Non-finite input returns `0.0`. pub fn wrap_degrees(deg: f64) -> f64 { if !deg.is_finite() { @@ -86,6 +172,35 @@ mod tests { assert!((n.length() - 1.0).abs() < 1e-12); } + #[test] + fn points_are_inside_or_outside_an_l_shaped_room() { + // An L: a 4 × 2 rectangle with a 2 × 2 square on its left half. + let l = [ + Vec2::new(0.0, 0.0), + Vec2::new(4.0, 0.0), + Vec2::new(4.0, 2.0), + Vec2::new(2.0, 2.0), + Vec2::new(2.0, 4.0), + Vec2::new(0.0, 4.0), + ]; + assert!(point_in_polygon(Vec2::new(1.0, 3.0), &l)); + assert!(point_in_polygon(Vec2::new(3.0, 1.0), &l)); + assert!(!point_in_polygon(Vec2::new(3.0, 3.0), &l)); + assert!(!point_in_polygon(Vec2::new(-1.0, 1.0), &l)); + assert!(!point_in_polygon(Vec2::new(1.0, 1.0), &l[..2])); + // The nearest edge of (1, 3) is x = 0 or y = 4 or x = 2: 1 m. + assert_eq!(distance_to_boundary(Vec2::new(1.0, 3.0), &l), 1.0); + assert_eq!(distance_to_boundary(Vec2::new(3.5, 1.0), &l), 0.5); + } + + #[test] + fn segment_distances_clamp_to_the_ends() { + let (a, b) = (Vec2::new(0.0, 0.0), Vec2::new(2.0, 0.0)); + assert_eq!(distance_to_segment(Vec2::new(1.0, 3.0), a, b), 3.0); + assert_eq!(distance_to_segment(Vec2::new(5.0, 4.0), a, b), 5.0); + assert_eq!(distance_to_segment(Vec2::new(0.0, 1.0), a, a), 1.0); + } + #[test] fn degrees_wrap_into_the_half_open_range() { assert_eq!(wrap_degrees(360.0), 0.0); diff --git a/crates/photometry/src/ldt.rs b/crates/photometry/src/ldt.rs index ca21b4b..658e62f 100644 --- a/crates/photometry/src/ldt.rs +++ b/crates/photometry/src/ldt.rs @@ -92,10 +92,10 @@ pub fn read_eulumdat(bytes: &[u8]) -> Result { "luminous area height at C180", "luminous area height at C270", "downward flux fraction", - "light output ratio", ] { fields.number(what)?; } + let (_, declared_lor) = fields.number("light output ratio")?; let (line, conversion) = fields.number("intensity conversion factor")?; if conversion <= 0.0 { return Err(invalid(line, "intensity conversion factor", "is not positive")); @@ -170,7 +170,9 @@ pub fn read_eulumdat(bytes: &[u8]) -> Result { other => invalid(c_first_line, "C angles", other.reason()), })?; let mut photometry = Photometry::new(table, flux_lm, lamps < 0.0, symmetry, Format::Eulumdat) - .map_err(|e| error(flux_line, ReadErrorKind::Table(e)))?; + .map_err(|e| error(flux_line, ReadErrorKind::Table(e)))? + // A ratio of 0 is how some writers leave it out. + .with_declared_lor((declared_lor > 0.0).then_some(declared_lor / 100.0)); photometry.manufacturer = company; photometry.name = name; photometry.catalogue = catalogue; diff --git a/crates/photometry/src/lib.rs b/crates/photometry/src/lib.rs index 4bacf11..0f5da97 100644 --- a/crates/photometry/src/lib.rs +++ b/crates/photometry/src/lib.rs @@ -24,7 +24,7 @@ mod text; pub use ies::{IesEdition, read_ies}; pub use ldt::read_eulumdat; -pub use model::{EmergencyOutput, Format, Photometry, Tilt, read}; +pub use model::{EmergencyOutput, Format, Photometry, PhotometryWarning, Tilt, read}; pub use symmetry::Symmetry; pub use text::{MAX_FILE_BYTES, ReadError, ReadErrorKind}; @@ -190,7 +190,8 @@ impl IntensityTable { let c = wrap_degrees(c); // Find the C planes either side of `c`, wrapping from the last plane back to C0 (= C360). let n_c = self.c_angles.len(); - let lower = self.c_angles.iter().rposition(|a| *a <= c).unwrap_or(0); + // The last plane at or below `c` (binary search: angles increase strictly). + let lower = self.c_angles.partition_point(|a| *a <= c).saturating_sub(1); let (upper, c_lo, c_hi) = match self.c_angles.get(lower + 1) { Some(hi) => (lower + 1, self.c_angles.get(lower).copied().unwrap_or(0.0), *hi), None => (0, self.c_angles.get(lower).copied().unwrap_or(0.0), 360.0), @@ -240,7 +241,7 @@ fn bracket(angles: &[f64], x: f64) -> (usize, usize, f64) { if x >= *last_a { return (last, last, 0.0); } - let hi = angles.iter().position(|a| *a > x).unwrap_or(last); + let hi = angles.partition_point(|a| *a <= x).min(last); let lo = hi.saturating_sub(1); let (a_lo, a_hi) = (angles.get(lo).copied().unwrap_or(0.0), angles.get(hi).copied().unwrap_or(0.0)); let t = if a_hi > a_lo { (x - a_lo) / (a_hi - a_lo) } else { 0.0 }; diff --git a/crates/photometry/src/model.rs b/crates/photometry/src/model.rs index dfc6486..0cc9095 100644 --- a/crates/photometry/src/model.rs +++ b/crates/photometry/src/model.rs @@ -77,6 +77,8 @@ pub struct Photometry { symmetry: Symmetry, format: Format, tilt: Option, + /// The light output ratio the file declares (EULUMDAT line 23), as a fraction. + declared_lor: Option, /// The manufacturer, as the file names it (may be empty). pub manufacturer: String, /// The luminaire's name or description (may be empty). @@ -104,12 +106,18 @@ impl Photometry { symmetry, format, tilt: None, + declared_lor: None, manufacturer: String::new(), name: String::new(), catalogue: String::new(), }) } + pub(crate) fn with_declared_lor(mut self, lor: Option) -> Self { + self.declared_lor = lor; + self + } + pub(crate) fn with_tilt(mut self, tilt: Option) -> Self { self.tilt = tilt; self @@ -149,6 +157,22 @@ impl Photometry { self.table.integrated_flux() / 1000.0 } + /// What looks wrong in the file without making it unreadable. A check using this photometry + /// should show these to the designer. + pub fn warnings(&self) -> Vec { + let mut warnings = Vec::new(); + let lor = self.light_output_ratio(); + if !self.absolute && lor > 1.0 + LOR_SLACK { + warnings.push(PhotometryWarning::MoreLightThanLamps { light_output_ratio: lor }); + } + if let Some(declared) = self.declared_lor + && (declared - lor).abs() > LOR_SLACK.max(0.05 * declared) + { + warnings.push(PhotometryWarning::DeclaredRatioDiffers { declared, integrated: lor }); + } + warnings + } + /// Intensity in candela towards `(c, gamma)` (degrees) at a flux of `flux_lm`. pub fn candela(&self, c: f64, gamma: f64, flux_lm: f64) -> f64 { self.table.intensity(c, gamma) * flux_lm / 1000.0 @@ -171,6 +195,36 @@ impl Photometry { } } +/// The integration's own error stays far below this (`IntensityTable::integrated_flux`). +const LOR_SLACK: f64 = 0.01; + +/// Something implausible in a photometric file. +#[derive(Debug, Clone, Copy, PartialEq)] +pub enum PhotometryWarning { + /// Relative photometry whose table emits more than its lamps (a light output ratio above 1). + MoreLightThanLamps { light_output_ratio: f64 }, + /// The light output ratio the file declares differs from its table's by more than 5 %. + DeclaredRatioDiffers { declared: f64, integrated: f64 }, +} + +impl std::fmt::Display for PhotometryWarning { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + Self::MoreLightThanLamps { light_output_ratio } => write!( + f, + "the table emits {:.1} % of the lamp flux: more light than the lamps give; is the flux or the multiplier wrong?", + light_output_ratio * 100.0 + ), + Self::DeclaredRatioDiffers { declared, integrated } => write!( + f, + "the file declares a light output ratio of {:.1} % but its table integrates to {:.1} %", + declared * 100.0, + integrated * 100.0 + ), + } + } +} + /// Reads an IES LM-63 or EULUMDAT file, telling them apart by content: IES files have a `TILT=` /// line, EULUMDAT files never do. pub fn read(bytes: &[u8]) -> Result { @@ -210,6 +264,29 @@ mod tests { assert_eq!(photometry(2000.0).candela(10.0, 20.0, 500.0), 500.0); } + #[test] + fn implausible_files_are_warned_about() { + // Isotropic 1000 cd/klm emits 4π × 1000 lm per klm: a ratio of 12.6. + let p = photometry(1000.0); + assert!(matches!(p.warnings()[..], [PhotometryWarning::MoreLightThanLamps { .. }])); + // 1000/(4π) cd/klm all round emits exactly the lamps' flux: no warning at a ratio of 1. + let exact = Photometry::new( + IntensityTable::isotropic(1000.0 / (4.0 * std::f64::consts::PI)).unwrap(), + 1000.0, + false, + Symmetry::Rotational, + Format::Eulumdat, + ) + .unwrap(); + assert_eq!(exact.warnings(), vec![]); + // Declared 100 %, integrated 100 %: fine; declared 80 %: warned. + assert_eq!(exact.clone().with_declared_lor(Some(1.0)).warnings(), vec![]); + assert!(matches!( + exact.with_declared_lor(Some(0.8)).warnings()[..], + [PhotometryWarning::DeclaredRatioDiffers { .. }] + )); + } + #[test] fn rated_flux_must_be_positive() { for bad in [0.0, -1.0, f64::NAN] { diff --git a/docs/gaps.md b/docs/gaps.md index c7f1683..7da338b 100644 --- a/docs/gaps.md +++ b/docs/gaps.md @@ -1,22 +1,19 @@ # Gaps: what luxpath cannot do yet -> **Last updated:** 2026-10-21 · **Change:** M1 closed rows 1–2; added rows 12–15 +> **Last updated:** 2026-10-21 · **Change:** M2 closed rows 3, 8, 9 and 13; added row 16 Every known shortfall, ranked by user impact. This is the work list: pick from the top, and remove or update a row in the same PR that closes it. **[alpha blocker]** rows fail the alpha gate (ROADMAP.md). | # | Gap | Why it ranks here | Milestone | Est. h | |---|---|---|---|---| -| 3 | **[alpha blocker]** No grids or multi-luminaire sums | Requirements are about areas and routes, not points | M2 | 8–12 | | 4 | **[alpha blocker]** Requirement values not verified against EN 1838:2024 text | Values in `standards-parity.md` come from secondary sources | M4 | 3 + buying the standard | | 5 | **[alpha blocker]** No scene format or `check` command | Users can't describe a floor | M3 | 15–20 | | 6 | **[alpha blocker]** No EN 1838 checks | The project's purpose | M4 | 25–40 | | 7 | **[alpha blocker]** No independent validation | Trust: nobody should rely on unvalidated numbers | M6 | 10–20 + human | -| 8 | Only downward-pointing luminaires | Wall-mounted luminaires and exit signs are common | M2 | 4–6 | -| 9 | No vertical illuminance | Needed for safety equipment checks | M2 | 2–4 | +| 16 | **[alpha blocker]** Luminaires are points; no obstruction | A wrong pass: direct light from a point ignores the luminaire's size (fine beyond about five times its size) and walls between it and the floor; a corridor around a corner reads lit when it is not | M4 | 6–10 | | 10 | No reports beyond the terminal | Designers hand over documents | M5 | 15–25 | | 11 | No CI on a hosted repository yet | The workflow exists but the repository isn't published | — | 1 | | 12 | **[alpha blocker]** Readers tried on only two real files | Synthetic files test the formats as we read them; manufacturers' files (kept in `plan/`, never committed) test them as written. Wanted: a few dozen IES and LDT files, including emergency luminaires | M6 | 2–4 + files from the owner | -| 13 | No plausibility warnings on photometry | A relative file emitting more than its lamps (the LM-63 example file integrates to 211 %), or a declared LDT light output ratio far from the table's, should be reported, not silently used | M2 | 2–3 | | 14 | Format details from our reading, not the specification texts | The LM-63 factor that changed meaning between editions, and the C0 orientation of each format, are unverified like the EN 1838 values (gap 4); LM-63-2019 is sold by the IES | M6 | 2 + buying LM-63-2019 | | 15 | Refused file variants | IES type A and B photometry, `TILT=`, and EULUMDAT files with several lamp sets are refused with a located error | M7 | 4–8 | \ No newline at end of file diff --git a/docs/roadmap.md b/docs/roadmap.md index f2dd6aa..4dcd7ad 100644 --- a/docs/roadmap.md +++ b/docs/roadmap.md @@ -1,6 +1,6 @@ # Milestones -> **Last updated:** 2026-10-21 · **Change:** M1 done +> **Last updated:** 2026-10-21 · **Change:** M1 and M2 done Estimates are agent-assisted hours (a human directing an agent, reviewing every change) and are rough until M2 calibrates them. Each milestone ends with its exit criteria met, `cargo xtask ci` green, and the parity entries @@ -10,7 +10,7 @@ updated. |---|---|---|---|---|---| | M0 | Foundations | Workspace, lints, CI, `xtask` gates (ci, layers, parity), docs, intensity table, direct point illuminance, demo CLI | All gates green; one hand-checked calculation end to end | 5–8 | **done** | | M1 | Photometry | IES LM-63 (1995/2002/2019) and EULUMDAT readers; symmetry expansion; tilt data; emergency-mode flux; validation errors with line numbers; property tests and fuzzing | Synthetic fixtures for every symmetry type read identically from IES and LDT; fuzzing finds no panic in 10 min | 20–35 | **done** | -| M2 | Engine | Calculation grids (areas, route bands), sums over luminaires, tilted/wall-mounted luminaires, vertical illuminance, min/avg/max/uniformity statistics | Analytic test cases (isotropic, Lambertian) match to 1e-6; performance: 100 luminaires × 10k points < 100 ms | 20–30 | planned | +| M2 | Engine | Calculation grids (areas, route bands), sums over luminaires, tilted/wall-mounted luminaires, vertical illuminance, min/avg/max/uniformity statistics | Analytic test cases (isotropic, Lambertian) match to 1e-6; performance: 100 luminaires × 10k points < 100 ms | 20–30 | **done** | | M3 | Scene and CLI | TOML scene format (versioned), 2D geometry (polygons, routes with width), `luxpath check scene.toml` with exit codes | Three example scenes in `examples/`; malformed scenes give located errors | 20–30 | planned | | M4 | EN 1838:2024 | Every requirement in `docs/standards-parity.md`: escape routes over the full width, uniformity, open areas, high-risk areas, points of emphasis, glare limits, declared duration/response | Each P0 requirement tested on both sides of its threshold; values verified against the standard's text | 25–40 | planned | | M5 | Reports | Text, JSON (stable schema), SVG plan with false colours and failing points; PDF later | A report a designer can read without knowing the tool | 15–25 | planned | diff --git a/docs/status.md b/docs/status.md index 0f038e9..e963729 100644 --- a/docs/status.md +++ b/docs/status.md @@ -2,9 +2,18 @@ > **Last updated:** 2026-10-21 -**Milestone:** M1 done (photometry readers). **Direction under study:** EN 1838 checks on IFC models +**Milestone:** M2 done (calculation engine); M1 done (photometry readers). **Direction under study:** EN 1838 checks on IFC models (docs/research.md §6–7). The IFC inspector (`luxpath inspect-ifc`, an early piece of M7) is done. +## M2 (2026-10-21) +`luxpath-calc` places luminaires with a heading and a tilt (90° = wall-mounted; IES tilt factors applied, a tilt +outside the file's table refused), computes direct illuminance on any surface (horizontal, vertical, inclined), +builds grids over areas (less a border) and along routes (over their full width), sums every luminaire and gives +min, mean, max, max/min and min/mean; one undefined point makes the statistics unverifiable. Isotropic and +Lambertian sources match their closed forms within 1e-6; 100 luminaires × 10,000 points take about 55 ms. +`luxpath photometry` warns about implausible files. The engine still treats luminaires as points and ignores walls +(gap 16, an alpha blocker before M4). + ## M1 (2026-10-21) `luxpath-photometry` reads IES LM-63 (1986–2019, type C, relative and absolute, tilt tables) and EULUMDAT (every symmetry indicator, Latin-1, decimal commas) into one full C-γ table; errors name their line. Synthetic files of @@ -20,11 +29,11 @@ FireExit flags. Luminaires collapse into few types (410 → 5), so a photometry - **A. Evidence:** inspect one real project IFC from a practising office (privately, in `plan/`). - **B. IFC scene:** read luminaire positions and space outlines (geometry) and define the mapping file (type name or product code → photometric file, emergency yes/no). -- **C. M2 engine:** grids, sums over luminaires, tilted luminaires, vertical illuminance (with M1's photometry). +- **C. M3 scene:** the TOML scene format and `luxpath check` (then obstruction, gap 16, before M4). - **D. Real files:** photometric files of luminaires the owner's practice uses, privately in `plan/` (gap 12). ## Blockers -- None for M2. **Before M4:** the EN 1838:2024 text (Italian edition UNI EN 1838:2025) is needed to verify values (gap 4). +- None for M3. **Before M4:** the EN 1838:2024 text (Italian edition UNI EN 1838:2025) is needed to verify values (gap 4). ## Open questions (owner) - Practitioner questions in `docs/research.md` §5. diff --git a/parity/features.toml b/parity/features.toml index ba6afa3..12328bb 100644 --- a/parity/features.toml +++ b/parity/features.toml @@ -67,7 +67,7 @@ tier = "P0" milestone = "M1" status = "partial" evidence = ["crates/photometry/src/ies.rs::a_quadrant_file_reads_into_a_full_table_in_cd_per_klm", "crates/photometry/src/ies.rs::absolute_photometry_takes_its_flux_from_the_table", "crates/photometry/src/ies.rs::older_editions_apply_the_ballast_lamp_factor_and_newer_ones_do_not", "crates/photometry/src/ies.rs::tilt_tables_are_read", "crates/photometry/src/ies.rs::errors_name_their_line", "crates/photometry/tests/fuzz.rs::fuzz_briefly"] -notes = "Editions 1986–2019, type C, relative and absolute photometry, TILT=INCLUDE, located errors, fuzzed. Missing: type A and B photometry (refused, not guessed), TILT=; the edition-specific factors are from our reading of the format, not yet checked against the LM-63 texts; no real manufacturer file tried yet." +notes = "Editions 1986–2019, type C, relative and absolute photometry, TILT=INCLUDE, located errors, fuzzed, warnings for files emitting more than their lamps. Missing: type A and B photometry (refused, not guessed), TILT=; the edition-specific factors are from our reading of the format, not yet checked against the LM-63 texts; no real manufacturer file tried yet." [[feature]] id = "photometry.ldt-read" @@ -77,7 +77,7 @@ tier = "P0" milestone = "M1" status = "partial" evidence = ["crates/photometry/src/ldt.rs::a_c0_c180_file_reads_into_a_full_table", "crates/photometry/src/ldt.rs::every_symmetry_indicator_selects_its_planes", "crates/photometry/src/ldt.rs::decimal_commas_and_latin1_are_accepted", "crates/photometry/src/ldt.rs::errors_name_their_line", "crates/photometry/tests/fuzz.rs::fuzz_briefly"] -notes = "All five symmetry indicators, Latin-1, decimal commas, located errors, fuzzed. Missing: files with more than one lamp set (refused), a check of the declared light output ratio against the table; no real manufacturer file tried yet." +notes = "All five symmetry indicators, Latin-1, decimal commas, located errors, fuzzed; warnings when the declared light output ratio and the table disagree. Missing: files with more than one lamp set (refused); no real manufacturer file tried yet beyond one public sample." [[feature]] id = "photometry.emergency-flux" @@ -107,9 +107,9 @@ title = "Direct horizontal illuminance at a point (inverse-square cosine law)" area = "calc" tier = "P0" milestone = "M0" -status = "partial" -evidence = ["crates/calc/src/lib.rs::off_axis_points_follow_the_cosine_law", "crates/calc/src/lib.rs::rotation_selects_the_c_plane"] -notes = "Downward-pointing luminaires rotated about z. Missing: tilted luminaires, vertical and inclined planes, obstruction." +status = "shipped" +evidence = ["crates/calc/src/tests.rs::off_axis_points_follow_the_cosine_law", "crates/calc/src/tests.rs::a_lambertian_source_matches_its_closed_form", "crates/calc/src/tests.rs::rotation_selects_the_c_plane"] +notes = "On any surface normal, isotropic and Lambertian sources within 1e-6 of their closed forms. Missing: obstruction (walls between luminaire and point), area sources (luminaires are points)." [[feature]] id = "calc.grids" @@ -117,7 +117,9 @@ title = "Calculation grids over areas and along routes, summed over all luminair area = "calc" tier = "P0" milestone = "M2" -status = "planned" +status = "shipped" +evidence = ["crates/calc/src/grid.rs::an_area_grid_takes_cell_centres_inside_less_the_border", "crates/calc/src/grid.rs::a_route_grid_covers_the_full_width_and_each_corner_once", "crates/calc/src/tests.rs::grids_sum_every_luminaire", "crates/calc/src/tests.rs::an_area_under_one_luminaire_matches_the_closed_form_point_by_point", "crates/calc/src/stats.rs::statistics_by_hand", "crates/calc/tests/performance.rs::a_hundred_luminaires_over_ten_thousand_points_within_budget"] +notes = "Area grids (polygon less a border), route grids over the full width, explicit points; min, mean, max, max/min and min/mean; any undefined point makes the statistics unverifiable. 100 luminaires × 10,000 points in about 55 ms. The spacings and borders EN 1838 asks for are M4's (from the standard's text)." [[feature]] id = "calc.tilt" @@ -125,7 +127,9 @@ title = "Tilted and wall-mounted luminaires" area = "calc" tier = "P1" milestone = "M2" -status = "planned" +status = "shipped" +evidence = ["crates/calc/src/tests.rs::a_tilted_luminaire_aims_its_light_towards_c0", "crates/calc/src/tests.rs::tilt_factors_from_the_file_scale_the_light_and_unknown_tilts_are_errors"] +notes = "Heading and tilt (about the luminaire's C90-C270 axis, towards C0); IES tilt factors applied, a tilt outside the file's table refused. Missing: a third rotation (roll about C0) for luminaires mounted on sloped ceilings." [[feature]] id = "calc.vertical-illuminance" @@ -133,7 +137,9 @@ title = "Vertical illuminance (e.g. at fire-fighting equipment)" area = "calc" tier = "P0" milestone = "M2" -status = "planned" +status = "shipped" +evidence = ["crates/calc/src/tests.rs::vertical_illuminance_follows_the_cosine_law_on_the_wall"] +notes = "Any surface normal, so vertical and inclined planes alike. Which points and normals EN 1838 asks for at equipment is M4's." [[feature]] id = "calc.first-reflection" diff --git a/xtask/src/main.rs b/xtask/src/main.rs index 8d04418..96b944c 100644 --- a/xtask/src/main.rs +++ b/xtask/src/main.rs @@ -1,6 +1,7 @@ //! `cargo xtask `: repository automation. //! -//! - `ci`: every gate CI runs (fmt, clippy -D warnings, tests, wasm check, layers, parity). +//! - `ci`: every gate CI runs (fmt, clippy -D warnings, tests, the release-mode performance budget, +//! wasm check, layers, parity). //! - `layers`: checks that crates only depend on crates in lower layers. //! - `parity`: validates `parity/features.toml` and prints progress by milestone. //! - `fuzz [seconds]`: fuzzes the photometric file readers (default 600 s, M1's exit criterion). @@ -43,6 +44,8 @@ fn ci() -> Result<(), String> { cargo(&["fmt", "--all", "--check"])?; cargo(&["clippy", "--workspace", "--all-targets", "--", "-D", "warnings"])?; cargo(&["test", "--workspace"])?; + // M2's budget: timing is meaningful only in release. + cargo(&["test", "--release", "-p", "luxpath-calc", "--test", "performance", "--", "--ignored"])?; if wasm_target_installed() { let mut args = vec!["check", "--target", "wasm32-unknown-unknown"]; for krate in WASM_CRATES {