earthsci_ast/

units.rs

//! Unit parsing and dimensional analysis.
//!
//! This module implements a hand-rolled dimensional algebra for ESM
//! expressions. Two layers are exposed:
//!
//! 1. [`parse_unit`] parses a unit string (e.g. `"kg*m/s^2"`, `"J/(mol*K)"`)
//!    into a [`Unit`] carrying a `HashMap<Dimension, i32>` and a scale
//!    factor. Used for variable/parameter unit metadata.
//!
//! 2. [`Unit::propagate`] walks an [`Expr`] AST and returns the resulting
//!    [`Unit`] for the whole expression. This is the expression-level
//!    dimensional propagation that lets a caller verify, for example, that
//!    `D(h)/dt` has the same dimensions as a `v` of units `m/s`.
//!
//! The propagation rules match the Python implementation
//! (`earthsci_ast/units.py::_get_expression_dimension`) and the Julia
//! implementation (`EarthSciAST.jl/src/units.jl::
//! get_expression_dimensions`).

use crate::types::{Equation, Expr, ExpressionNode};
use std::collections::{BTreeMap, HashMap};
use thiserror::Error;

/// A dimension exponent.
///
/// Exponents are RATIONAL, not integer: `1/s^0.5` — the intensity of an SDE
/// noise term — is a legitimate unit that appears in the shared corpus
/// (`tests/fixtures/sde/*.esm`), and `sqrt` halves whatever it is given. An
/// integer exponent could represent neither `sqrt(m^3)` nor `s^0.5`, and since
/// a dimensional mismatch is a hard error, legitimate files would be
/// rejected.
///
/// Always stored in lowest terms with a positive denominator, so `PartialEq`
/// is exact dimensional equality (½ and 2/4 compare equal) and `HashMap` lookups
/// are sound.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
pub struct Rational {
    num: i32,
    den: i32,
}

impl Rational {
    /// A whole number exponent.
    #[must_use]
    pub fn int(n: i32) -> Self {
        Rational { num: n, den: 1 }
    }

    /// `num/den`, reduced. A zero denominator is treated as `0/1`.
    #[must_use]
    pub fn new(num: i32, den: i32) -> Self {
        if den == 0 {
            return Rational { num: 0, den: 1 };
        }
        let sign = if den < 0 { -1 } else { 1 };
        let (num, den) = (num * sign, den * sign);
        let g = gcd(num.abs(), den);
        if g == 0 {
            return Rational { num: 0, den: 1 };
        }
        Rational {
            num: num / g,
            den: den / g,
        }
    }

    /// Convert a decimal literal (`0.5`, `-1.5`, `2`) to an exact rational.
    ///
    /// Goes through the DECIMAL representation rather than the binary float, so
    /// `0.5` is exactly ½ and `0.1` is exactly 1/10 — a continued-fraction
    /// expansion of the `f64` would introduce a spurious huge denominator.
    /// Returns `None` for a non-finite value or one needing more precision than
    /// an `i32` numerator can hold.
    #[must_use]
    pub fn from_f64(x: f64) -> Option<Self> {
        if !x.is_finite() {
            return None;
        }
        if x.fract() == 0.0 && x.abs() < i32::MAX as f64 {
            return Some(Rational::int(x as i32));
        }
        // Up to 6 decimal places is far more than any real unit exponent needs.
        for places in 1..=6u32 {
            let den = 10i32.checked_pow(places)?;
            let scaled = x * f64::from(den);
            if scaled.fract().abs() < 1e-9 && scaled.abs() < i32::MAX as f64 {
                return Some(Rational::new(scaled.round() as i32, den));
            }
        }
        None
    }

    fn add(self, other: Self) -> Self {
        Rational::new(
            self.num * other.den + other.num * self.den,
            self.den * other.den,
        )
    }

    fn sub(self, other: Self) -> Self {
        Rational::new(
            self.num * other.den - other.num * self.den,
            self.den * other.den,
        )
    }

    fn mul(self, other: Self) -> Self {
        Rational::new(self.num * other.num, self.den * other.den)
    }

    fn is_zero(self) -> bool {
        self.num == 0
    }

    /// `num/den` from 64-bit parts, reduced; `None` when the reduced value does
    /// not fit the `i32` representation (or `den` is zero).
    fn checked_from_i64(num: i64, den: i64) -> Option<Self> {
        if den == 0 {
            return None;
        }
        let sign = if den < 0 { -1 } else { 1 };
        let (num, den) = (num.checked_mul(sign)?, den.checked_mul(sign)?);
        let g = gcd_i64(num.abs(), den).max(1);
        Some(Rational {
            num: i32::try_from(num / g).ok()?,
            den: i32::try_from(den / g).ok()?,
        })
    }

    /// Exact sum, or `None` on overflow of the `i32` representation.
    fn checked_add(self, other: Self) -> Option<Self> {
        let (a, b, c, d) = self.parts_with(other);
        Rational::checked_from_i64(a * d + c * b, b * d)
    }

    /// Exact difference, or `None` on overflow.
    fn checked_sub(self, other: Self) -> Option<Self> {
        let (a, b, c, d) = self.parts_with(other);
        Rational::checked_from_i64(a * d - c * b, b * d)
    }

    /// Exact product, or `None` on overflow.
    fn checked_mul(self, other: Self) -> Option<Self> {
        let (a, b, c, d) = self.parts_with(other);
        Rational::checked_from_i64(a * c, b * d)
    }

    /// Exact quotient, or `None` on overflow or a zero divisor.
    fn checked_div(self, other: Self) -> Option<Self> {
        let (a, b, c, d) = self.parts_with(other);
        Rational::checked_from_i64(a * d, b * c)
    }

    /// Both operands' parts widened to `i64`, where a product of two cannot
    /// overflow.
    fn parts_with(self, other: Self) -> (i64, i64, i64, i64) {
        (
            i64::from(self.num),
            i64::from(self.den),
            i64::from(other.num),
            i64::from(other.den),
        )
    }

    /// As an `f64`, for scaling a unit's magnitude by this exponent.
    fn as_f64(self) -> f64 {
        f64::from(self.num) / f64::from(self.den)
    }
}

impl std::fmt::Display for Rational {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        if self.den == 1 {
            write!(f, "{}", self.num)
        } else {
            write!(f, "{}/{}", self.num, self.den)
        }
    }
}

fn gcd(a: i32, b: i32) -> i32 {
    if b == 0 { a } else { gcd(b, a % b) }
}

fn gcd_i64(a: i64, b: i64) -> i64 {
    if b == 0 { a } else { gcd_i64(b, a % b) }
}

/// One factor of an [`ExactScale`]: a prime number, or π.
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
pub enum ScaleFactor {
    /// A prime number.
    Prime(u64),
    /// The circle constant, which only the degree of arc carries (`deg` is π/180 rad).
    Pi,
}

/// The EXACT scale of a unit relative to SI (esm-spec §4.8.1 "Scales are EXACT").
///
/// Stored as a product of prime powers and a power of π, each with a rational
/// exponent: `mi` is 2⁴·3²·5⁻³·11·127. Multiplying units adds exponents, dividing
/// subtracts them and a rational power multiplies them, so a composite unit's scale
/// never rounds and `sqrt(km)` stays exact. Equality is exact, which is what lets a
/// checker decide `m + km` or `m/s = mi/h` without a floating-point tolerance. The
/// floating-point [`Unit::scale`] is kept for NUMERIC conversion only.
#[derive(Debug, Clone, PartialEq, Eq, Default)]
pub struct ExactScale {
    factors: BTreeMap<ScaleFactor, Rational>,
}

impl ExactScale {
    /// The scale 1.
    #[must_use]
    pub fn one() -> Self {
        Self::default()
    }

    /// A positive whole number.
    ///
    /// # Panics
    ///
    /// Panics on zero, which is not a scale.
    #[must_use]
    pub fn integer(n: u64) -> Self {
        assert!(n > 0, "a unit scale must be positive");
        let mut out = Self::one();
        let mut rest = n;
        let mut p = 2u64;
        while p * p <= rest {
            while rest.is_multiple_of(p) {
                out.bump(ScaleFactor::Prime(p), Rational::int(1));
                rest /= p;
            }
            p += if p == 2 { 1 } else { 2 };
        }
        if rest > 1 {
            out.bump(ScaleFactor::Prime(rest), Rational::int(1));
        }
        out
    }

    /// `num/den` for positive whole numbers.
    #[must_use]
    pub fn ratio(num: u64, den: u64) -> Self {
        Self::integer(num).divide(&Self::integer(den))
    }

    /// `10^k`.
    #[must_use]
    pub fn power_of_ten(k: i32) -> Self {
        let mut out = Self::one();
        out.bump(ScaleFactor::Prime(2), Rational::int(k));
        out.bump(ScaleFactor::Prime(5), Rational::int(k));
        out
    }

    /// π.
    #[must_use]
    pub fn pi() -> Self {
        let mut out = Self::one();
        out.bump(ScaleFactor::Pi, Rational::int(1));
        out
    }

    /// The exact value of a positive decimal literal as the registry writes it
    /// (`"0.3048"`, `"133.322387415"`, `"2.6867e20"`).
    ///
    /// # Panics
    ///
    /// Panics on anything that is not such a literal: the registry is a fixed
    /// table, so a typo in it is a bug to fail on loudly.
    #[must_use]
    pub fn decimal(s: &str) -> Self {
        let (mantissa, exp) = match s.split_once(['e', 'E']) {
            Some((m, e)) => (m, e.parse::<i32>().expect("registry decimal exponent")),
            None => (s, 0),
        };
        let (int_part, frac_part) = mantissa.split_once('.').unwrap_or((mantissa, ""));
        let digits = format!("{int_part}{frac_part}");
        let n: u64 = digits
            .trim_start_matches('0')
            .parse()
            .expect("registry decimal mantissa");
        let places = i32::try_from(frac_part.len()).expect("registry decimal places");
        Self::integer(n).multiply(&Self::power_of_ten(exp - places))
    }

    fn bump(&mut self, factor: ScaleFactor, by: Rational) {
        let e = self.factors.entry(factor).or_insert(Rational::int(0));
        *e = e.add(by);
        if e.is_zero() {
            self.factors.remove(&factor);
        }
    }

    /// The scale of a product of units.
    #[must_use]
    pub fn multiply(&self, other: &Self) -> Self {
        let mut out = self.clone();
        for (f, e) in &other.factors {
            out.bump(*f, *e);
        }
        out
    }

    /// The scale of a quotient of units.
    #[must_use]
    pub fn divide(&self, other: &Self) -> Self {
        let mut out = self.clone();
        for (f, e) in &other.factors {
            out.bump(*f, e.mul(Rational::int(-1)));
        }
        out
    }

    /// The scale of a unit raised to a rational power.
    #[must_use]
    pub fn power(&self, exponent: Rational) -> Self {
        let mut out = Self::one();
        for (f, e) in &self.factors {
            out.bump(*f, e.mul(exponent));
        }
        out
    }

    /// Whether this is exactly 1.
    #[must_use]
    pub fn is_one(&self) -> bool {
        self.factors.is_empty()
    }

    /// The nearest `f64`, for diagnostics and tests. Two scales are never
    /// compared through this.
    #[must_use]
    pub fn to_f64(&self) -> f64 {
        self.factors.iter().fold(1.0, |acc, (f, e)| {
            let base = match f {
                ScaleFactor::Prime(p) => *p as f64,
                ScaleFactor::Pi => std::f64::consts::PI,
            };
            acc * base.powf(e.as_f64())
        })
    }

    /// `p/q`, `p/q*pi` or `p/q*pi^k` in lowest terms (`/q` omitted when it is 1),
    /// the spelling `tests/conformance/unit_registry` pins. `None` when an exponent
    /// is not whole (`sqrt(km)`) or the numbers do not fit in 128 bits.
    #[must_use]
    pub fn ratio_string(&self) -> Option<String> {
        let (mut num, mut den, mut pi) = (1u128, 1u128, 0i32);
        for (f, e) in &self.factors {
            if e.den != 1 {
                return None;
            }
            match f {
                ScaleFactor::Pi => pi = e.num,
                ScaleFactor::Prime(p) => {
                    let pow = u128::from(*p).checked_pow(e.num.unsigned_abs())?;
                    if e.num > 0 {
                        num = num.checked_mul(pow)?;
                    } else {
                        den = den.checked_mul(pow)?;
                    }
                }
            }
        }
        let mut s = num.to_string();
        if den != 1 {
            s = format!("{s}/{den}");
        }
        Some(match pi {
            0 => s,
            1 => format!("{s}*pi"),
            k => format!("{s}*pi^{k}"),
        })
    }
}

impl std::fmt::Display for ExactScale {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self.ratio_string() {
            Some(s) => write!(f, "{s}"),
            None => write!(f, "{:e}", self.to_f64()),
        }
    }
}

/// Represents a physical unit with dimensions
#[derive(Debug, Clone, PartialEq)]
pub struct Unit {
    /// Base dimensions with their (rational) powers
    dimensions: HashMap<Dimension, Rational>,
    /// Scale factor for NUMERIC unit conversions. It never decides a
    /// dimensional verdict; see [`ExactScale`].
    scale: f64,
    /// The exact scale every scale agreement is decided on (esm-spec §4.8.1).
    exact: ExactScale,
}

/// Base physical dimensions — the eight canonical axes shared across the
/// bindings (`m kg s mol K A cd rad`).
///
/// `Angle` is a full axis rather than a synonym for dimensionless, matching the
/// Go reference. The trigonometric functions therefore accept an angle OR a
/// dimensionless argument (see `propagate_transcendental_dim`); `exp`/`log` still
/// demand strict dimensionlessness.
#[derive(Debug, Clone, PartialEq, Eq, Hash)]
pub enum Dimension {
    Mass,        // M   — kg
    Length,      // L   — m
    Time,        // T   — s
    Current,     // I   — A  (electric current)
    Temperature, // Θ   — K  (thermodynamic temperature)
    Amount,      // N   — mol (amount of substance)
    Luminosity,  // J   — cd (luminous intensity)
    Angle,       // rad (plane angle)
}

/// Error types for unit operations
#[derive(Debug, Clone, PartialEq, Error)]
pub enum UnitError {
    #[error("Parse error: {0}")]
    ParseError(String),
    #[error("Dimension mismatch: {0}")]
    DimensionMismatch(String),
    #[error("Unknown unit: {0}")]
    UnknownUnit(String),
}

/// Severity of a dimensional finding — the cross-binding policy that decides
/// whether a finding blocks validation.
///
/// The split is drawn at *provability*, and it is decided AT THE POINT the
/// finding is raised (mirroring the TypeScript reference, whose `UnitWarning`
/// carries the same classification):
///
/// * [`UnitSeverity::Error`] — a PROVABLE dimensional inconsistency: every
///   operand's dimension was determined, and they are incompatible. That is a
///   defect in the FILE, so it is promoted by `structural.rs` to a
///   `unit_inconsistency` structural error and makes `is_valid` false. The
///   shared corpus requires this: `tests/invalid/expected_errors.json` pins the
///   `units_*.esm` fixtures as `is_valid: false` with a structural error, so a
///   binding that keeps them as warnings ACCEPTS files the corpus pins invalid.
///
/// * [`UnitSeverity::Analysis`] — the checker could not DETERMINE a dimension
///   (unknown variable, unparseable unit string, symbolic/non-literal exponent,
///   wrong arity, a broadcast missing its `fn`). This reports what the checker
///   could not conclude, NOT a defect in the file, so it stays a non-blocking
///   warning and the affected subexpression propagates as [`Dim::Unknown`].
///
/// An operator with NO dimensional rule (a rewrite-target sugar op such as
/// `grad`/`div`/`laplacian`, or any unregistered user op) is the one
/// undeterminable case that emits NO finding at all: per esm-spec §4.8.4 the
/// checker reports UNKNOWN and skips silently, matching the other four bindings
/// (Julia/TS/Python/Go). See [`propagate_operator_dim`].
///
/// An unknown *variable* stays `Analysis` here only because it is already a
/// hard `undefined_variable` structural error — the file is still rejected; we
/// just do not double-report it.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum UnitSeverity {
    /// Provable dimensional inconsistency — blocks validation.
    Error,
    /// Undeterminable dimension — non-blocking warning.
    Analysis,
}

/// The cross-binding `code` values a dimensional finding can carry, as surfaced
/// on [`crate::validate::UnitWarning::code`]. Mirrors Go's `UnitFinding*`
/// constants (`pkg/earthsci-ast-go/pkg/esm/units.go`) and TypeScript's
/// `UnitWarning['code']` union (`pkg/earthsci-ast-ts/src/units.ts`).
///
/// * [`UNIT_FINDING_DIMENSIONAL_MISMATCH`] — a PROVABLE inconsistency: adding
///   metres to kilograms, an equation whose sides cannot agree. A defect in the
///   FILE; this binding promotes it to a `unit_inconsistency` structural error.
/// * [`UNIT_FINDING_UNPARSEABLE`] — a declared unit string that does not denote
///   a real unit. Also a defect in the FILE; promoted to `unit_parse_error` at
///   the declaration's own pointer.
/// * [`UNIT_FINDING_ANALYSIS`] — the checker cannot DETERMINE a dimension. A
///   limit of the ANALYSIS, not a defect, so it stays a non-blocking warning.
///
/// These strings are a cross-binding contract; they must not change.
pub const UNIT_FINDING_DIMENSIONAL_MISMATCH: &str = crate::diagnostic::codes::DIMENSIONAL_MISMATCH;
/// See [`UNIT_FINDING_DIMENSIONAL_MISMATCH`].
pub const UNIT_FINDING_UNPARSEABLE: &str = crate::diagnostic::codes::UNPARSEABLE_UNIT;
/// See [`UNIT_FINDING_DIMENSIONAL_MISMATCH`].
pub const UNIT_FINDING_ANALYSIS: &str = crate::diagnostic::codes::ANALYSIS;

/// A single dimensional finding produced while propagating an expression.
#[derive(Debug, Clone, PartialEq)]
pub struct UnitFinding {
    /// Whether this finding blocks validation. See [`UnitSeverity`].
    pub severity: UnitSeverity,
    /// Human-readable description of the finding.
    pub message: String,
}

impl UnitFinding {
    fn error(message: impl Into<String>) -> Self {
        UnitFinding {
            severity: UnitSeverity::Error,
            message: message.into(),
        }
    }

    fn analysis(message: impl Into<String>) -> Self {
        UnitFinding {
            severity: UnitSeverity::Analysis,
            message: message.into(),
        }
    }

    /// True when this finding is a provable inconsistency.
    pub fn is_error(&self) -> bool {
        self.severity == UnitSeverity::Error
    }

    /// The cross-binding finding code implied by this finding's severity, taken
    /// at the raise site rather than recovered from the message prose.
    ///
    /// An unparseable declared unit never reaches a [`UnitFinding`] in this
    /// binding — it is caught at the declaration and reported directly as a
    /// `unit_parse_error` structural error — so only the two severities map
    /// here. See [`UNIT_FINDING_DIMENSIONAL_MISMATCH`].
    pub fn code(&self) -> &'static str {
        match self.severity {
            UnitSeverity::Error => UNIT_FINDING_DIMENSIONAL_MISMATCH,
            UnitSeverity::Analysis => UNIT_FINDING_ANALYSIS,
        }
    }
}

/// The dimension of a subexpression: either determined, or not.
///
/// `Unknown` is what makes multi-finding collection sound. Propagation NEVER
/// aborts an expression at the first finding — it records the finding, yields
/// `Unknown` for that subtree, and keeps walking its siblings. An
/// undeterminable operand therefore cannot SUPPRESS a provable mismatch
/// elsewhere in the same expression (the old `Result`-based propagation bailed
/// at the first `Err`, so a single unknown variable hid every real error after
/// it), and it equally cannot MANUFACTURE one: `Unknown` compares against
/// nothing, so no finding is raised for it.
#[derive(Debug, Clone, PartialEq)]
enum Dim {
    /// Dimension determined.
    Known(Unit),
    /// Dimension could not be determined. Most undeterminable causes record an
    /// `Analysis` finding; an op with no dimensional rule records none (it is
    /// reported as unknown and skipped, per esm-spec §4.8.4).
    Unknown,
}

impl Dim {
    fn known(&self) -> Option<&Unit> {
        match self {
            Dim::Known(u) => Some(u),
            Dim::Unknown => None,
        }
    }
}

impl Unit {
    /// The multiplicative factor taking this unit to its SI base (`atm` → 101325).
    ///
    /// NOTE: the representation is purely multiplicative — there is no affine
    /// OFFSET — so `degC` and `K` are indistinguishable here. A check that must
    /// separate them (the `default_units` identity check) cannot do so by
    /// dimension and scale alone.
    pub fn scale(&self) -> f64 {
        self.scale
    }

    /// The exact scale this unit's agreements are decided on (esm-spec §4.8.1).
    pub fn exact_scale(&self) -> &ExactScale {
        &self.exact
    }

    /// Whether two units agree in dimension AND exact scale — what `+`, `-`,
    /// comparisons and the two sides of an equation require (esm-spec §4.8.3).
    pub fn same_unit(&self, other: &Unit) -> bool {
        self.is_compatible(other) && self.exact == other.exact
    }

    /// This unit with its exact scale set. Every registry entry whose scale is
    /// not 1 is built through this.
    fn with_exact(mut self, exact: ExactScale) -> Self {
        self.exact = exact;
        self
    }

    /// Whether two units measure the same physical quantity, ignoring scale
    /// (`atm` and `Pa` are the same dimension at different scales).
    pub fn same_dimensions(&self, other: &Unit) -> bool {
        self.dimensions == other.dimensions
    }

    /// Create a dimensionless unit
    pub fn dimensionless() -> Self {
        Unit {
            dimensions: HashMap::new(),
            scale: 1.0,
            exact: ExactScale::one(),
        }
    }

    /// Create a unit with a single dimension raised to an integer power.
    ///
    /// `scale` is the NUMERIC conversion factor only; the unit's
    /// [`ExactScale`] is 1, so a scaled unit built here agrees exactly with its
    /// unscaled base (esm-spec §4.8.3). Use [`parse_unit`] for a scaled unit.
    pub fn base(dimension: Dimension, power: i32, scale: f64) -> Self {
        Unit::base_rational(dimension, Rational::int(power), scale)
    }

    /// Create a unit with a single dimension raised to a RATIONAL power (e.g.
    /// `s^-1/2`, the dimension of an SDE noise intensity). As with [`Unit::base`],
    /// the exact scale is 1 whatever `scale` is.
    pub fn base_rational(dimension: Dimension, power: Rational, scale: f64) -> Self {
        let mut dimensions = HashMap::new();
        if !power.is_zero() {
            dimensions.insert(dimension, power);
        }
        Unit {
            dimensions,
            scale,
            exact: ExactScale::one(),
        }
    }

    /// Check if two units have compatible dimensions
    pub fn is_compatible(&self, other: &Unit) -> bool {
        self.dimensions == other.dimensions
    }

    /// Check if this unit is dimensionless
    pub fn is_dimensionless(&self) -> bool {
        self.dimensions.is_empty()
    }

    /// True when this unit is dimensionless **at scale 1** — a PURE NUMBER.
    ///
    /// This is what esm-spec §4.8.3's "the argument MUST be dimensionless"
    /// means for a strict transcendental (issue #409): `ppm` and `percent` are
    /// dimensionless too, and `log(x [ppm])` has two defensible readings —
    /// the log of the ppm NUMBER, or the log of the mole fraction — that differ
    /// by `ln(1e-6) = 13.8155…`. Dimension alone cannot tell them apart, so the
    /// checker refuses rather than picking one.
    pub fn is_pure_number(&self) -> bool {
        self.dimensions.is_empty() && self.exact.is_one()
    }

    /// True when this unit is a plane angle to the FIRST power (`rad`, `deg`,
    /// `mrad`) — the shape a circular trig function's argument may carry.
    ///
    /// The exponent matters because the argument is CONVERTED to radians
    /// (issue #409): `rad^2` is not an angle any conversion turns into one, and
    /// accepting it would multiply by `scale` where `scale^2` was meant.
    pub fn is_plane_angle(&self) -> bool {
        self.dimensions.len() == 1
            && self
                .dimensions
                .get(&Dimension::Angle)
                .is_some_and(|p| *p == Rational::int(1))
    }

    /// Multiply two units
    pub fn multiply(&self, other: &Unit) -> Unit {
        let mut unit = self.combine(other, Rational::add, self.scale * other.scale);
        unit.exact = self.exact.multiply(&other.exact);
        unit
    }

    /// Divide two units
    pub fn divide(&self, other: &Unit) -> Unit {
        let mut unit = self.combine(other, Rational::sub, self.scale / other.scale);
        unit.exact = self.exact.divide(&other.exact);
        unit
    }

    /// Merge `other`'s exponents into a copy of `self`'s with `op`, dropping any
    /// axis whose exponent cancels to zero.
    fn combine(&self, other: &Unit, op: fn(Rational, Rational) -> Rational, scale: f64) -> Unit {
        let mut dimensions = self.dimensions.clone();
        for (dim, power) in &other.dimensions {
            let entry = dimensions.entry(dim.clone()).or_insert(Rational::int(0));
            *entry = op(*entry, *power);
            if entry.is_zero() {
                dimensions.remove(dim);
            }
        }
        Unit {
            dimensions,
            scale,
            exact: ExactScale::one(),
        }
    }

    /// Raise unit to an integer power.
    pub fn power(&self, exponent: i32) -> Unit {
        self.power_rational(Rational::int(exponent))
    }

    /// Raise unit to a RATIONAL power. `sqrt` is `power_rational(1/2)`, and a
    /// literal `^0.5` in an expression lands here too.
    pub fn power_rational(&self, exponent: Rational) -> Unit {
        let dimensions = self
            .dimensions
            .iter()
            .map(|(dim, power)| (dim.clone(), power.mul(exponent)))
            .filter(|(_, power)| !power.is_zero())
            .collect();

        Unit {
            dimensions,
            scale: self.scale.powf(exponent.as_f64()),
            exact: self.exact.power(exponent),
        }
    }

    /// Propagate units through an expression AST.
    ///
    /// Given a variable→unit environment, returns the [`Unit`] of the whole
    /// expression by recursively walking its AST. Returns a
    /// [`UnitError::DimensionMismatch`] when the expression is dimensionally
    /// inconsistent (e.g. adding `m` to `s`, using a non-dimensionless
    /// exponent, or passing a dimensional argument to `exp`/`log`/`sin`).
    ///
    /// The `t` identifier is treated as the implicit time variable with units
    /// of seconds if not otherwise present in `env`.
    ///
    /// This is the single-outcome (`Result`) view, kept for callers that only
    /// need "does this expression have a unit". Validation uses
    /// [`check_expression_dimensions`], which reports EVERY finding rather than
    /// only the first.
    ///
    /// # Arguments
    ///
    /// * `expr` - Expression to analyse.
    /// * `env`  - Map from variable name to its [`Unit`].
    pub fn propagate(expr: &Expr, env: &HashMap<String, Unit>) -> Result<Unit, UnitError> {
        Self::propagate_in(expr, env)
    }

    /// [`Unit::propagate`] against any [`UnitEnv`].
    pub(crate) fn propagate_in(expr: &Expr, env: &dyn UnitEnv) -> Result<Unit, UnitError> {
        let mut findings = Vec::new();
        let dim = propagate_dim(expr, env, &mut findings);
        // A provable inconsistency outranks an undeterminable dimension.
        if let Some(err) = findings.iter().find(|f| f.is_error()) {
            return Err(UnitError::DimensionMismatch(err.message.clone()));
        }
        match dim {
            Dim::Known(unit) => Ok(unit),
            Dim::Unknown => Err(UnitError::UnknownUnit(
                findings
                    .first()
                    .map(|f| f.message.clone())
                    .unwrap_or_else(|| "undetermined dimension".to_string()),
            )),
        }
    }
}

/// Propagate the dimension of `expr`, recording every finding in `findings`.
///
/// Never aborts: an undeterminable subtree yields [`Dim::Unknown`] and the walk
/// continues, so all provable mismatches in the expression are reported.
fn propagate_dim(expr: &Expr, env: &dyn UnitEnv, findings: &mut Vec<UnitFinding>) -> Dim {
    match expr {
        // A BARE NUMERIC LITERAL has an INDETERMINATE dimension, not a
        // dimensionless one: nothing in the AST says whether `273.15` is a pure
        // number or a temperature offset, or whether `6.022e23` is Avogadro's
        // constant. Calling it dimensionless manufactured a mismatch on every
        // line of the valid corpus that uses an implicit-unit constant
        // (`T - 273.15`, `1 - phi`).
        //
        // This matters BECAUSE these findings are now hard errors: a checker
        // that fails the build must not fabricate a dimension it cannot know.
        // It costs nothing on the pinned invalid corpus — every fixture there
        // states its inconsistency between DECLARED quantities (`length + mass`,
        // `ln(mass)`, `m^kg`), never via a literal. Literals still behave
        // correctly wherever their meaning IS determined: additively they are
        // neutral and adopt their sibling's dimension, an all-literal expression
        // is dimensionless, and an exponent is read BY VALUE.
        Expr::Number(_) | Expr::Integer(_) => Dim::Unknown,
        Expr::Variable(name) => {
            if let Some(unit) = env.unit_of(name) {
                Dim::Known(unit.clone())
            } else if name == "t" {
                // Implicit time variable.
                Dim::Known(Unit::base(Dimension::Time, 1, 1.0))
            } else {
                // Dimension unknown — but NOT reported here, because the reason
                // is always already reported at the variable itself, and this
                // path fires once per REFERENCE:
                //   * genuinely undeclared  ⇒ hard `undefined_variable` error;
                //   * unparseable unit      ⇒ `build_unit_env` warns once;
                //   * no `units` field      ⇒ a legitimate choice, not a defect.
                // Re-reporting per reference would only duplicate those.
                Dim::Unknown
            }
        }
        Expr::Operator(op) => propagate_operator_dim(op, env, findings),
    }
}

/// Propagate every argument, returning their dimensions positionally.
fn propagate_args(
    op: &ExpressionNode,
    env: &dyn UnitEnv,
    findings: &mut Vec<UnitFinding>,
) -> Vec<Dim> {
    op.args
        .iter()
        .map(|a| propagate_dim(a, env, findings))
        .collect()
}

/// Require exactly `n` arguments. Arity is enforced by the operator registry
/// and the schema; here a wrong arity simply means we cannot analyse the node,
/// so it is an `Analysis` finding rather than a dimensional defect.
fn require_arity(op: &ExpressionNode, n: usize, findings: &mut Vec<UnitFinding>) -> bool {
    if op.args.len() == n {
        return true;
    }
    findings.push(UnitFinding::analysis(format!(
        "'{}' expects {} argument(s) but has {}; not dimension-checked",
        op.op,
        n,
        op.args.len()
    )));
    false
}

/// Like [`require_arity`], but a LOWER bound: the node is well-formed when it has
/// at least `n` operands. Used for the ops whose operand list is open-ended at
/// this end — a rewrite-target `D` carrying trailing auxiliary boundary operands
/// (esm-spec §4.2 "Arity of `D`").
fn require_min_arity(op: &ExpressionNode, n: usize, findings: &mut Vec<UnitFinding>) -> bool {
    if op.args.len() >= n {
        return true;
    }
    findings.push(UnitFinding::analysis(format!(
        "'{}' expects at least {} argument(s) but has {}; not dimension-checked",
        op.op,
        n,
        op.args.len()
    )));
    false
}

/// Dispatch one operator node to its family's propagation helper. Kept as a
/// thin match so each family's dimensional rules live in one named function;
/// the rules mirror the Python (`units.py::_get_expression_dimension`) and
/// Julia (`units.jl::get_expression_dimensions`) implementations.
fn propagate_operator_dim(
    op: &ExpressionNode,
    env: &dyn UnitEnv,
    findings: &mut Vec<UnitFinding>,
) -> Dim {
    match op.op.as_str() {
        "+" | "-" => propagate_additive_dim(op, env, findings),
        "*" | "/" => propagate_multiplicative_dim(op, env, findings),
        "^" | "power" | "pow" => propagate_power_dim(op, env, findings),
        "D" | "ic" => propagate_calculus_dim(op, env, findings),
        "exp" | "log" | "log10" | "ln" | "sin" | "cos" | "tan" | "asin" | "acos" | "atan"
        | "sinh" | "cosh" | "tanh" | "asinh" | "acosh" | "atanh" => {
            propagate_transcendental_dim(op, env, findings)
        }
        "sqrt" => propagate_sqrt_dim(op, env, findings),
        // abs/sign/floor/ceil preserve dimensions; `Pre` is an initial-value
        // marker — same dimensions as its argument.
        "abs" | "floor" | "ceil" | "round" | "sign" | "Pre" => {
            if !require_arity(op, 1, findings) {
                return Dim::Unknown;
            }
            propagate_dim(&op.args[0], env, findings)
        }
        "min" | "max" => propagate_matching_dim(op, env, findings),
        "atan2" => {
            if !require_arity(op, 2, findings) {
                return Dim::Unknown;
            }
            // Both arguments must share dimensions (their RATIO is what is fed
            // to the arctangent). The result is an ANGLE — `atan2(y, x)` is the
            // canonical spelling of a bearing.
            propagate_matching_dim(op, env, findings);
            Dim::Known(radian())
        }
        "ifelse" => propagate_ifelse_dim(op, env, findings),
        ">" | "<" | ">=" | "<=" | "==" | "!=" => {
            if !require_arity(op, 2, findings) {
                return Dim::Known(Unit::dimensionless());
            }
            // Operands must be comparable; the flag itself is dimensionless.
            propagate_matching_dim(op, env, findings);
            Dim::Known(Unit::dimensionless())
        }
        "and" | "or" | "not" => {
            // A boolean connective's result is a dimensionless boolean. Its
            // operands carry no unit requirement of their own, but they are
            // walked so a mismatch inside one (`not(x [m] > z [kg])`) is
            // reported (esm-spec §4.8.3).
            propagate_args(op, env, findings);
            Dim::Known(Unit::dimensionless())
        }
        // Array operators: propagate the element dimension. Shape and
        // indexing are orthogonal to dimension (see gt-t5c / gt-vt3 — shapes
        // are a separate concern from unit checking).
        "faq" | "makearray" | "index" | "reshape" | "transpose" | "concat" | "broadcast" => {
            propagate_array_dim(op, env, findings)
        }
        // A `const` that DECLARES its units has that unit (esm-spec §4.8.5);
        // without `units` it is undeterminable, like a bare literal. An
        // unresolvable string is reported at the containing expression field by
        // `structural.rs` (see [`unresolvable_const_units`]), not here.
        "const" => match op.units.as_deref().map(parse_unit) {
            Some(Ok(unit)) => Dim::Known(unit),
            _ => Dim::Unknown,
        },
        // No dimensional rule for this operator — an unregistered user op, or a
        // rewrite-target sugar op (`grad`/`div`/`laplacian`/`fn`/`table_lookup`/
        // `godunov_hamiltonian`/…) whose dimension is UNDETERMINABLE until a
        // discretization rule lowers it (esm-spec §4.8.4). The checker reports
        // UNKNOWN and SKIPS: it emits NO finding, matching the other four
        // bindings (Julia/TS/Python/Go), which return an unknown dimension and
        // let their callers skip the check. The node's dimension is `Unknown` —
        // NOT a silent `dimensionless`, which would manufacture false mismatches
        // in the enclosing expression. `findings` is left untouched so nothing
        // is singled out; `_` treats `grad` exactly like any other no-rule op.
        _ => Dim::Unknown,
    }
}

/// Report a provable mismatch among the operands that DID resolve, and return
/// the shared dimension. Used by `+`/`-`, `min`/`max`, comparisons and `atan2`.
///
/// Two operands are only ever compared when BOTH dimensions were determined —
/// an undeterminable operand is skipped, never assumed dimensionless, so it can
/// neither hide nor manufacture a mismatch (esm-spec §4.8.3). A numeric literal
/// is undeterminable, so it adopts the dimension of what it is combined with
/// (`T - 273.15` is a temperature), and `x + 0.5 * y` has the dimension of `x`.
/// If NO operand resolves (`1 + 2`), the result is undeterminable, never
/// dimensionless (§4.8.4).
fn propagate_matching_dim(
    op: &ExpressionNode,
    env: &dyn UnitEnv,
    findings: &mut Vec<UnitFinding>,
) -> Dim {
    let mut first: Option<Unit> = None;
    for arg in &op.args {
        let dim = propagate_dim(arg, env, findings);
        let Some(unit) = dim.known() else {
            continue;
        };
        match &first {
            None => first = Some(unit.clone()),
            Some(f) if !f.is_compatible(unit) => {
                findings.push(UnitFinding::error(format!(
                    "Incompatible dimensions in '{}': {} vs {}",
                    op.op,
                    describe(f),
                    describe(unit)
                )));
            }
            Some(f) if f.exact != unit.exact => {
                findings.push(UnitFinding::error(format!(
                    "Incompatible scales in '{}': {} vs {}",
                    op.op,
                    describe_scaled(f),
                    describe_scaled(unit)
                )));
            }
            _ => {}
        }
    }
    match first {
        Some(unit) => Dim::Known(unit),
        None => Dim::Unknown,
    }
}

/// `+` / `-`: every operand must share dimensions; the result carries them.
/// A unary `+` or `-` carries its single operand's dimension unchanged, so
/// `+(2)` is as undeterminable as `2`.
fn propagate_additive_dim(
    op: &ExpressionNode,
    env: &dyn UnitEnv,
    findings: &mut Vec<UnitFinding>,
) -> Dim {
    propagate_matching_dim(op, env, findings)
}

/// `*` / `/`: dimensions multiply / divide, no compatibility requirement.
fn propagate_multiplicative_dim(
    op: &ExpressionNode,
    env: &dyn UnitEnv,
    findings: &mut Vec<UnitFinding>,
) -> Dim {
    let dims = propagate_args(op, env, findings);
    if op.op == "*" {
        let mut result = Unit::dimensionless();
        for d in &dims {
            match d.known() {
                Some(u) => result = result.multiply(u),
                None => return Dim::Unknown,
            }
        }
        return Dim::Known(result);
    }
    if !require_arity(op, 2, findings) {
        return Dim::Unknown;
    }
    match (dims[0].known(), dims[1].known()) {
        (Some(num), Some(den)) => Dim::Known(num.divide(den)),
        _ => Dim::Unknown,
    }
}

/// `^` / `power` / `pow`: the exponent must be dimensionless; a dimensional
/// base additionally requires a literal integer exponent (only integer powers
/// of dimensions are representable).
fn propagate_power_dim(
    op: &ExpressionNode,
    env: &dyn UnitEnv,
    findings: &mut Vec<UnitFinding>,
) -> Dim {
    if !require_arity(op, 2, findings) {
        return Dim::Unknown;
    }
    let base = propagate_dim(&op.args[0], env, findings);
    let exp = propagate_dim(&op.args[1], env, findings);

    // A dimensional exponent is provably wrong regardless of the base.
    if let Some(e) = exp.known()
        && !e.is_dimensionless()
    {
        findings.push(UnitFinding::error(format!(
            "Exponent in '{}' must be dimensionless, got {}",
            op.op,
            describe(e)
        )));
        return Dim::Unknown;
    }

    let Some(base_unit) = base.known() else {
        return Dim::Unknown;
    };
    // A dimensionless base of scale 1 stays exactly that under any exponent, so
    // a non-literal exponent is only a problem for a DIMENSIONAL or SCALED base
    // (`x^2` with `x` in `%` has scale 1/10000).
    if base_unit.is_dimensionless() && base_unit.exact.is_one() {
        return Dim::Known(Unit::dimensionless());
    }

    // The exponent is read BY VALUE, not by dimension — a literal's dimension is
    // indeterminate, but its VALUE is exactly what a power needs. A constant
    // expression over literals (`3/2`, `-(1)`) has an exact value too, and is
    // folded to it (esm-spec §4.8.3 `^`). Exponents are RATIONAL, so a
    // fractional power of a dimensional quantity (`x^0.5`) propagates exactly.
    if let Some(r) = constant_exponent(&op.args[1]) {
        return Dim::Known(base_unit.power_rational(r));
    }
    match &op.args[1] {
        Expr::Integer(_) | Expr::Number(_) => {
            let n = literal_f64(&op.args[1]).unwrap_or(f64::NAN);
            findings.push(UnitFinding::analysis(format!(
                "Exponent {n} is not a representable rational, so the dimension of a \
                 dimensional base is unknown"
            )));
        }
        // Symbolic exponent (`x^n` for a parameter `n`): the resulting
        // dimension depends on a value we do not have. Undeterminable, NOT a
        // mismatch (§4.8.4).
        _ => {
            findings.push(UnitFinding::analysis(format!(
                "Exponent of '{}' is not a constant expression, so the dimension of a \
                 dimensional base is unknown",
                op.op
            )));
        }
    }
    Dim::Unknown
}

/// A numeric literal's value. It may arrive as EITHER `Expr::Number` (JSON
/// `2.0`) or `Expr::Integer` (JSON `2`).
fn literal_f64(expr: &Expr) -> Option<f64> {
    match expr {
        Expr::Integer(i) => Some(*i as f64),
        Expr::Number(n) => Some(*n),
        _ => None,
    }
}

/// The exact rational value of a CONSTANT-EXPRESSION exponent (esm-spec §4.8.3
/// `^`), or `None` when the exponent is symbolic.
///
/// Folds: a numeric literal; a `const` that declares no `units`, or a
/// dimensionless unit at exact scale 1 (its value IS the pure number it
/// denotes); and `+`, `*` (one or more operands), `-` (unary negation, or a left
/// fold of two or more) and `/` (exactly two operands, nonzero divisor) over
/// constant expressions, nested to any depth. Anything else — a variable, an
/// expression containing one, a `const` at a scale other than 1 (`%`: its value
/// and the quantity it denotes disagree), any other operator, an overflow of
/// the rational representation — is symbolic, which leaves the power of a
/// dimensional base undeterminable (§4.8.4).
pub fn constant_exponent(expr: &Expr) -> Option<Rational> {
    if let Some(n) = literal_f64(expr) {
        return Rational::from_f64(n);
    }
    let Expr::Operator(node) = expr else {
        return None;
    };
    if node.op == "const" {
        if let Some(units) = &node.units {
            let unit = parse_unit(units).ok()?;
            if !unit.is_pure_number() {
                return None;
            }
        }
        return node
            .value
            .as_ref()
            .and_then(serde_json::Value::as_f64)
            .and_then(Rational::from_f64);
    }
    let values = node
        .args
        .iter()
        .map(constant_exponent)
        .collect::<Option<Vec<Rational>>>()?;
    let (first, rest) = values.split_first()?;
    match node.op.as_str() {
        "+" => rest.iter().try_fold(*first, |acc, v| acc.checked_add(*v)),
        "*" => rest.iter().try_fold(*first, |acc, v| acc.checked_mul(*v)),
        "-" if rest.is_empty() => Rational::int(0).checked_sub(*first),
        "-" => rest.iter().try_fold(*first, |acc, v| acc.checked_sub(*v)),
        "/" if rest.len() == 1 => first.checked_div(rest[0]),
        _ => None,
    }
}

/// Calculus-family operators: the time/coordinate derivative `D` divides its
/// argument's dimensions by the `wrt` variable's, and the initial-condition
/// marker `ic` (v0.8.0) is dimension-preserving — the initial value of a field
/// carries the same units as the field itself.
fn propagate_calculus_dim(
    op: &ExpressionNode,
    env: &dyn UnitEnv,
    findings: &mut Vec<UnitFinding>,
) -> Dim {
    // esm-spec §4.2 "Arity of `D`": a REWRITE-TARGET `D` (spatial `wrt`) MAY
    // carry trailing auxiliary operands after `args[0]` — the per-face
    // boundary/halo values a discretization rule binds and consumes. They carry
    // no evaluator semantics, so the dimensional rule reads `args[0]` alone and
    // MUST NOT report them as an arity defect. The STRUCTURAL time derivative
    // (and `ic`) stay strictly unary.
    let well_formed = if crate::op_registry::is_rewrite_target_derivative(op) {
        require_min_arity(op, 1, findings)
    } else {
        require_arity(op, 1, findings)
    };
    if !well_formed {
        return Dim::Unknown;
    }
    let arg = propagate_dim(&op.args[0], env, findings);
    if op.op == "ic" {
        return arg;
    }
    let wrt = crate::op_registry::derivative_wrt(op);
    // The independent variable's unit is only known if it was DECLARED. We do
    // not assume an undeclared `t` means seconds: that would give `D(x)/dt` a
    // dimension of `[x]/second`, and every dimensionless toy model
    // (`D(x)/dt = -x`, the idiom used across the shared corpus) would then look
    // like a mismatch of `time^-1` against `dimensionless`. An undeclared
    // independent variable leaves the derivative INDETERMINATE here; the
    // equation-level time-ratio rule in `check_equation_dimensions` still
    // catches derivative equations that no time unit could reconcile.
    let (Some(arg_unit), Some(wrt_unit)) = (arg.known(), env.unit_of(wrt)) else {
        return Dim::Unknown;
    };
    Dim::Known(arg_unit.divide(wrt_unit))
}

/// True when `expr` is `D(x)` taken with respect to an UNDECLARED independent
/// variable, which makes its dimension indeterminate. Mirrors the TypeScript
/// reference (`units.ts::derivativeOfUndeclaredTime`).
fn derivative_of_undeclared_time<'a>(expr: &'a Expr, env: &dyn UnitEnv) -> Option<&'a Expr> {
    let Expr::Operator(op) = expr else {
        return None;
    };
    if op.op != "D" || op.args.len() != 1 {
        return None;
    }
    let wrt = crate::op_registry::derivative_wrt(op);
    (!env.contains(wrt)).then(|| &op.args[0])
}

/// The weaker-but-still-provable rule for `d(state)/dt = rhs` when the time
/// unit is unknown: whatever time unit `t` carries, it can only ever contribute
/// powers of TIME to the ratio. So if `dims(state) / dims(rhs)` contains any
/// NON-time dimension, no choice of time unit could reconcile the two sides and
/// the equation is provably inconsistent. A ratio that is a pure power of time
/// (including time^0) is accepted. Mirrors `units.ts::derivativeTimeMismatch`.
fn derivative_time_mismatch(state: &Unit, rhs: &Unit) -> Option<String> {
    let ratio = state.divide(rhs);
    let unreconcilable = ratio
        .dimensions
        .iter()
        .any(|(d, p)| *d != Dimension::Time && !p.is_zero());
    unreconcilable.then(|| {
        format!(
            "No time unit can reconcile d({})/dt with {} (their ratio {} is not a power of time)",
            describe(state),
            describe(rhs),
            describe(&ratio)
        )
    })
}

/// Transcendental and trigonometric functions: argument must be dimensionless,
/// result is dimensionless.
///
/// `sqrt` is NOT in this family — it halves its argument's dimensions (see
/// [`propagate_sqrt_dim`]).
///
/// The TRIGONOMETRIC members additionally accept a plane ANGLE at ANY scale:
/// `rad` is a real dimension axis here, so `sin(theta)` with `theta` in radians
/// — or in `deg` — is correct, not an error, and the flatten path converts the
/// scaled spelling to radians before anything evaluates it. Everything else
/// demands a PURE NUMBER: dimensionless AND at scale 1 (issue #409).
fn propagate_transcendental_dim(
    op: &ExpressionNode,
    env: &dyn UnitEnv,
    findings: &mut Vec<UnitFinding>,
) -> Dim {
    // `rad` is a real axis here, so the circular functions are NOT symmetric —
    // and getting this wrong is silently destructive in both directions:
    //
    // * FORWARD circular (`sin`/`cos`/`tan`) map an ANGLE to a ratio: they
    //   accept an angle OR a dimensionless argument, and return DIMENSIONLESS.
    //   (`sin(kg)` is still an error.)
    // * INVERSE circular (`asin`/`acos`/`atan`) map a ratio to an ANGLE: they
    //   require a DIMENSIONLESS argument and RETURN AN ANGLE. Asserting a
    //   dimensionless result while `rad` is an axis makes every
    //   `zenith: "rad"` computed by `acos(...)` a guaranteed false mismatch.
    // * The HYPERBOLIC family takes and returns pure numbers — a hyperbolic
    //   "angle" is an area, not a plane angle, so `rad` is not involved.
    // * `exp`/`log` demand strict dimensionlessness: their Taylor series adds
    //   powers of the argument.
    let is_inverse_circular = matches!(op.op.as_str(), "asin" | "acos" | "atan");
    let result = if is_inverse_circular {
        radian()
    } else {
        Unit::dimensionless()
    };

    if !require_arity(op, 1, findings) {
        // The result is fixed by the operator whatever the argument turns out
        // to be.
        return Dim::Known(result);
    }

    // Only the forward circular functions accept an angle.
    let angle_ok = matches!(op.op.as_str(), "sin" | "cos" | "tan");

    let arg = propagate_dim(&op.args[0], env, findings);
    if let Some(u) = arg.known() {
        // Three outcomes, and the middle one is the point of issue #409.
        //
        // * A PURE NUMBER (dimensionless at scale 1) is always admissible.
        // * An ANGLE is admissible to `sin`/`cos`/`tan` WHATEVER its scale,
        //   because `deg` → `rad` is exact and has no second reading: the
        //   evaluation path converts it (`normalize_angle_arguments`).
        // * Dimensionless at a scale OTHER than 1 (`ppm`, `percent`) is
        //   REFUSED. It passes any dimension-only test, which is how
        //   `log(x [ppm])` evaluated `log` of the ppm number with no
        //   diagnostic; the two readings differ by `ln(1e-6)` and nothing in
        //   the document says which was meant.
        if angle_ok && u.is_plane_angle() {
            // Converted, not refused.
        } else if u.is_pure_number() {
            // Nothing to say.
        } else if u.is_dimensionless() {
            findings.push(UnitFinding::error(scaled_dimensionless_message(&op.op, u)));
        } else {
            findings.push(UnitFinding::error(format!(
                "Argument to '{}' must be dimensionless{}, got {}",
                op.op,
                if angle_ok { " (or an angle)" } else { "" },
                describe(u)
            )));
        }
    }
    Dim::Known(result)
}

/// The factor that brings a quantity in `unit` to RADIANS, when `unit` is a
/// plane angle at a scale other than 1 (esm-spec §4.8.3, issue #409).
///
/// `None` for anything else — a pure number, a `rad`, a dimensional unit — so a
/// document that declares no scaled angle is rewritten not at all.
#[must_use]
pub fn angle_normalization_factor(unit: &Unit) -> Option<f64> {
    (unit.is_plane_angle() && !unit.exact.is_one()).then_some(unit.scale)
}

/// Rewrite every `sin`/`cos`/`tan` whose argument is an angle at a scale other
/// than 1 so the argument reaches the evaluator in RADIANS.
///
/// esm-spec §4.8.3 names angles as *the* exception to the dimensionless-argument
/// rule, and `deg` is a registry unit at scale π/180 — so `sin(theta)` with
/// `theta` in `deg` is a CONFORMING document, and before this it evaluated
/// `sin(90)` = 0.894 where 1 was meant, with no diagnostic (issue #409). The
/// conversion is exact and has no second reading, which is why this half
/// converts where the `ppm` half refuses.
///
/// Returns `None` when nothing changed, so the common case clones nothing.
/// The rewrite runs on the EVALUATION path only (the flatten funnel); the
/// checker still sees the authored spelling, and so still reports the argument's
/// declared unit rather than a literal-poisoned `unknown`.
#[must_use]
pub fn normalize_angle_arguments(expr: &Expr, env: &HashMap<String, Unit>) -> Option<Expr> {
    let Expr::Operator(node) = expr else {
        return None;
    };

    // Rewrite the children first, so a nested `sin(deg)` is converted too —
    // EVERY expression-bearing child (`ExpressionNode::for_each_child`), not a
    // hand list: a trig call in a `makearray` region value, a `table_lookup`
    // axis input, a `faq` filter or an aggregate key is as much an evaluated
    // argument as one in `args`. The node is cloned only when something in it
    // changes.
    let mut children: Vec<Option<Expr>> = Vec::new();
    node.for_each_child(&mut |child| children.push(normalize_angle_arguments(child, env)));

    // Then this node: a circular trig whose (possibly rewritten) argument is a
    // scaled angle gets the factor folded in. `args` are the first children
    // visited, so `children[0]` is the rewritten `args[0]`.
    let mut factor = None;
    if matches!(node.op.as_str(), "sin" | "cos" | "tan") && node.args.len() == 1 {
        let arg = children[0].as_ref().unwrap_or(&node.args[0]);
        let mut findings = Vec::new();
        factor = propagate_dim(arg, env, &mut findings)
            .known()
            .and_then(angle_normalization_factor);
    }

    if factor.is_none() && children.iter().all(Option::is_none) {
        return None;
    }
    let mut out = (**node).clone();
    let mut rewritten = children.into_iter();
    out.for_each_child_mut(&mut |child| {
        if let Some(Some(next)) = rewritten.next() {
            *child = next;
        }
    });
    if let Some(factor) = factor {
        let scaled = ExpressionNode {
            op: "*".to_string(),
            args: vec![out.args[0].clone(), Expr::Number(factor)],
            ..Default::default()
        };
        out.args[0] = Expr::operator(scaled);
    }
    Some(Expr::operator(out))
}

/// The radian — the unit of the `Angle` axis, returned by the inverse circular
/// functions.
fn radian() -> Unit {
    Unit::base(Dimension::Angle, 1, 1.0)
}

/// The registry spellings of the dimensionless-but-scaled units, most common
/// first, used to NAME a scale in a diagnostic (`… (ppm); divide by 1 ppm …`).
///
/// A fixed, ordered table rather than a reverse sweep of the registry: the
/// registry is a `HashMap`, so a sweep would pick `ppmv` or `ppm` depending on
/// iteration order and the five bindings would print different messages for the
/// same document. The same table, in the same order, is in every binding.
const SCALED_DIMENSIONLESS_SPELLINGS: &[(&str, i32)] =
    &[("percent", -2), ("ppm", -6), ("ppb", -9), ("ppt", -12)];

/// The canonical registry spelling of a dimensionless unit at `exact`, when one
/// exists (`1/100` → `percent`, `1e-6` → `ppm`).
#[must_use]
pub fn scaled_dimensionless_spelling(exact: &ExactScale) -> Option<&'static str> {
    SCALED_DIMENSIONLESS_SPELLINGS
        .iter()
        .find(|(_, k)| ExactScale::power_of_ten(*k) == *exact)
        .map(|(name, _)| *name)
}

/// The §4.8.3 refusal for a dimensionless-but-SCALED argument to an op that
/// requires a pure number (esm-spec §9.6.6 `unit_transcendental_scaled_argument`).
///
/// It names the REPAIR, not only the refusal: the author states the reading by
/// dividing by a quantity carrying the scale, which costs one node and records
/// the decision in the document. Normalizing silently instead would change the
/// numbers of every document that already passes a `percent` or a `ppm` into
/// `exp`/`log`.
fn scaled_dimensionless_message(op: &str, unit: &Unit) -> String {
    match scaled_dimensionless_spelling(&unit.exact) {
        Some(name) => format!(
            "Argument to '{op}' must be dimensionless at scale 1, but is \
             dimensionless at scale {} ({name}); divide by 1 {name}, or by the \
             scale you mean, to state which reading is intended",
            unit.exact
        ),
        None => format!(
            "Argument to '{op}' must be dimensionless at scale 1, but is \
             dimensionless at scale {}; divide by a quantity carrying that \
             scale to state which reading is intended",
            unit.exact
        ),
    }
}

/// Square root HALVES its argument's dimensions — it is not a transcendental
/// function and does NOT require a dimensionless argument. Because exponents are
/// rational, `sqrt(m^3)` is exactly `m^3/2` rather than "not representable".
fn propagate_sqrt_dim(
    op: &ExpressionNode,
    env: &dyn UnitEnv,
    findings: &mut Vec<UnitFinding>,
) -> Dim {
    if !require_arity(op, 1, findings) {
        return Dim::Unknown;
    }
    let arg = propagate_dim(&op.args[0], env, findings);
    match arg.known() {
        Some(unit) => Dim::Known(unit.power_rational(Rational::new(1, 2))),
        None => Dim::Unknown,
    }
}

/// `ifelse`: the two branches follow the `+` rule (esm-spec §4.8.3). When both
/// are determinable they must share dimension and scale; when only one is, the
/// result is its unit (`ifelse(c, x, 0.5 * y)` has the unit of `x`); when
/// neither is (`ifelse(c, 1, 2)`), the result is undeterminable.
fn propagate_ifelse_dim(
    op: &ExpressionNode,
    env: &dyn UnitEnv,
    findings: &mut Vec<UnitFinding>,
) -> Dim {
    if !require_arity(op, 3, findings) {
        return Dim::Unknown;
    }
    // The condition (arg 0) need not be dimensionless — comparison ops already
    // produce a dimensionless Boolean, and we don't want to reject bare scalars
    // used as truthiness flags. It is still walked so findings inside it (e.g.
    // a mismatched comparison) are reported.
    propagate_dim(&op.args[0], env, findings);
    let t = propagate_dim(&op.args[1], env, findings);
    let f = propagate_dim(&op.args[2], env, findings);
    match (t.known(), f.known()) {
        (Some(a), Some(b)) if !a.is_compatible(b) => {
            findings.push(UnitFinding::error(format!(
                "'ifelse' branches must share dimensions: {} vs {}",
                describe(a),
                describe(b)
            )));
            Dim::Unknown
        }
        (Some(a), Some(b)) if a.exact != b.exact => {
            findings.push(UnitFinding::error(format!(
                "'ifelse' branches must share a scale: {} vs {}",
                describe_scaled(a),
                describe_scaled(b)
            )));
            Dim::Unknown
        }
        (Some(a), _) | (None, Some(a)) => Dim::Known(a.clone()),
        (None, None) => Dim::Unknown,
    }
}

/// Array operators carry their ELEMENT's unit (esm-spec §4.8.3, "Array
/// operators"). Shape and indexing are orthogonal to dimension (gt-t5c /
/// gt-vt3); what an array holds is not. The same rules serve the checker and
/// the evaluation path's angle conversion, so the two cannot disagree.
///
/// * `index`, `reshape`, `transpose` — the unit of the source array
///   (`args[0]`); the subscripts and shape operands are not part of it, but
///   are walked for findings like any other operand.
/// * `concat` — its operands agree like those of `+`.
/// * `makearray` — its region `values` agree like the operands of `+`: an
///   undeterminable region is skipped, two determinable ones that disagree
///   are a provable mismatch.
/// * `broadcast` — the unit of `{op: fn, args}`.
/// * `faq` — decided by its ⊕ (see [`faq_unit`]); `filter`, `expr` and any
///   operands are always walked for findings.
///
/// Every operand of every array operator is walked, so a mismatch inside a
/// subscript (`index(lat, a [m] + b [kg])`) is reported.
///
/// A node with nothing to read a unit from is UNDETERMINABLE, never
/// dimensionless (§4.8.4).
fn propagate_array_dim(
    op: &ExpressionNode,
    env: &dyn UnitEnv,
    findings: &mut Vec<UnitFinding>,
) -> Dim {
    match op.op.as_str() {
        "faq" => {
            // The `filter` is a boolean over the node's ranges. Its value never
            // reaches the element unit, but it is an expression like any other,
            // so a comparison in it is checked by the §4.8.3 rules exactly as
            // the same comparison in the body is (EarthSciAST#414).
            if let Some(filter) = &op.filter {
                propagate_dim(filter, env, findings);
            }
            let body = op
                .expr
                .as_ref()
                .map(|body| propagate_dim(body, env, findings));
            propagate_args(op, env, findings);
            match body {
                Some(body) => faq_unit(op, body),
                None => Dim::Unknown,
            }
        }
        "makearray" => {
            propagate_args(op, env, findings);
            let mut first: Option<Unit> = None;
            for value in op.values.iter().flatten() {
                let dim = propagate_dim(value, env, findings);
                let Some(unit) = dim.known() else {
                    continue;
                };
                match &first {
                    None => first = Some(unit.clone()),
                    Some(f) if !f.same_unit(unit) => {
                        findings.push(UnitFinding::error(format!(
                            "'makearray' regions must share units: {} vs {}",
                            describe_scaled(f),
                            describe_scaled(unit)
                        )));
                    }
                    _ => {}
                }
            }
            first.map_or(Dim::Unknown, Dim::Known)
        }
        "concat" => propagate_matching_dim(op, env, findings),
        "broadcast" => {
            // Elementwise map over arrays with `fn` naming the scalar
            // operator. Construct a synthetic scalar node and recurse so we
            // reuse the same dimensional rules. (That synthetic node is the
            // same identity the EVALUATORS now implement: `broadcast(fn=F, A)`
            // means `{op: F, args: A}` applied element-wise — see
            // `op_registry::check_broadcast_fn`.)
            //
            // A missing or unusable `fn` is a HARD `invalid_broadcast_fn`
            // structural error raised by `structural.rs`; what is reported here
            // is only the consequence for the DIMENSION, which the checker
            // genuinely cannot determine.
            let Some(fn_name) = op.broadcast_fn.as_deref() else {
                propagate_args(op, env, findings);
                findings.push(UnitFinding::analysis(
                    "'broadcast' is missing its 'fn'; its dimension is unknown".to_string(),
                ));
                return Dim::Unknown;
            };
            let synthetic = ExpressionNode {
                op: fn_name.to_string(),
                args: op.args.clone(),
                ..ExpressionNode::default()
            };
            propagate_operator_dim(&synthetic, env, findings)
        }
        // "index" | "reshape" | "transpose": the source array's unit.
        _ => propagate_args(op, env, findings)
            .into_iter()
            .next()
            .unwrap_or(Dim::Unknown),
    }
}

/// The unit of an `faq` node whose body has unit `body` (esm-spec §4.8.3).
///
/// The ⊗ of the semiring is written by the author inside `expr`, so it is
/// already typed; `join` clauses are key equalities and contribute no factor.
/// What remains is the ⊕ reduction over the contracted indices:
///
/// * `+`, `max`, `min` (`sum_product`, `max_product`, `min_sum`, `max_sum`,
///   or `reduce` `+`/`max`/`min`) — the body's unit;
/// * `*` (`reduce: "*"`) — the body's unit only when the body is a pure
///   number, otherwise undeterminable: the power is the range length;
/// * `∨` (`bool_and_or`) — dimensionless, a boolean;
/// * value invention (`distinct`, or a `key`) — undeterminable: the result
///   is ids, not a quantity. So is an unknown ⊕ spelling, which is rejected
///   elsewhere.
fn faq_unit(op: &ExpressionNode, body: Dim) -> Dim {
    if op.distinct == Some(true) || op.key.is_some() {
        return Dim::Unknown;
    }
    match crate::faq::effective_reduce_kind(op.semiring.as_deref(), op.reduce.as_deref()) {
        Ok(
            crate::faq::ReduceKind::Sum | crate::faq::ReduceKind::Max | crate::faq::ReduceKind::Min,
        ) => body,
        Ok(crate::faq::ReduceKind::Product) => {
            if body.known().is_some_and(Unit::is_pure_number) {
                body
            } else {
                Dim::Unknown
            }
        }
        Ok(crate::faq::ReduceKind::Or | crate::faq::ReduceKind::And) => {
            Dim::Known(Unit::dimensionless())
        }
        Err(_) => Dim::Unknown,
    }
}

/// Render a unit's dimensions for a diagnostic message.
fn describe(unit: &Unit) -> String {
    if unit.is_dimensionless() {
        return "dimensionless".to_string();
    }
    let mut parts: Vec<String> = unit
        .dimensions
        .iter()
        .map(|(d, p)| {
            let name = match d {
                Dimension::Mass => "mass",
                Dimension::Length => "length",
                Dimension::Time => "time",
                Dimension::Current => "current",
                Dimension::Temperature => "temperature",
                Dimension::Amount => "amount",
                Dimension::Luminosity => "luminosity",
                Dimension::Angle => "angle",
            };
            if *p == Rational::int(1) {
                name.to_string()
            } else {
                format!("{name}^{p}")
            }
        })
        .collect();
    parts.sort();
    parts.join("*")
}

/// [`describe`], followed by the exact scale when it is not 1 — what a scale
/// disagreement has to show, since the dimensions alone are identical.
fn describe_scaled(unit: &Unit) -> String {
    if unit.exact.is_one() {
        describe(unit)
    } else {
        format!("{} (scale {})", describe(unit), unit.exact)
    }
}

/// Build a `HashMap<String, Unit>` environment from model variable metadata,
/// together with a list of warnings for any variables whose declared unit
/// string could not be parsed.
///
/// A variable is entered into the environment ONLY when its dimension is
/// actually known. Both an *unparseable* unit string and a *missing* one leave
/// the variable OUT of the map, so [`propagate_dim`] yields [`Dim::Unknown`] for
/// it and every expression containing it is skipped for dimensional checking.
///
/// Neither case is coerced to dimensionless. Doing so would both HIDE genuine
/// mismatches (when the variable is used consistently) and MANUFACTURE false
/// ones (when the real unit was, e.g., `m`) — and since a provable mismatch is
/// now a hard validation error, a fabricated dimension would fail the build on
/// a legitimate file. A model that declares no units at all (a common idiom in
/// the shared corpus) is therefore simply not dimension-checked, rather than
/// being treated as an all-dimensionless model in which `D(u)/dt = A*sin(w*t)`
/// looks inconsistent.
///
/// The two cases are NOT alike, and only one of them is a defect:
///
/// * A MISSING `units` field is a legitimate choice — the variable is simply
///   not dimension-checked.
/// * An UNPARSEABLE `units` field is a HARD ERROR (`unit_parse_error`,
///   esm-spec §4.8.4). A unit string either resolves against the shared
///   registry or it is a defect in the file: silently treating a typo
///   (`1/time`, `m/s2`) as "dimension unknown" disables every dimensional check
///   downstream, which is precisely the class of bug the checker exists to
///   catch. The failures are returned to the caller, which reports each one at
///   the JSON-Pointer of the offending DECLARATION.
///
/// Note this is a *different* case from a genuinely UNDETERMINABLE dimension (a
/// symbolic exponent), which stays a non-blocking warning: there the checker
/// cannot conclude, whereas here the author wrote something that is not a unit.
pub fn build_unit_env(
    variables: &indexmap::IndexMap<String, crate::ModelVariable>,
) -> (HashMap<String, Unit>, Vec<UnitParseFailure>) {
    let mut env = HashMap::new();
    let mut failures = Vec::new();
    // Many declarations share one unit string; each is parsed once.
    let mut parsed: HashMap<&str, Option<Unit>> = HashMap::new();
    for (name, var) in variables {
        let Some(declared) = &var.units else {
            // No declared units — dimension unknown, not dimensionless.
            continue;
        };
        let unit = parsed
            .entry(declared.as_str())
            .or_insert_with(|| parse_unit(declared).ok());
        match unit {
            Some(unit) => {
                env.insert(name.clone(), unit.clone());
            }
            None => failures.push(UnitParseFailure {
                name: name.clone(),
                units: declared.clone(),
            }),
        }
    }
    // `variables` is a HashMap, so iteration order is nondeterministic; sort so
    // a multi-failure file reports its errors in a stable order.
    failures.sort_by(|a, b| a.name.cmp(&b.name));
    (env, failures)
}

/// Where the dimensional checker looks up the unit of a name. A name it does
/// not know is UNDETERMINABLE, never dimensionless (esm-spec §4.8.4).
pub(crate) trait UnitEnv {
    /// The unit `name` carries, or `None` when it is undeterminable.
    fn unit_of(&self, name: &str) -> Option<&Unit>;

    /// Whether `name` has a known unit.
    fn contains(&self, name: &str) -> bool {
        self.unit_of(name).is_some()
    }
}

impl<S: std::hash::BuildHasher> UnitEnv for HashMap<String, Unit, S> {
    fn unit_of(&self, name: &str) -> Option<&Unit> {
        self.get(name)
    }
}

/// The units of every `<System>.<member>` scoped reference (esm-spec §4.9.2)
/// into the document's top-level models and reaction systems, at any depth,
/// each parsed once. Members that declare no units, or units that do not
/// parse (reported at their own declaration), are left out.
///
/// It depends only on the document, so validation builds it once and layers
/// each component's own names over it ([`ComponentUnitEnv`]).
pub(crate) fn document_scoped_units(esm_file: &crate::EsmFile) -> HashMap<String, Unit> {
    let mut parsed: HashMap<String, Option<Unit>> = HashMap::new();
    crate::structural::toplevel_reference_units(esm_file)
        .into_iter()
        .filter_map(|(name, units)| {
            parsed
                .entry(units)
                .or_insert_with_key(|units| parse_unit(units).ok())
                .clone()
                .map(|unit| (name, unit))
        })
        .collect()
}

/// The unit environment one component's expressions are checked against: its
/// own declarations, then the `<sub>.<member>` paths relative to its own
/// subsystems, then the document's top-level scoped references
/// ([`document_scoped_units`]), in that precedence.
///
/// A mounted component is, after resolution, identical to one written inline
/// (esm-spec §4.7), so a variable read across a `{ref}` mount edge carries its
/// declared units into the mounting component exactly as a local one does
/// (§4.8.4, EarthSciAST#414). A relative path occupies its name even when its
/// target declares no units (or units that do not parse): it then shadows a
/// top-level member of the same dotted name and is undeterminable, so a check
/// that reads it is skipped rather than judged against the wrong declaration.
pub(crate) struct ComponentUnitEnv<'d> {
    own: HashMap<String, Unit>,
    relative: HashMap<String, Option<Unit>>,
    document: &'d HashMap<String, Unit>,
}

impl<'d> ComponentUnitEnv<'d> {
    /// Layer a component's own unit environment `own` and its `subsystems`
    /// over the document's top-level scoped units.
    pub(crate) fn new(
        own: HashMap<String, Unit>,
        subsystems: Option<&indexmap::IndexMap<String, serde_json::Value>>,
        document: &'d HashMap<String, Unit>,
    ) -> Self {
        let mut parsed: HashMap<String, Option<Unit>> = HashMap::new();
        let relative = subsystems
            .map(crate::structural::subsystem_reference_units)
            .unwrap_or_default()
            .into_iter()
            .filter(|(name, _)| !own.contains_key(name))
            .map(|(name, units)| {
                let unit = units.and_then(|units| {
                    parsed
                        .entry(units)
                        .or_insert_with_key(|units| parse_unit(units).ok())
                        .clone()
                });
                (name, unit)
            })
            .collect();
        ComponentUnitEnv {
            own,
            relative,
            document,
        }
    }

    /// Every name with a known unit, own declarations winning, as one map.
    pub(crate) fn into_map(self) -> HashMap<String, Unit> {
        let mut env = self.own;
        for (name, unit) in self.document {
            if !env.contains_key(name) && !self.relative.contains_key(name) {
                env.insert(name.clone(), unit.clone());
            }
        }
        for (name, unit) in self.relative {
            if let Some(unit) = unit {
                env.entry(name).or_insert(unit);
            }
        }
        env
    }
}

impl UnitEnv for ComponentUnitEnv<'_> {
    fn unit_of(&self, name: &str) -> Option<&Unit> {
        if let Some(unit) = self.own.get(name) {
            return Some(unit);
        }
        match self.relative.get(name) {
            Some(unit) => unit.as_ref(),
            None => self.document.get(name),
        }
    }
}

/// The unit environment a model's equations are checked against, as one map:
/// its own declarations ([`build_unit_env`]) plus every scoped reference its
/// expressions can make, with the precedence of [`ComponentUnitEnv`]. The
/// parse failures returned are the model's own.
pub fn model_unit_env(
    esm_file: &crate::EsmFile,
    model: &crate::Model,
) -> (HashMap<String, Unit>, Vec<UnitParseFailure>) {
    let document = document_scoped_units(esm_file);
    let (own, failures) = build_unit_env(&model.variables);
    let env = ComponentUnitEnv::new(own, model.subsystems.as_ref(), &document).into_map();
    (env, failures)
}

/// A declared unit string that denotes no real unit.
///
/// Carries the NAME of the declaration alongside the offending string so the
/// caller can report `unit_parse_error` at the declaration's JSON-Pointer
/// rather than at the model.
#[derive(Debug, Clone, PartialEq, Eq)]
pub struct UnitParseFailure {
    /// The declaration's name (variable, species, or parameter).
    pub name: String,
    /// The unit string exactly as the author wrote it.
    pub units: String,
}

/// Validate that an equation's LHS and RHS have matching dimensions.
///
/// Returns `Ok(())` when both sides propagate to dimensionally-equal units,
/// and a [`UnitError`] with a descriptive message otherwise. Propagation
/// errors from either side are surfaced verbatim.
///
/// This is the single-outcome view. Validation uses
/// [`check_equation_dimensions`], which reports EVERY finding and separates
/// provable mismatches from undeterminable dimensions.
///
/// A time derivative `D(x)/Dt` whose `t` declares no units has an
/// undeterminable left-hand side (esm-spec §4.8.4), so the full comparison is
/// skipped and the weaker rule applies: if no time unit can reconcile the two
/// sides, the result is [`UnitError::DimensionMismatch`]; if one can, it is
/// [`UnitError::UnknownUnit`] ("skipped"), not `Ok(())`.
pub fn validate_equation_dimensions(
    eq: &Equation,
    env: &HashMap<String, Unit>,
) -> Result<(), UnitError> {
    if let Some(err) = check_equation_dimensions(eq, env)
        .iter()
        .find(|f| f.is_error())
    {
        return Err(UnitError::DimensionMismatch(err.message.clone()));
    }
    // Single-outcome contract: distinguish "checked, consistent" (`Ok`) from
    // "could not be checked" (`UnknownUnit`) — an equation with an indeterminate
    // side was SKIPPED, and reporting `Ok` would claim we verified it.
    let mut sink = Vec::new();
    let lhs = propagate_dim(&eq.lhs, env, &mut sink);
    let rhs = propagate_dim(&eq.rhs, env, &mut sink);
    if lhs.known().is_none() || rhs.known().is_none() {
        return Err(UnitError::UnknownUnit(
            "equation has an indeterminate side; skipped for dimensional checking".to_string(),
        ));
    }
    Ok(())
}

/// Dimension-check an equation, returning EVERY finding.
///
/// Both sides are always walked — a finding on the LHS never short-circuits the
/// RHS — so one undeterminable subexpression cannot hide a provable mismatch
/// elsewhere in the same equation. The LHS/RHS comparison itself is only made
/// when BOTH sides resolved: comparing against an unknown dimension could only
/// ever produce a false mismatch.
pub fn check_equation_dimensions(eq: &Equation, env: &HashMap<String, Unit>) -> Vec<UnitFinding> {
    equation_findings(eq, env)
}

/// [`check_equation_dimensions`] against any [`UnitEnv`].
pub(crate) fn equation_findings(eq: &Equation, env: &dyn UnitEnv) -> Vec<UnitFinding> {
    let mut findings = Vec::new();
    let lhs = propagate_dim(&eq.lhs, env, &mut findings);
    let rhs = propagate_dim(&eq.rhs, env, &mut findings);

    // `D(x)/dt` with an undeclared `t` is indeterminate, so the plain LHS-vs-RHS
    // comparison below cannot see it. Apply the weaker time-ratio rule instead.
    if let Some(state) = derivative_of_undeclared_time(&eq.lhs, env) {
        let state_dim = propagate_dim(state, env, &mut Vec::new());
        if let (Some(s), Some(r)) = (state_dim.known(), rhs.known())
            && let Some(message) = derivative_time_mismatch(s, r)
        {
            findings.push(UnitFinding::error(message));
        }
        return findings;
    }

    if let (Some(l), Some(r)) = (lhs.known(), rhs.known())
        && !l.same_unit(r)
    {
        findings.push(UnitFinding::error(format!(
            "Left-hand side has units of {} but right-hand side has units of {}",
            describe_scaled(l),
            describe_scaled(r)
        )));
    }
    findings
}

/// Dimension-check a standalone expression (an observed variable's defining
/// expression), returning EVERY finding.
///
/// When `declared` is `Some`, the expression's dimension is additionally
/// compared against the variable's declared units — the observed-variable
/// analogue of the equation LHS/RHS check.
pub fn check_expression_dimensions(
    expr: &Expr,
    declared: Option<&Unit>,
    env: &HashMap<String, Unit>,
) -> Vec<UnitFinding> {
    expression_findings(expr, declared, env)
}

/// [`check_expression_dimensions`] against any [`UnitEnv`].
pub(crate) fn expression_findings(
    expr: &Expr,
    declared: Option<&Unit>,
    env: &dyn UnitEnv,
) -> Vec<UnitFinding> {
    let mut findings = Vec::new();
    let actual = propagate_dim(expr, env, &mut findings);

    // If the expression is ALREADY internally inconsistent, its overall
    // dimension is whatever the first resolved operand happened to be — so
    // comparing it against the declaration would report the SAME defect a
    // second time at the same path. Report the internal mismatch only.
    if findings.iter().any(UnitFinding::is_error) {
        return findings;
    }

    if let (Some(declared), Some(actual)) = (declared, actual.known())
        && !declared.same_unit(actual)
    {
        findings.push(UnitFinding::error(format!(
            "Declared units of {} do not match the expression's units of {}",
            describe_scaled(declared),
            describe_scaled(actual)
        )));
    }
    findings
}

/// The dimensional analyser's verdict on one equation, in the form a report
/// prints: checked and consistent, provably inconsistent, or not checked.
#[derive(Debug, Clone, PartialEq)]
pub enum EquationDimensionVerdict {
    /// Both sides resolved and agree. Carries the propagated dimension, spelled
    /// over the canonical axes (`length*time^-1`).
    Consistent(String),
    /// A provable mismatch. Carries every error finding's message.
    Mismatch(Vec<String>),
    /// A side's dimension could not be determined, so the comparison was
    /// skipped (esm-spec §4.8.4). Carries the side and the operand to blame.
    NotChecked(String),
}

/// Classify one equation the way [`check_equation_dimensions`] judges it, and
/// say WHY when it could not be judged.
///
/// A `Consistent` verdict is only returned when the comparison actually ran, so
/// an equation skipped because of an indeterminate operand — a bare literal in a
/// product, a `const` node that declares no `units`, a variable with no declared
/// units — is never reported as fine.
pub fn equation_dimension_verdict(
    eq: &Equation,
    env: &HashMap<String, Unit>,
) -> EquationDimensionVerdict {
    let errors: Vec<String> = check_equation_dimensions(eq, env)
        .into_iter()
        .filter(UnitFinding::is_error)
        .map(|f| f.message)
        .collect();
    if !errors.is_empty() {
        return EquationDimensionVerdict::Mismatch(errors);
    }
    let rhs = propagate_dim(&eq.rhs, env, &mut Vec::new());
    // `D(x)` against an undeclared time is judged by the time-ratio rule, which
    // needs the state's dimension rather than the derivative's.
    let (lhs_label, lhs_expr) = match derivative_of_undeclared_time(&eq.lhs, env) {
        Some(state) => ("differentiated state", state),
        None => ("left-hand side", &eq.lhs),
    };
    let lhs = propagate_dim(lhs_expr, env, &mut Vec::new());
    match (lhs.known(), rhs.known()) {
        (Some(_), Some(r)) => EquationDimensionVerdict::Consistent(describe(r)),
        (None, _) => EquationDimensionVerdict::NotChecked(format!(
            "{lhs_label}: {}",
            indeterminate_operand(lhs_expr, env)
        )),
        (_, None) => EquationDimensionVerdict::NotChecked(format!(
            "right-hand side: {}",
            indeterminate_operand(&eq.rhs, env)
        )),
    }
}

/// Name the innermost operand that makes `expr`'s dimension indeterminate.
/// Called only on an expression whose dimension did not resolve.
fn indeterminate_operand(expr: &Expr, env: &dyn UnitEnv) -> String {
    match expr {
        Expr::Number(n) => format!("operand `{n}` has an indeterminate dimension"),
        Expr::Integer(n) => format!("operand `{n}` has an indeterminate dimension"),
        Expr::Variable(name) => format!("variable `{name}` has no declared units"),
        Expr::Operator(op) => {
            // A bare literal is the least informative operand to blame — it is
            // indeterminate everywhere — so a NAMED operand that is equally
            // indeterminate is reported in preference to it.
            let mut literal: Option<String> = None;
            for arg in &op.args {
                if propagate_dim(arg, env, &mut Vec::new()).known().is_some() {
                    continue;
                }
                if matches!(arg, Expr::Number(_) | Expr::Integer(_)) {
                    literal.get_or_insert_with(|| indeterminate_operand(arg, env));
                    continue;
                }
                return indeterminate_operand(arg, env);
            }
            literal.unwrap_or_else(|| format!("`{}` node has no determinable dimension", op.op))
        }
    }
}

/// Every declared `const` unit string in `expr` that does not resolve against
/// the registry (esm-spec §4.8.5 item 2), in walk order.
pub fn unresolvable_const_units(expr: &Expr) -> Vec<String> {
    fn walk(expr: &Expr, out: &mut Vec<String>) {
        if let Expr::Operator(node) = expr {
            if node.op == "const"
                && let Some(units) = &node.units
                && parse_unit(units).is_err()
            {
                out.push(units.clone());
            }
            node.for_each_child(&mut |child| walk(child, out));
        }
    }
    let mut out = Vec::new();
    walk(expr, &mut out);
    out
}

/// Reject declared `units` on any expression node in a document declaring
/// `esm` < 1.2.0 (esm-spec §4.8.5 item 6), naming the first offending node.
/// Runs on the raw JSON before schema validation, like the `solver` gate.
pub fn reject_const_units_pre_v12(
    view: &serde_json::Value,
) -> Result<(), crate::diagnostic::DiagnosticError> {
    reject_const_units_pre_v12_found(view, || {
        crate::json_visit::find_value_path(view, &mut |v| {
            v.as_object().is_some_and(|o| {
                crate::json_visit::small_get(o, "op").is_some()
                    && crate::json_visit::small_get(o, "units").is_some()
            })
        })
    })
}

/// [`reject_const_units_pre_v12`] with the search for the first expression node
/// declaring `units` supplied by the caller, and run only when the declared
/// version is below 1.2.0.
pub(crate) fn reject_const_units_pre_v12_found(
    view: &serde_json::Value,
    offending: impl FnOnce() -> Option<String>,
) -> Result<(), crate::diagnostic::DiagnosticError> {
    let Some(esm) = view.get("esm").and_then(|v| v.as_str()) else {
        return Ok(());
    };
    let Some((major, minor, _)) = crate::diagnostic::parse_semver(esm) else {
        return Ok(());
    };
    if (major, minor) >= (1, 2) {
        return Ok(());
    }
    match offending() {
        None => Ok(()),
        Some(path) => Err(crate::diagnostic::err(
            crate::diagnostic::codes::CONST_UNITS_VERSION_TOO_OLD,
            format!(
                "declared `units` on an expression node require esm >= 1.2.0; file declares \
                 {esm}. Offending path: {path}"
            ),
        )),
    }
}

/// Parse a unit string into a Unit struct
///
/// Supports common scientific unit notations:
/// - base units (`"m"`, `"s"`, `"kg"`, `"mol"`, `"K"`, `"Pa"`, …)
/// - products (`"kg*m"`, `"N*m"`)
/// - quotients (`"m/s"`, `"mol/L"`)
/// - parenthesised groups (`"J/(mol*K)"`)
/// - integer powers (`"m^2"`, `"s^-1"`)
/// - dimensionless (`""`, `"1"`, `"dimensionless"`)
pub fn parse_unit(unit_str: &str) -> Result<Unit, UnitError> {
    parse_normalized(&normalize_unit_str(unit_str), unit_str)
}

/// Rewrite the ordinary earth-science spellings of a unit into the ASCII
/// grammar, BEFORE parsing.
///
/// These are not exotic: `W/m²`, `cm³`, `J/(kg·K)`, `kg⋅m/s`, `μg/m^3`, `°C` and
/// `Pa*m**3` all appear in the shared corpus and all failed to parse, leaving
/// their variables dimensionless-unknown and silently unchecked.
///
/// * superscript digits `⁰¹²³⁴⁵⁶⁷⁸⁹` and superscript minus `⁻` → `^n`
/// * middot `·` (U+00B7) and dot-operator `⋅` (U+22C5) → `*`
/// * `µ` (U+00B5) / `μ` (U+03BC) → `u`  (so `μg` is the `u` + `g` prefix form)
/// * `°C` → `degC`, `Ω` → `Ohm`
/// * `**` → `^`  (the Python spelling of a power)
fn normalize_unit_str(s: &str) -> String {
    // `°C` and `**` are two-character sequences, so handle them first.
    let s = s.replace("°C", "degC").replace("**", "^");

    let mut out = String::with_capacity(s.len());
    let mut in_superscript = false;
    for c in s.chars() {
        let sup = match c {
            '⁰' => Some('0'),
            '¹' => Some('1'),
            '²' => Some('2'),
            '³' => Some('3'),
            '⁴' => Some('4'),
            '⁵' => Some('5'),
            '⁶' => Some('6'),
            '⁷' => Some('7'),
            '⁸' => Some('8'),
            '⁹' => Some('9'),
            '⁻' => Some('-'),
            _ => None,
        };
        match sup {
            // A run of superscripts is one exponent: emit a single `^` before it.
            Some(ascii) => {
                if !in_superscript {
                    out.push('^');
                    in_superscript = true;
                }
                out.push(ascii);
            }
            None => {
                in_superscript = false;
                match c {
                    '·' | '⋅' | '×' => out.push('*'),
                    'µ' | 'μ' => out.push('u'),
                    'Ω' => out.push_str("Ohm"),
                    _ => out.push(c),
                }
            }
        }
    }
    out
}

/// Parse an already-normalized unit string. `original` is carried only for
/// error messages, so a diagnostic quotes what the author actually wrote.
fn parse_normalized(s: &str, original: &str) -> Result<Unit, UnitError> {
    let s = s.trim();
    if s.is_empty() || s == "1" || s == "dimensionless" {
        return Ok(Unit::dimensionless());
    }

    // Strip a single layer of surrounding parentheses when they balance.
    if s.starts_with('(') && s.ends_with(')') && parens_balance(&s[1..s.len() - 1]) {
        return parse_normalized(&s[1..s.len() - 1], original);
    }

    let base_units = get_base_units();
    if let Some(unit) = base_units.get(s) {
        return Ok(unit.clone());
    }
    // An SI-PREFIXED symbol (`mg`, `um`, `nm`, `mL`, `umol`, `dm`). Tried only
    // after the exact table, so a name that merely LOOKS prefixed (`min`, `mol`,
    // `cd`, `day`) keeps its own meaning. This is what makes the `µ`/`μ` → `u`
    // normalization above mean anything: `μg/m^3` reaches here as `ug/m^3`.
    if let Some(unit) = parse_si_prefixed(s, base_units) {
        return Ok(unit);
    }

    // Product / quotient. `*` and `/` share one precedence level and bind
    // LEFT-associatively, so we split at the LAST top-level operator and
    // recurse on the left. Splitting on `/` as a separate, looser level (as
    // this parser used to) pulls every `*` factor appearing after the last `/`
    // into the DENOMINATOR: `J/mol*K` parsed as `J/(mol*K)`, silently negating
    // K's exponent, and `kg/m^3*s` came out as `kg/(m^3*s)`. Equal-precedence
    // left association is the reading of the Go reference parser and of
    // pint/Unitful.
    //
    // WHITESPACE between two terms is an implicit multiplication (`ppb^-1 s^-1`),
    // so it is a top-level operator too.
    if let Some((idx, len, sym)) = find_last_top_level_operator(s) {
        let (left, right) = (&s[..idx], &s[idx + len..]);
        if left.trim().is_empty() || right.trim().is_empty() {
            return Err(UnitError::ParseError(format!(
                "Dangling '{sym}' in unit \"{original}\""
            )));
        }
        let l = parse_normalized(left, original)?;
        let r = parse_normalized(right, original)?;
        return Ok(match sym {
            '/' => l.divide(&r),
            _ => l.multiply(&r),
        });
    }

    // Power: right-most '^'. The exponent is RATIONAL — `s^0.5` and `m^(1/2)`
    // are legitimate (an SDE noise intensity carries `1/s^0.5`).
    if let Some(idx) = s.rfind('^') {
        let (base_s, pow_s) = (&s[..idx], &s[idx + 1..]);
        let base_unit = parse_normalized(base_s, original)?;
        let power = parse_exponent(pow_s).ok_or_else(|| {
            UnitError::ParseError(format!(
                "Invalid exponent \"{pow_s}\" in unit \"{original}\""
            ))
        })?;
        return Ok(base_unit.power_rational(power));
    }

    // A numeric atom is admissible ONLY when its value is exactly 1 -- the
    // leading `1` of `1/s`. Any other number is a SCALING FACTOR, and a unit string
    // denotes a UNIT, not a quantity. This is not pedantry: `(m/s)^-1/3` -- an author
    // reaching for a rational exponent -- parses under this grammar as `((m/s)^-1)/3`,
    // and the five bindings gave it three different meanings (magnitude dropped,
    // magnitude retained, whole string rejected). Rejecting the scaling factor makes
    // `^(-1/3)` the only spelling of that unit, so the same string cannot silently
    // mean two things.
    //
    // Before this rule Rust RETAINED the magnitude here, so `"1000/s"` came out
    // as 1/s scaled by 1000 while Julia read the same string as a bare 1/s.
    if let Ok(v) = s.parse::<f64>() {
        if v != 1.0 {
            return Err(UnitError::ParseError(format!(
                "unit \"{original}\": a number other than 1 is a scaling factor, not a unit (esm-spec 4.8.2); a rational exponent is spelled ^(p/q)"
            )));
        }
        return Ok(Unit {
            dimensions: HashMap::new(),
            scale: 1.0,
            exact: ExactScale::one(),
        });
    }

    Err(UnitError::UnknownUnit(original.to_string()))
}

/// Decompose an SI-prefixed symbol into prefix × base (`mg` → milli × gram).
///
/// Only the symbols in [`PREFIXABLE`] take a prefix, which keeps the
/// decomposition from inventing units out of arbitrary identifiers — `not_a_unit`
/// must still fail to parse, and a count like `units` must not be read as
/// micro-something. Prefixes are tried LONGEST-first so `da` (deca) is not
/// mistaken for `d` (deci).
fn parse_si_prefixed(s: &str, base_units: &HashMap<String, Unit>) -> Option<Unit> {
    /// (symbol, factor, power of ten), longest symbol first.
    const PREFIXES: [(&str, f64, i32); 20] = [
        ("da", 1e1, 1),
        ("Y", 1e24, 24),
        ("Z", 1e21, 21),
        ("E", 1e18, 18),
        ("P", 1e15, 15),
        ("T", 1e12, 12),
        ("G", 1e9, 9),
        ("M", 1e6, 6),
        ("k", 1e3, 3),
        ("h", 1e2, 2),
        ("d", 1e-1, -1),
        ("c", 1e-2, -2),
        ("m", 1e-3, -3),
        ("u", 1e-6, -6),
        ("n", 1e-9, -9),
        ("p", 1e-12, -12),
        ("f", 1e-15, -15),
        ("a", 1e-18, -18),
        ("z", 1e-21, -21),
        ("y", 1e-24, -24),
    ];
    /// Symbols that admit an SI prefix.
    const PREFIXABLE: [&str; 13] = [
        "m", "g", "s", "mol", "K", "A", "cd", "L", "N", "Pa", "J", "W", "rad",
    ];

    for (prefix, factor, power) in PREFIXES {
        let Some(base) = s.strip_prefix(prefix) else {
            continue;
        };
        if !PREFIXABLE.contains(&base) {
            continue;
        }
        let mut unit = base_units.get(base)?.clone();
        unit.scale *= factor;
        unit.exact = unit.exact.multiply(&ExactScale::power_of_ten(power));
        return Some(unit);
    }
    None
}

/// Parse a unit exponent: an integer (`2`, `-1`), a decimal (`0.5`, `-1.5`), or
/// a parenthesised fraction (`(1/2)`).
fn parse_exponent(s: &str) -> Option<Rational> {
    let s = s.trim();
    let inner = s
        .strip_prefix('(')
        .and_then(|rest| rest.strip_suffix(')'))
        .unwrap_or(s)
        .trim();
    if let Some((num, den)) = inner.split_once('/') {
        return Some(Rational::new(
            num.trim().parse().ok()?,
            den.trim().parse().ok()?,
        ));
    }
    if let Ok(n) = inner.parse::<i32>() {
        return Some(Rational::int(n));
    }
    Rational::from_f64(inner.parse::<f64>().ok()?)
}

/// Find the *last* top-level (outside any parenthesised group) `*` or `/`,
/// returning its byte index and which operator it is.
///
/// `parse_unit` splits here and recurses on the *left* substring, which makes
/// the two operators share a precedence level and associate left-to-right:
/// `a/b/c` is `(a/b)/c`, `a*b/c` is `(a*b)/c`, and `a/b*c` is `(a/b)*c`. This
/// gives the conventional reading of `W/m^2/K` as `W/(m^2*K)` while keeping
/// `J/mol*K` as `(J/mol)*K`, and matches the sibling parsers (Go's recursive
/// descent, Python/pint, Julia's Unitful `uparse`).
/// Returns `(byte index, byte length, operator)`. The length matters because an
/// implicit multiplication is a RUN of whitespace, which the caller must skip
/// whole.
fn find_last_top_level_operator(s: &str) -> Option<(usize, usize, char)> {
    let bytes = s.as_bytes();
    let mut depth = 0i32;
    let mut last = None;
    let mut i = 0usize;
    while i < bytes.len() {
        let c = bytes[i] as char;
        match c {
            '(' => depth += 1,
            ')' => depth -= 1,
            '*' | '/' if depth == 0 => last = Some((i, 1, c)),
            ' ' | '\t' if depth == 0 => {
                // A run of whitespace BETWEEN two terms is an implicit `*`.
                // Whitespace merely padding an explicit operator (`m / s`) is
                // not: the operator arm above already claimed that position, and
                // a run adjacent to one must not be treated as a second operator.
                let start = i;
                while i < bytes.len() && matches!(bytes[i] as char, ' ' | '\t') {
                    i += 1;
                }
                let before = s[..start].trim_end().chars().last();
                let after = s[i..].chars().next();
                let padding = matches!(before, Some('*' | '/' | '^') | None)
                    || matches!(after, Some('*' | '/' | '^') | None);
                if !padding {
                    last = Some((start, i - start, '*'));
                }
                continue;
            }
            _ => {}
        }
        i += 1;
    }
    last
}

fn parens_balance(s: &str) -> bool {
    let mut depth = 0i32;
    for c in s.chars() {
        match c {
            '(' => depth += 1,
            ')' => {
                depth -= 1;
                if depth < 0 {
                    return false;
                }
            }
            _ => {}
        }
    }
    depth == 0
}

/// Get commonly used base units
/// The cached base-unit registry. `parse_unit` recurses (parentheses, `/`,
/// `*` splits) and runs inside per-equation validation loops, so rebuilding
/// the ~50-entry map on every call was pure waste; the table is immutable
/// and deterministic, so one process-wide build is safe.
fn get_base_units() -> &'static HashMap<String, Unit> {
    static BASE_UNITS: std::sync::LazyLock<HashMap<String, Unit>> =
        std::sync::LazyLock::new(build_base_units);
    &BASE_UNITS
}

fn build_base_units() -> HashMap<String, Unit> {
    let mut units = HashMap::new();

    // Length units
    units.insert("m".to_string(), Unit::base(Dimension::Length, 1, 1.0));
    units.insert(
        "cm".to_string(),
        Unit::base(Dimension::Length, 1, 0.01).with_exact(ExactScale::power_of_ten(-2)),
    );
    units.insert(
        "km".to_string(),
        Unit::base(Dimension::Length, 1, 1000.0).with_exact(ExactScale::power_of_ten(3)),
    );
    units.insert(
        "mm".to_string(),
        Unit::base(Dimension::Length, 1, 0.001).with_exact(ExactScale::power_of_ten(-3)),
    );
    units.insert("meter".to_string(), Unit::base(Dimension::Length, 1, 1.0));
    units.insert("meters".to_string(), Unit::base(Dimension::Length, 1, 1.0));
    // The international foot, exact by definition since 1959: 1 ft = 0.3048 m.
    // Emission inventories are written in it — the EPA FF10 point-source format
    // stores STKHGT and STKDIAM in feet — and a format for air-quality models
    // that cannot spell the unit its own input files use forces every such
    // column to be declared in a unit it is not stored in.
    units.insert(
        "ft".to_string(),
        Unit::base(Dimension::Length, 1, 0.3048).with_exact(ExactScale::decimal("0.3048")),
    );
    // The international mile, exact by definition since the same 1959
    // agreement: 1 mi = 5280 ft = 1609.344 m. It is the unit the US onroad
    // transportation inventory is written in end to end — EPA MOVES stores
    // `link.linkLength` in miles, `link.linkAvgSpeed` and
    // `driveSchedule.averageSpeed` in `mi/h`, and its whole activity model is
    // built on vehicle-MILES travelled — so a table that has `ft` and not `mi`
    // can spell a stack height and not a road.
    //
    // `mi/h` composes from this entry and `h`; `mph` is deliberately NOT a
    // separate name, because it is a fused spelling of a compound the grammar
    // already builds, and every fused alias is a second way to say one thing.
    // Nor are `mile`/`miles` here: `foot`/`feet` are pinned as REJECTS by
    // tests/conformance/unit_registry, so the imperial family is symbol-only
    // by contract and this follows it.
    //
    // The foot is READ BACK OUT of the table rather than retyped as 1609.344,
    // so this entry cannot drift away from the one that defines it.
    units.insert(
        "mi".to_string(),
        Unit::base(Dimension::Length, 1, 5280.0 * units["ft"].scale)
            .with_exact(ExactScale::integer(5280).multiply(&units["ft"].exact)),
    );

    // Time units
    units.insert("s".to_string(), Unit::base(Dimension::Time, 1, 1.0));
    units.insert(
        "min".to_string(),
        Unit::base(Dimension::Time, 1, 60.0).with_exact(ExactScale::integer(60)),
    );
    units.insert(
        "h".to_string(),
        Unit::base(Dimension::Time, 1, 3600.0).with_exact(ExactScale::integer(3600)),
    );
    units.insert(
        "hr".to_string(),
        Unit::base(Dimension::Time, 1, 3600.0).with_exact(ExactScale::integer(3600)),
    );
    units.insert(
        "hour".to_string(),
        Unit::base(Dimension::Time, 1, 3600.0).with_exact(ExactScale::integer(3600)),
    );
    // The canonical spelling of the day is `day`. A bare `d` is deliberately NOT
    // a unit (esm-spec §4.8.1): a one-letter symbol reads as the deci- prefix or
    // as a differential, so admitting it would make `d` ambiguous at every site.
    units.insert(
        "day".to_string(),
        Unit::base(Dimension::Time, 1, 86400.0).with_exact(ExactScale::integer(86400)),
    );

    // Mass units
    units.insert("kg".to_string(), Unit::base(Dimension::Mass, 1, 1.0));
    units.insert(
        "g".to_string(),
        Unit::base(Dimension::Mass, 1, 0.001).with_exact(ExactScale::power_of_ten(-3)),
    );
    // The international avoirdupois pound, exact by definition since 1959:
    // 1 lb = 0.45359237 kg. US emission rates are tabulated in it — MOVES's
    // NONROAD brake-specific fuel consumption is `lb/(hp*h)` and its gasoline
    // density constant CMFGAS is 6.237 lb/gal — and it is the unit `short_ton`
    // below is DEFINED in, so having the ton and not the pound left the
    // derived unit spellable and the base one not.
    units.insert(
        "lb".to_string(),
        Unit::base(Dimension::Mass, 1, 0.45359237).with_exact(ExactScale::decimal("0.45359237")),
    );
    // The two tons, both spelled UNAMBIGUOUSLY and neither spelled `ton`. A bare
    // `ton` is three different masses (short 907.18474 kg, metric 1000 kg, long
    // 1016.0469088 kg), and a table whose job is to make a declared unit mean
    // ONE thing cannot hold a name that means three. `t` is excluded for the
    // same reason `d` is. `short_ton` is exactly 2000 international pounds —
    // what a US emissions inventory means by "tons", and exactly InMAP's
    // 907184740000 ug/short-ton emission-conversion constant.
    // The pound is READ BACK OUT of the table rather than retyped as
    // 907.18474, for the same reason `mi` reads the foot back out.
    units.insert(
        "short_ton".to_string(),
        Unit::base(Dimension::Mass, 1, 2000.0 * units["lb"].scale)
            .with_exact(ExactScale::integer(2000).multiply(&units["lb"].exact)),
    );
    units.insert(
        "tonne".to_string(),
        Unit::base(Dimension::Mass, 1, 1000.0).with_exact(ExactScale::power_of_ten(3)),
    );

    // Amount units
    units.insert("mol".to_string(), Unit::base(Dimension::Amount, 1, 1.0));

    // Temperature units. `degC` (and its `°C` spelling, normalized to `degC`
    // before parsing) shares the TEMPERATURE dimension with kelvin: the two
    // differ by an affine OFFSET, which a multiplicative dimension algebra
    // cannot express, and dimensional analysis only cares about the axis.
    units.insert("K".to_string(), Unit::base(Dimension::Temperature, 1, 1.0));
    units.insert(
        "degC".to_string(),
        Unit::base(Dimension::Temperature, 1, 1.0),
    );
    units.insert(
        "Celsius".to_string(),
        Unit::base(Dimension::Temperature, 1, 1.0),
    );

    // Current / luminous intensity
    units.insert("A".to_string(), Unit::base(Dimension::Current, 1, 1.0));
    units.insert("cd".to_string(), Unit::base(Dimension::Luminosity, 1, 1.0));

    // Plane angle — a real axis (the eighth), not a synonym for dimensionless.
    // The trig functions accept it; see `propagate_transcendental_dim`.
    units.insert("rad".to_string(), Unit::base(Dimension::Angle, 1, 1.0));
    // Degrees of arc — the same axis, scaled. (`deg` is the short spelling; the
    // TEMPERATURE `degC`/`degF` are separate entries and are matched first, so
    // there is no collision.)
    let degree = Unit::base(Dimension::Angle, 1, std::f64::consts::PI / 180.0)
        .with_exact(ExactScale::pi().divide(&ExactScale::integer(180)));
    units.insert("degrees".to_string(), degree.clone());
    units.insert("degree".to_string(), degree.clone());
    units.insert("deg".to_string(), degree);

    // Solid angle — the steradian, `sr = rad²` (esm-spec §4.8.1). `rad` is the
    // Angle axis, so a solid angle is that axis SQUARED, not a ninth axis, just
    // as an area is `m²` and not a new axis. Spherical-mesh cell areas (a patch
    // of the unit sphere: MPAS Voronoi cells, Girard triangles) declare
    // `units: "sr"`, and this makes those declarations resolve.
    units.insert("sr".to_string(), Unit::base(Dimension::Angle, 2, 1.0));

    // Volume (L = dm³ = 10⁻³ m³)
    units.insert(
        "L".to_string(),
        Unit::base(Dimension::Length, 3, 0.001).with_exact(ExactScale::power_of_ten(-3)),
    );
    // The US liquid gallon, exact by definition: 231 in³ = 3.785 411 784 L
    // (NIST SP 811 App. B). It is spelled `gal` and not `gallon` for the same
    // reason `ft` has no `foot`, and it is the US liquid gallon and not the
    // imperial one — those differ by 20%, which is exactly the kind of silent
    // scale error tests/conformance/unit_registry exists to pin. US fuel data
    // is per gallon: MOVES stores `fueltype.fuelDensity` in g/gal, its
    // refuelling spill rate in g/gal, and its dioxin and metal emission rates
    // in g/gal.
    units.insert(
        "gal".to_string(),
        Unit::base(Dimension::Length, 3, 0.003_785_411_784)
            .with_exact(ExactScale::decimal("0.003785411784")),
    );

    // Frequency: Hz = s⁻¹.
    units.insert("Hz".to_string(), Unit::base(Dimension::Time, -1, 1.0));

    // Fahrenheit shares the TEMPERATURE axis with kelvin, like `degC`, at a
    // scale of exactly 5/9 K (esm-spec §4.8.1). Its affine offset is not
    // modelled, and not checked.
    units.insert(
        "degF".to_string(),
        Unit::base(Dimension::Temperature, 1, 5.0 / 9.0).with_exact(ExactScale::ratio(5, 9)),
    );

    // Year — the JULIAN year, 365.25 days = 31_557_600 s exactly. That is the
    // astronomical/climate convention and it is what Julia, Python, Go and
    // TypeScript all carry; Rust alone had 3.1536e7 (a 365-day year), which is
    // 0.0685% short. A scale error is INVISIBLE to dimensional analysis — the
    // dimension is [time] either way — so nothing in the corpus could see it
    // until `tests/conformance/unit_registry` started pinning scales, which is
    // exactly what that fixture is for.
    units.insert(
        "year".to_string(),
        Unit::base(Dimension::Time, 1, 31_557_600.0).with_exact(ExactScale::integer(31_557_600)),
    );
    units.insert(
        "yr".to_string(),
        Unit::base(Dimension::Time, 1, 31_557_600.0).with_exact(ExactScale::integer(31_557_600)),
    );

    // Derived units
    // Velocity: m/s
    units.insert(
        "m/s".to_string(),
        Unit::base(Dimension::Length, 1, 1.0).divide(&Unit::base(Dimension::Time, 1, 1.0)),
    );

    // Concentration: mol/L (mol/m^3)
    let liter = Unit::base(Dimension::Length, 3, 0.001).with_exact(ExactScale::power_of_ten(-3));
    units.insert(
        "mol/L".to_string(),
        Unit::base(Dimension::Amount, 1, 1.0).divide(&liter),
    );
    units.insert(
        "M".to_string(),
        Unit::base(Dimension::Amount, 1, 1.0).divide(&liter),
    );

    // Force: kg*m/s^2 (Newton)
    let newton = Unit::base(Dimension::Mass, 1, 1.0)
        .multiply(&Unit::base(Dimension::Length, 1, 1.0))
        .divide(&Unit::base(Dimension::Time, 2, 1.0));
    units.insert("N".to_string(), newton.clone());

    // Pressure: kg/(m*s^2) (Pascal)
    let pascal = Unit::base(Dimension::Mass, 1, 1.0)
        .divide(&Unit::base(Dimension::Length, 1, 1.0))
        .divide(&Unit::base(Dimension::Time, 2, 1.0));
    units.insert("Pa".to_string(), pascal);

    // Energy: kg*m^2/s^2 (Joule)
    let joule = Unit::base(Dimension::Mass, 1, 1.0)
        .multiply(&Unit::base(Dimension::Length, 2, 1.0))
        .divide(&Unit::base(Dimension::Time, 2, 1.0));
    units.insert("J".to_string(), joule.clone());
    let mut kj = joule.clone();
    kj.scale *= 1000.0;
    kj.exact = kj.exact.multiply(&ExactScale::power_of_ten(3));
    units.insert("kJ".to_string(), kj);
    let mut cal = joule.clone();
    cal.scale *= 4.184;
    cal.exact = cal.exact.multiply(&ExactScale::decimal("4.184"));
    units.insert("cal".to_string(), cal);
    let mut kcal = joule.clone();
    kcal.scale *= 4184.0;
    kcal.exact = kcal.exact.multiply(&ExactScale::integer(4184));
    units.insert("kcal".to_string(), kcal);

    // Power: kg*m^2/s^3 (Watt)
    let watt = Unit::base(Dimension::Mass, 1, 1.0)
        .multiply(&Unit::base(Dimension::Length, 2, 1.0))
        .divide(&Unit::base(Dimension::Time, 3, 1.0));
    units.insert("W".to_string(), watt.clone());

    // Mechanical (imperial) horsepower — 550 ft·lbf/s, which is
    // 745.699 871 582 270 2 W: NIST SP 811 App. B gives 7.456 999 E+02 W. The
    // foot and the pound are READ BACK OUT of the table rather than retyped as
    // literals here, so this entry cannot drift away from the two entries that
    // define it. NOT the metric horsepower (PS, 735.49875 W), which is a
    // different unit by 1.4%.
    //
    // Engine ratings are the axis MOVES's NONROAD model bins on:
    // `nrsourceusetype.hpAvg` is horsepower and every `nremissionrate` row is
    // `g/(hp*h)`, so without this entry the rate tables that ARE the model have
    // no honest declaration.
    let foot = units["ft"].scale;
    let pound = units["lb"].scale;
    let mut horsepower = watt.clone();
    horsepower.scale *= 550.0 * foot * pound * 9.80665;
    horsepower.exact = horsepower.exact.multiply(
        &ExactScale::integer(550)
            .multiply(&units["ft"].exact)
            .multiply(&units["lb"].exact)
            .multiply(&ExactScale::decimal("9.80665")),
    );
    units.insert("hp".to_string(), horsepower);

    // Electromagnetic family, all derived from the ampere.
    //
    // `C` is the COULOMB (A·s), NOT Celsius — degrees Celsius is spelled
    // `degC`, and the `°C` form normalizes to it before parsing. A unit string
    // carries DIMENSIONS ONLY, so a bare `C` can never mean a chemical species
    // tag either.
    let ampere = Unit::base(Dimension::Current, 1, 1.0);
    let coulomb = ampere.multiply(&Unit::base(Dimension::Time, 1, 1.0));
    units.insert("C".to_string(), coulomb.clone());

    // Volt = W/A = kg*m^2/(s^3*A)
    let volt = watt.divide(&ampere);
    units.insert("V".to_string(), volt.clone());

    // Ohm = V/A. Reached from the `Ω` spelling too, which normalizes to `Ohm`.
    units.insert("Ohm".to_string(), volt.divide(&ampere));

    // Farad = C/V (so `F/m`, the permittivity of the corpus, is F ÷ metre).
    units.insert("F".to_string(), coulomb.divide(&volt));

    // Tesla = kg/(s^2*A). The weber appears in the corpus as the composite
    // `V*s`, which the parser builds from `V` and `s`, so it needs no entry.
    units.insert(
        "T".to_string(),
        Unit::base(Dimension::Mass, 1, 1.0)
            .divide(&Unit::base(Dimension::Time, 2, 1.0))
            .divide(&ampere),
    );

    // Non-SI energies, as multiples of the joule.
    let joule_scaled = |scale: f64, exact: ExactScale| {
        let mut u = joule.clone();
        u.scale *= scale;
        u.exact = u.exact.multiply(&exact);
        u
    };
    units.insert(
        "erg".to_string(),
        joule_scaled(1e-7, ExactScale::power_of_ten(-7)),
    );
    units.insert(
        "BTU".to_string(),
        joule_scaled(1_055.055_852_62, ExactScale::decimal("1055.05585262")),
    );
    // The watt-hour, exactly 3600 J, and the kilowatt-hour it is the base of.
    // esm-spec §4.8.1's Energy row lists BOTH; Julia, Python, Go and TypeScript
    // carried both while this table had only `kWh`, so `units: "Wh"` resolved in
    // four bindings and was a hard error in this one. `W*h` composes either way,
    // which is why the gap was silent.
    units.insert(
        "Wh".to_string(),
        joule_scaled(3600.0, ExactScale::integer(3600)),
    );
    units.insert(
        "kWh".to_string(),
        joule_scaled(3.6e6, ExactScale::integer(3_600_000)),
    );

    // ESM-specific units standard (docs/units-standard.md).
    // Mole-fraction family: dimensionless with scale factors.
    // ppmv/ppbv/pptv are volume-mixing-ratio aliases of ppm/ppb/ppt under
    // the ideal-gas approximation — identical dimension and scale.
    let dimensionless_scaled = |scale: f64, exact: ExactScale| Unit {
        dimensions: HashMap::new(),
        scale,
        exact,
    };
    units.insert(
        "ppm".to_string(),
        dimensionless_scaled(1e-6, ExactScale::power_of_ten(-6)),
    );
    units.insert(
        "ppmv".to_string(),
        dimensionless_scaled(1e-6, ExactScale::power_of_ten(-6)),
    );
    units.insert(
        "ppb".to_string(),
        dimensionless_scaled(1e-9, ExactScale::power_of_ten(-9)),
    );
    units.insert(
        "ppbv".to_string(),
        dimensionless_scaled(1e-9, ExactScale::power_of_ten(-9)),
    );
    units.insert(
        "ppt".to_string(),
        dimensionless_scaled(1e-12, ExactScale::power_of_ten(-12)),
    );
    units.insert(
        "pptv".to_string(),
        dimensionless_scaled(1e-12, ExactScale::power_of_ten(-12)),
    );
    // Percent — dimensionless with a scale, like the mole-fraction family.
    units.insert(
        "%".to_string(),
        dimensionless_scaled(0.01, ExactScale::power_of_ten(-2)),
    );
    units.insert(
        "percent".to_string(),
        dimensionless_scaled(0.01, ExactScale::power_of_ten(-2)),
    );

    // COUNTS carry no dimension axis: they count entities, they do not measure
    // one. Composites therefore fall out of the parser with the right dimension
    // on their own — `molec/cm^3` is `[length]^-3`, `individuals/km^2` is
    // `[length]^-2`.
    for count in [
        "molec",
        "molecule",
        "count",
        "individuals",
        "vehicles",
        "units",
    ] {
        units.insert(count.to_string(), Unit::dimensionless());
    }

    // Practical salinity — dimensionless by definition (PSS-78).
    units.insert("psu".to_string(), Unit::dimensionless());

    // Pressure: atmosphere, and the micro-atmosphere used for seawater pCO2.
    units.insert(
        "atm".to_string(),
        pascal_scaled(&units, 101_325.0, ExactScale::integer(101_325)),
    );
    units.insert(
        "uatm".to_string(),
        pascal_scaled(&units, 0.101_325, ExactScale::decimal("0.101325")),
    );
    // Non-SI pressures, as multiples of the pascal.
    units.insert(
        "bar".to_string(),
        pascal_scaled(&units, 1e5, ExactScale::power_of_ten(5)),
    );
    // Torr is EXACTLY atm/760 (esm-spec §4.8.1) — not a rounded decimal, and not
    // mmHg, which is the conventional 133.322387415 Pa.
    units.insert(
        "Torr".to_string(),
        pascal_scaled(&units, 101_325.0 / 760.0, ExactScale::ratio(101_325, 760)),
    );
    units.insert(
        "mmHg".to_string(),
        pascal_scaled(
            &units,
            133.322_387_415,
            ExactScale::decimal("133.322387415"),
        ),
    );
    // Inch of mercury — exactly 25.4 mmHg, the conventional value (NIST SP 811).
    // US barometric datasets store pressure in inHg; without this entry such a
    // column has no honest declaration, because a unit string carries no numeric
    // scale factor, so `25.4 mmHg` cannot be spelled either.
    units.insert(
        "inHg".to_string(),
        pascal_scaled(
            &units,
            3_386.388_640_341,
            ExactScale::decimal("3386.388640341"),
        ),
    );
    // psi is lbf/in² = pound × standard gravity ÷ inch², every factor exact
    // (esm-spec §4.8.1), folded in that order for the numeric scale.
    let psi_scale = units["lb"].scale * 9.80665 / (0.0254 * 0.0254);
    let psi_exact = units["lb"]
        .exact
        .multiply(&ExactScale::decimal("9.80665"))
        .divide(&ExactScale::decimal("0.0254").power(Rational::int(2)));
    units.insert(
        "psi".to_string(),
        pascal_scaled(&units, psi_scale, psi_exact),
    );

    // Dobson unit: areal number density of ozone molecules.
    // Since `molec` is dimensionless, Dobson resolves to `[length]^-2` with
    // scale 2.6867e20 molec/m^2.
    units.insert(
        "Dobson".to_string(),
        Unit::base(Dimension::Length, -2, 2.6867e20).with_exact(ExactScale::decimal("2.6867e20")),
    );
    units.insert(
        "DU".to_string(),
        Unit::base(Dimension::Length, -2, 2.6867e20).with_exact(ExactScale::decimal("2.6867e20")),
    );

    units
}

/// A pressure unit expressed as a multiple of the pascal already in `units`.
fn pascal_scaled(units: &HashMap<String, Unit>, scale: f64, exact: ExactScale) -> Unit {
    let mut u = units["Pa"].clone();
    u.scale *= scale;
    u.exact = u.exact.multiply(&exact);
    u
}

/// Check dimensional consistency of an equation
///
/// # Arguments
///
/// * `lhs_unit` - Units of the left-hand side
/// * `rhs_unit` - Units of the right-hand side
///
/// # Returns
///
/// * `Result<(), UnitError>` - Ok if consistent, error otherwise
pub fn check_dimensional_consistency(lhs_unit: &Unit, rhs_unit: &Unit) -> Result<(), UnitError> {
    if lhs_unit.is_compatible(rhs_unit) {
        Ok(())
    } else {
        Err(UnitError::DimensionMismatch(
            "Left and right sides have incompatible dimensions".to_string(),
        ))
    }
}

/// Convert between compatible units
///
/// # Arguments
///
/// * `value` - Value to convert
/// * `from_unit` - Source unit
/// * `to_unit` - Target unit
///
/// # Returns
///
/// * `Result<f64, UnitError>` - Converted value or error
pub fn convert_units(value: f64, from_unit: &Unit, to_unit: &Unit) -> Result<f64, UnitError> {
    if !from_unit.is_compatible(to_unit) {
        return Err(UnitError::DimensionMismatch(
            "Units have incompatible dimensions".to_string(),
        ));
    }

    Ok(value * from_unit.scale / to_unit.scale)
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::types::ExpressionNode;

    fn op(name: &str, args: Vec<Expr>) -> Expr {
        Expr::operator(ExpressionNode {
            op: name.to_string(),
            args,
            ..ExpressionNode::default()
        })
    }

    fn op_with_wrt(name: &str, args: Vec<Expr>, wrt: &str) -> Expr {
        Expr::operator(ExpressionNode {
            op: name.to_string(),
            args,
            wrt: Some(wrt.to_string()),
            ..ExpressionNode::default()
        })
    }

    fn env_of(pairs: &[(&str, &str)]) -> HashMap<String, Unit> {
        pairs
            .iter()
            .map(|(k, v)| ((*k).to_string(), parse_unit(v).unwrap()))
            .collect()
    }

    #[test]
    fn test_parse_dimensionless() {
        assert!(parse_unit("").unwrap().is_dimensionless());
        assert!(parse_unit("1").unwrap().is_dimensionless());
        assert!(parse_unit("dimensionless").unwrap().is_dimensionless());
    }

    /// Exponents are RATIONAL. `1/s^0.5` is the intensity of an SDE noise term
    /// and appears in the shared corpus (`tests/fixtures/sde/*.esm`); an
    /// integer-only grammar could not parse it at all.
    #[test]
    fn test_rational_exponents() {
        let half = Rational::new(-1, 2);
        assert_eq!(
            parse_unit("1/s^0.5")
                .unwrap()
                .dimensions
                .get(&Dimension::Time),
            Some(&half)
        );
        // The fraction and decimal spellings agree.
        assert_eq!(
            parse_unit("s^(-1/2)")
                .unwrap()
                .dimensions
                .get(&Dimension::Time),
            Some(&half)
        );
        // `**` is the Python spelling of a power.
        assert_eq!(
            parse_unit("m**3")
                .unwrap()
                .dimensions
                .get(&Dimension::Length),
            Some(&Rational::int(3))
        );
        // Reduced to lowest terms, so 2/4 IS 1/2.
        assert_eq!(Rational::new(2, 4), Rational::new(1, 2));
        assert_eq!(Rational::from_f64(0.5), Some(Rational::new(1, 2)));
    }

    /// Ordinary earth-science spellings must parse: superscripts, the middot and
    /// dot-operator, micro, degree-Celsius and ohm.
    #[test]
    fn test_unicode_normalization() {
        // W/m² — superscript exponent.
        let u = parse_unit("W/m²").unwrap();
        assert_eq!(u.dimensions.get(&Dimension::Mass), Some(&Rational::int(1)));
        assert_eq!(u.dimensions.get(&Dimension::Time), Some(&Rational::int(-3)));
        // Length cancels: W carries m^2, divided by m^2.
        assert_eq!(u.dimensions.get(&Dimension::Length), None);

        assert_eq!(
            parse_unit("cm³")
                .unwrap()
                .dimensions
                .get(&Dimension::Length),
            Some(&Rational::int(3))
        );
        // Superscript minus.
        assert_eq!(
            parse_unit("m⁻³")
                .unwrap()
                .dimensions
                .get(&Dimension::Length),
            Some(&Rational::int(-3))
        );
        // Middot (U+00B7) and dot-operator (U+22C5) are multiplication.
        assert_eq!(
            parse_unit("J/(kg·K)").unwrap(),
            parse_unit("J/(kg*K)").unwrap()
        );
        assert_eq!(parse_unit("kg⋅m/s").unwrap(), parse_unit("kg*m/s").unwrap());
        // Micro, in both its spellings.
        assert_eq!(parse_unit("μg/m^3").unwrap(), parse_unit("ug/m^3").unwrap());
        assert_eq!(parse_unit("µg").unwrap(), parse_unit("ug").unwrap());
        // °C shares the temperature axis with K (they differ by an affine
        // offset, which dimensional analysis does not model).
        assert!(
            parse_unit("°C")
                .unwrap()
                .is_compatible(&parse_unit("K").unwrap())
        );
    }

    /// Whitespace between two terms is an implicit multiplication.
    #[test]
    fn test_implicit_multiplication() {
        assert_eq!(
            parse_unit("ppb^-1 s^-1").unwrap(),
            parse_unit("ppb^-1*s^-1").unwrap()
        );
        // But whitespace merely PADDING an explicit operator is not a second
        // operator.
        assert_eq!(parse_unit("m / s").unwrap(), parse_unit("m/s").unwrap());
        assert_eq!(parse_unit("kg * m").unwrap(), parse_unit("kg*m").unwrap());
    }

    /// Inches of mercury — exactly 25.4 mmHg BY DEFINITION, so the two must agree
    /// to the last bit rather than merely within a tolerance. `inHg` is how US
    /// barometric data is stored, and before this entry the registry had `mmHg`
    /// and no route to it: unit strings carry no numeric scale factor, so
    /// `25.4 mmHg` was not spellable either, and such a column had to declare a
    /// WRONG unit or none at all.
    #[test]
    fn test_inches_of_mercury() {
        let inhg = parse_unit("inHg").unwrap();
        let mmhg = parse_unit("mmHg").unwrap();
        assert!(
            inhg.is_compatible(&parse_unit("Pa").unwrap()),
            "inHg must carry the pressure dimension"
        );
        assert_eq!(
            inhg.scale,
            mmhg.scale * 25.4,
            "inHg must be EXACTLY 25.4 mmHg, not merely close to it"
        );
    }

    /// The international mile — EXACTLY 5280 international feet by the same
    /// 1959 agreement that fixes `ft`, so it is asserted against the neighbouring
    /// entry bit-for-bit rather than within a tolerance: the two cannot be
    /// allowed to drift apart, and 5280 × 0.3048 is 1609.344 with no rounding in
    /// binary64. Before this entry the table could spell a stack height in feet
    /// and not a road in miles, which is the unit the entire US onroad
    /// transportation inventory is written in.
    #[test]
    fn test_international_mile() {
        let mi = parse_unit("mi").unwrap();
        let ft = parse_unit("ft").unwrap();
        assert!(
            mi.is_compatible(&parse_unit("m").unwrap()),
            "mi must carry the length dimension"
        );
        assert_eq!(mi.scale, 1609.344, "the international mile is 1609.344 m");
        assert_eq!(
            mi.scale,
            ft.scale * 5280.0,
            "mi must be EXACTLY 5280 ft, not merely close to it"
        );
        // `mph` is NOT a name; the compound is what spells a speed, and it
        // composes from the two entries without a fused alias.
        let mph = parse_unit("mi/h").unwrap();
        assert!(mph.is_compatible(&parse_unit("m/s").unwrap()));
        // 1 mi/h is EXACTLY 0.44704 m/s — the constant MOVES's own SQL writes.
        assert_eq!(mph.scale, 0.44704);
        assert!(
            parse_unit("mph").is_err(),
            "`mph` must stay unresolvable: a fused spelling of a compound the \
             grammar already builds is a second way to say one thing"
        );
    }

    /// The international avoirdupois pound — EXACTLY 0.45359237 kg, and exactly
    /// a two-thousandth of the `short_ton` already in the table, which is what
    /// the assertion pins. Having the ton and not the pound meant the DERIVED
    /// unit was spellable and the one it is defined in was not.
    #[test]
    fn test_avoirdupois_pound() {
        let lb = parse_unit("lb").unwrap();
        assert!(lb.is_compatible(&parse_unit("kg").unwrap()));
        assert_eq!(lb.scale, 0.45359237);
        assert_eq!(
            parse_unit("short_ton").unwrap().scale,
            lb.scale * 2000.0,
            "short_ton is DEFINED as 2000 lb; the two must agree to the last bit"
        );
    }

    /// Mechanical horsepower — 550 ft·lbf/s = 745.6998715822702 W (NIST SP 811
    /// gives 7.456 999 E+02 W). Asserted BOTH ways: against the literal, and
    /// against the foot / pound / standard-gravity product it is built from, so
    /// a future edit to `ft` or `lb` that forgot this entry would fail here.
    /// It is emphatically not the metric horsepower (PS, 735.49875 W).
    #[test]
    fn test_mechanical_horsepower() {
        let hp = parse_unit("hp").unwrap();
        assert!(hp.is_compatible(&parse_unit("W").unwrap()));
        assert_eq!(hp.scale, 745.699_871_582_270_2);
        assert_eq!(
            hp.scale,
            550.0 * parse_unit("ft").unwrap().scale * parse_unit("lb").unwrap().scale * 9.80665,
            "hp must be exactly 550 ft*lbf/s in terms of the table's own ft and lb"
        );
        assert!(
            (hp.scale - 735.498_75).abs() > 1.0,
            "hp must be the MECHANICAL horsepower, not the metric one"
        );
    }

    /// The US liquid gallon — exactly 231 in³ = 3.785411784 L. Pinned against
    /// the litre already in the table, and pinned as NOT the imperial gallon,
    /// which is 20% larger: a dimension-only check cannot tell the two apart.
    #[test]
    fn test_us_liquid_gallon() {
        let gal = parse_unit("gal").unwrap();
        let l = parse_unit("L").unwrap();
        assert_eq!(
            gal.dimensions.get(&Dimension::Length),
            Some(&Rational::int(3)),
            "a gallon is a volume, i.e. [length]^3"
        );
        assert_eq!(gal.scale, 0.003_785_411_784);
        assert!(
            (gal.scale / l.scale - 3.785_411_784).abs() < 1e-12,
            "the US liquid gallon is 3.785411784 L"
        );
        assert!(
            (gal.scale - 0.004_546_09).abs() > 1e-6,
            "must be the US liquid gallon, not the imperial one"
        );
    }

    /// The compounds the four new entries exist to make spellable. Each is a
    /// real MOVES column: `mi/h` is `link.linkAvgSpeed`, `g/(hp*h)` is
    /// `nremissionrate.meanBaseRate`, `lb/(hp*h)` is the brake-specific fuel
    /// consumption carrier, and `g/gal` is `fueltype.fuelDensity`. None of them
    /// needs a table entry of its own — they fall out of the grammar — which is
    /// the point of adding the BASE units rather than the compounds.
    #[test]
    fn test_moves_compounds_derive() {
        assert!(
            parse_unit("mi/h")
                .unwrap()
                .is_compatible(&parse_unit("m/s").unwrap())
        );
        assert!(
            parse_unit("g/(hp*h)")
                .unwrap()
                .is_compatible(&parse_unit("kg/J").unwrap())
        );
        assert!(
            parse_unit("lb/(hp*h)")
                .unwrap()
                .is_compatible(&parse_unit("kg/J").unwrap())
        );
        assert!(
            parse_unit("g/gal")
                .unwrap()
                .is_compatible(&parse_unit("kg/m^3").unwrap())
        );
        assert!(
            parse_unit("g/mi")
                .unwrap()
                .is_compatible(&parse_unit("kg/m").unwrap())
        );
        assert!(
            parse_unit("kJ/gal")
                .unwrap()
                .is_compatible(&parse_unit("J/m^3").unwrap())
        );
    }

    /// SI prefixes resolve against the prefixable symbols — which is what makes
    /// the `µ` → `u` normalization useful — but must not invent units out of
    /// arbitrary identifiers, nor shadow a name that merely looks prefixed.
    #[test]
    fn test_si_prefixes() {
        assert_eq!(
            parse_unit("mg").unwrap().dimensions.get(&Dimension::Mass),
            Some(&Rational::int(1))
        );
        assert!((parse_unit("mg").unwrap().scale - 1e-6).abs() < 1e-18);
        assert!((parse_unit("um").unwrap().scale - 1e-6).abs() < 1e-18);
        assert!((parse_unit("nm").unwrap().scale - 1e-9).abs() < 1e-21);
        assert_eq!(parse_unit("mL").unwrap(), {
            let mut l = parse_unit("L").unwrap();
            l.scale *= 1e-3;
            l.exact = l.exact.multiply(&ExactScale::power_of_ten(-3));
            l
        });
        assert_eq!(
            parse_unit("umol")
                .unwrap()
                .dimensions
                .get(&Dimension::Amount),
            Some(&Rational::int(1))
        );
        // Exact names win over a prefix reading: `min` is minutes, not milli-in.
        assert!((parse_unit("min").unwrap().scale - 60.0).abs() < 1e-9);
        assert_eq!(
            parse_unit("mol")
                .unwrap()
                .dimensions
                .get(&Dimension::Amount),
            Some(&Rational::int(1))
        );
        // A prefix on a non-prefixable identifier is NOT a unit.
        assert!(parse_unit("not_a_unit").is_err());
    }

    /// Counts carry no dimension axis, so composites come out with the right
    /// dimension on their own.
    #[test]
    fn test_counts_are_dimensionless() {
        for c in [
            "molec",
            "molecule",
            "count",
            "individuals",
            "vehicles",
            "units",
        ] {
            assert!(
                parse_unit(c).unwrap().is_dimensionless(),
                "{c} must be dimensionless"
            );
        }
        // individuals/km^2 is an inverse area.
        assert_eq!(
            parse_unit("individuals/km^2")
                .unwrap()
                .dimensions
                .get(&Dimension::Length),
            Some(&Rational::int(-2))
        );
        assert!(parse_unit("psu").unwrap().is_dimensionless());
        assert!(parse_unit("%").unwrap().is_dimensionless());
        assert!((parse_unit("%").unwrap().scale - 0.01).abs() < 1e-12);
        // uatm is a pressure.
        assert!(
            parse_unit("uatm")
                .unwrap()
                .is_compatible(&parse_unit("Pa").unwrap())
        );
    }

    /// `sqrt` HALVES a dimension — it is not a transcendental and does not
    /// require a dimensionless argument. With rational exponents `sqrt(m^3)` is
    /// exactly `m^3/2` rather than "not representable".
    #[test]
    fn test_sqrt_halves_dimensions() {
        let env = env_of(&[("area", "m^2"), ("vol", "m^3")]);

        let e = op("sqrt", vec![Expr::Variable("area".into())]);
        let u = Unit::propagate(&e, &env).unwrap();
        assert_eq!(
            u.dimensions.get(&Dimension::Length),
            Some(&Rational::int(1))
        );
        assert!(check_expression_dimensions(&e, None, &env).is_empty());

        let e = op("sqrt", vec![Expr::Variable("vol".into())]);
        let u = Unit::propagate(&e, &env).unwrap();
        assert_eq!(
            u.dimensions.get(&Dimension::Length),
            Some(&Rational::new(3, 2))
        );
        assert!(
            !check_expression_dimensions(&e, None, &env)
                .iter()
                .any(UnitFinding::is_error),
            "sqrt of an odd power is representable, not an error"
        );
    }

    /// `rad` is a real axis, so a trig function accepts an ANGLE — but `exp` and
    /// `log` still demand strict dimensionlessness.
    #[test]
    fn test_trig_accepts_angle_but_exp_does_not() {
        let env = env_of(&[("theta", "rad"), ("m", "m")]);

        let e = op("sin", vec![Expr::Variable("theta".into())]);
        assert!(
            !check_expression_dimensions(&e, None, &env)
                .iter()
                .any(UnitFinding::is_error),
            "sin(theta) with theta in radians is correct"
        );

        // A LENGTH is still wrong for sin.
        let e = op("sin", vec![Expr::Variable("m".into())]);
        assert!(
            check_expression_dimensions(&e, None, &env)
                .iter()
                .any(UnitFinding::is_error)
        );

        // exp of an angle is a dimensional argument.
        let e = op("exp", vec![Expr::Variable("theta".into())]);
        assert!(
            check_expression_dimensions(&e, None, &env)
                .iter()
                .any(UnitFinding::is_error),
            "exp requires strict dimensionlessness"
        );
    }

    /// esm-spec §4.8.3 / issue #409: "dimensionless" means dimensionless AT
    /// SCALE 1 for every op whose argument the table requires to be
    /// dimensionless, and an ANGLE is admitted to a circular function at ANY
    /// scale because the evaluation path converts it.
    ///
    /// Before this, `ppm` and `percent` satisfied every dimension-only test, so
    /// `log(c [ppm])` was accepted and evaluated the log of the ppm NUMBER with
    /// no diagnostic — and the two readings differ by `ln(1e-6)`.
    #[test]
    fn test_transcendental_argument_must_be_at_scale_one() {
        let env = env_of(&[
            ("c", "ppm"),
            ("p", "percent"),
            ("n", "1"),
            ("theta", "rad"),
            ("lat", "deg"),
            ("solid", "sr"),
        ]);
        let errors = |e: &Expr| {
            check_expression_dimensions(e, None, &env)
                .iter()
                .filter(|f| f.is_error())
                .map(|f| f.message.clone())
                .collect::<Vec<_>>()
        };

        // The strict transcendentals and the inverse circular functions refuse a
        // dimensionless-but-SCALED argument, and the message names the repair.
        for name in [
            "log", "ln", "log10", "exp", "sinh", "cosh", "tanh", "asinh", "acosh", "atanh", "asin",
            "acos", "atan",
        ] {
            let found = errors(&op(name, vec![Expr::Variable("c".into())]));
            assert_eq!(found.len(), 1, "{name}(c [ppm]) must be refused");
            assert!(
                found[0].contains("dimensionless at scale 1")
                    && found[0].contains("(ppm)")
                    && found[0].contains("divide by 1 ppm"),
                "{name}: the diagnostic must name the repair, got {}",
                found[0]
            );
            assert!(
                errors(&op(name, vec![Expr::Variable("n".into())])).is_empty(),
                "{name}(n [1]) is a pure number and stays accepted"
            );
        }
        assert!(
            errors(&op("exp", vec![Expr::Variable("p".into())]))[0].contains("divide by 1 percent"),
            "the percent spelling is named too"
        );

        // A circular function CONVERTS an angle rather than refusing it, at any
        // scale — but a dimensionless-but-scaled argument is refused there too.
        for name in ["sin", "cos", "tan"] {
            assert!(
                errors(&op(name, vec![Expr::Variable("lat".into())])).is_empty(),
                "{name}(lat [deg]) is a conforming document"
            );
            assert!(
                errors(&op(name, vec![Expr::Variable("theta".into())])).is_empty(),
                "{name}(theta [rad]) is a conforming document"
            );
            assert_eq!(
                errors(&op(name, vec![Expr::Variable("p".into())])).len(),
                1,
                "{name}(p [percent]) leaves the reading unstated"
            );
            // `sr` is rad^2. No conversion turns a solid angle into a plane one,
            // so accepting it would multiply by `scale` where `scale^2` was meant.
            assert_eq!(
                errors(&op(name, vec![Expr::Variable("solid".into())])).len(),
                1,
                "{name}(solid [sr]) is not a plane angle"
            );
        }
    }

    /// The repair the diagnostic names — dividing by a quantity whose DECLARED
    /// units carry the scale — reaches scale 1 exactly and is accepted.
    #[test]
    fn test_scaled_dimensionless_repair_is_accepted() {
        let env = env_of(&[("c", "ppm"), ("one_ppm", "ppm")]);
        let quotient = op(
            "/",
            vec![Expr::Variable("c".into()), Expr::Variable("one_ppm".into())],
        );
        assert!(
            !check_expression_dimensions(&op("log", vec![quotient]), None, &env)
                .iter()
                .any(UnitFinding::is_error),
            "log(c / one_ppm) states the reading and must be accepted"
        );
    }

    /// The angle-normalization factor is `None` for everything a rewrite must
    /// leave alone, so a document declaring no scaled angle is not walked.
    #[test]
    fn test_angle_normalization_factor() {
        assert!(angle_normalization_factor(&parse_unit("rad").unwrap()).is_none());
        assert!(angle_normalization_factor(&parse_unit("1").unwrap()).is_none());
        assert!(angle_normalization_factor(&parse_unit("ppm").unwrap()).is_none());
        assert!(angle_normalization_factor(&parse_unit("m").unwrap()).is_none());
        assert!(angle_normalization_factor(&parse_unit("sr").unwrap()).is_none());
        let deg = angle_normalization_factor(&parse_unit("deg").unwrap()).unwrap();
        assert_eq!(deg, std::f64::consts::PI / 180.0);
        // 90 deg is exactly a quarter turn under this factor.
        assert_eq!((90.0 * deg).sin(), 1.0);
    }

    /// The rewrite folds the factor into the ARGUMENT, leaves `rad` alone, and
    /// reaches a `sin` nested inside another argument.
    #[test]
    fn test_normalize_angle_arguments_rewrite() {
        let env = env_of(&[("lat", "deg"), ("theta", "rad")]);
        let untouched = op("sin", vec![Expr::Variable("theta".into())]);
        assert!(
            normalize_angle_arguments(&untouched, &env).is_none(),
            "a `rad` argument is at scale 1 and must not be rewritten"
        );

        let nested = op(
            "+",
            vec![
                Expr::Number(1.0),
                op("cos", vec![Expr::Variable("lat".into())]),
            ],
        );
        let rewritten =
            normalize_angle_arguments(&nested, &env).expect("the `deg` cos is rewritten");
        let Expr::Operator(sum) = &rewritten else {
            panic!("expected the sum back")
        };
        let Expr::Operator(cos) = &sum.args[1] else {
            panic!("expected the cos back")
        };
        let Expr::Operator(product) = &cos.args[0] else {
            panic!("expected the argument to be a product")
        };
        assert_eq!(product.op, "*");
        assert_eq!(product.args[0], Expr::Variable("lat".into()));
        assert_eq!(
            product.args[1],
            Expr::Number(std::f64::consts::PI / 180.0),
            "the factor is the declared scale, applied ONCE"
        );
    }

    /// The rewrite reaches a circular function in EVERY expression-bearing
    /// field, not only `args`: a `makearray` region value, a `table_lookup`
    /// axis input, a `faq` filter and an aggregate key.
    #[test]
    fn test_normalize_angle_arguments_reaches_every_child() {
        let env = env_of(&[("lat", "deg")]);
        let sin = || op("sin", vec![Expr::Variable("lat".into())]);
        let folded = |e: &Expr| {
            let Expr::Operator(n) = e else { return false };
            n.op == "sin"
                && matches!(&n.args[0], Expr::Operator(p) if p.op == "*"
                    && p.args[1] == Expr::Number(std::f64::consts::PI / 180.0))
        };
        let nodes = [
            ExpressionNode {
                op: "makearray".into(),
                values: Some(vec![sin()]),
                ..Default::default()
            },
            ExpressionNode {
                op: "table_lookup".into(),
                axes: Some([("x".to_string(), sin())].into_iter().collect()),
                ..Default::default()
            },
            ExpressionNode {
                op: "faq".into(),
                filter: Some(Box::new(sin())),
                key: Some(Box::new(sin())),
                ..Default::default()
            },
        ];
        for node in nodes {
            let op_name = node.op.clone();
            let rewritten = normalize_angle_arguments(&Expr::operator(node), &env)
                .unwrap_or_else(|| panic!("{op_name}: the nested `deg` sin was not reached"));
            let Expr::Operator(out) = &rewritten else {
                panic!("{op_name}: expected the node back")
            };
            let mut children = 0;
            out.for_each_child(&mut |c| {
                assert!(folded(c), "{op_name}: {c:?} was not folded");
                children += 1;
            });
            assert!(children > 0, "{op_name}");
        }
    }

    /// The fold is IDEMPOTENT. It has to be: the array runtime runs the same
    /// rewrite on a `FlattenedSystem` that `flatten` already folded
    /// (`simulate_array::compile::normalize_model_angle_arguments`), and a
    /// second application would evaluate `sin(90 · 0.01745²)` with no error.
    /// What makes it idempotent is the rewrite's own output: the product
    /// carries a BARE LITERAL, so its unit is undeterminable (§4.8.4) and
    /// `angle_normalization_factor` is never reached a second time.
    #[test]
    fn test_normalize_angle_arguments_is_idempotent() {
        let env = env_of(&[("lat", "deg")]);
        let once = normalize_angle_arguments(&op("sin", vec![Expr::Variable("lat".into())]), &env)
            .expect("the `deg` argument is folded on the first pass");
        assert!(
            normalize_angle_arguments(&once, &env).is_none(),
            "a second pass must find nothing to rewrite, got {once:?}"
        );
    }

    /// The derivative rule leaves the TIME exponent free and requires only that
    /// the non-time dimensions reconcile — so `x` in metres with an RHS in
    /// `m/s^2` is accepted (some time unit reconciles it), while `m` against
    /// `kg` never can be.
    #[test]
    fn test_derivative_time_exponent_is_free() {
        let env = env_of(&[("x", "m"), ("accel", "m/s^2"), ("mass", "kg")]);

        let eq = Equation {
            comment: None,
            lhs: op_with_wrt("D", vec![Expr::Variable("x".into())], "t"),
            rhs: Expr::Variable("accel".into()),
        };
        assert!(
            !check_equation_dimensions(&eq, &env)
                .iter()
                .any(UnitFinding::is_error),
            "ratio m/(m/s^2) = s^2 is a power of time ⇒ reconcilable"
        );

        let eq = Equation {
            comment: None,
            lhs: op_with_wrt("D", vec![Expr::Variable("x".into())], "t"),
            rhs: Expr::Variable("mass".into()),
        };
        assert!(
            check_equation_dimensions(&eq, &env)
                .iter()
                .any(UnitFinding::is_error),
            "ratio m/kg is not a power of time ⇒ provable mismatch"
        );
    }

    /// Every registry entry carries an exact scale, and it is the same number as
    /// the floating-point one used for numeric conversion. An entry that forgot
    /// its exact form would silently default to 1 and fail here.
    #[test]
    fn every_registry_scale_has_a_matching_exact_form() {
        for (name, unit) in build_base_units() {
            let exact = unit.exact_scale().to_f64();
            assert!(
                (exact - unit.scale()).abs() <= 1e-12 * unit.scale().abs(),
                "{name}: exact scale {exact} disagrees with the float {}",
                unit.scale()
            );
        }
    }

    /// The three entries the bindings disagreed on, and the scales whose exact
    /// form is least obvious (esm-spec §4.8.1).
    #[test]
    fn reconciled_registry_scales_are_exact() {
        let exact = |s: &str| parse_unit(s).unwrap().exact_scale().ratio_string().unwrap();
        assert_eq!(exact("Torr"), "20265/152");
        assert_eq!(exact("psi"), "8896443230521/1290320000");
        assert_eq!(exact("degF"), "5/9");
        assert_eq!(exact("deg"), "1/180*pi");
        assert_eq!(exact("mi"), "201168/125");
        assert_eq!(exact("hp"), "37284993579113511/50000000000000");
        assert_eq!(exact("DU"), "268670000000000000000");
        assert_eq!(parse_unit("psi").unwrap().scale(), 6_894.757_293_168_361);
        assert_eq!(parse_unit("Torr").unwrap().scale(), 101_325.0 / 760.0);
        assert_eq!(parse_unit("degF").unwrap().scale(), 5.0 / 9.0);
    }

    #[test]
    fn adding_metres_to_kilometres_is_a_scale_mismatch() {
        let env = env_of(&[("x", "m"), ("y", "km"), ("z", "m")]);
        let var = |n: &str| Expr::Variable(n.into());
        let findings = check_expression_dimensions(&op("+", vec![var("x"), var("y")]), None, &env);
        assert!(
            findings
                .iter()
                .any(|f| f.is_error() && f.message.contains("scale")),
            "m + km must be a scale mismatch: {findings:?}"
        );
        let findings = check_expression_dimensions(&op("+", vec![var("x"), var("z")]), None, &env);
        assert!(!findings.iter().any(UnitFinding::is_error), "{findings:?}");
    }

    #[test]
    fn an_equation_between_m_per_s_and_mi_per_h_is_a_scale_mismatch() {
        let env = env_of(&[
            ("speed_ms", "m/s"),
            ("speed_mph", "mi/h"),
            ("ms_per_mph", "m*h/(mi*s)"),
        ]);
        let var = |n: &str| Expr::Variable(n.into());
        let bare = Equation {
            comment: None,
            lhs: var("speed_ms"),
            rhs: var("speed_mph"),
        };
        assert!(
            check_equation_dimensions(&bare, &env)
                .iter()
                .any(UnitFinding::is_error),
            "m/s = mi/h agrees in dimension and not in scale"
        );
        let converted = Equation {
            comment: None,
            lhs: var("speed_ms"),
            rhs: op("*", vec![var("speed_mph"), var("ms_per_mph")]),
        };
        assert!(
            !check_equation_dimensions(&converted, &env)
                .iter()
                .any(UnitFinding::is_error),
            "mi/h times m*h/(mi*s) is exactly m/s"
        );
    }

    #[test]
    fn ifelse_branches_must_agree_in_scale() {
        let env = env_of(&[("flag", "1"), ("x", "m"), ("y", "km")]);
        let var = |n: &str| Expr::Variable(n.into());
        let e = op("ifelse", vec![var("flag"), var("x"), var("y")]);
        assert!(
            check_expression_dimensions(&e, None, &env)
                .iter()
                .any(UnitFinding::is_error)
        );
    }

    #[test]
    fn an_observed_expression_must_agree_with_its_declared_scale() {
        let env = env_of(&[("y", "km")]);
        let declared = parse_unit("m").unwrap();
        assert!(
            check_expression_dimensions(&Expr::Variable("y".into()), Some(&declared), &env)
                .iter()
                .any(UnitFinding::is_error)
        );
    }

    /// A `const` that declares its units has that unit, a `const` without units
    /// stays undeterminable, and an unresolvable declared unit is listed for the
    /// structural layer (esm-spec §4.8.5).
    #[test]
    fn a_const_with_declared_units_has_that_unit() {
        let konst = |units: Option<&str>| {
            Expr::Operator(std::sync::Arc::new(ExpressionNode {
                op: "const".into(),
                value: Some(serde_json::json!(0.44704)),
                units: units.map(str::to_string),
                ..ExpressionNode::default()
            }))
        };
        let env = env_of(&[
            ("speed_mph", "mi/h"),
            ("speed_ms", "m/s"),
            ("speed_kg", "kg"),
        ]);
        let rhs = |units| op("*", vec![Expr::Variable("speed_mph".into()), konst(units)]);
        let eq = |lhs: &str, units| Equation {
            comment: None,
            lhs: Expr::Variable(lhs.into()),
            rhs: rhs(units),
        };
        assert!(
            !check_equation_dimensions(&eq("speed_ms", Some("m*h/(mi*s)")), &env)
                .iter()
                .any(UnitFinding::is_error),
            "mi/h times a const declared m*h/(mi*s) is exactly m/s"
        );
        assert!(
            check_equation_dimensions(&eq("speed_kg", Some("m*h/(mi*s)")), &env)
                .iter()
                .any(UnitFinding::is_error),
            "declared kg against m/s is a provable mismatch"
        );
        assert!(
            !check_equation_dimensions(&eq("speed_kg", None), &env)
                .iter()
                .any(UnitFinding::is_error),
            "a const without units is undeterminable, so nothing is checked"
        );
        assert_eq!(
            unresolvable_const_units(&rhs(Some("mph"))),
            vec!["mph".to_string()]
        );
        assert!(unresolvable_const_units(&rhs(Some("m*h/(mi*s)"))).is_empty());
    }

    #[test]
    fn const_units_are_gated_at_esm_1_2_0() {
        let doc = |esm: &str| {
            serde_json::json!({
                "esm": esm,
                "models": {"M": {"equations": [{"lhs": "x", "rhs":
                    {"op": "const", "args": [], "value": 1.0, "units": "m"}}]}}
            })
        };
        let err = reject_const_units_pre_v12(&doc("1.1.0")).expect_err("1.1.0 must be rejected");
        assert!(
            err.to_string().contains("const_units_version_too_old"),
            "{err}"
        );
        assert!(
            err.to_string().contains("/models/M/equations/0/rhs"),
            "{err}"
        );
        assert!(reject_const_units_pre_v12(&doc("1.2.0")).is_ok());
    }

    #[test]
    fn scaled_powers_and_roots_stay_exact() {
        assert_eq!(
            parse_unit("km^2")
                .unwrap()
                .power_rational(Rational::new(1, 2))
                .exact_scale(),
            parse_unit("km").unwrap().exact_scale()
        );
        let env = env_of(&[("p", "%"), ("q", "1")]);
        let eq = Equation {
            comment: None,
            lhs: Expr::Variable("q".into()),
            rhs: op("^", vec![Expr::Variable("p".into()), Expr::Integer(2)]),
        };
        assert!(
            check_equation_dimensions(&eq, &env)
                .iter()
                .any(UnitFinding::is_error),
            "a percentage squared has scale 1/10000, not 1"
        );
    }

    #[test]
    fn test_parse_base_units() {
        assert_eq!(
            parse_unit("m").unwrap().dimensions.get(&Dimension::Length),
            Some(&Rational::int(1))
        );
        assert_eq!(
            parse_unit("s").unwrap().dimensions.get(&Dimension::Time),
            Some(&Rational::int(1))
        );
        assert_eq!(
            parse_unit("kg").unwrap().dimensions.get(&Dimension::Mass),
            Some(&Rational::int(1))
        );
    }

    #[test]
    fn test_parse_compound_units() {
        let u = parse_unit("m/s").unwrap();
        assert_eq!(
            u.dimensions.get(&Dimension::Length),
            Some(&Rational::int(1))
        );
        assert_eq!(u.dimensions.get(&Dimension::Time), Some(&Rational::int(-1)));

        let u = parse_unit("mol/L").unwrap();
        assert_eq!(
            u.dimensions.get(&Dimension::Amount),
            Some(&Rational::int(1))
        );
        assert_eq!(
            u.dimensions.get(&Dimension::Length),
            Some(&Rational::int(-3))
        );
    }

    #[test]
    fn test_parse_parenthesised_units() {
        // J/(mol*K) == kg*m^2/(s^2*mol*K)
        let u = parse_unit("J/(mol*K)").unwrap();
        assert_eq!(u.dimensions.get(&Dimension::Mass), Some(&Rational::int(1)));
        assert_eq!(
            u.dimensions.get(&Dimension::Length),
            Some(&Rational::int(2))
        );
        assert_eq!(u.dimensions.get(&Dimension::Time), Some(&Rational::int(-2)));
        assert_eq!(
            u.dimensions.get(&Dimension::Amount),
            Some(&Rational::int(-1))
        );
        assert_eq!(
            u.dimensions.get(&Dimension::Temperature),
            Some(&Rational::int(-1))
        );
    }

    #[test]
    fn test_parse_pascal_and_kg_per_m3() {
        let pa = parse_unit("Pa").unwrap();
        assert_eq!(pa.dimensions.get(&Dimension::Mass), Some(&Rational::int(1)));
        assert_eq!(
            pa.dimensions.get(&Dimension::Length),
            Some(&Rational::int(-1))
        );
        assert_eq!(
            pa.dimensions.get(&Dimension::Time),
            Some(&Rational::int(-2))
        );

        let rho = parse_unit("kg/m^3").unwrap();
        assert_eq!(
            rho.dimensions.get(&Dimension::Mass),
            Some(&Rational::int(1))
        );
        assert_eq!(
            rho.dimensions.get(&Dimension::Length),
            Some(&Rational::int(-3))
        );
    }

    #[test]
    fn test_parse_negative_power() {
        let inv_s = parse_unit("s^-1").unwrap();
        assert_eq!(
            inv_s.dimensions.get(&Dimension::Time),
            Some(&Rational::int(-1))
        );
    }

    #[test]
    fn test_division_is_left_associative() {
        // `W/m^2/K` must parse as `(W/m^2)/K = W/(m^2*K)`, not the
        // right-associative `W/(m^2/K) = W*K/m^2`. W = kg*m^2/s^3, so the
        // correct result is kg*s^-3*K^-1: Temperature power -1 (a +1 would
        // betray the old right-associative bug).
        let u = parse_unit("W/m^2/K").unwrap();
        assert_eq!(u.dimensions.get(&Dimension::Mass), Some(&Rational::int(1)));
        assert_eq!(u.dimensions.get(&Dimension::Time), Some(&Rational::int(-3)));
        assert_eq!(
            u.dimensions.get(&Dimension::Temperature),
            Some(&Rational::int(-1))
        );
        // Length cancels: m^2 (from W) divided by m^2 leaves power 0 (absent).
        assert_eq!(u.dimensions.get(&Dimension::Length), None);
    }

    #[test]
    fn test_find_last_top_level_operator() {
        // `a/b/c` splits at the *last* top-level operator so parse_unit recurses
        // on the left `a/b`, yielding left-associative `(a/b)/c`.
        assert_eq!(find_last_top_level_operator("a/b/c"), Some((3, 1, '/')));
        assert_eq!(find_last_top_level_operator("a*b*c"), Some((3, 1, '*')));
        // `*` and `/` share one precedence level, so the LAST of either wins:
        // `a/b*c` is `(a/b)*c`, not `a/(b*c)`.
        assert_eq!(find_last_top_level_operator("a/b*c"), Some((3, 1, '*')));
        // Separators inside parentheses are hidden from the top-level scan.
        assert_eq!(find_last_top_level_operator("a/(b/c)"), Some((1, 1, '/')));
        assert_eq!(find_last_top_level_operator("kg"), None);
    }

    /// `*` and `/` are the SAME precedence level and associate left-to-right.
    /// Splitting `/` first (looser) pulled every factor after the last `/` into
    /// the denominator, so `J/mol*K` parsed as `J/(mol*K)` — silently negating
    /// K's exponent.
    #[test]
    fn test_mixed_product_quotient_is_left_associative() {
        let u = parse_unit("J/mol*K").unwrap();
        // (J/mol)*K ⇒ temperature appears with power +1, not -1.
        assert_eq!(
            u.dimensions.get(&Dimension::Temperature),
            Some(&Rational::int(1))
        );
        assert_eq!(
            u.dimensions.get(&Dimension::Amount),
            Some(&Rational::int(-1))
        );

        let u = parse_unit("kg/m^3*s").unwrap();
        assert_eq!(u.dimensions.get(&Dimension::Time), Some(&Rational::int(1)));
        assert_eq!(
            u.dimensions.get(&Dimension::Length),
            Some(&Rational::int(-3))
        );

        // Repeated division stays left-associative: m/s/s == m/s^2.
        let u = parse_unit("m/s/s").unwrap();
        assert_eq!(u.dimensions.get(&Dimension::Time), Some(&Rational::int(-2)));
        assert_eq!(
            u.dimensions.get(&Dimension::Length),
            Some(&Rational::int(1))
        );

        // Parenthesised denominators are unaffected.
        let u = parse_unit("J/(mol*K)").unwrap();
        assert_eq!(
            u.dimensions.get(&Dimension::Temperature),
            Some(&Rational::int(-1))
        );
        assert_eq!(
            u.dimensions.get(&Dimension::Amount),
            Some(&Rational::int(-1))
        );
    }

    /// A dangling operator is malformed, not silently dimensionless. `"m/"` used
    /// to parse as metres and `"/s"` as s^-1.
    #[test]
    fn test_dangling_operator_is_rejected() {
        assert!(matches!(parse_unit("m/"), Err(UnitError::ParseError(_))));
        assert!(matches!(parse_unit("/s"), Err(UnitError::ParseError(_))));
        assert!(matches!(parse_unit("kg*"), Err(UnitError::ParseError(_))));
    }

    #[test]
    fn test_parse_mixed_product_quotient() {
        // Mixed `kg*m/s^2` is the Newton: mass^1 length^1 time^-2.
        let u = parse_unit("kg*m/s^2").unwrap();
        assert_eq!(u.dimensions.get(&Dimension::Mass), Some(&Rational::int(1)));
        assert_eq!(
            u.dimensions.get(&Dimension::Length),
            Some(&Rational::int(1))
        );
        assert_eq!(u.dimensions.get(&Dimension::Time), Some(&Rational::int(-2)));
    }

    #[test]
    fn test_unit_compatibility() {
        let m = parse_unit("m").unwrap();
        let cm = parse_unit("cm").unwrap();
        let s = parse_unit("s").unwrap();
        assert!(m.is_compatible(&cm));
        assert!(!m.is_compatible(&s));
    }

    #[test]
    fn test_unit_arithmetic() {
        let m = parse_unit("m").unwrap();
        let s = parse_unit("s").unwrap();
        let v = m.divide(&s);
        assert_eq!(
            v.dimensions.get(&Dimension::Length),
            Some(&Rational::int(1))
        );
        assert_eq!(v.dimensions.get(&Dimension::Time), Some(&Rational::int(-1)));
        assert_eq!(
            m.power(2).dimensions.get(&Dimension::Length),
            Some(&Rational::int(2))
        );
    }

    #[test]
    fn test_dimensional_consistency() {
        let m = parse_unit("m").unwrap();
        let cm = parse_unit("cm").unwrap();
        let s = parse_unit("s").unwrap();
        assert!(check_dimensional_consistency(&m, &cm).is_ok());
        assert!(check_dimensional_consistency(&m, &s).is_err());
    }

    #[test]
    fn test_unit_conversion() {
        let m = parse_unit("m").unwrap();
        let cm = parse_unit("cm").unwrap();
        assert!((convert_units(1.0, &m, &cm).unwrap() - 100.0).abs() < 1e-10);
        assert!((convert_units(100.0, &cm, &m).unwrap() - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_unknown_unit() {
        let result = parse_unit("unknown_unit");
        assert!(matches!(result, Err(UnitError::UnknownUnit(_))));
    }

    // ---- propagate ---------------------------------------------------

    #[test]
    fn propagate_literal_is_indeterminate() {
        // A bare numeric literal carries NO knowable dimension: `273.15` may be
        // a pure number or a temperature offset. It is therefore indeterminate,
        // not dimensionless — a checker whose findings are hard errors must not
        // fabricate a dimension it cannot know.
        let env = HashMap::new();
        assert!(matches!(
            Unit::propagate(&Expr::Number(2.5), &env),
            Err(UnitError::UnknownUnit(_))
        ));

        // Additively, though, a literal is NEUTRAL: it adopts its sibling's
        // dimension, so `T - 273.15` is a temperature and raises no finding.
        let env = env_of(&[("T", "K")]);
        let expr = op("-", vec![Expr::Variable("T".into()), Expr::Number(273.15)]);
        let u = Unit::propagate(&expr, &env).unwrap();
        assert_eq!(
            u.dimensions.get(&Dimension::Temperature),
            Some(&Rational::int(1))
        );
        assert!(check_expression_dimensions(&expr, None, &env).is_empty());

        // But a sum with NO determinable operand is undeterminable, never
        // dimensionless (esm-spec §4.8.3, §4.8.4): for an integer as well as a
        // float literal, in `min`/`max` as in `+`/`-`, and under a unary `+`.
        let literal_pairs = [
            vec![Expr::Number(1.0), Expr::Number(2.0)],
            vec![Expr::Integer(1), Expr::Integer(2)],
        ];
        for args in literal_pairs {
            for name in ["+", "-", "min", "max"] {
                assert!(
                    matches!(
                        Unit::propagate(&op(name, args.clone()), &HashMap::new()),
                        Err(UnitError::UnknownUnit(_))
                    ),
                    "{name} over literals must be undeterminable"
                );
            }
        }
        for lit in [Expr::Integer(2), Expr::Number(2.5)] {
            assert!(matches!(
                Unit::propagate(&op("+", vec![lit]), &HashMap::new()),
                Err(UnitError::UnknownUnit(_))
            ));
        }
    }

    /// A sum whose operands are partly undeterminable has the unit of the
    /// determinable ones (esm-spec §4.8.3), and a unary `+` carries its
    /// operand's unit.
    #[test]
    fn propagate_sum_with_undeterminable_operand_takes_known_unit() {
        let env = env_of(&[("x", "m"), ("y", "m")]);
        let half_y = op("*", vec![Expr::Number(0.5), Expr::Variable("y".into())]);
        for name in ["+", "-", "min", "max"] {
            let expr = op(name, vec![Expr::Variable("x".into()), half_y.clone()]);
            let u = Unit::propagate(&expr, &env).unwrap();
            assert_eq!(
                u.dimensions.get(&Dimension::Length),
                Some(&Rational::int(1)),
                "{name}"
            );
        }
        let plus_x = op("+", vec![Expr::Variable("x".into())]);
        let u = Unit::propagate(&plus_x, &env).unwrap();
        assert_eq!(
            u.dimensions.get(&Dimension::Length),
            Some(&Rational::int(1))
        );
    }

    /// `ifelse` takes the unit of its determinable branch, whichever branch that
    /// is, and with no determinable branch it is undeterminable (esm-spec
    /// §4.8.3).
    #[test]
    fn ifelse_takes_the_unit_of_its_determinable_branch() {
        let env = env_of(&[("x", "m"), ("y", "m"), ("c", "1")]);
        let cond = op(">", vec![Expr::Variable("c".into()), Expr::Integer(0)]);
        let x = Expr::Variable("x".into());
        let half_y = op("*", vec![Expr::Number(0.5), Expr::Variable("y".into())]);
        let branch_pairs = [
            (x.clone(), half_y.clone()),
            (half_y.clone(), x.clone()),
            (Expr::Integer(2), x.clone()),
        ];
        for (then_branch, else_branch) in branch_pairs {
            let expr = op("ifelse", vec![cond.clone(), then_branch, else_branch]);
            let u = Unit::propagate(&expr, &env).unwrap();
            assert_eq!(
                u.dimensions.get(&Dimension::Length),
                Some(&Rational::int(1))
            );
        }
        let literal_branches = op("ifelse", vec![cond, Expr::Integer(1), Expr::Number(2.0)]);
        assert!(matches!(
            Unit::propagate(&literal_branches, &env),
            Err(UnitError::UnknownUnit(_))
        ));
    }

    /// A boolean connective places no requirement on its operands' units, but a
    /// mismatch inside an operand is still reported (esm-spec §4.8.3).
    #[test]
    fn connectives_report_a_mismatch_inside_their_operands() {
        let env = env_of(&[("x", "m"), ("z", "kg"), ("c", "1")]);
        let mismatch = op(
            ">",
            vec![Expr::Variable("x".into()), Expr::Variable("z".into())],
        );
        let ok = op(">", vec![Expr::Variable("c".into()), Expr::Integer(0)]);
        for expr in [
            op("not", vec![mismatch.clone()]),
            op("and", vec![mismatch.clone(), ok.clone()]),
            op("or", vec![ok.clone(), mismatch.clone()]),
        ] {
            assert!(
                !check_expression_dimensions(&expr, None, &env).is_empty(),
                "{expr:?}"
            );
        }
        let dimensional_operands = op(
            "and",
            vec![Expr::Variable("x".into()), Expr::Variable("z".into())],
        );
        assert!(check_expression_dimensions(&dimensional_operands, None, &env).is_empty());
    }

    /// An exponent is read BY VALUE, so `L^2` still yields an area even though
    /// the literal `2` has no dimension of its own.
    #[test]
    fn propagate_literal_exponent_is_read_by_value() {
        let env = env_of(&[("L", "m")]);
        let expr = op("^", vec![Expr::Variable("L".into()), Expr::Integer(2)]);
        let u = Unit::propagate(&expr, &env).unwrap();
        assert_eq!(
            u.dimensions.get(&Dimension::Length),
            Some(&Rational::int(2))
        );
    }

    /// An `Expr::Integer` exponent is an i64 and is read at full width. A value
    /// above `i32::MAX` has no representable rational dimension exponent, so it
    /// reports `UnknownUnit` naming the ORIGINAL value — it must not truncate
    /// to a bogus i32 and propagate a confidently wrong dimension.
    #[test]
    fn propagate_huge_integer_exponent_does_not_wrap() {
        let env = env_of(&[("L", "m")]);
        let exp = i64::from(i32::MAX) + 6;
        let expr = op("^", vec![Expr::Variable("L".into()), Expr::Integer(exp)]);
        let err = Unit::propagate(&expr, &env).unwrap_err();
        match err {
            UnitError::UnknownUnit(msg) => assert!(
                msg.contains(&exp.to_string()),
                "finding must name the untruncated exponent, got: {msg}"
            ),
            other => panic!("expected UnknownUnit, got {other:?}"),
        }
    }

    #[test]
    fn propagate_variable_lookup() {
        let env = env_of(&[("v", "m/s")]);
        let u = Unit::propagate(&Expr::Variable("v".into()), &env).unwrap();
        assert_eq!(
            u.dimensions.get(&Dimension::Length),
            Some(&Rational::int(1))
        );
        assert_eq!(u.dimensions.get(&Dimension::Time), Some(&Rational::int(-1)));
    }

    #[test]
    fn propagate_time_variable_default() {
        let env = HashMap::new();
        let u = Unit::propagate(&Expr::Variable("t".into()), &env).unwrap();
        assert_eq!(u.dimensions.get(&Dimension::Time), Some(&Rational::int(1)));
    }

    #[test]
    fn propagate_unknown_variable_errors() {
        let env = HashMap::new();
        let err = Unit::propagate(&Expr::Variable("nope".into()), &env).unwrap_err();
        assert!(matches!(err, UnitError::UnknownUnit(_)));
    }

    #[test]
    fn propagate_addition_matches() {
        // Both in metres: `+` needs its operands to agree in scale as well as
        // dimension (esm-spec §4.8.3).
        let env = env_of(&[("h1", "m"), ("h2", "m")]);
        let e = op(
            "+",
            vec![Expr::Variable("h1".into()), Expr::Variable("h2".into())],
        );
        let u = Unit::propagate(&e, &env).unwrap();
        assert_eq!(
            u.dimensions.get(&Dimension::Length),
            Some(&Rational::int(1))
        );
    }

    #[test]
    fn propagate_addition_mismatch_errors() {
        let env = env_of(&[("h", "m"), ("t", "s")]);
        let e = op(
            "+",
            vec![Expr::Variable("h".into()), Expr::Variable("t".into())],
        );
        let err = Unit::propagate(&e, &env).unwrap_err();
        assert!(matches!(err, UnitError::DimensionMismatch(_)));
    }

    #[test]
    fn propagate_multiplication_combines_dims() {
        let env = env_of(&[("rho", "kg/m^3"), ("v", "m/s")]);
        let e = op(
            "*",
            vec![Expr::Variable("rho".into()), Expr::Variable("v".into())],
        );
        let u = Unit::propagate(&e, &env).unwrap();
        assert_eq!(u.dimensions.get(&Dimension::Mass), Some(&Rational::int(1)));
        assert_eq!(
            u.dimensions.get(&Dimension::Length),
            Some(&Rational::int(-2))
        );
        assert_eq!(u.dimensions.get(&Dimension::Time), Some(&Rational::int(-1)));
    }

    #[test]
    fn propagate_division_subtracts_dims() {
        let env = env_of(&[("m1", "kg"), ("V", "m^3")]);
        let e = op(
            "/",
            vec![Expr::Variable("m1".into()), Expr::Variable("V".into())],
        );
        let u = Unit::propagate(&e, &env).unwrap();
        assert_eq!(u.dimensions.get(&Dimension::Mass), Some(&Rational::int(1)));
        assert_eq!(
            u.dimensions.get(&Dimension::Length),
            Some(&Rational::int(-3))
        );
    }

    #[test]
    fn propagate_power_integer() {
        let env = env_of(&[("s", "m")]);
        let e = op("^", vec![Expr::Variable("s".into()), Expr::Number(3.0)]);
        let u = Unit::propagate(&e, &env).unwrap();
        assert_eq!(
            u.dimensions.get(&Dimension::Length),
            Some(&Rational::int(3))
        );
    }

    #[test]
    fn propagate_power_rejects_dimensional_exponent() {
        // h ^ t is invalid: exponent must be dimensionless.
        let env = env_of(&[("h", "m"), ("t", "s")]);
        let e = op(
            "^",
            vec![Expr::Variable("h".into()), Expr::Variable("t".into())],
        );
        let err = Unit::propagate(&e, &env).unwrap_err();
        assert!(matches!(err, UnitError::DimensionMismatch(_)));
    }

    #[test]
    fn propagate_derivative_divides_by_wrt() {
        // Canonical bead example: D(h)/dt should have dimensions of v (m/s).
        // `t` must be DECLARED for the derivative to be determinate — an
        // undeclared independent variable leaves it unknown (see
        // `derivative_of_undeclared_time`).
        let env = env_of(&[("h", "m"), ("v", "m/s"), ("t", "s")]);
        let dh = op_with_wrt("D", vec![Expr::Variable("h".into())], "t");
        let dh_dim = Unit::propagate(&dh, &env).unwrap();
        let v_dim = Unit::propagate(&Expr::Variable("v".into()), &env).unwrap();
        assert!(dh_dim.is_compatible(&v_dim));
    }

    /// With an UNDECLARED `t` the derivative's dimension is unknown, so no time
    /// unit is assumed and the plain LHS/RHS comparison is skipped.
    #[test]
    fn propagate_derivative_of_undeclared_time_is_unknown() {
        let env = env_of(&[("h", "m")]);
        let dh = op_with_wrt("D", vec![Expr::Variable("h".into())], "t");
        assert!(matches!(
            Unit::propagate(&dh, &env),
            Err(UnitError::UnknownUnit(_))
        ));
    }

    #[test]
    fn propagate_exp_requires_dimensionless_arg() {
        let env = env_of(&[("h", "m")]);
        let e = op("exp", vec![Expr::Variable("h".into())]);
        let err = Unit::propagate(&e, &env).unwrap_err();
        assert!(matches!(err, UnitError::DimensionMismatch(_)));
    }

    #[test]
    fn propagate_sqrt_halves_even_powers() {
        let env = env_of(&[("a", "m^2")]);
        let e = op("sqrt", vec![Expr::Variable("a".into())]);
        let u = Unit::propagate(&e, &env).unwrap();
        assert_eq!(
            u.dimensions.get(&Dimension::Length),
            Some(&Rational::int(1))
        );
    }

    #[test]
    fn propagate_ifelse_branches_must_match() {
        let env = env_of(&[("h", "m"), ("x", "s")]);
        let e = op(
            "ifelse",
            vec![
                Expr::Number(1.0),
                Expr::Variable("h".into()),
                Expr::Variable("x".into()),
            ],
        );
        assert!(matches!(
            Unit::propagate(&e, &env).unwrap_err(),
            UnitError::DimensionMismatch(_)
        ));
    }

    #[test]
    fn propagate_comparison_requires_matching_dims() {
        let env = env_of(&[("h", "m"), ("x", "s")]);
        let e = op(
            ">",
            vec![Expr::Variable("h".into()), Expr::Variable("x".into())],
        );
        assert!(matches!(
            Unit::propagate(&e, &env).unwrap_err(),
            UnitError::DimensionMismatch(_)
        ));
    }

    /// The spatial-calculus sugar ops (`grad` / `div` / `laplacian` / `curl` /
    /// `∇` / `integral`) are ORDINARY open-tier rewrite targets with NO
    /// privileged dimensional rule (esm-spec §4.2 / §4.8.3). Their dimension is
    /// UNDETERMINABLE until a discretization rule lowers them, so propagation
    /// yields an unknown dimension via the same no-rule path as any unregistered
    /// user op — never the old bespoke "divide by a metre denominator"
    /// behaviour. Per esm-spec §4.8.4 an unlowered op is reported as UNKNOWN and
    /// SKIPPED: it emits NO finding at all, matching the other four bindings
    /// (Julia/TS/Python/Go). `grad` is not singled out.
    #[test]
    fn propagate_spatial_sugar_ops_are_undeterminable() {
        let env = env_of(&[("c", "mol/m^3")]);
        for op_name in ["grad", "div", "laplacian", "curl", "∇", "integral"] {
            let e = op(op_name, vec![Expr::Variable("c".into())]);
            // No dimensional rule ⇒ undeterminable, so `propagate` cannot
            // conclude a unit.
            assert!(
                matches!(Unit::propagate(&e, &env), Err(UnitError::UnknownUnit(_))),
                "{op_name} must have an undeterminable dimension"
            );
            // Report-unknown-and-skip: NO finding at all — not a promotable
            // `Error`, and not even a non-blocking `Analysis`. The checker just
            // declines to fabricate a dimension it cannot know.
            let findings = check_expression_dimensions(&e, None, &env);
            assert!(
                findings.is_empty(),
                "{op_name} must emit no finding (report unknown + skip): {findings:?}"
            );
        }
    }

    #[test]
    fn propagate_faq_body() {
        // aggregate { expr: x * scale } with x in m/s and a dimensionless
        // `scale` -> result m/s. (A bare numeric literal would leave the product
        // INDETERMINATE — see `propagate_literal_is_indeterminate` — so the
        // element dimension is exercised with a declared factor.)
        let env = env_of(&[("x", "m/s"), ("scale", "1")]);
        let body = op(
            "*",
            vec![Expr::Variable("x".into()), Expr::Variable("scale".into())],
        );
        let node = Expr::operator(ExpressionNode {
            op: "faq".to_string(),
            args: vec![],
            expr: Some(Box::new(body)),
            ..ExpressionNode::default()
        });
        let u = Unit::propagate(&node, &env).unwrap();
        assert_eq!(
            u.dimensions.get(&Dimension::Length),
            Some(&Rational::int(1))
        );
        assert_eq!(u.dimensions.get(&Dimension::Time), Some(&Rational::int(-1)));
    }

    #[test]
    fn propagate_index_inherits_from_source() {
        let env = env_of(&[("arr", "kg")]);
        let node = op(
            "index",
            vec![Expr::Variable("arr".into()), Expr::Number(0.0)],
        );
        let u = Unit::propagate(&node, &env).unwrap();
        assert_eq!(u.dimensions.get(&Dimension::Mass), Some(&Rational::int(1)));
    }

    #[test]
    fn propagate_broadcast_dispatches_to_fn() {
        // broadcast(fn='+', a, b) with matching dims works.
        let env = env_of(&[("a", "m"), ("b", "m")]);
        let node = Expr::operator(ExpressionNode {
            op: "broadcast".to_string(),
            broadcast_fn: Some("+".to_string()),
            args: vec![Expr::Variable("a".into()), Expr::Variable("b".into())],
            ..ExpressionNode::default()
        });
        let u = Unit::propagate(&node, &env).unwrap();
        assert_eq!(
            u.dimensions.get(&Dimension::Length),
            Some(&Rational::int(1))
        );
    }

    #[test]
    fn propagate_unknown_operator_reports_unknown_and_skips() {
        // An operator with no dimensional rule means the CHECKER cannot
        // conclude anything — it is not a defect in the file. Per esm-spec
        // §4.8.4 it reports UNKNOWN and SKIPS: an unknown dimension and NO
        // finding at all (matching Julia/TS/Python/Go), never a promotable
        // `Error` and never a silent `dimensionless`.
        let env = HashMap::new();
        let e = op("zorp", vec![Expr::Number(1.0)]);
        let findings = check_expression_dimensions(&e, None, &env);
        assert!(findings.is_empty(), "expected no finding, got {findings:?}");
        assert!(matches!(
            Unit::propagate(&e, &env),
            Err(UnitError::UnknownUnit(_))
        ));
    }

    /// An undeterminable operand must not SUPPRESS a provable mismatch
    /// elsewhere in the same expression: propagation yields `Unknown` for it and
    /// keeps walking instead of bailing at the first one.
    #[test]
    fn undeterminable_operand_does_not_hide_a_real_mismatch() {
        let env = env_of(&[("len", "m"), ("mass", "kg")]);
        // (zorp(1) + len) + mass — the unknown op is undeterminable (reported as
        // unknown and skipped, so it emits NO finding of its own), but the `len`
        // vs `mass` mismatch after it is still proven.
        let expr = op(
            "+",
            vec![
                op("zorp", vec![Expr::Number(1.0)]),
                Expr::Variable("len".into()),
                Expr::Variable("mass".into()),
            ],
        );
        let findings = check_expression_dimensions(&expr, None, &env);
        assert!(
            findings.iter().any(UnitFinding::is_error),
            "the m-vs-kg mismatch must still be proven: {findings:?}"
        );
        // The no-rule op adds nothing: the only finding is the real mismatch.
        assert!(
            findings.iter().all(UnitFinding::is_error),
            "the undeterminable op must not add a finding of its own: {findings:?}"
        );
    }

    #[test]
    fn validate_equation_dimensions_passes() {
        // `t` is declared, so D(h)/dt is determinate (m/s) and matches v.
        let env = env_of(&[("h", "m"), ("v", "m/s"), ("t", "s")]);
        let eq = Equation {
            comment: None,
            lhs: op_with_wrt("D", vec![Expr::Variable("h".into())], "t"),
            rhs: Expr::Variable("v".into()),
        };
        validate_equation_dimensions(&eq, &env).unwrap();
    }

    #[test]
    fn validate_equation_dimensions_detects_mismatch() {
        // With `t` DECLARED the derivative is determinate, so the plain
        // LHS-vs-RHS comparison proves the mismatch: D(h)/dt is m/s, h is m.
        let env = env_of(&[("h", "m"), ("v", "m/s"), ("t", "s")]);
        let eq = Equation {
            comment: None,
            lhs: op_with_wrt("D", vec![Expr::Variable("h".into())], "t"),
            rhs: Expr::Variable("h".into()),
        };
        assert!(validate_equation_dimensions(&eq, &env).is_err());

        // The consistent equation is accepted.
        let eq = Equation {
            comment: None,
            lhs: op_with_wrt("D", vec![Expr::Variable("h".into())], "t"),
            rhs: Expr::Variable("v".into()),
        };
        assert!(validate_equation_dimensions(&eq, &env).is_ok());
    }

    /// With an UNDECLARED `t`, the weaker time-ratio rule still proves a
    /// derivative equation that NO time unit could reconcile — while accepting
    /// one that some time unit could.
    #[test]
    fn derivative_time_ratio_rule() {
        let env = env_of(&[("h", "m"), ("v", "m/s"), ("mass", "kg")]);

        // d(m)/dt = kg: ratio m/kg is not a power of time ⇒ provable mismatch.
        let eq = Equation {
            comment: None,
            lhs: op_with_wrt("D", vec![Expr::Variable("h".into())], "t"),
            rhs: Expr::Variable("mass".into()),
        };
        let findings = check_equation_dimensions(&eq, &env);
        assert!(findings.iter().any(UnitFinding::is_error), "{findings:?}");

        // d(m)/dt = m/s: ratio is exactly `s` ⇒ reconcilable, no finding.
        let eq = Equation {
            comment: None,
            lhs: op_with_wrt("D", vec![Expr::Variable("h".into())], "t"),
            rhs: Expr::Variable("v".into()),
        };
        assert!(
            !check_equation_dimensions(&eq, &env)
                .iter()
                .any(UnitFinding::is_error)
        );

        // The dimensionless toy idiom `D(x)/dt = -x` (ratio time^0) is accepted:
        // it is the idiom used across the shared valid corpus.
        let env = env_of(&[("x", "dimensionless")]);
        let eq = Equation {
            comment: None,
            lhs: op_with_wrt("D", vec![Expr::Variable("x".into())], "t"),
            rhs: op("-", vec![Expr::Variable("x".into())]),
        };
        assert!(
            !check_equation_dimensions(&eq, &env)
                .iter()
                .any(UnitFinding::is_error)
        );
    }

    /// Build a bare state [`crate::ModelVariable`] carrying only a declared
    /// unit string (all other metadata omitted).
    fn state_var_with_units(units: Option<&str>) -> crate::ModelVariable {
        crate::test_support::var(crate::VariableType::Unknown, units)
    }

    #[test]
    fn build_unit_env_unparseable_unit_is_unknown_not_dimensionless() {
        let mut variables = indexmap::IndexMap::new();
        // `x` has an unparseable unit; `y` has a real unit of metres.
        variables.insert("x".to_string(), state_var_with_units(Some("not_a_unit")));
        variables.insert("y".to_string(), state_var_with_units(Some("m")));

        let (env, failures) = build_unit_env(&variables);

        // The unparseable variable is OMITTED (unknown), NOT coerced to a
        // dimensionless unit; the parseable one is present.
        assert!(
            !env.contains_key("x"),
            "unparseable unit must not be in env"
        );
        assert!(env.contains_key("y"));

        // The failure is REPORTED, naming the offending variable and the string
        // the author actually wrote, so the caller can raise a hard
        // `unit_parse_error` at `/models/<M>/variables/x` (esm-spec §4.8.4).
        assert_eq!(
            failures,
            vec![UnitParseFailure {
                name: "x".to_string(),
                units: "not_a_unit".to_string(),
            }]
        );

        // An equation `y = x` referencing the unknown-unit variable propagates
        // to UnknownUnit — a non-blocking, skip-for-dim-checking outcome — NOT
        // a false DimensionMismatch (which a dimensionless coercion of `x`
        // against `y`'s metres would have produced).
        let eq = Equation {
            comment: None,
            lhs: Expr::Variable("y".into()),
            rhs: Expr::Variable("x".into()),
        };
        let result = validate_equation_dimensions(&eq, &env);
        assert!(
            matches!(result, Err(UnitError::UnknownUnit(_))),
            "expected UnknownUnit (skipped), not a false mismatch, got {result:?}"
        );
    }

    fn konst(value: f64, units: Option<&str>) -> Expr {
        Expr::operator(ExpressionNode {
            op: "const".to_string(),
            value: Some(serde_json::json!(value)),
            units: units.map(str::to_string),
            ..ExpressionNode::default()
        })
    }

    #[test]
    fn constant_exponents_fold_to_exact_rationals() {
        let int = Expr::Integer;
        let cases = [
            (op("/", vec![int(3), int(2)]), Rational::new(3, 2)),
            (op("-", vec![int(1)]), Rational::int(-1)),
            (op("-", vec![int(5), int(3)]), Rational::int(2)),
            (
                op(
                    "/",
                    vec![op("+", vec![int(1), int(2)]), op("*", vec![int(4), int(2)])],
                ),
                Rational::new(3, 8),
            ),
            (konst(1.5, None), Rational::new(3, 2)),
            (konst(0.5, Some("1")), Rational::new(1, 2)),
            (Expr::Number(0.25), Rational::new(1, 4)),
        ];
        for (expr, expected) in cases {
            assert_eq!(constant_exponent(&expr), Some(expected), "{expr:?}");
        }
    }

    #[test]
    fn symbolic_exponents_do_not_fold() {
        let int = Expr::Integer;
        let var = |n: &str| Expr::Variable(n.into());
        for expr in [
            var("n"),
            op("+", vec![var("n"), int(1)]),
            konst(50.0, Some("%")),
            konst(2.0, Some("m")),
            op("/", vec![int(1), int(0)]),
            op("/", vec![int(1), int(2), int(3)]),
            op("sqrt", vec![int(4)]),
            op(
                "*",
                vec![int(i64::from(i32::MAX)), int(i64::from(i32::MAX))],
            ),
        ] {
            assert_eq!(constant_exponent(&expr), None, "{expr:?}");
        }
    }

    #[test]
    fn a_folded_exponent_is_checked_and_a_symbolic_one_is_skipped() {
        let env = env_of(&[("T", "K"), ("n", "1"), ("T_32", "K^(3/2)"), ("T_2", "K^2")]);
        let var = |n: &str| Expr::Variable(n.into());
        let three_halves = op("/", vec![Expr::Integer(3), Expr::Integer(2)]);
        let pow = op("^", vec![var("T"), three_halves]);
        let ok = check_expression_dimensions(&pow, env.get("T_32"), &env);
        assert!(!ok.iter().any(UnitFinding::is_error), "{ok:?}");
        let bad = check_expression_dimensions(&pow, env.get("T_2"), &env);
        assert!(bad.iter().any(UnitFinding::is_error), "{bad:?}");
        let sym = op("^", vec![var("T"), var("n")]);
        let skipped = check_expression_dimensions(&sym, env.get("T_2"), &env);
        assert!(!skipped.iter().any(UnitFinding::is_error), "{skipped:?}");
    }

    #[test]
    fn an_faq_filter_is_dimension_checked() {
        let env = env_of(&[("x", "m"), ("y", "m"), ("z", "kg")]);
        let var = |n: &str| Expr::Variable(n.into());
        let faq = |filter: Expr| {
            Expr::operator(ExpressionNode {
                op: "faq".to_string(),
                filter: Some(Box::new(filter)),
                expr: Some(Box::new(konst(1.0, Some("1")))),
                ..ExpressionNode::default()
            })
        };
        let bad = check_expression_dimensions(&faq(op("<=", vec![var("x"), var("z")])), None, &env);
        assert!(bad.iter().any(UnitFinding::is_error), "{bad:?}");
        let ok = check_expression_dimensions(&faq(op("<=", vec![var("x"), var("y")])), None, &env);
        assert!(!ok.iter().any(UnitFinding::is_error), "{ok:?}");
    }

    fn faq_node(body: Expr, semiring: Option<&str>, reduce: Option<&str>) -> ExpressionNode {
        ExpressionNode {
            op: "faq".to_string(),
            expr: Some(Box::new(body)),
            semiring: semiring.map(str::to_string),
            reduce: reduce.map(str::to_string),
            ..ExpressionNode::default()
        }
    }

    #[test]
    fn array_operators_carry_their_element_unit() {
        let env = env_of(&[("v", "h"), ("w", "1"), ("a", "m"), ("b", "kg")]);
        let var = |n: &str| Expr::Variable(n.into());
        let unit = |e: &Expr| propagate_dim(e, &env, &mut Vec::new());
        let h = parse_unit("h").unwrap();
        let idx = op("index", vec![var("v"), var("i")]);
        assert_eq!(unit(&idx).known(), Some(&h));
        assert!(unit(&op("index", vec![])).known().is_none());
        // faq by semiring.
        for (sr, rd) in [
            (Some("sum_product"), None),
            (Some("max_product"), None),
            (Some("min_sum"), None),
            (None, Some("max")),
            (None, None),
        ] {
            let f = Expr::operator(faq_node(idx.clone(), sr, rd));
            assert_eq!(unit(&f).known(), Some(&h), "{sr:?} {rd:?}");
        }
        let prod = Expr::operator(faq_node(idx.clone(), None, Some("*")));
        assert!(unit(&prod).known().is_none());
        let wprod = Expr::operator(faq_node(
            op("index", vec![var("w"), var("i")]),
            None,
            Some("*"),
        ));
        assert!(unit(&wprod).known().is_some_and(Unit::is_pure_number));
        let any = Expr::operator(faq_node(idx.clone(), Some("bool_and_or"), None));
        assert!(unit(&any).known().is_some_and(Unit::is_dimensionless));
        let mut distinct = faq_node(idx.clone(), None, None);
        distinct.distinct = Some(true);
        assert!(unit(&Expr::operator(distinct)).known().is_none());
        // makearray: a literal region is skipped; disagreeing regions are a mismatch.
        let mk = |values: Vec<Expr>| {
            Expr::operator(ExpressionNode {
                op: "makearray".to_string(),
                values: Some(values),
                ..ExpressionNode::default()
            })
        };
        assert_eq!(
            unit(&mk(vec![var("v"), Expr::Number(0.0)])).known(),
            Some(&h)
        );
        assert!(unit(&mk(vec![])).known().is_none());
        let mut f = Vec::new();
        propagate_dim(&mk(vec![var("a"), var("b")]), &env, &mut f);
        assert!(f.iter().any(UnitFinding::is_error), "{f:?}");
        // concat agrees like +.
        let mut f = Vec::new();
        propagate_dim(&op("concat", vec![var("a"), var("b")]), &env, &mut f);
        assert!(f.iter().any(UnitFinding::is_error), "{f:?}");
    }

    #[test]
    fn array_operators_check_every_operand() {
        let env = env_of(&[("v", "h"), ("a", "m"), ("b", "kg")]);
        let var = |n: &str| Expr::Variable(n.into());
        let h = parse_unit("h").unwrap();
        let bad = || op("+", vec![var("a"), var("b")]);
        // A mismatch in a subscript or shape operand is reported, while the
        // result is still the source array's unit.
        for name in ["index", "reshape", "transpose"] {
            let mut f = Vec::new();
            let dim = propagate_dim(&op(name, vec![var("v"), bad()]), &env, &mut f);
            assert_eq!(dim.known(), Some(&h), "{name}");
            assert!(f.iter().any(UnitFinding::is_error), "{name}: {f:?}");
            // A dimensional subscript is not itself an error.
            let mut f = Vec::new();
            propagate_dim(&op(name, vec![var("v"), var("a")]), &env, &mut f);
            assert!(!f.iter().any(UnitFinding::is_error), "{name}: {f:?}");
        }
        // faq, makearray and a broadcast with no `fn` walk their operands too.
        let mut fq = faq_node(var("v"), None, None);
        fq.args = vec![bad()];
        let mut mk = ExpressionNode {
            op: "makearray".to_string(),
            values: Some(vec![var("v")]),
            ..ExpressionNode::default()
        };
        mk.args = vec![bad()];
        let bc = ExpressionNode {
            op: "broadcast".to_string(),
            args: vec![bad()],
            ..ExpressionNode::default()
        };
        for node in [fq, mk, bc] {
            let name = node.op.clone();
            let mut f = Vec::new();
            propagate_dim(&Expr::operator(node), &env, &mut f);
            assert!(f.iter().any(UnitFinding::is_error), "{name}: {f:?}");
        }
    }
}