ogeom-offset 0.9.9

//! Fitted skins, fillings and projections measured between their samples:
//! each result sampled densely and compared with the geometry it stands
//! for, independently of the points it was fitted to.
#![allow(clippy::unwrap_used, clippy::expect_used, reason = "test code")]

use ogeom_core::Tolerances;
use ogeom_geom::{Curve3d as _, Surface as _, SurfaceGeometry};
use ogeom_math::{Circle, Direction, Frame, Plane, Point};
use ogeom_topo::{EdgeRepr, Filter, Model, NodeData, Shape, ShapeType, explore};

const T: Tolerances = Tolerances::millimetres();
const PI: f64 = core::f64::consts::PI;

/// The spline surfaces of every face of `shape`.
fn spline_surfaces(model: &Model, shape: &Shape) -> Vec<ogeom_geom::BSplineSurface> {
    let mut out = Vec::new();
    for face in explore(model, shape, Filter::OfType(ShapeType::Face)).unwrap() {
        let NodeData::Face(data) = model.node(&face).unwrap().data() else {
            continue;
        };
        if let Some(SurfaceGeometry::BSpline(patch)) = model.geometry().surface(data.surface) {
            out.push(patch.clone());
        }
    }
    out
}

/// The parameter `k / n` of the way across `range`.
fn across(range: (f64, f64), k: usize, n: usize) -> f64 {
    #[allow(clippy::cast_precision_loss)]
    let f = k as f64 / n as f64;
    range.0 + (range.1 - range.0) * f
}

/// The distance from `p` to the polyline `line`.
fn to_polyline(p: Point, line: &[Point]) -> f64 {
    let mut best = f64::INFINITY;
    for w in line.windows(2) {
        let d = w[1] - w[0];
        let len2 = d.dot(d);
        let t = if len2 > 0.0 {
            ((p - w[0]).dot(d) / len2).clamp(0.0, 1.0)
        } else {
            0.0
        };
        best = best.min(p.distance(w[0] + d * t));
    }
    best
}

fn edge_of(model: &mut Model, curve: ogeom_geom::Curve) -> Shape {
    let domain = curve.domain();
    ogeom_algo::make_edge(model, curve, domain, T)
        .unwrap()
        .shape
}

/// A cubic spline spine riding `periods` periods of a sine of `amplitude`
/// along `length` of x: its control points on the sine, twenty a period,
/// on uniform knots. Also returns the spine densely sampled.
fn sine_spine(model: &mut Model, amplitude: f64, length: f64, periods: f64) -> (Shape, Vec<Point>) {
    let wave = |x: f64| Point::new(x, amplitude * (2.0 * PI * periods * x / length).sin(), 0.0);
    #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
    let n = (20.0 * periods) as usize;
    let control: Vec<Point> = (0..=n).map(|k| wave(across((0.0, length), k, n))).collect();
    let mut knots = vec![0.0; 4];
    knots.extend((1..n - 2).map(|k| across((0.0, 1.0), k, n - 2)));
    knots.extend([1.0; 4]);
    let curve = ogeom_geom::Curve::BSpline(
        ogeom_geom::BSplineCurve::new(ogeom_math::KnotVector::new(knots, 3).unwrap(), control, T)
            .unwrap(),
    );
    let (lo, hi) = curve.domain();
    let dense: Vec<Point> = (0..=20_000)
        .map(|k| curve.point_at(across((lo, hi), k, 20_000), T).unwrap())
        .collect();
    (edge_of(model, curve), dense)
}

#[test]
fn a_skinned_pipe_along_a_wavy_spine_holds_its_radius_between_stations() {
    let mut model = Model::new();
    let (spine, dense) = sine_spine(&mut model, 1.0, 40.0, 8.0);
    let (r, tolerance) = (0.3, 1e-3);
    let pipe = ogeom_offset::make_pipe_skinned(&mut model, &spine, r, tolerance, T)
        .unwrap()
        .shape;
    let diagnosis = ogeom_algo::check(&model, &pipe, T).unwrap();
    assert!(diagnosis.is_valid(), "{diagnosis}");
    let walls = spline_surfaces(&model, &pipe);
    assert_eq!(walls.len(), 1);
    let wall = &walls[0];
    let (ud, vd) = wall.domain();
    let mut worst = 0.0_f64;
    for j in 0..=400 {
        for i in 0..=40 {
            let p = wall
                .point_at(across(ud, i, 40), across(vd, j, 400), T)
                .unwrap();
            // The spine runs along x: only its stretch within reach counts.
            let lo = dense.partition_point(|q| q.x < p.x - 1.0);
            let hi = dense
                .partition_point(|q| q.x <= p.x + 1.0)
                .min(dense.len() - 1);
            worst = worst.max((to_polyline(p, &dense[lo.saturating_sub(1)..=hi]) - r).abs());
        }
    }
    eprintln!("pipe wall off its radius by {worst}");
    assert!(worst <= tolerance, "the wall strays {worst} from the tube");
}

/// The distance from `p` to the circle of `radius` about the origin in the
/// XY plane.
fn to_flat_circle(p: Point, radius: f64) -> f64 {
    (p.x.hypot(p.y) - radius).hypot(p.z)
}

#[test]
fn a_filling_holds_its_curved_border_between_samples() {
    // A half disc's arc below the x axis, closed by three lines through a
    // height of one.
    let mut model = Model::new();
    let arc = ogeom_geom::CircleCurve::new(Circle::new(Frame::WORLD, 1.0, T).unwrap());
    let corners = [
        Point::new(1.0, 0.0, 0.0),
        Point::new(1.0, 0.0, 1.0),
        Point::new(-1.0, 0.0, 1.0),
        Point::new(-1.0, 0.0, 0.0),
    ];
    let bottom = ogeom_algo::make_edge(&mut model, arc.into(), (PI, 2.0 * PI), T)
        .unwrap()
        .shape;
    let line = |model: &mut Model, a: Point, b: Point| {
        ogeom_algo::make_edge(
            model,
            ogeom_geom::LineCurve::segment(a, b, T).unwrap().into(),
            (0.0, a.distance(b)),
            T,
        )
        .unwrap()
        .shape
    };
    let right = line(&mut model, corners[0], corners[1]);
    let top = line(&mut model, corners[1], corners[2]);
    let left = line(&mut model, corners[2], corners[3]);
    let tolerance = 1e-4;
    let filled =
        ogeom_offset::make_filling(&mut model, &[bottom, right, top, left], 4, tolerance, T)
            .unwrap()
            .shape;
    let patches = spline_surfaces(&model, &filled);
    let patch = &patches[0];
    let (ud, vd) = patch.domain();
    let mut worst = 0.0_f64;
    for i in 0..=2000 {
        let p = patch.point_at(across(ud, i, 2000), vd.0, T).unwrap();
        worst = worst.max(to_flat_circle(p, 1.0));
    }
    eprintln!("filling border off its arc by {worst}");
    assert!(worst <= tolerance, "the border strays {worst} from its arc");
}

#[test]
fn a_skinned_loft_through_thin_ellipses_holds_its_end_sections() {
    let mut model = Model::new();
    let (a, b) = (10.0, 0.2);
    let ellipse = |z: f64| {
        let frame = Frame::new(Point::new(0.0, 0.0, z), Direction::Z, Direction::X, T).unwrap();
        ogeom_geom::EllipseCurve::new(ogeom_math::Ellipse::new(frame, a, b, T).unwrap())
    };
    let sections: Vec<Shape> = [0.0, 5.0, 10.0]
        .iter()
        .map(|&z| {
            let edge = edge_of(&mut model, ellipse(z).into());
            ogeom_algo::make_wire(&mut model, &[edge], T).unwrap().shape
        })
        .collect();
    let tolerance = 1e-4;
    let loft = ogeom_offset::make_loft_skinned(&mut model, &sections, tolerance, T)
        .unwrap()
        .shape;
    let diagnosis = ogeom_algo::check(&model, &loft, T).unwrap();
    assert!(diagnosis.is_valid(), "{diagnosis}");
    // The ellipse at its height, densely enough that the polyline's sag is
    // far below the tolerance.
    let ring = |z: f64| -> Vec<Point> {
        let curve = ellipse(z);
        (0..=40_000)
            .map(|k| {
                curve
                    .point_at(across((0.0, 2.0 * PI), k, 40_000), T)
                    .unwrap()
            })
            .collect()
    };
    let rings = [ring(0.0), ring(10.0)];
    let walls = spline_surfaces(&model, &loft);
    let wall = &walls[0];
    let (ud, vd) = wall.domain();
    let mut worst = 0.0_f64;
    for v in [vd.0, vd.1] {
        for i in 0..=1000 {
            let p = wall.point_at(across(ud, i, 1000), v, T).unwrap();
            let line = if p.z < 5.0 { &rings[0] } else { &rings[1] };
            worst = worst.max(to_polyline(p, line));
        }
    }
    eprintln!("loft end rings off their ellipses by {worst}");
    assert!(
        worst <= tolerance,
        "an end ring strays {worst} from its ellipse"
    );
}

/// A gently curved spline spine from the origin, setting off along +z.
fn bent_spine(model: &mut Model) -> Shape {
    let curve = ogeom_geom::Curve::BSpline(
        ogeom_geom::BSplineCurve::new(
            ogeom_math::KnotVector::new(vec![0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0], 3).unwrap(),
            vec![
                Point::new(0.0, 0.0, 0.0),
                Point::new(0.0, 0.0, 3.0),
                Point::new(1.0, 0.0, 6.0),
                Point::new(3.0, 0.0, 9.0),
            ],
            T,
        )
        .unwrap(),
    );
    let edge = edge_of(model, curve);
    ogeom_algo::make_wire(model, &[edge], T).unwrap().shape
}

#[test]
fn a_pipe_shell_holds_its_profile_arc_between_samples() {
    // A half disc in the XY plane: the arc of radius one over +y, closed
    // by its diameter.
    let mut model = Model::new();
    let (a, b) = (Point::new(-1.0, 0.0, 0.0), Point::new(1.0, 0.0, 0.0));
    let (va, vb) = (
        ogeom_algo::make_vertex(&mut model, a).shape,
        ogeom_algo::make_vertex(&mut model, b).shape,
    );
    let arc = ogeom_algo::make_edge_between(
        &mut model,
        ogeom_geom::CircleCurve::new(Circle::new(Frame::WORLD, 1.0, T).unwrap()).into(),
        (0.0, PI),
        &vb,
        &va,
        T,
    )
    .unwrap()
    .shape;
    let diameter = ogeom_algo::make_edge_between(
        &mut model,
        ogeom_geom::LineCurve::segment(a, b, T).unwrap().into(),
        (0.0, 2.0),
        &va,
        &vb,
        T,
    )
    .unwrap()
    .shape;
    let wire = ogeom_algo::make_wire(&mut model, &[arc, diameter], T)
        .unwrap()
        .shape;
    let plane: SurfaceGeometry = ogeom_geom::PlaneSurface::new(Plane::new(Frame::WORLD)).into();
    let profile = ogeom_algo::make_face(&mut model, plane, &[wire], T)
        .unwrap()
        .shape;
    let spine = bent_spine(&mut model);
    let tolerance = 1e-5;
    let pipe = ogeom_offset::make_pipe_shell(&mut model, &profile, &spine, false, tolerance, T)
        .unwrap()
        .shape;
    let diagnosis = ogeom_algo::check(&model, &pipe, T).unwrap();
    assert!(diagnosis.is_valid(), "{diagnosis}");
    // Every spline border lying in the profile's plane off its diameter is
    // the arc.
    let mut worst = 0.0_f64;
    let mut seen = 0;
    for patch in spline_surfaces(&model, &pipe) {
        let (ud, vd) = patch.domain();
        for k in 0..=1000 {
            for (u, v) in [
                (across(ud, k, 1000), vd.0),
                (across(ud, k, 1000), vd.1),
                (ud.0, across(vd, k, 1000)),
                (ud.1, across(vd, k, 1000)),
            ] {
                let p = patch.point_at(u, v, T).unwrap();
                if p.z.abs() < 1e-3 && p.y > 1e-2 {
                    seen += 1;
                    worst = worst.max(to_flat_circle(p, 1.0));
                }
            }
        }
    }
    assert!(seen > 500, "the arc's strip has a border on the profile");
    eprintln!("pipe shell border off its arc by {worst}");
    assert!(worst <= tolerance, "the border strays {worst} from its arc");
}

#[test]
fn a_projected_circle_holds_its_stated_tolerance_between_stations() {
    // A circle of radius 3 at height 10 over a ball of radius 5: its foot
    // on the ball is the circle of radius 15 / sqrt(109) at height
    // 50 / sqrt(109).
    let mut model = Model::new();
    let ball = ogeom_algo::make_sphere(&mut model, Frame::WORLD, 5.0, T)
        .unwrap()
        .shape;
    let frame = Frame::new(Point::new(0.0, 0.0, 10.0), Direction::Z, Direction::X, T).unwrap();
    let circle = edge_of(
        &mut model,
        ogeom_geom::CircleCurve::new(Circle::new(frame, 3.0, T).unwrap()).into(),
    );
    let wire = ogeom_algo::make_wire(&mut model, &[circle], T)
        .unwrap()
        .shape;
    let tolerance = 1e-3;
    let (landed, _) =
        ogeom_offset::normal_projection(&mut model, &ball, &wire, 8, tolerance, T).unwrap();
    assert!(!landed.is_empty());
    let (radius, height) = (15.0 / 109.0_f64.sqrt(), 50.0 / 109.0_f64.sqrt());
    for stretch in &landed {
        let data = model.node(&stretch.edge).unwrap().data().as_edge().unwrap();
        let EdgeRepr::Curve3d { curve, range, .. } = data.curve3d().unwrap() else {
            unreachable!()
        };
        let geometry = model.geometry().curve(*curve).unwrap();
        let mut worst = 0.0_f64;
        for k in 0..=2000 {
            let p = geometry.point_at(across(*range, k, 2000), T).unwrap();
            worst = worst.max((p.x.hypot(p.y) - radius).hypot(p.z - height));
        }
        eprintln!(
            "projection off its foot by {worst}, stating {}",
            stretch.tolerance
        );
        assert!(
            worst <= stretch.tolerance.max(1e-9) && stretch.tolerance <= tolerance,
            "the projection strays {worst} and states {}",
            stretch.tolerance
        );
    }
    let diagnosis = ogeom_algo::check(&model, &landed[0].edge, T).unwrap();
    assert!(diagnosis.is_valid(), "{diagnosis}");
}

#[test]
fn a_pipe_through_thin_sections_holds_its_end_section() {
    // An ellipse 10 by 0.2 square to a quarter arc of radius 20 at each of
    // its ends, the arc turning in the ellipse's long direction.
    let mut model = Model::new();
    let ellipse = |model: &mut Model, frame: Frame| {
        let curve =
            ogeom_geom::EllipseCurve::new(ogeom_math::Ellipse::new(frame, 10.0, 0.2, T).unwrap());
        let edge = edge_of(model, curve.into());
        ogeom_algo::make_wire(model, &[edge], T).unwrap().shape
    };
    let start = ellipse(
        &mut model,
        Frame::new(Point::ORIGIN, -Direction::Y, Direction::X, T).unwrap(),
    );
    let end = ellipse(
        &mut model,
        Frame::new(Point::new(20.0, -20.0, 0.0), Direction::X, -Direction::Y, T).unwrap(),
    );
    let arc = Circle::new(
        Frame::new(Point::new(20.0, 0.0, 0.0), Direction::Z, Direction::X, T).unwrap(),
        20.0,
        T,
    )
    .unwrap();
    let spine = ogeom_algo::make_edge(
        &mut model,
        ogeom_geom::CircleCurve::new(arc).into(),
        (PI, 1.5 * PI),
        T,
    )
    .unwrap()
    .shape;
    let spine = ogeom_algo::make_wire(&mut model, &[spine], T)
        .unwrap()
        .shape;
    let tolerance = 1e-3;
    let pipe =
        ogeom_offset::make_pipe_sections(&mut model, &[start, end], &spine, false, tolerance, T)
            .unwrap()
            .shape;
    let diagnosis = ogeom_algo::check(&model, &pipe, T).unwrap();
    assert!(diagnosis.is_valid(), "{diagnosis}");
    let ring: Vec<Point> = (0..=40_000)
        .map(|k| {
            let a = 2.0 * PI * f64::from(k) / 40_000.0;
            Point::new(10.0 * a.cos(), 0.0, 0.2 * a.sin())
        })
        .collect();
    let mut worst = 0.0_f64;
    let mut seen = 0;
    for patch in spline_surfaces(&model, &pipe) {
        let (ud, vd) = patch.domain();
        for k in 0..=2000 {
            for (u, v) in [
                (across(ud, k, 2000), vd.0),
                (across(ud, k, 2000), vd.1),
                (ud.0, across(vd, k, 2000)),
                (ud.1, across(vd, k, 2000)),
            ] {
                let p = patch.point_at(u, v, T).unwrap();
                if p.y.abs() < 1e-6 {
                    seen += 1;
                    worst = worst.max(to_polyline(p, &ring));
                }
            }
        }
    }
    assert!(seen > 1000, "the pipe has a border on the start section");
    eprintln!("pipe sections' start ring off its ellipse by {worst}");
    assert!(worst <= tolerance, "the start ring strays {worst}");
}

/// The distance from `p` to `curve` over `domain`: the nearest of `coarse`
/// samples of it (parameter and point), then a golden-section search on the
/// samples either side.
fn to_curve(p: Point, curve: &ogeom_geom::Curve, coarse: &[(f64, Point)]) -> f64 {
    let k = (0..coarse.len())
        .min_by(|a, b| {
            p.distance(coarse[*a].1)
                .total_cmp(&p.distance(coarse[*b].1))
        })
        .unwrap();
    let (mut lo, mut hi) = (
        coarse[k.saturating_sub(1)].0,
        coarse[(k + 1).min(coarse.len() - 1)].0,
    );
    let off = |t: f64| p.distance(curve.point_at(t, T).unwrap());
    let g = (5.0_f64.sqrt() - 1.0) / 2.0;
    for _ in 0..80 {
        let (a, b) = (hi - g * (hi - lo), lo + g * (hi - lo));
        if off(a) < off(b) {
            hi = b;
        } else {
            lo = a;
        }
    }
    off(f64::midpoint(lo, hi))
}

/// A disc of radius `r` about `frame`'s origin, square to its `z`.
fn disc(model: &mut Model, frame: Frame, r: f64) -> Shape {
    let circle = edge_of(
        model,
        ogeom_geom::CircleCurve::new(Circle::new(frame, r, T).unwrap()).into(),
    );
    let wire = ogeom_algo::make_wire(model, &[circle], T).unwrap().shape;
    let plane: SurfaceGeometry = ogeom_geom::PlaneSurface::new(Plane::new(frame)).into();
    ogeom_algo::make_face(model, plane, &[wire], T)
        .unwrap()
        .shape
}

/// How far the spline walls of a disc of radius `r` swept along `spine`
/// stand off the tube of that radius round it, sampled densely.
fn off_the_tube(model: &Model, pipe: &Shape, spine: &ogeom_geom::Curve, r: f64) -> f64 {
    let domain = spine.domain();
    let coarse: Vec<(f64, Point)> = (0..=2000)
        .map(|k| {
            let t = across(domain, k, 2000);
            (t, spine.point_at(t, T).unwrap())
        })
        .collect();
    let mut worst = 0.0_f64;
    let walls = spline_surfaces(model, pipe);
    assert!(!walls.is_empty(), "the pipe has a fitted wall");
    for wall in walls {
        let (ud, vd) = wall.domain();
        for j in 0..=400 {
            for i in 0..=40 {
                let p = wall
                    .point_at(across(ud, i, 40), across(vd, j, 400), T)
                    .unwrap();
                worst = worst.max((to_curve(p, spine, &coarse) - r).abs());
            }
        }
    }
    worst
}

#[test]
fn a_pipe_shell_round_a_tight_bend_holds_its_radius_between_stations() {
    let spine = tight_bend();
    let (r, tolerance) = (0.3, 1e-4);
    for frenet in [false, true] {
        let mut model = Model::new();
        let edge = edge_of(&mut model, spine.clone());
        let wire = ogeom_algo::make_wire(&mut model, &[edge], T).unwrap().shape;
        let profile = disc(&mut model, Frame::WORLD, r);
        let pipe = ogeom_offset::make_pipe_shell(&mut model, &profile, &wire, frenet, tolerance, T)
            .unwrap()
            .shape;
        let diagnosis = ogeom_algo::check(&model, &pipe, T).unwrap();
        assert!(diagnosis.is_valid(), "{diagnosis}");
        let worst = off_the_tube(&model, &pipe, &spine, r);
        eprintln!("pipe shell round a bend (Frenet {frenet}) off its radius by {worst}");
        assert!(worst <= tolerance, "the wall strays {worst} from the tube");
    }
}

#[test]
fn a_pipe_shell_under_a_law_holds_its_radius_between_stations() {
    use ogeom_geom::Transformable as _;
    let spine = tight_bend();
    let r = 0.3;
    for (guided, tolerance) in [(true, 1e-4), (false, 1e-5)] {
        let mut model = Model::new();
        let edge = edge_of(&mut model, spine.clone());
        let wire = ogeom_algo::make_wire(&mut model, &[edge], T).unwrap().shape;
        // The spine moved one along y: it crosses every plane square to
        // the spine one along y from the spine.
        let beside = spine
            .transformed(
                &ogeom_math::Transform::translation(ogeom_math::Vector::new(0.0, 1.0, 0.0)),
                T,
            )
            .unwrap();
        let guide = edge_of(&mut model, beside);
        let law = if guided {
            ogeom_offset::PipeLaw::Auxiliary { guide: &guide }
        } else {
            ogeom_offset::PipeLaw::Binormal(Direction::Y)
        };
        let profile = disc(&mut model, Frame::WORLD, r);
        let pipe =
            ogeom_offset::make_pipe_shell_law(&mut model, &profile, &wire, law, tolerance, T)
                .unwrap()
                .shape;
        let diagnosis = ogeom_algo::check(&model, &pipe, T).unwrap();
        assert!(diagnosis.is_valid(), "{diagnosis}");
        let worst = off_the_tube(&model, &pipe, &spine, r);
        eprintln!("pipe shell under a law (guided {guided}) off its radius by {worst}");
        assert!(worst <= tolerance, "the wall strays {worst} from the tube");
    }
}

/// A cubic spline running up z, turning through a quarter within one span
/// of its parameter and running on along x: a pipe's stations stand evenly
/// in that parameter, a few across the bend.
fn tight_bend() -> ogeom_geom::Curve {
    let control = vec![
        Point::new(0.0, 0.0, 0.0),
        Point::new(0.0, 0.0, 3.0),
        Point::new(0.0, 0.0, 6.0),
        Point::new(0.0, 0.0, 9.0),
        Point::new(1.0, 0.0, 9.0),
        Point::new(4.0, 0.0, 9.0),
        Point::new(7.0, 0.0, 9.0),
        Point::new(10.0, 0.0, 9.0),
    ];
    ogeom_geom::Curve::BSpline(
        ogeom_geom::BSplineCurve::new(
            ogeom_math::KnotVector::new(
                vec![0.0, 0.0, 0.0, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0],
                3,
            )
            .unwrap(),
            control,
            T,
        )
        .unwrap(),
    )
}

#[test]
fn a_closed_pipe_shell_holds_its_radius_between_stations() {
    // An ellipse 5 by 2: its stations stand evenly in its angle, sparsest
    // in length round the tight ends.
    let spine: ogeom_geom::Curve =
        ogeom_geom::EllipseCurve::new(ogeom_math::Ellipse::new(Frame::WORLD, 5.0, 2.0, T).unwrap())
            .into();
    let (r, tolerance) = (0.3, 1e-4);
    let mut model = Model::new();
    let edge = edge_of(&mut model, spine.clone());
    let wire = ogeom_algo::make_wire(&mut model, &[edge], T).unwrap().shape;
    let start = Frame::new(Point::new(5.0, 0.0, 0.0), Direction::Y, Direction::X, T).unwrap();
    let profile = disc(&mut model, start, r);
    let pipe = ogeom_offset::make_pipe_shell(&mut model, &profile, &wire, false, tolerance, T)
        .unwrap()
        .shape;
    let diagnosis = ogeom_algo::check(&model, &pipe, T).unwrap();
    assert!(diagnosis.is_valid(), "{diagnosis}");
    let worst = off_the_tube(&model, &pipe, &spine, r);
    eprintln!("closed pipe shell off its radius by {worst}");
    assert!(worst <= tolerance, "the wall strays {worst} from the tube");
}

#[test]
fn a_wide_helical_sweep_holds_its_profile_between_stations() {
    // A circle of radius one, a thousand out from the z axis in the XZ
    // plane, screwed once round at a pitch of three: its stations, 48 a
    // quarter turn, stand about thirty apart.
    let (out, r, pitch) = (1000.0, 1.0, 3.0);
    let mut model = Model::new();
    let frame = Frame::new(Point::new(out, 0.0, 0.0), -Direction::Y, Direction::X, T).unwrap();
    let profile = disc(&mut model, frame, r);
    let axis = ogeom_math::Axis {
        location: Point::ORIGIN,
        direction: Direction::Z,
    };
    let sweep =
        ogeom_offset::make_helical_sweep(&mut model, &profile, axis, pitch, 1.0, false, 0.0, T)
            .unwrap()
            .shape;
    let diagnosis = ogeom_algo::check(&model, &sweep, T).unwrap();
    assert!(diagnosis.is_valid(), "{diagnosis}");
    // Each wall point read in the half plane through the axis it stands
    // in, screwed back to the start by the turn it has made (the one of
    // its angle's whole turns that lands it nearest the profile): there it
    // is on the profile's circle.
    let mut worst = 0.0_f64;
    for wall in spline_surfaces(&model, &sweep) {
        let (ud, vd) = wall.domain();
        for j in 0..=200 {
            for i in 0..=40 {
                let q = wall
                    .point_at(across(ud, i, 40), across(vd, j, 200), T)
                    .unwrap();
                let angle = q.y.atan2(q.x);
                let off = [angle, angle + 2.0 * PI]
                    .iter()
                    .map(|turned| {
                        let back = q.z - pitch * turned / (2.0 * PI);
                        ((q.x.hypot(q.y) - out).hypot(back) - r).abs()
                    })
                    .fold(f64::INFINITY, f64::min);
                worst = worst.max(off);
            }
        }
    }
    // The sweep's own target.
    let tolerance = T.confusion() * 100.0;
    eprintln!("wide helical sweep off its profile by {worst}");
    assert!(worst <= tolerance, "the wall strays {worst} from the screw");
}

/// The distance from `p` to the triangle `a`, `b`, `c`.
fn to_triangle(p: Point, a: Point, b: Point, c: Point) -> f64 {
    let n = (b - a).cross(c - a);
    let m = n.magnitude();
    if m > 0.0 {
        let n = n / m;
        let h = (p - a).dot(n);
        let f = p - n * h;
        if (b - a).cross(f - a).dot(n) >= 0.0
            && (c - b).cross(f - b).dot(n) >= 0.0
            && (a - c).cross(f - c).dot(n) >= 0.0
        {
            return h.abs();
        }
    }
    to_polyline(p, &[a, b, c, a])
}

#[test]
fn a_skinned_loft_closes_a_wavy_end_on_its_cone() {
    // Rings of radius 10 waving 0.8 in and out twelve times round, the
    // last also rising and falling 2 twice round: no plane caps it, so
    // the loft closes it with the cone from the ring to its centroid.
    let ring = |z: f64, lift: f64| {
        move |t: f64| {
            let a = 2.0 * PI * t;
            let r = 10.0 + 0.8 * (12.0 * a).sin();
            Point::new(r * a.cos(), r * a.sin(), z + lift * (2.0 * a).sin())
        }
    };
    let mut model = Model::new();
    let ts: Vec<f64> = (0..=800).map(|k| across((0.0, 1.0), k, 800)).collect();
    let mut sections = Vec::new();
    let mut end = None;
    for (z, lift) in [(0.0, 0.0), (5.0, 0.0), (10.0, 2.0)] {
        let at = ring(z, lift);
        let curve: ogeom_geom::Curve =
            ogeom_geom::fit::fit_curve_sampled(|t| Ok(at(t)), &ts, true, 3, 1e-7, T)
                .unwrap()
                .curve
                .into();
        end = Some(curve.clone());
        let edge = edge_of(&mut model, curve);
        sections.push(ogeom_algo::make_wire(&mut model, &[edge], T).unwrap().shape);
    }
    let tolerance = 1e-4;
    let loft = ogeom_offset::make_loft_skinned(&mut model, &sections, tolerance, T)
        .unwrap()
        .shape;
    let diagnosis = ogeom_algo::check(&model, &loft, T).unwrap();
    assert!(diagnosis.is_valid(), "{diagnosis}");
    // The end ring densely: its chords sag well under the tolerance.
    let end = end.unwrap();
    let domain = end.domain();
    let count = 8000;
    let rim: Vec<Point> = (0..=count)
        .map(|k| end.point_at(across(domain, k, count), T).unwrap())
        .collect();
    let mut worst = 0.0_f64;
    let mut seen = 0;
    for patch in spline_surfaces(&model, &loft) {
        let (ud, vd) = patch.domain();
        let apex = patch.point_at(ud.0, vd.1, T).unwrap();
        if apex.distance(patch.point_at(f64::midpoint(ud.0, ud.1), vd.1, T).unwrap()) > 1e-9 {
            continue;
        }
        seen += 1;
        for j in 0..=20 {
            for i in 0..=400 {
                let q = patch
                    .point_at(across(ud, i, 400), across(vd, j, 20), T)
                    .unwrap();
                // The fan from the apex to the rim near the angle q stands
                // at: the rim runs round with its angle.
                let turn = q.y.atan2(q.x).rem_euclid(2.0 * PI) / (2.0 * PI);
                #[allow(
                    clippy::cast_possible_truncation,
                    clippy::cast_sign_loss,
                    clippy::cast_precision_loss
                )]
                let at = (turn * count as f64) as usize;
                let mut best = f64::INFINITY;
                for k in (at + count - 100)..(at + count + 100) {
                    let (a, b) = (rim[k % count], rim[(k + 1) % count]);
                    best = best.min(to_triangle(q, apex, a, b));
                }
                worst = worst.max(best);
            }
        }
    }
    assert_eq!(seen, 1, "one end closes on a cone");
    eprintln!("wavy end's cap off its cone by {worst}");
    assert!(worst <= tolerance, "the cap strays {worst} from its cone");
}

/// A quarter arc of radius `r` about the origin from `(r, 0)` to `(0, r)`,
/// then a straight leg of 20 up `+y`: a corner between a curved run and a
/// straight one, turning in the arc's plane.
fn arc_then_leg(model: &mut Model, r: f64) -> Shape {
    let (a, b, c) = (
        Point::new(r, 0.0, 0.0),
        Point::new(0.0, r, 0.0),
        Point::new(0.0, r + 20.0, 0.0),
    );
    let va = ogeom_algo::make_vertex(model, a).shape;
    let vb = ogeom_algo::make_vertex(model, b).shape;
    let vc = ogeom_algo::make_vertex(model, c).shape;
    let arc = ogeom_algo::make_edge_between(
        model,
        ogeom_geom::CircleCurve::new(Circle::new(Frame::WORLD, r, T).unwrap()).into(),
        (0.0, PI / 2.0),
        &va,
        &vb,
        T,
    )
    .unwrap()
    .shape;
    let leg = ogeom_algo::make_edge_between(
        model,
        ogeom_geom::LineCurve::segment(b, c, T).unwrap().into(),
        (0.0, 20.0),
        &vb,
        &vc,
        T,
    )
    .unwrap()
    .shape;
    ogeom_algo::make_wire(model, &[arc, leg], T).unwrap().shape
}

#[test]
fn a_pipe_shell_holds_its_walls_between_rows_at_a_curved_corner() {
    // The arc's wall runs up to where each generator crosses the leg's,
    // past the arc's end along its straight extension on the outside of the
    // turn; the leg's runs back to the same crossing. Each wall point lies
    // on the arc's tube (its torus, then the extension's cylinder along -x
    // past x = 0) or on the leg's cylinder about the y axis.
    let (r, w, tolerance) = (20.0, 2.0, 1e-4);
    let off_round = |p: Point| -> f64 {
        let arc = if p.x >= 0.0 {
            (p.x.hypot(p.y) - r).hypot(p.z)
        } else {
            (p.y - r).hypot(p.z)
        };
        (arc - w).abs().min((p.x.hypot(p.z) - w).abs())
    };
    // The square's walls: the larger of the two offsets across the
    // section from its centre, against the half side.
    let off_square = |p: Point| -> f64 {
        let arc = if p.x >= 0.0 {
            (p.x.hypot(p.y) - r).abs().max(p.z.abs())
        } else {
            (p.y - r).abs().max(p.z.abs())
        };
        (arc - w).abs().min((p.x.abs().max(p.z.abs()) - w).abs())
    };
    for round in [true, false] {
        let mut model = Model::new();
        let spine = arc_then_leg(&mut model, r);
        let start = Frame::new(Point::new(r, 0.0, 0.0), Direction::Y, Direction::X, T).unwrap();
        let profile = if round {
            disc(&mut model, start, w)
        } else {
            let corners: Vec<Point> = [(-w, -w), (w, -w), (w, w), (-w, w)]
                .iter()
                .map(|(a, b)| Point::new(r + a, 0.0, *b))
                .collect();
            let wire = ogeom_algo::make_polygon(&mut model, &corners, true, T)
                .unwrap()
                .shape;
            let plane: SurfaceGeometry = ogeom_geom::PlaneSurface::new(Plane::new(start)).into();
            ogeom_algo::make_face(&mut model, plane, &[wire], T)
                .unwrap()
                .shape
        };
        let pipe = ogeom_offset::make_pipe_shell(&mut model, &profile, &spine, false, tolerance, T)
            .unwrap()
            .shape;
        let diagnosis = ogeom_algo::check(&model, &pipe, T).unwrap();
        assert!(diagnosis.is_valid(), "{diagnosis}");
        let mut worst = 0.0_f64;
        for wall in spline_surfaces(&model, &pipe) {
            let (ud, vd) = wall.domain();
            for j in 0..=400 {
                for i in 0..=40 {
                    let p = wall
                        .point_at(across(ud, i, 40), across(vd, j, 400), T)
                        .unwrap();
                    worst = worst.max(if round { off_round(p) } else { off_square(p) });
                }
            }
        }
        eprintln!("pipe shell at a curved corner (round {round}) off its walls by {worst}");
        assert!(worst <= tolerance, "a wall strays {worst} at the corner");
    }
}