SOLIDWORKS Bridge

//! THE SIMULATOR'S SOLID, BY COLUMNS (0.11.0) — so an exam part's mass comes
//! out of the simulator to the key, partial cuts included.
//!
//! The arithmetic model (mock_model.rs) keeps a box and a volume per body:
//! exact for a boss, a cut wholly inside the material and a revolve, and only
//! approximate for a cut that is partly outside (the CSWA Tool Block's R41
//! slot, its under-notch). This model keeps every feature as a PRIMITIVE —
//! a prism (a sketch's closed profiles between two heights, optionally
//! drafted) or a full revolve — added or subtracted in order (CSG), and
//! integrates the part by COLUMNS: along one model axis every primitive's
//! extent is an exact set of intervals (the profile sliced by a line —
//! lines and arcs intersected in closed form), and the columns sample the
//! other two axes on a grid. The error is the grid's, at internal edges only
//! (the grid is aligned to the part's box): about 0.1 % on the Tool Block,
//! well inside the exam's 1 %.
//!
//! What it cannot represent (fillets, chamfers, shells, circular patterns,
//! partial revolves, a revolve whose axis is not the column axis) makes it
//! decline, saying why; the arithmetic model answers then.

use serde::{Deserialize, Serialize};

use crate::modeling::{self, PlaneFrame, V3};
use crate::script::Entity;

/// One feature's solid, added or removed.
#[derive(Clone, Debug, Serialize, Deserialize)]
pub struct Prim {
    pub add: bool,
    /// The feature that made it (mirror and patterns copy by feature).
    pub feature: String,
    pub kind: Kind,
}

#[derive(Clone, Debug, Serialize, Deserialize)]
pub enum Kind {
    /// The closed profiles of `entities` (plane coordinates of `frame`; nested
    /// profiles alternate solid / hole) between heights `lo` and `hi` along the
    /// frame's normal. `draft`: (k, h0) — the profile is offset INWARD by
    /// k·|h − h0| at height h (k < 0: outward); rectangles and circles only.
    Prism { frame: PlaneFrame, lo: f64, hi: f64, entities: Vec<Entity>, draft: Option<(f64, f64)> },
    /// A FULL revolve of `entities` (in `frame` at height `h`) about the line
    /// entity `axis`.
    Revolve { frame: PlaneFrame, h: f64, entities: Vec<Entity>, axis: usize },
}

/// The model axis (0 X, 1 Y, 2 Z) a unit vector of the default frames runs
/// along, and its sign.
fn axis_of(v: V3) -> (usize, f64) {
    let i = (0..3).max_by(|a, b| v[*a].abs().partial_cmp(&v[*b].abs()).unwrap()).unwrap_or(0);
    (i, v[i].signum())
}

/// Where the line `p + s·e` (plane coordinates) crosses the boundary of the
/// profiles: the sorted parameters `s`. Between consecutive pairs the line is
/// inside (even-odd: a hole inside a profile is outside).
pub fn crossings(entities: &[Entity], p: (f64, f64), e: (f64, f64), shrink: f64) -> Vec<f64> {
    let mut out = Vec::new();
    let seg = |a: (f64, f64), b: (f64, f64), out: &mut Vec<f64>| {
        let d = (b.0 - a.0, b.1 - a.1);
        let den = e.0 * d.1 - e.1 * d.0;
        if den.abs() < 1e-18 {
            return;
        }
        let w = (a.0 - p.0, a.1 - p.1);
        let s = (w.0 * d.1 - w.1 * d.0) / den;
        let t = (w.0 * e.1 - w.1 * e.0) / den;
        if (0.0..1.0).contains(&t) {
            out.push(s);
        }
    };
    let circle = |c: (f64, f64), r: f64, sweep: Option<(f64, f64)>, out: &mut Vec<f64>| {
        if r <= 0.0 {
            return;
        }
        let w = (p.0 - c.0, p.1 - c.1);
        let b = w.0 * e.0 + w.1 * e.1;
        let cc = w.0 * w.0 + w.1 * w.1 - r * r;
        let disc = b * b - cc;
        if disc <= 0.0 {
            return;
        }
        for s in [-b - disc.sqrt(), -b + disc.sqrt()] {
            let q = (p.0 + s * e.0 - c.0, p.1 + s * e.1 - c.1);
            match sweep {
                None => out.push(s),
                Some((t0, sw)) => {
                    let phi = q.1.atan2(q.0);
                    // Half-open at the END of the arc, as segments are.
                    let tau = std::f64::consts::TAU;
                    let rel = if sw >= 0.0 { (phi - t0).rem_euclid(tau) } else { (t0 - phi).rem_euclid(tau) };
                    if rel < sw.abs() {
                        out.push(s);
                    }
                }
            }
        }
    };
    for en in entities {
        match *en {
            Entity::Rect { x0, y0, x1, y1 } => {
                let (a0, b0, a1, b1) = (x0.min(x1) + shrink, y0.min(y1) + shrink, x0.max(x1) - shrink, y0.max(y1) - shrink);
                if a0 >= a1 || b0 >= b1 {
                    continue;
                }
                for (a, b) in [((a0, b0), (a1, b0)), ((a1, b0), (a1, b1)), ((a1, b1), (a0, b1)), ((a0, b1), (a0, b0))] {
                    seg(a, b, &mut out);
                }
            }
            Entity::Circle { cx, cy, r } => circle((cx, cy), r - shrink, None, &mut out),
            Entity::Line { x0, y0, x1, y1, construction: false } => seg((x0, y0), (x1, y1), &mut out),
            Entity::Arc { construction: false, .. } => {
                if let Some((cx, cy, r, t0, sw)) = crate::script::arc_geom(en) {
                    circle((cx, cy), r, Some((t0, sw)), &mut out);
                }
            }
            _ => {}
        }
    }
    out.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
    out
}

/// Is `p` inside the closed profile that `members` (indices into `entities`)
/// make? Even-odd along a ray in a direction no drawn vertex lines up with
/// (0.11.4, the simulator's sketch regions).
pub fn point_in(entities: &[Entity], members: &[usize], p: (f64, f64)) -> bool {
    let sub: Vec<Entity> = members.iter().filter_map(|i| entities.get(*i).cloned()).collect();
    let a: f64 = 0.618_033_988_75;
    crossings(&sub, p, (a.cos(), a.sin()), 0.0).iter().filter(|s| **s > 0.0).count() % 2 == 1
}

/// Points along a profile's boundary, `n` per entity.
fn boundary_samples(entities: &[Entity], members: &[usize], n: usize) -> Vec<(f64, f64)> {
    let mut out = Vec::new();
    let seg = |a: (f64, f64), b: (f64, f64), out: &mut Vec<(f64, f64)>| {
        for k in 0..n {
            let t = (k as f64 + 0.5) / n as f64;
            out.push((a.0 + (b.0 - a.0) * t, a.1 + (b.1 - a.1) * t));
        }
    };
    for i in members {
        match entities.get(*i).cloned() {
            Some(Entity::Rect { x0, y0, x1, y1 }) => {
                for (a, b) in [((x0, y0), (x1, y0)), ((x1, y0), (x1, y1)), ((x1, y1), (x0, y1)), ((x0, y1), (x0, y0))] {
                    seg(a, b, &mut out);
                }
            }
            Some(Entity::Circle { cx, cy, r }) => {
                for k in 0..4 * n {
                    let t = std::f64::consts::TAU * (k as f64 + 0.5) / (4 * n) as f64;
                    out.push((cx + r * t.cos(), cy + r * t.sin()));
                }
            }
            Some(Entity::Line { x0, y0, x1, y1, .. }) => seg((x0, y0), (x1, y1), &mut out),
            Some(e @ Entity::Arc { .. }) => {
                if let Some((cx, cy, r, t0, sw)) = crate::script::arc_geom(&e) {
                    for k in 0..n {
                        let t = t0 + sw * (k as f64 + 0.5) / n as f64;
                        out.push((cx + r * t.cos(), cy + r * t.sin()));
                    }
                }
            }
            None => {}
        }
    }
    out
}

/// Do two closed profiles' boundaries CROSS — one partly inside the other
/// (sampled: 64 points per entity, either way round)?
pub fn profiles_cross(entities: &[Entity], a: &[usize], b: &[usize]) -> bool {
    let mixed = |samples: &[(f64, f64)], other: &[usize]| {
        let ins = samples.iter().filter(|p| point_in(entities, other, **p)).count();
        ins > 0 && ins < samples.len()
    };
    mixed(&boundary_samples(entities, a, 64), b) || mixed(&boundary_samples(entities, b, 64), a)
}

/// Pairs of crossings as intervals (an odd count drops the last: a tangency).
fn pairs(c: &[f64]) -> Vec<(f64, f64)> {
    c.chunks_exact(2).map(|w| (w[0], w[1])).filter(|(a, b)| b > a).collect()
}

fn union(a: &[(f64, f64)], b: &[(f64, f64)]) -> Vec<(f64, f64)> {
    let mut all: Vec<(f64, f64)> = a.iter().chain(b.iter()).copied().collect();
    all.sort_by(|x, y| x.0.partial_cmp(&y.0).unwrap_or(std::cmp::Ordering::Equal));
    let mut out: Vec<(f64, f64)> = Vec::new();
    for (s, e) in all {
        match out.last_mut() {
            Some(l) if s <= l.1 => l.1 = l.1.max(e),
            _ => out.push((s, e)),
        }
    }
    out
}

fn minus(a: &[(f64, f64)], b: &[(f64, f64)]) -> Vec<(f64, f64)> {
    let mut out = Vec::new();
    for &(s, e) in a {
        let mut cur = vec![(s, e)];
        for &(bs, be) in b {
            let mut next = Vec::new();
            for (cs, ce) in cur {
                if be <= cs || bs >= ce {
                    next.push((cs, ce));
                } else {
                    if bs > cs {
                        next.push((cs, bs));
                    }
                    if be < ce {
                        next.push((be, ce));
                    }
                }
            }
            cur = next;
        }
        out.extend(cur);
    }
    out
}

fn intersect(a: &[(f64, f64)], lo: f64, hi: f64) -> Vec<(f64, f64)> {
    a.iter().filter_map(|&(s, e)| {
        let (s, e) = (s.max(lo), e.min(hi));
        (e > s).then_some((s, e))
    }).collect()
}

/// A rectangle or circle alone (no holes): what a draft is modelled on.
fn single_shape(entities: &[Entity]) -> Option<&Entity> {
    let solid: Vec<&Entity> = entities.iter().filter(|e| !e.is_construction()).collect();
    match solid.as_slice() {
        [e @ (Entity::Rect { .. } | Entity::Circle { .. })] => Some(e),
        _ => None,
    }
}

/// Can this primitive be represented? (Why not, if not.)
pub fn supported(kind: &Kind) -> Result<(), String> {
    match kind {
        Kind::Prism { draft: Some(_), entities, .. } if single_shape(entities).is_none() => {
            Err("a draft on a profile other than one rectangle or one circle".into())
        }
        _ => Ok(()),
    }
}

/// 8-point Gauss-Legendre nodes and weights on [-1, 1].
const GAUSS: [(f64, f64); 8] = [
    (-0.960_289_856_497_536_2, 0.101_228_536_290_376_26),
    (-0.796_666_477_413_626_7, 0.222_381_034_453_374_47),
    (-0.525_532_409_916_329_0, 0.313_706_645_877_887_3),
    (-0.183_434_642_495_649_8, 0.362_683_783_378_362_0),
    (0.183_434_642_495_649_8, 0.362_683_783_378_362_0),
    (0.525_532_409_916_329_0, 0.313_706_645_877_887_3),
    (0.796_666_477_413_626_7, 0.222_381_034_453_374_47),
    (0.960_289_856_497_536_2, 0.101_228_536_290_376_26),
];

/// Quadrature nodes on [a, b]: Gauss-Legendre after the substitution
/// x = mid + half·sin φ, which smooths the square-root behaviour a circle's
/// chord length has at its ends (a band ends exactly there) — and integrates
/// a straight edge's linear length to ~1e-11.
fn nodes(a: f64, b: f64) -> Vec<(f64, f64)> {
    let (mid, half) = ((a + b) / 2.0, (b - a) / 2.0);
    let h = std::f64::consts::FRAC_PI_2;
    GAUSS.iter().map(|(x, w)| {
        let phi = x * h;
        (mid + half * phi.sin(), w * h * half * phi.cos())
    }).collect()
}

/// The coordinates along one plane direction where a profile has a vertex or
/// an extreme: where the slices of it change shape.
fn profile_coords(entities: &[Entity], along_u: bool) -> Vec<f64> {
    let pick = |x: f64, y: f64| if along_u { x } else { y };
    let mut out = Vec::new();
    for e in entities.iter().filter(|e| !e.is_construction()) {
        match *e {
            Entity::Rect { x0, y0, x1, y1 } => out.extend([pick(x0, y0), pick(x1, y1)]),
            Entity::Circle { cx, cy, r } => out.extend([pick(cx, cy) - r, pick(cx, cy), pick(cx, cy) + r]),
            Entity::Line { x0, y0, x1, y1, .. } => out.extend([pick(x0, y0), pick(x1, y1)]),
            Entity::Arc { x0, y0, x1, y1, .. } => {
                out.extend([pick(x0, y0), pick(x1, y1)]);
                let (lo, hi) = if along_u { crate::script::extent_along(e, 1.0, 0.0) } else { crate::script::extent_along(e, 0.0, 1.0) };
                out.extend([lo, hi]);
            }
        }
    }
    out
}

/// The column axis: a revolve's axis if there is one (every revolve must
/// share it), else the axis most prisms are extruded along.
fn column_axis(prims: &[Prim]) -> Result<usize, String> {
    let mut col: Option<usize> = None;
    for p in prims {
        if let Kind::Revolve { frame, entities, axis, .. } = &p.kind {
            let ((ax, ay), (bx, by)) = entities.get(*axis).and_then(|e| e.ends()).ok_or("a revolve without its axis line")?;
            let d3 = modeling::sub(frame.point(bx, by, 0.0), frame.point(ax, ay, 0.0));
            let (i, _) = axis_of(d3);
            if d3[i].abs() / modeling::norm(d3) < 1.0 - 1e-9 {
                return Err("a revolve axis that is not along a model axis".into());
            }
            if col.map(|c| c != i).unwrap_or(false) {
                return Err("revolves about different axes".into());
            }
            col = Some(i);
        }
    }
    Ok(col.unwrap_or_else(|| {
        let mut count = [0usize; 3];
        for p in prims {
            if let Kind::Prism { frame, .. } = &p.kind {
                count[axis_of(frame.n).0] += 1;
            }
        }
        (0..3).max_by_key(|i| (count[*i], 3 - *i)).unwrap_or(2)
    }))
}

/// The plane coordinate (and sign) of model axis `m` in `frame`, if `m` lies
/// in the plane: (true for u / false for v, sign).
fn in_plane(frame: &PlaneFrame, m: usize) -> Option<(bool, f64)> {
    let (xa, xs) = axis_of(frame.x);
    let (ya, ys) = axis_of(frame.y);
    if xa == m {
        Some((true, xs))
    } else if ya == m {
        Some((false, ys))
    } else {
        None
    }
}

/// Breakpoints of a revolve along model axis `m` (≠ the column axis): where
/// the distance from its axis crosses a profile radius.
fn revolve_radii(frame: &PlaneFrame, h: f64, entities: &[Entity], axis: usize) -> Option<(V3, Vec<f64>)> {
    let ((ax, ay), (bx, by)) = entities.get(axis)?.ends()?;
    let l = (bx - ax).hypot(by - ay);
    let nrm = (-(by - ay) / l, (bx - ax) / l);
    let o = nrm.0 * ax + nrm.1 * ay;
    let mut radii = Vec::new();
    for (i, e) in entities.iter().enumerate() {
        if i == axis || e.is_construction() {
            continue;
        }
        let (lo, hi) = crate::script::extent_along(e, nrm.0, nrm.1);
        radii.extend([(lo - o).abs(), (hi - o).abs()]);
        if let Some(((x0, y0), (x1, y1))) = e.ends() {
            radii.extend([(nrm.0 * x0 + nrm.1 * y0 - o).abs(), (nrm.0 * x1 + nrm.1 * y1 - o).abs()]);
        }
    }
    Some((frame.point(ax, ay, h), radii))
}

/// THE INTEGRATION. Along the column axis every primitive is an exact set of
/// intervals. Across it, the part is cut into BANDS at every coordinate where
/// a profile has a vertex or an extreme (and at every primitive's own ends),
/// and each band is integrated by 8-point Gauss-Legendre — so a polygon's
/// straight edges carry no grid error at all, and arcs converge fast. Returns
/// the volume (m³) and the centre of mass.
pub fn integrate(prims: &[Prim], bounds: (V3, V3), _n: usize) -> Result<(f64, V3), String> {
    if prims.is_empty() {
        return Err("nothing to integrate".into());
    }
    let c = column_axis(prims)?;
    let (r_ax, q_ax) = match c {
        0 => (1, 2),
        1 => (0, 2),
        _ => (0, 1),
    };
    let (lo, hi) = bounds;
    // Breakpoints along q, for the whole part.
    let mut qb: Vec<f64> = vec![lo[q_ax], hi[q_ax]];
    for p in prims {
        match &p.kind {
            Kind::Prism { frame, lo: h0, hi: h1, entities, .. } => {
                if axis_of(frame.n).0 == q_ax {
                    qb.extend([*h0, *h1]);
                } else if let Some((u, sg)) = in_plane(frame, q_ax) {
                    qb.extend(profile_coords(entities, u).into_iter().map(|x| sg * x));
                }
            }
            Kind::Revolve { frame, h, entities, axis } => {
                if let Some((a3, radii)) = revolve_radii(frame, *h, entities, *axis) {
                    for r in radii {
                        qb.extend([a3[q_ax] - r, a3[q_ax] + r]);
                    }
                    qb.push(a3[q_ax]);
                }
            }
        }
    }
    let clean = |mut v: Vec<f64>, a: f64, b: f64| -> Vec<f64> {
        v.retain(|x| x.is_finite() && *x >= a && *x <= b);
        v.sort_by(|x, y| x.partial_cmp(y).unwrap_or(std::cmp::Ordering::Equal));
        v.dedup_by(|x, y| (*x - *y).abs() < 1e-12);
        v
    };
    let qb = clean(qb, lo[q_ax], hi[q_ax]);
    let continuous = prims.iter().any(|p| match &p.kind {
        Kind::Revolve { .. } => true,
        Kind::Prism { frame, draft, .. } => draft.is_some() || axis_of(frame.n).0 == q_ax,
    });
    let mut vol = 0.0;
    let mut mom = [0.0; 3];
    for w in qb.windows(2) {
        let (q0, q1) = (w[0], w[1]);
        if q1 - q0 < 1e-15 {
            continue;
        }
        for (q, wq) in nodes(q0, q1) {
            // Breakpoints along r, on this row.
            let mut rb: Vec<f64> = vec![lo[r_ax], hi[r_ax]];
            for p in prims {
                match &p.kind {
                    Kind::Prism { frame, lo: h0, hi: h1, entities, .. } => {
                        let pa = axis_of(frame.n).0;
                        if pa == r_ax {
                            rb.extend([*h0, *h1]);
                        } else if pa == c {
                            // The row's line in the plane: q fixed, r varying.
                            if let (Some((uq, sq)), Some((_, sr))) = (in_plane(frame, q_ax), in_plane(frame, r_ax)) {
                                let (pt, dir) = if uq { ((sq * q, 0.0), (0.0, 1.0)) } else { ((0.0, sq * q), (1.0, 0.0)) };
                                rb.extend(crossings(entities, pt, dir, 0.0).into_iter().map(|s| sr * s));
                            }
                        } else if let Some((u, sg)) = in_plane(frame, r_ax) {
                            rb.extend(profile_coords(entities, u).into_iter().map(|x| sg * x));
                        }
                    }
                    Kind::Revolve { frame, h, entities, axis } => {
                        if let Some((a3, radii)) = revolve_radii(frame, *h, entities, *axis) {
                            let dq = q - a3[q_ax];
                            for r in radii {
                                if r * r > dq * dq {
                                    let d = (r * r - dq * dq).sqrt();
                                    rb.extend([a3[r_ax] - d, a3[r_ax] + d]);
                                }
                            }
                        }
                    }
                }
            }
            let rb = clean(rb, lo[r_ax], hi[r_ax]);
            for v in rb.windows(2) {
                let (r0, r1) = (v[0], v[1]);
                if r1 - r0 < 1e-15 {
                    continue;
                }
                let at: Vec<(f64, f64)> = if continuous { nodes(r0, r1) } else { vec![((r0 + r1) / 2.0, r1 - r0)] };
                for (r, wr) in at {
                    let mut point = [0.0; 3];
                    point[r_ax] = r;
                    point[q_ax] = q;
                    let mut set: Vec<(f64, f64)> = Vec::new();
                    for p in prims {
                        let iv = column(&p.kind, c, point);
                        set = if p.add { union(&set, &iv) } else { minus(&set, &iv) };
                    }
                    let cell = wq * wr;
                    for (a, b) in set {
                        let len = b - a;
                        vol += len * cell;
                        mom[c] += (b * b - a * a) / 2.0 * cell;
                        mom[r_ax] += len * r * cell;
                        mom[q_ax] += len * q * cell;
                    }
                }
            }
        }
    }
    if vol <= 0.0 {
        return Ok((0.0, [0.0; 3]));
    }
    Ok((vol, modeling::scale(mom, 1.0 / vol)))
}

/// One primitive's extent along the column through `point` (axis `c`).
fn column(kind: &Kind, c: usize, point: V3) -> Vec<(f64, f64)> {
    match kind {
        Kind::Prism { frame, lo, hi, entities, draft } => {
            let (pa, _) = axis_of(frame.n);
            if pa == c {
                // In-plane point; along the column, the prism's own heights.
                let (u, v) = (modeling::dot(point, frame.x), modeling::dot(point, frame.y));
                if let (Some((k, h0)), Some(shape)) = (draft, single_shape(entities)) {
                    // How far from the sketch plane this point stays inside.
                    let room = match *shape {
                        Entity::Circle { cx, cy, r } => r - (u - cx).hypot(v - cy),
                        Entity::Rect { x0, y0, x1, y1 } => {
                            ((x1 - x0).abs() / 2.0 - (u - (x0 + x1) / 2.0).abs()).min((y1 - y0).abs() / 2.0 - (v - (y0 + y1) / 2.0).abs())
                        }
                        _ => return Vec::new(),
                    };
                    if *k > 0.0 {
                        if room <= 0.0 {
                            return Vec::new();
                        }
                        let m = room / k;
                        return intersect(&[(h0 - m, h0 + m)], *lo, *hi);
                    }
                    let m = (-room / -k).max(0.0);
                    return minus(&[(*lo, *hi)], &[(h0 - m, h0 + m)]);
                }
                let xs = crossings(entities, (u, v), (1.0, 0.0), 0.0);
                let inside = xs.iter().filter(|s| **s > 0.0).count() % 2 == 1;
                if inside { vec![(*lo, *hi)] } else { Vec::new() }
            } else {
                // The column runs IN the sketch plane: inside the prism's
                // heights, the profile sliced by the column's line.
                let h = point[pa];
                if h < *lo || h > *hi {
                    return Vec::new();
                }
                let (xa, xs) = axis_of(frame.x);
                let (ya, ys) = axis_of(frame.y);
                let (p, e, sign) = if xa == c {
                    ((0.0, ys * point[ya]), (1.0, 0.0), xs)
                } else if ya == c {
                    ((xs * point[xa], 0.0), (0.0, 1.0), ys)
                } else {
                    return Vec::new();
                };
                let (shrink, ents): (f64, std::borrow::Cow<[Entity]>) = match (draft, single_shape(entities)) {
                    (Some((k, h0)), Some(_)) => (k * (h - h0).abs(), std::borrow::Cow::Borrowed(entities.as_slice())),
                    _ => (0.0, std::borrow::Cow::Borrowed(entities.as_slice())),
                };
                pairs(&crossings(&ents, p, e, shrink))
                    .into_iter()
                    .map(|(a, b)| if sign > 0.0 { (a, b) } else { (-b, -a) })
                    .collect()
            }
        }
        Kind::Revolve { frame, h, entities, axis } => {
            let Some(((ax, ay), (bx, by))) = entities.get(*axis).and_then(|e| e.ends()) else { return Vec::new() };
            let l = (bx - ax).hypot(by - ay);
            let dir = ((bx - ax) / l, (by - ay) / l);
            let nrm = (-dir.1, dir.0);
            let a3 = frame.point(ax, ay, *h);
            let d3 = modeling::scale(modeling::sub(frame.point(bx, by, *h), a3), 1.0 / l);
            // The column's distance from the axis.
            let mut off = modeling::sub(point, a3);
            off[c] = 0.0;
            let rho = modeling::norm(off);
            // Which side of the axis the profile lies on.
            let side = entities.iter().enumerate().filter(|(i, e)| *i != *axis && !e.is_construction()).map(|(_, e)| {
                let (lo, hi) = crate::script::extent_along(e, nrm.0, nrm.1);
                let o = nrm.0 * ax + nrm.1 * ay;
                (lo - o) + (hi - o)
            }).sum::<f64>().signum();
            let p = (ax + nrm.0 * side * rho, ay + nrm.1 * side * rho);
            let profile: Vec<Entity> = entities.iter().enumerate().filter(|(i, _)| *i != *axis || !entities[*i].is_construction()).map(|(_, e)| e.clone()).collect();
            let base = a3[c];
            let sgn = d3[c].signum();
            pairs(&crossings(&profile, p, dir, 0.0))
                .into_iter()
                .map(|(s0, s1)| {
                    let (u0, u1) = (base + sgn * s0, base + sgn * s1);
                    (u0.min(u1), u0.max(u1))
                })
                .collect()
        }
    }
}

/// A primitive reflected across the model plane `coord` along `axis`.
pub fn reflect(kind: &Kind, axis: usize, coord: f64) -> Kind {
    let flip = |es: &[Entity], along_u: bool, about: f64| -> Vec<Entity> {
        let m = |x: f64| 2.0 * about - x;
        es.iter()
            .map(|e| match *e {
                Entity::Rect { x0, y0, x1, y1 } => {
                    if along_u { Entity::Rect { x0: m(x1), y0, x1: m(x0), y1 } } else { Entity::Rect { x0, y0: m(y1), x1, y1: m(y0) } }
                }
                Entity::Circle { cx, cy, r } => if along_u { Entity::Circle { cx: m(cx), cy, r } } else { Entity::Circle { cx, cy: m(cy), r } },
                Entity::Line { x0, y0, x1, y1, construction } => {
                    if along_u { Entity::Line { x0: m(x0), y0, x1: m(x1), y1, construction } } else { Entity::Line { x0, y0: m(y0), x1, y1: m(y1), construction } }
                }
                Entity::Arc { cx, cy, x0, y0, x1, y1, ccw, construction } => {
                    if along_u {
                        Entity::Arc { cx: m(cx), cy, x0: m(x0), y0, x1: m(x1), y1, ccw: !ccw, construction }
                    } else {
                        Entity::Arc { cx, cy: m(cy), x0, y0: m(y0), x1, y1: m(y1), ccw: !ccw, construction }
                    }
                }
            })
            .collect()
    };
    let in_plane = |frame: &PlaneFrame, es: &[Entity]| -> Option<Vec<Entity>> {
        let (xa, xs) = axis_of(frame.x);
        let (ya, ys) = axis_of(frame.y);
        if xa == axis {
            Some(flip(es, true, xs * coord))
        } else if ya == axis {
            Some(flip(es, false, ys * coord))
        } else {
            None
        }
    };
    match kind {
        Kind::Prism { frame, lo, hi, entities, draft } => {
            if axis_of(frame.n).0 == axis {
                Kind::Prism { frame: *frame, lo: 2.0 * coord - hi, hi: 2.0 * coord - lo, entities: entities.clone(), draft: draft.map(|(k, h0)| (k, 2.0 * coord - h0)) }
            } else {
                Kind::Prism { frame: *frame, lo: *lo, hi: *hi, entities: in_plane(frame, entities).unwrap_or_else(|| entities.clone()), draft: *draft }
            }
        }
        Kind::Revolve { frame, h, entities, axis: a } => {
            if axis_of(frame.n).0 == axis {
                Kind::Revolve { frame: *frame, h: 2.0 * coord - h, entities: entities.clone(), axis: *a }
            } else {
                Kind::Revolve { frame: *frame, h: *h, entities: in_plane(frame, entities).unwrap_or_else(|| entities.clone()), axis: *a }
            }
        }
    }
}

/// A primitive moved by `d` along model axis `axis`.
pub fn translate(kind: &Kind, axis: usize, d: f64) -> Kind {
    // A translation is two reflections: about 0, then about d/2.
    reflect(&reflect(kind, axis, 0.0), axis, d / 2.0)
}

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

    fn front() -> PlaneFrame {
        PlaneFrame::defaults()[0]
    }

    fn sq(x0: f64, y0: f64, x1: f64, y1: f64) -> Entity {
        Entity::Rect { x0, y0, x1, y1 }
    }

    #[test]
    fn a_partial_cut_is_integrated_where_the_box_model_guesses() {
        // A 40 x 20 x 10 block; a Ø30 circle centred on its right edge cut through.
        let block = Prim { add: true, feature: "B".into(), kind: Kind::Prism { frame: front(), lo: 0.0, hi: 0.010, entities: vec![sq(0.0, 0.0, 0.040, 0.020)], draft: None } };
        let cut = Prim { add: false, feature: "C".into(), kind: Kind::Prism { frame: front(), lo: 0.0, hi: 0.010, entities: vec![Entity::Circle { cx: 0.040, cy: 0.010, r: 0.015 }], draft: None } };
        let (v, _) = integrate(&[block.clone(), cut], ([0.0, 0.0, 0.0], [0.040, 0.020, 0.010]), 800).unwrap();
        // The disc ∩ rectangle: the half disc clipped to |y| ≤ 10 of a r 15 circle.
        let r: f64 = 15.0;
        let seg = |h: f64| r * r * (h / r).asin() + h * (r * r - h * h).sqrt(); // ∫_{-h}^{h} sqrt(r²−y²) dy
        let removed = seg(10.0) * 10.0; // (half of the band area: x from 25 to 40) × depth
        let want = 40.0 * 20.0 * 10.0 - removed;
        assert!((v * 1e9 - want).abs() / want < 1e-6, "{} vs {want}", v * 1e9);
        // The block alone: exact (edges on the grid), centroid at its middle.
        let (v, c) = integrate(&[block], ([0.0, 0.0, 0.0], [0.040, 0.020, 0.010]), 200).unwrap();
        assert!((v * 1e9 - 8000.0).abs() < 1e-6 && (c[0] - 0.020).abs() < 1e-12 && (c[2] - 0.005).abs() < 1e-12);
    }

    #[test]
    fn a_prism_across_the_columns_is_sliced_exactly() {
        // A cylinder along Y (Top Plane) through a block extruded along Z: the
        // columns run along Z, the cylinder is sliced by each one.
        let top = PlaneFrame::defaults()[1];
        let block = Prim { add: true, feature: "B".into(), kind: Kind::Prism { frame: front(), lo: 0.0, hi: 0.060, entities: vec![sq(0.0, 0.0, 0.100, 0.020)], draft: None } };
        let hole = Prim { add: false, feature: "H".into(), kind: Kind::Prism { frame: top, lo: 0.0, hi: 0.020, entities: vec![Entity::Circle { cx: 0.025, cy: -0.030, r: 0.010 }], draft: None } };
        let (v, c) = integrate(&[block, hole], ([0.0, 0.0, 0.0], [0.100, 0.020, 0.060]), 600).unwrap();
        let want = 120000.0 - std::f64::consts::PI * 100.0 * 20.0;
        assert!((v * 1e9 - want).abs() / want < 1e-7, "{} vs {want}", v * 1e9);
        assert!((c[0] * 1e3 - 51.38).abs() < 0.005 && (c[1] * 1e3 - 10.0).abs() < 1e-9 && (c[2] * 1e3 - 30.0).abs() < 0.01, "{c:?}");
    }

    #[test]
    fn interval_sets_union_and_subtract() {
        assert_eq!(union(&[(0.0, 2.0)], &[(1.0, 3.0), (5.0, 6.0)]), vec![(0.0, 3.0), (5.0, 6.0)]);
        assert_eq!(minus(&[(0.0, 10.0)], &[(2.0, 3.0), (5.0, 12.0)]), vec![(0.0, 2.0), (3.0, 5.0)]);
        assert_eq!(intersect(&[(0.0, 10.0)], 4.0, 20.0), vec![(4.0, 10.0)]);
        // A ring: crossings of a line through two concentric circles alternate.
        let ring = [Entity::Circle { cx: 0.0, cy: 0.0, r: 2.0 }, Entity::Circle { cx: 0.0, cy: 0.0, r: 1.0 }];
        assert_eq!(pairs(&crossings(&ring, (-5.0, 0.0), (1.0, 0.0), 0.0)), vec![(3.0, 4.0), (6.0, 7.0)]);
    }
}