hpr_design/

fins.rs

//! Fin sets and tube fins: planforms, cross-sections, tabs, and mass properties from geometry.
//!
//! **Fin frame.** A fin set's frame has its origin on the body axis at the station of the root
//! chord's leading edge. Fin `k` of `N` lies in the plane through the axis at roll angle
//! `φ_k = φ_0 + 2πk/N` from `x_B` toward `y_B`. In a fin's own coordinates, `x` runs aft from the
//! root leading edge, the span `h` runs outward from the body surface (radius `R_b + h`), and the
//! thickness coordinate `τ` is normal to the fin plane. Planform definitions follow S. Niskanen,
//! *OpenRocket technical documentation* v13.05, §3.2.2, pp. 25–29: a trapezoidal fin has root chord
//! `c_r`, tip chord `c_t` parallel to the body, span `s`, and sweep `x_t`, the axial distance from
//! the root leading edge to the tip leading edge; an elliptical fin has chord
//! `c(h) = c_r √(1 − (h/s)²)` centered on the root chord; a freeform fin is a polygon.
//!
//! **Mass.** Each fin is a plate whose local thickness varies along each chord with the
//! cross-section. At span `h`, a chord from `a` to `b` (length `c`) contributes the moments
//! `M_k = ∫ x^k t(x) dx` and `T = ∫ t(x)³/12 dx`:
//!
//! ```text
//! square:   t(x) = t
//! rounded:  semicircular leading and trailing edges of radius a_r = min(t, c)/2
//! airfoil:  t(x) = 10 t P(ξ),  P = 0.2969√ξ − 0.1260ξ − 0.3516ξ² + 0.2843ξ³ − 0.1015ξ⁴
//! ```
//!
//! The airfoil is the NACA four-digit symmetric thickness distribution (I. H. Abbott and A. E.
//! von Doenhoff, *Theory of Wing Sections*, Dover, 1959, eq. 6.2) with its maximum thickness,
//! `1.0003 t` at `ξ = 0.2998`, scaled to the fin thickness. OpenRocket's documentation uses the
//! cross-section only for drag (§3.4.4, pp. 49–50); hpr-design also counts the volume it removes.
//! The fin's moments about the fin-set frame are then
//!
//! ```text
//! m = ρ ∫ M_0 dh,  ∫ r dm = ρ ∫ (R_b + h) M_0 dh,  ∫ r² dm = ρ ∫ (R_b + h)² M_0 dh,
//! ∫ x dm = ρ ∫ M_1 dh,  ∫ x² dm = ρ ∫ M_2 dh,  ∫ r x dm = ρ ∫ (R_b + h) M_1 dh,  ∫ τ² dm = ρ ∫ T dh
//! ```
//!
//! and with the fin at `φ = 0` (points at `(r, τ, −x)` in body axes)
//!
//! ```text
//! I_xx = ∫(τ² + x²) dm,  I_yy = ∫(r² + x²) dm,  I_zz = ∫(r² + τ²) dm,  I_xz = ∫ r x dm,  I_xy = I_yz = 0
//! ```
//!
//! about the origin. A tab is a square-section slab below the root (`−h_tab ≤ h ≤ 0`). The flat
//! root is taken to sit on the body at radius `R_b`; the sliver between a flat root and the curved
//! tube, `t²/8R_b` deep, is ignored. On a nose cone or a transition the surface is not level
//! along the root: `R_b` is the body's radius at the root leading edge, `h` is measured from it,
//! and a freeform fin's root runs through its root points on the surface (OpenRocket's reading of
//! a fin on a nose cone, [ADR-166][adr-166]). **Cant** `δ` turns each fin (with its tab) by `δ` about its own
//! outward span axis through the root mid-chord, right-handed, so a positive cant turns fin 0's
//! leading edge toward `−y_B`. See `docs/physics/mass.md`.
//!
//! **Fillets.** A fillet of radius `r` runs the root chord `c_r` on each face of each fin. Its
//! section is the region between the fin's mid-plane, the body's circle of radius `R_b`, and a
//! circle of radius `r` tangent to both: OpenRocket 24.12's reading, which leaves the fin's
//! thickness out. Its mass and center of mass match OpenRocket's to 1e-15 on nine probe designs
//! ([ADR-096][adr-096], the decision on fillets). With the fillet circle's center at
//! `(c, r)`, `c = √(R_b² + 2 R_b r)` and `θ = atan(r / c)`, the section is the triangle
//! `(0, 0), (c, 0), (c, r)` less the body's sector of angle `θ` and the fillet's of `π/2 − θ`:
//!
//! ```text
//! A = c r / 2 − R_b² θ / 2 − r² (π/2 − θ) / 2
//! ```
//!
//! which tends to `r² (1 − π/4)` on a flat body. Its moments come the same way, and the fillets
//! are a prism of that section along the root, in the fillet's own material.
//!
//! [adr-166]: https://github.com/nrdptel/hpr-sim/blob/main/docs/decisions/0166-a-fin-root-that-follows-the-body.md
//! [adr-096]: https://github.com/nrdptel/hpr-sim/blob/main/docs/DECISIONS.md#adr-096-fin-fillets-and-an-automatic-radius-inside-a-nose-cone-read-as-openrocket-reads-them-2026-09-28

use std::f64::consts::PI;

use hpr_core::quadrature::{Tolerance, integrate};
use hpr_core::{DMat3, DQuat, DVec3};
use serde::{Deserialize, Serialize};

use crate::error::DesignError;
use crate::mass::MassProperties;
use crate::material::Material;
use crate::shapes::{Profile, check_dimension};

/// The outline of one fin.
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, schemars::JsonSchema)]
#[serde(tag = "kind", rename_all = "snake_case", deny_unknown_fields)]
#[non_exhaustive]
pub enum FinPlanform {
    /// A trapezoid with its tip chord parallel to the root.
    Trapezoidal {
        /// Root chord, m.
        root_chord_m: f64,
        /// Tip chord, m (zero for a pointed fin).
        tip_chord_m: f64,
        /// Span from the body surface to the tip, m.
        span_m: f64,
        /// Axial distance from the root leading edge aft to the tip leading edge, m.
        sweep_m: f64,
    },
    /// Half an ellipse on the root chord.
    Elliptical {
        /// Root chord, m.
        root_chord_m: f64,
        /// Span, m.
        span_m: f64,
    },
    /// A polygon given as `[x, h]` points from the root leading edge, which must be `[0, 0]`,
    /// around to the root trailing edge `[c_r, h_r]` with `c_r > 0` and `h_r ≥ 0`, closed along
    /// the root; `x` runs aft and `h` outward, in meters. On a body tube the root is level,
    /// `h_r = 0`; on a nose cone or a transition it follows the surface, through `root_m`.
    Freeform {
        /// The outline, m.
        points_m: Vec<[f64; 2]>,
        /// The root's points between its trailing and leading edges, aft to fore, m: none for a
        /// straight root. The root is straight between them.
        #[serde(default, skip_serializing_if = "Vec::is_empty")]
        root_m: Vec<[f64; 2]>,
    },
}

/// The shape of a fin's section along its chord.
#[derive(
    Debug, Clone, Copy, PartialEq, Eq, Hash, Default, Serialize, Deserialize, schemars::JsonSchema,
)]
#[serde(rename_all = "snake_case")]
#[non_exhaustive]
pub enum FinCrossSection {
    /// Constant thickness, square edges.
    #[default]
    Square,
    /// Semicircular leading and trailing edges.
    Rounded,
    /// The NACA four-digit symmetric thickness distribution.
    Airfoil,
}

/// A rectangular tab below a fin's root, reaching into the body.
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize, schemars::JsonSchema)]
#[serde(deny_unknown_fields)]
pub struct FinTab {
    /// Depth below the root, m.
    pub height_m: f64,
    /// Length along the root, m.
    pub length_m: f64,
    /// Distance from the root leading edge aft to the tab's leading edge, m.
    pub offset_m: f64,
}

/// Fillets along each fin's root: a concave joint on both faces, running the root chord.
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, schemars::JsonSchema)]
#[serde(deny_unknown_fields)]
pub struct FinFillet {
    /// Radius of the fillet's concave face, m.
    pub radius_m: f64,
    /// Material (bulk).
    pub material: Material,
}

/// A set of identical fins spaced evenly around the body.
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, schemars::JsonSchema)]
#[serde(deny_unknown_fields)]
pub struct FinSet {
    /// Number of fins, 1 to 64 (`FinSet::MAX_COUNT`).
    pub count: u32,
    /// Outline of each fin.
    pub planform: FinPlanform,
    /// Maximum thickness, m.
    pub thickness_m: f64,
    /// Section shape.
    #[serde(default)]
    pub cross_section: FinCrossSection,
    /// Optional tab below each root.
    #[serde(default)]
    pub tab: Option<FinTab>,
    /// Optional fillets along each root.
    #[serde(default, skip_serializing_if = "Option::is_none")]
    pub fillet: Option<FinFillet>,
    /// Cant angle, rad.
    #[serde(default)]
    pub cant_rad: f64,
    /// Roll angle of the first fin from `x_B` toward `y_B`, rad.
    #[serde(default)]
    pub base_angle_rad: f64,
    /// Material (bulk).
    pub material: Material,
}

/// Area and centroid of a planform.
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
pub struct PlanformGeometry {
    /// Area of one fin, m².
    pub area_m2: f64,
    /// Centroid, m aft of the root leading edge.
    pub centroid_x_m: f64,
    /// Centroid, m outward from the root.
    pub centroid_span_m: f64,
}

/// NACA four-digit thickness polynomial `P(ξ)` (Abbott and von Doenhoff, eq. 6.2); the mass code
/// uses its moments, and the tests integrate it directly.
#[cfg(test)]
fn naca_polynomial(xi: f64) -> f64 {
    0.2969 * xi.sqrt() - 0.1260 * xi - 0.3516 * xi * xi + 0.2843 * xi.powi(3) - 0.1015 * xi.powi(4)
}

/// `10 ∫₀¹ P dξ`, `10 ∫₀¹ ξ P dξ` and `10 ∫₀¹ ξ² P dξ`, term by term.
const AIRFOIL_PHI: [f64; 3] = [
    10.0 * (0.2969 * 2.0 / 3.0 - 0.1260 / 2.0 - 0.3516 / 3.0 + 0.2843 / 4.0 - 0.1015 / 5.0),
    10.0 * (0.2969 * 2.0 / 5.0 - 0.1260 / 3.0 - 0.3516 / 4.0 + 0.2843 / 5.0 - 0.1015 / 6.0),
    10.0 * (0.2969 * 2.0 / 7.0 - 0.1260 / 4.0 - 0.3516 / 5.0 + 0.2843 / 6.0 - 0.1015 / 7.0),
];

/// `1000 ∫₀¹ P³ dξ`, evaluated with mpmath at 30 digits; `airfoil_constants_are_the_integrals`
/// recomputes it.
const AIRFOIL_PSI: f64 = 0.472_889_488_894_451_523_708_489_652_762;

/// Chordwise moments `[M_0, M_1, M_2, T]` of a chord from `a` to `b` of thickness `t`.
fn chord_moments(section: FinCrossSection, a: f64, b: f64, t: f64) -> [f64; 4] {
    let c = b - a;
    if c <= 0.0 {
        return [0.0; 4];
    }
    match section {
        FinCrossSection::Square => [
            t * c,
            0.5 * t * (b * b - a * a),
            t * (b.powi(3) - a.powi(3)) / 3.0,
            t.powi(3) * c / 12.0,
        ],
        FinCrossSection::Airfoil => {
            let [p0, p1, p2] = AIRFOIL_PHI;
            [
                t * c * p0,
                t * (a * c * p0 + c * c * p1),
                t * (a * a * c * p0 + 2.0 * a * c * c * p1 + c.powi(3) * p2),
                t.powi(3) * c * AIRFOIL_PSI / 12.0,
            ]
        }
        FinCrossSection::Rounded => {
            // A slab of the middle thickness minus what the two rounded edges remove. With
            // v = a_r − u the depth from the edge's center, the removed thickness is
            // δ = 2a_r − 2√(a_r² − v²), whose moments from the edge are D_0, D_1, D_2.
            let mid = t.min(c);
            let r = 0.5 * mid;
            let d0 = r * r * (2.0 - PI / 2.0);
            let d1 = r * d0 - r.powi(3) / 3.0;
            let d2 = r * r * d0 - PI * r.powi(4) / 8.0;
            let e0 = 8.0 * r.powi(4) - 1.5 * PI * r.powi(4);
            let m0 = mid * c - 2.0 * d0;
            [
                m0,
                m0 * 0.5 * (a + b),
                mid * (b.powi(3) - a.powi(3)) / 3.0
                    - (a * a * d0 + 2.0 * a * d1 + d2)
                    - (b * b * d0 - 2.0 * b * d1 + d2),
                (mid.powi(3) * c - 2.0 * e0) / 12.0,
            ]
        }
    }
}

impl FinPlanform {
    /// Root chord, m.
    pub fn root_chord_m(&self) -> f64 {
        match self {
            Self::Trapezoidal { root_chord_m, .. } | Self::Elliptical { root_chord_m, .. } => {
                *root_chord_m
            }
            Self::Freeform { points_m, .. } => points_m.last().map_or(0.0, |p| p[0]),
        }
    }

    /// Span, m.
    pub fn span_m(&self) -> f64 {
        match self {
            Self::Trapezoidal { span_m, .. } | Self::Elliptical { span_m, .. } => *span_m,
            Self::Freeform { points_m, root_m } => {
                points_m.iter().chain(root_m).fold(0.0, |m, p| m.max(p[1]))
            }
        }
    }

    /// Checks dimensions, and that a freeform outline is a simple polygon, nowhere below the root
    /// leading edge, that starts and ends on the root.
    ///
    /// # Errors
    ///
    /// [`DesignError::Domain`] or [`DesignError::Geometry`].
    pub fn validate(&self) -> Result<(), DesignError> {
        match self {
            Self::Trapezoidal {
                root_chord_m,
                tip_chord_m,
                span_m,
                sweep_m,
            } => {
                check_dimension("fin root chord", *root_chord_m, false)?;
                check_dimension("fin tip chord", *tip_chord_m, true)?;
                check_dimension("fin span", *span_m, false)?;
                if !sweep_m.is_finite() {
                    return Err(DesignError::Domain {
                        what: "fin sweep",
                        value: *sweep_m,
                    });
                }
            }
            Self::Elliptical {
                root_chord_m,
                span_m,
            } => {
                check_dimension("fin root chord", *root_chord_m, false)?;
                check_dimension("fin span", *span_m, false)?;
            }
            Self::Freeform { points_m, root_m } => validate_outline(points_m, root_m)?,
        }
        Ok(())
    }

    /// The chord intervals `(a, b)` at span `h` (m outward from the root), sorted along `x` (m aft
    /// of the root leading edge), replacing the contents of `out`. A span at or beyond the tip
    /// gives none. The planform should be valid ([`FinPlanform::validate`]).
    pub fn chords_at(&self, h: f64, out: &mut Vec<(f64, f64)>) {
        out.clear();
        match self {
            Self::Trapezoidal {
                root_chord_m,
                tip_chord_m,
                span_m,
                sweep_m,
            } => {
                if (0.0..*span_m).contains(&h) {
                    let f = h / span_m;
                    let lead = sweep_m * f;
                    out.push((lead, lead + root_chord_m + (tip_chord_m - root_chord_m) * f));
                }
            }
            Self::Elliptical {
                root_chord_m,
                span_m,
            } => {
                if (0.0..*span_m).contains(&h) {
                    let f = h / span_m;
                    let chord = root_chord_m * (1.0 - f * f).max(0.0).sqrt();
                    let lead = 0.5 * (root_chord_m - chord);
                    out.push((lead, lead + chord));
                }
            }
            Self::Freeform { points_m, root_m } => {
                let mut crossings: Vec<f64> = Vec::new();
                let polygon = || points_m.iter().chain(root_m);
                for (&p, &q) in polygon().zip(polygon().cycle().skip(1)) {
                    let (lo, hi) = (p[1].min(q[1]), p[1].max(q[1]));
                    if p[1] != q[1] && h >= lo && h < hi {
                        crossings.push(p[0] + (h - p[1]) * (q[0] - p[0]) / (q[1] - p[1]));
                    }
                }
                crossings.sort_by(f64::total_cmp);
                let (pairs, _) = crossings.as_chunks::<2>();
                out.extend(pairs.iter().map(|&[a, b]| (a, b)));
            }
        }
    }

    /// Span stations where the chord intervals change form: 0, the span, and freeform vertices.
    fn breakpoints(&self) -> Vec<f64> {
        let mut points = vec![0.0, self.span_m()];
        if let Self::Freeform { points_m, root_m } = self {
            points.extend(points_m.iter().chain(root_m).map(|p| p[1]));
        }
        points.sort_by(f64::total_cmp);
        points.dedup();
        points
    }

    /// Area and centroid of the planform.
    ///
    /// # Errors
    ///
    /// As [`FinPlanform::validate`], or [`DesignError::Numerics`] if an integral fails.
    pub fn geometry(&self) -> Result<PlanformGeometry, DesignError> {
        self.validate()?;
        let scale = self.root_chord_m() + self.span_m();
        let mut chords = Vec::new();
        let mut total = [0.0; 3];
        for pair in self.breakpoints().windows(2) {
            let piece = integrate(
                |h| {
                    self.chords_at(h * scale, &mut chords);
                    chords.iter().fold([0.0; 3], |acc, &(a, b)| {
                        let (a, b) = (a / scale, b / scale);
                        [
                            acc[0] + (b - a),
                            acc[1] + 0.5 * (b * b - a * a),
                            acc[2] + h * (b - a),
                        ]
                    })
                },
                pair[0] / scale,
                pair[1] / scale,
                PLATE,
            )?;
            for (sum, value) in total.iter_mut().zip(piece.value) {
                *sum += value;
            }
        }
        let area = total[0] * scale * scale;
        if area <= 0.0 {
            return Err(DesignError::Geometry("the fin has no area".to_owned()));
        }
        Ok(PlanformGeometry {
            area_m2: area,
            centroid_x_m: scale * total[1] / total[0],
            centroid_span_m: scale * total[2] / total[0],
        })
    }
}

/// Quadrature settings for plate integrands scaled to order one.
const PLATE: Tolerance = Tolerance {
    relative: 1e-12,
    absolute: 1e-16,
    max_intervals: 4000,
};

/// Checks a freeform outline and its root points: finite points, at least three in the outline,
/// the first at the root leading edge and the last aft of it, none below it, the root's points
/// running fore between the root's ends, positive area, and no two non-adjacent edges touching.
fn validate_outline(outline: &[[f64; 2]], root: &[[f64; 2]]) -> Result<(), DesignError> {
    if outline.len() < 3 {
        return Err(DesignError::Geometry(format!(
            "a freeform fin needs at least 3 points, got {}",
            outline.len()
        )));
    }
    let points: Vec<[f64; 2]> = outline.iter().chain(root).copied().collect();
    for p in &points {
        for (what, value) in [("freeform fin x", p[0]), ("freeform fin span", p[1])] {
            if !value.is_finite() {
                return Err(DesignError::Domain { what, value });
            }
        }
        if p[1] < 0.0 {
            return Err(DesignError::Domain {
                what: "freeform fin span",
                value: p[1],
            });
        }
    }
    let (first, last) = (outline[0], outline[outline.len() - 1]);
    if first != [0.0, 0.0] || last[0] <= 0.0 {
        return Err(DesignError::Geometry(
            "a freeform fin must start at the root leading edge [0, 0] and end at the root \
             trailing edge [c, h] with c > 0"
                .to_owned(),
        ));
    }
    // The root runs fore from the trailing edge to the leading edge, one point after another.
    let mut aft = last[0];
    for p in root {
        if !(p[0] < aft && p[0] > 0.0) {
            return Err(DesignError::Geometry(format!(
                "a freeform fin's root points must run fore, strictly between its trailing edge \
                 at x = {} m and its leading edge; one is at x = {} m",
                last[0], p[0]
            )));
        }
        aft = p[0];
    }
    let points = points.as_slice();
    let n = points.len();
    let edge = |i: usize| (points[i], points[(i + 1) % n]);
    for i in 0..n {
        for j in (i + 1)..n {
            // Adjacent edges share a vertex; the closing edge is adjacent to the first.
            if j == i + 1 || (i == 0 && j == n - 1) {
                continue;
            }
            let (p, q) = edge(i);
            let (r, s) = edge(j);
            if segments_touch(p, q, r, s) {
                return Err(DesignError::Geometry(format!(
                    "the freeform fin outline crosses itself (edges {i} and {j})"
                )));
            }
        }
    }
    let twice_area: f64 = (0..n)
        .map(|i| {
            let (p, q) = edge(i);
            p[0] * q[1] - q[0] * p[1]
        })
        .sum();
    if twice_area == 0.0 {
        return Err(DesignError::Geometry("the fin has no area".to_owned()));
    }
    // From the root leading edge out along the outline and back along the root, a fin above its
    // root turns clockwise in `[x, h]`: every outline over a level root does. One that turns the
    // other way runs below its root, inside the body.
    if twice_area > 0.0 {
        return Err(DesignError::Geometry(
            "the freeform fin outline runs below its root".to_owned(),
        ));
    }
    Ok(())
}

/// Whether the closed segments `pq` and `rs` share a point.
fn segments_touch(p: [f64; 2], q: [f64; 2], r: [f64; 2], s: [f64; 2]) -> bool {
    let cross = |o: [f64; 2], a: [f64; 2], b: [f64; 2]| {
        (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])
    };
    let within = |a: [f64; 2], b: [f64; 2], c: [f64; 2]| {
        c[0] >= a[0].min(b[0])
            && c[0] <= a[0].max(b[0])
            && c[1] >= a[1].min(b[1])
            && c[1] <= a[1].max(b[1])
    };
    let d1 = cross(r, s, p);
    let d2 = cross(r, s, q);
    let d3 = cross(p, q, r);
    let d4 = cross(p, q, s);
    if ((d1 > 0.0 && d2 < 0.0) || (d1 < 0.0 && d2 > 0.0))
        && ((d3 > 0.0 && d4 < 0.0) || (d3 < 0.0 && d4 > 0.0))
    {
        return true;
    }
    (d1 == 0.0 && within(r, s, p))
        || (d2 == 0.0 && within(r, s, q))
        || (d3 == 0.0 && within(p, q, r))
        || (d4 == 0.0 && within(p, q, s))
}

/// Volume integrals of a plate per unit density, about the fin-set origin with the fin at roll 0:
/// `[V, ∫r, ∫r², ∫x, ∫x², ∫rx, ∫τ²]` (each `dV`).
type PlateIntegrals = [f64; 7];

impl FinSet {
    /// How far a root's point may stand off the surface of the nose cone or transition it sits
    /// on, m: a micron, well inside any build tolerance and far outside the rounding of a file's
    /// seventeen digits.
    pub const ROOT_ON_SURFACE_M: f64 = 1e-6;

    /// The most planform, as a share of the fin's, that a root's straight pieces may cut off or
    /// add across a curved surface: 0.1%, a tenth of the 1% to which `cargo xtask ork` holds a
    /// design's mass to OpenRocket's. On the cockpit of OpenRocket's *Pods--airframes and
    /// winglets*, OpenRocket's own 20 pieces come to 0.032%, the `.ork` reader's 64 to 0.0031%,
    /// and one straight piece, the chord, to 11.4% of its area (12.8% more than the fin's).
    pub const ROOT_SLIVER_SHARE: f64 = 1e-3;

    /// The root's points from leading to trailing edge: its ends, and a freeform fin's root
    /// points between them.
    fn root_points(&self) -> Vec<[f64; 2]> {
        let mut points = vec![[0.0, 0.0]];
        match &self.planform {
            FinPlanform::Freeform { points_m, root_m } => {
                points.extend(root_m.iter().rev());
                points.extend(points_m.last());
            }
            other => points.push([other.root_chord_m(), 0.0]),
        }
        points
    }

    /// Checks that the set can sit on a body tube: its root is level, ending at `h = 0` with no
    /// root points.
    ///
    /// # Errors
    ///
    /// [`DesignError::Geometry`] for a root that isn't level.
    pub fn check_level_root(&self) -> Result<(), DesignError> {
        if self.root_points().iter().any(|p| p[1] != 0.0) {
            return Err(DesignError::Geometry(
                "a fin set on a body tube needs a level root: a freeform outline that ends at \
                 h = 0, with every root point there too"
                    .to_owned(),
            ));
        }
        Ok(())
    }

    /// How the root sits on a nose cone's or transition's surface, with its leading edge `fore_m`
    /// aft of the body's forward end: the farthest any of its points stands off it, m, and the
    /// planform its straight pieces cut off or add across the curve, m². A point at `[x, h]` is on
    /// the surface when `h = r(x_LE + x) − r(x_LE)`. Each piece's sliver is the parabolic segment
    /// `(2/3) δ Δx`, with `δ` the surface's distance from the piece's midpoint, which is exact to
    /// leading order in `Δx` for a smooth surface.
    pub fn root_on_surface(&self, profile: &Profile, fore_m: f64) -> (f64, f64) {
        let base = profile.radius_m(fore_m);
        let miss = |[x, h]: [f64; 2]| h - (profile.radius_m(fore_m + x) - base);
        let points = self.root_points();
        // NaN-propagating, so a broken profile or point is refused by name below.
        let off = points
            .iter()
            .map(|&p| miss(p).abs())
            .fold(0.0, |a: f64, m| {
                if a.is_nan() || m.is_nan() {
                    f64::NAN
                } else {
                    a.max(m)
                }
            });
        let sliver = points
            .windows(2)
            .map(|pair| {
                let middle = [
                    0.5 * (pair[0][0] + pair[1][0]),
                    0.5 * (pair[0][1] + pair[1][1]),
                ];
                2.0 / 3.0 * miss(middle).abs() * (pair[1][0] - pair[0][0])
            })
            .sum();
        (off, sliver)
    }
    /// The radius of a nose cone's or transition's surface at the root leading edge, `fore_m` aft
    /// of the body's forward end, after checking that the set can sit there: the root within the
    /// body's length, its points within [`FinSet::ROOT_ON_SURFACE_M`] of the surface, its pieces'
    /// slivers within [`FinSet::ROOT_SLIVER_SHARE`] of the fin's area
    /// ([`FinSet::root_on_surface`]), and no tab or fillet, whose models take a body tube. The
    /// root's height is measured from that radius.
    ///
    /// # Errors
    ///
    /// [`DesignError::Geometry`] for a set that can't sit on the body as given.
    pub fn root_radius_on(&self, profile: &Profile, fore_m: f64) -> Result<f64, DesignError> {
        let length = profile.length_m();
        let chord = self.planform.root_chord_m();
        let slack = Self::ROOT_ON_SURFACE_M;
        if !(fore_m >= -slack && fore_m + chord <= length + slack) {
            return Err(DesignError::Geometry(format!(
                "a fin set's root runs from {fore_m} m to {} m along a body {length} m long; on a \
                 nose cone or a transition it must stay on it",
                fore_m + chord
            )));
        }
        if self.tab.is_some() || self.fillet.is_some() {
            return Err(DesignError::Geometry(
                "a fin set on a nose cone or a transition with a tab or a fillet: their models \
                 take a body tube"
                    .to_owned(),
            ));
        }
        let (off, sliver) = self.root_on_surface(profile, fore_m);
        // A NaN from a broken profile is refused too.
        if off.is_nan() || off > slack {
            return Err(DesignError::Geometry(format!(
                "a fin set's root stands {off} m off the surface of the nose cone or transition it \
                 sits on, at one of its points"
            )));
        }
        let share = sliver / self.planform.geometry()?.area_m2;
        if share.is_nan() || share > Self::ROOT_SLIVER_SHARE {
            return Err(DesignError::Geometry(format!(
                "a fin set's root cuts across the curved surface it sits on between its points, by \
                 {:.3}% of the fin's area, more than {}%: draw it through more points on the \
                 surface",
                100.0 * share,
                100.0 * Self::ROOT_SLIVER_SHARE
            )));
        }
        Ok(profile.radius_m(fore_m))
    }

    /// The most fins a set may have, the same bound as [`crate::PodSet::MAX_COUNT`].
    ///
    /// The set's mass properties build one body per fin and add them up, so a count is an
    /// allocation and a loop: four billion fins would ask for hundreds of gigabytes and never
    /// return, where a refusal is an error. 64 is far more than any rocket carries: the most fins
    /// in one set among the `.ork` designs hpr's checks read is 8 (a one-off count over 78
    /// designs, not a committed survey), the aerodynamics take
    /// 1 to 8 (the fin–fin interference factors of the OpenRocket technical documentation 13.05,
    /// eq. 3.54, from MIL-HDBK-762(MI) p. 5-24, and its tumbling efficiencies, Table 3.4, stop at
    /// 8), and the `.ork` reader leaves out any part counted more than 64 times. A set of 9 to 64
    /// fins is weighed; the aerodynamics refuse it with their own error.
    pub const MAX_COUNT: u32 = 64;

    /// Checks the set's count, dimensions and planform.
    ///
    /// # Errors
    ///
    /// [`DesignError::Domain`] for no fins or more than [`Self::MAX_COUNT`] (`what` is
    /// `"fin count (1 to 64)"`, the value the count), a bad dimension, cant or base angle, and
    /// [`DesignError::Geometry`] for a planform or tab that doesn't close or fit.
    pub fn validate(&self) -> Result<(), DesignError> {
        if self.count == 0 || self.count > Self::MAX_COUNT {
            return Err(DesignError::Domain {
                what: "fin count (1 to 64)",
                value: f64::from(self.count),
            });
        }
        self.planform.validate()?;
        check_dimension("fin thickness", self.thickness_m, false)?;
        for (what, value) in [
            ("fin cant", self.cant_rad),
            ("fin base angle", self.base_angle_rad),
        ] {
            if !value.is_finite() {
                return Err(DesignError::Domain { what, value });
            }
        }
        if let Some(fillet) = &self.fillet {
            check_dimension("fin fillet radius", fillet.radius_m, true)?;
        }
        if let Some(tab) = self.tab {
            check_dimension("fin tab height", tab.height_m, false)?;
            check_dimension("fin tab length", tab.length_m, false)?;
            if !tab.offset_m.is_finite() {
                return Err(DesignError::Domain {
                    what: "fin tab offset",
                    value: tab.offset_m,
                });
            }
            let root = self.planform.root_chord_m();
            if tab.offset_m < 0.0 || tab.offset_m + tab.length_m > root * (1.0 + 1e-12) {
                return Err(DesignError::Geometry(format!(
                    "a fin tab must lie along the root chord ({root} m)"
                )));
            }
        }
        Ok(())
    }

    /// The planform integrals of one fin, without its tab.
    fn fin_integrals(&self, body_radius_m: f64) -> Result<PlateIntegrals, DesignError> {
        let scale = self.planform.root_chord_m() + self.planform.span_m();
        let rb = body_radius_m / scale;
        let t = self.thickness_m / scale;
        let mut chords = Vec::new();
        let mut total = [0.0; 7];
        for pair in self.planform.breakpoints().windows(2) {
            let piece = integrate(
                |h| {
                    self.planform.chords_at(h * scale, &mut chords);
                    let mut m = [0.0; 4];
                    for &(a, b) in &chords {
                        let c = chord_moments(self.cross_section, a / scale, b / scale, t);
                        for k in 0..4 {
                            m[k] += c[k];
                        }
                    }
                    let r = rb + h;
                    [m[0], r * m[0], r * r * m[0], m[1], m[2], r * m[1], m[3]]
                },
                pair[0] / scale,
                pair[1] / scale,
                PLATE,
            )?;
            for (sum, value) in total.iter_mut().zip(piece.value) {
                *sum += value;
            }
        }
        // Undo the scaling: lengths to the powers 3, 4, 5, 4, 5, 5, 5.
        let powers = [3, 4, 5, 4, 5, 5, 5];
        Ok(std::array::from_fn(|k| total[k] * scale.powi(powers[k])))
    }

    /// The integrals of the tab slab `offset ≤ x ≤ offset + length`, `−height ≤ h ≤ 0`.
    fn tab_integrals(&self, body_radius_m: f64) -> PlateIntegrals {
        let Some(tab) = self.tab else {
            return [0.0; 7];
        };
        let t = self.thickness_m;
        let (x0, x1) = (tab.offset_m, tab.offset_m + tab.length_m);
        let (r0, r1) = (body_radius_m - tab.height_m, body_radius_m);
        let (h, l) = (tab.height_m, tab.length_m);
        let r_first = 0.5 * (r1 * r1 - r0 * r0);
        let x_first = 0.5 * (x1 * x1 - x0 * x0);
        [
            t * l * h,
            t * l * r_first,
            t * l * (r1.powi(3) - r0.powi(3)) / 3.0,
            t * h * x_first,
            t * h * (x1.powi(3) - x0.powi(3)) / 3.0,
            t * x_first * r_first,
            t.powi(3) / 12.0 * l * h,
        ]
    }

    /// Mass properties of one fin (and its tab) at roll angle 0 and no cant, in the fin-set
    /// frame, on a body of radius `body_radius_m`.
    ///
    /// # Errors
    ///
    /// As [`FinSet::validate`], plus [`DesignError::Domain`] for a negative body radius or a
    /// fillet more than [`FILLET_RATIO_MAX`] times the body radius, and
    /// [`DesignError::Numerics`] if an integral fails.
    pub fn single_fin(&self, body_radius_m: f64) -> Result<MassProperties, DesignError> {
        self.validate()?;
        check_dimension("fin body radius", body_radius_m, true)?;
        if let Some(tab) = self.tab
            && tab.height_m > body_radius_m
        {
            return Err(DesignError::Geometry(format!(
                "a fin tab must reach no deeper than the body radius ({body_radius_m} m)"
            )));
        }
        let density = self.material.bulk_kg_m3("fin set")?;
        let fin = self.fin_integrals(body_radius_m)?;
        let tab = self.tab_integrals(body_radius_m);
        let [v, r1, r2, x1, x2, rx, tau2] = std::array::from_fn(|k| fin[k] + tab[k]);
        if v <= 0.0 {
            return Err(DesignError::Geometry("the fin has no volume".to_owned()));
        }
        let mass = density * v;
        // Points sit at (r, τ, −x): I_xz = −∫ x_B z_B dm = ∫ r x dm.
        let about_origin = DMat3::from_cols(
            DVec3::new(tau2 + x2, 0.0, rx),
            DVec3::new(0.0, r2 + x2, 0.0),
            DVec3::new(rx, 0.0, r2 + tau2),
        ) * density;
        let cg = DVec3::new(r1 / v, 0.0, -x1 / v);
        let at_origin = MassProperties {
            mass_kg: mass,
            cg_m: cg,
            inertia_kg_m2: DMat3::ZERO,
        };
        // Move the tensor from the origin to the center: I_cg = I_o − m(|c|²E − c cᵀ).
        let inertia = about_origin - at_origin.inertia_about(DVec3::ZERO);
        let plate = MassProperties {
            inertia_kg_m2: inertia,
            ..at_origin
        };
        let fin = match self.fillets(body_radius_m)? {
            Some(fillets) => MassProperties::combine([&plate, &fillets]),
            None => plate,
        };
        if self.cant_rad == 0.0 {
            return Ok(fin);
        }
        let pivot = DVec3::new(body_radius_m, 0.0, -0.5 * self.planform.root_chord_m());
        Ok(fin
            .translated(-pivot)
            .rotated(DQuat::from_rotation_x(self.cant_rad))
            .translated(pivot))
    }

    /// Mass properties of the two fillets along one fin's root, at roll angle 0 and no cant, in
    /// the fin-set frame; `None` for a fin with no fillets or fillets of no radius. The pair is a
    /// prism of length `c_r` whose section is [`fillet_section`] on each side of the fin's plane,
    /// so with `A`, `S_x`, `S_xx`, `S_yy` one side's area and moments about the body axis:
    ///
    /// ```text
    /// m = 2ρA c_r,  I_xx = 2ρ(S_yy c_r + A c_r³/3),  I_yy = 2ρ(S_xx c_r + A c_r³/3),
    /// I_zz = 2ρ(S_xx + S_yy) c_r,  I_xz = ρ S_x c_r²
    /// ```
    ///
    /// about the origin, then moved to the center, `(S_x/A, 0, −c_r/2)`.
    fn fillets(&self, body_radius_m: f64) -> Result<Option<MassProperties>, DesignError> {
        let Some(fillet) = self.fillet.as_ref().filter(|f| f.radius_m > 0.0) else {
            return Ok(None);
        };
        let density = fillet.material.bulk_kg_m3("fin fillet")?;
        // No body, no fillet: the section is none ([`fillet_section`]).
        if body_radius_m == 0.0 {
            return Ok(None);
        }
        // The section is a triangle less two sectors, and far from real fillets the difference
        // loses digits to cancellation. At the bounds below it keeps 10 (a millionth of the
        // body's radius: 8.6e-11 off 60-digit arithmetic) and 12 (a thousand times: 7.6e-13).
        // A thinner fillet weighs under 1e-12 of the body's radius squared per meter of root:
        // none. A wider one is refused, not weighed on digits that are noise.
        let ratio = fillet.radius_m / body_radius_m;
        if ratio < 1e-6 {
            return Ok(None);
        }
        if ratio > FILLET_RATIO_MAX {
            return Err(DesignError::Domain {
                what: "fin fillet radius over the body radius (at most 1000)",
                value: ratio,
            });
        }
        let [area, sx, sxx, syy] = fillet_section(body_radius_m, fillet.radius_m);
        // Invariant between the bounds above: the section is positive and finite there.
        if !(area > 0.0 && [sx, sxx, syy].iter().all(|v| v.is_finite())) {
            return Err(DesignError::Domain {
                what: "fin fillet section area",
                value: area,
            });
        }
        let l = self.planform.root_chord_m();
        let mass = 2.0 * density * area * l;
        if mass == 0.0 {
            return Ok(None);
        }
        let ixz = density * sx * l * l;
        let about_origin = DMat3::from_cols(
            DVec3::new(2.0 * density * (syy * l + area * l.powi(3) / 3.0), 0.0, ixz),
            DVec3::new(0.0, 2.0 * density * (sxx * l + area * l.powi(3) / 3.0), 0.0),
            DVec3::new(ixz, 0.0, 2.0 * density * (sxx + syy) * l),
        );
        let at_origin = MassProperties {
            mass_kg: mass,
            cg_m: DVec3::new(sx / area, 0.0, -0.5 * l),
            inertia_kg_m2: DMat3::ZERO,
        };
        Ok(Some(MassProperties {
            inertia_kg_m2: about_origin - at_origin.inertia_about(DVec3::ZERO),
            ..at_origin
        }))
    }

    /// Mass properties of the whole set in its frame, on a body of radius `body_radius_m`.
    ///
    /// # Errors
    ///
    /// As [`FinSet::single_fin`].
    pub fn mass_properties(&self, body_radius_m: f64) -> Result<MassProperties, DesignError> {
        let fin = self.single_fin(body_radius_m)?;
        let fins: Vec<MassProperties> = (0..self.count)
            .map(|k| {
                fin.rolled(self.base_angle_rad + 2.0 * PI * f64::from(k) / f64::from(self.count))
            })
            .collect();
        Ok(MassProperties::combine(&fins))
    }
}

/// The widest fin fillet weighed, as a multiple of the body radius ([`FinSet::single_fin`]).
pub const FILLET_RATIO_MAX: f64 = 1e3;

/// One fillet's section on the `+y` side of a fin in the plane `y = 0`, with `x` outward along
/// the fin: its area and moments `[A, S_x, S_xx, S_yy]` about the body axis (`S_x = ∫ x dA`,
/// `S_xx = ∫ x² dA`, `S_yy = ∫ y² dA`), m², m³, m⁴, m⁴. The section is the triangle
/// `(0, 0), (c, 0), (c, r)` less two circular sectors, the body's about the axis from `0` to `θ`
/// and the fillet's about `(c, r)` from `θ − π` to `−π/2` (module docs). A body of no radius
/// gives a section of none.
pub(crate) fn fillet_section(body_radius_m: f64, radius_m: f64) -> [f64; 4] {
    let (rb, r) = (body_radius_m, radius_m);
    let c = (rb * rb + 2.0 * rb * r).sqrt();
    let theta = r.atan2(c);
    // The triangle: vertices (0, 0), (c, 0), (c, r).
    let triangle_area = 0.5 * c * r;
    let triangle = [
        triangle_area,
        triangle_area * 2.0 * c / 3.0,
        triangle_area * c * c / 2.0,
        triangle_area * r * r / 6.0,
    ];
    let body = sector([0.0, 0.0], rb, 0.0, theta);
    let joint = sector([c, r], r, theta - PI, -0.5 * PI);
    std::array::from_fn(|k| triangle[k] - body[k] - joint[k])
}

/// A circular sector about `center`, of radius `rho`, from angle `from` to `to` (`to ≥ from`):
/// `[A, S_x, S_xx, S_yy]` about the origin. With `Δ = to − from` and `σ = to + from`, the moments
/// about the center are `ρ³/3 · 2 sin(Δ/2) [cos(σ/2), sin(σ/2)]` and `ρ⁴/8 · (Δ ± sin Δ cos σ)`,
/// the second written as `(Δ − sin Δ) + sin Δ (1 ± cos σ)` so that the body's thin sector, whose
/// `Δ` and `σ` are both small, loses no digits to cancellation.
fn sector([x0, y0]: [f64; 2], rho: f64, from: f64, to: f64) -> [f64; 4] {
    let (span, sum) = (to - from, to + from);
    let area = 0.5 * rho * rho * span;
    let chord = 2.0 * (0.5 * span).sin() * rho.powi(3) / 3.0;
    let (sx, sy) = (chord * (0.5 * sum).cos(), chord * (0.5 * sum).sin());
    let (lean, sine) = (span_less_sine(span), span.sin());
    let (sxx, syy) = (
        rho.powi(4) / 8.0 * (lean + sine * 2.0 * (0.5 * sum).cos().powi(2)),
        rho.powi(4) / 8.0 * (lean + sine * 2.0 * (0.5 * sum).sin().powi(2)),
    );
    [
        area,
        sx + x0 * area,
        sxx + 2.0 * x0 * sx + x0 * x0 * area,
        syy + 2.0 * y0 * sy + y0 * y0 * area,
    ]
}

/// `x − sin x`, by its Taylor series below `x = 0.25`, where the difference would cancel: the
/// terms through `x¹⁵` leave a remainder under 1e-18 of the sum.
fn span_less_sine(x: f64) -> f64 {
    if x.abs() >= 0.25 {
        return x - x.sin();
    }
    // x³/3! − x⁵/5! + …, each term the last times −x²/((2k)(2k+1)).
    let mut term = x.powi(3) / 6.0;
    let mut sum = 0.0;
    for k in 2..=8 {
        sum += term;
        term *= -x * x / f64::from((2 * k) * (2 * k + 1));
    }
    sum
}

/// Tube fins: open tubes parallel to the body, touching it, spaced evenly around it.
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, schemars::JsonSchema)]
#[serde(deny_unknown_fields)]
pub struct TubeFinSet {
    /// Number of tubes, 1 to 64 (`TubeFinSet::MAX_COUNT`).
    pub count: u32,
    /// Length, m.
    pub length_m: f64,
    /// Outer radius, m.
    pub outer_radius_m: f64,
    /// Wall thickness, m.
    pub thickness_m: f64,
    /// Roll angle of the first tube's axis from `x_B` toward `y_B`, rad.
    #[serde(default)]
    pub base_angle_rad: f64,
    /// Material (bulk).
    pub material: Material,
}

impl TubeFinSet {
    /// The most tubes a set may have: as [`FinSet::MAX_COUNT`], a bound on the bodies its mass
    /// properties build one by one, far above any rocket's (the most tubes in one set among the
    /// `.ork` designs read is 6, a one-off count over 78 designs, not a committed survey).
    pub const MAX_COUNT: u32 = 64;

    /// The outer radius of `count` tubes that close the ring around a body of radius
    /// `body_radius_m`, each touching the body and its two neighbours, which is the radius
    /// OpenRocket 24.12 gives a tube fin set written `auto` ([ADR-098][adr-098]).
    ///
    /// The axes of `N` tubes of radius `r` sit on a circle of radius `R + r`, `2π/N` apart, so
    /// neighbouring axes are a chord `2(R + r) sin(π/N)` apart. Touching means that chord is `2r`:
    ///
    /// `r = R sin(π/N) / (1 − sin(π/N))`.
    ///
    /// For one or two tubes `sin(π/N)` is 0 or 1 and no finite ring closes; OpenRocket 24.12 gives
    /// those the body's radius, measured on its probes on 20 and 50 mm bodies, and so does this;
    /// so does a count of none, which [`Self::mass_properties`] refuses. On a 50 mm body, six
    /// tubes are 50 mm (`sin(π/6) = 1/2`), four are 120.7 mm and eight are 31.0 mm.
    ///
    /// [adr-098]: https://github.com/nrdptel/hpr-sim/blob/main/docs/DECISIONS.md#adr-098-a-tube-fin-sets-automatic-radius-read-as-openrocket-reads-it-2026-09-28
    pub fn closing_radius_m(body_radius_m: f64, count: u32) -> f64 {
        if count < 3 {
            return body_radius_m;
        }
        let s = (PI / f64::from(count)).sin();
        body_radius_m * s / (1.0 - s)
    }

    /// Mass properties in the set's frame (origin on the body axis at the tubes' forward end), on
    /// a body of radius `body_radius_m`. Each tube is a hollow cylinder,
    /// `I_a = m(R² + r²)/2` and `I_t = m((R² + r²)/4 + L²/12)`, with its axis at `R_b + R`.
    ///
    /// # Errors
    ///
    /// [`DesignError::Domain`] for a bad dimension, or a count of none or more than
    /// [`Self::MAX_COUNT`] (`what` is `"tube fin count (1 to 64)"`), [`DesignError::Geometry`] for
    /// a wall thicker than the radius, and material errors.
    pub fn mass_properties(&self, body_radius_m: f64) -> Result<MassProperties, DesignError> {
        if self.count == 0 || self.count > Self::MAX_COUNT {
            return Err(DesignError::Domain {
                what: "tube fin count (1 to 64)",
                value: f64::from(self.count),
            });
        }
        check_dimension("body radius", body_radius_m, true)?;
        if !self.base_angle_rad.is_finite() {
            return Err(DesignError::Domain {
                what: "tube fin base angle",
                value: self.base_angle_rad,
            });
        }
        let density = self.material.bulk_kg_m3("tube fin set")?;
        let tube = crate::parts::hollow_cylinder(
            "tube fin",
            density,
            self.length_m,
            self.outer_radius_m,
            self.thickness_m,
        )?;
        let one = tube.translated(DVec3::new(body_radius_m + self.outer_radius_m, 0.0, 0.0));
        let tubes: Vec<MassProperties> = (0..self.count)
            .map(|k| {
                one.rolled(self.base_angle_rad + 2.0 * PI * f64::from(k) / f64::from(self.count))
            })
            .collect();
        Ok(MassProperties::combine(&tubes))
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use hpr_core::quadrature::integrate_scalar;

    fn close(got: f64, want: f64, rel: f64, what: &str) {
        let err = if want == 0.0 {
            got.abs()
        } else {
            ((got - want) / want).abs()
        };
        assert!(err <= rel, "{what}: {got} vs {want} (relative {err:e})");
    }

    fn mat_close(a: DMat3, b: DMat3, rel: f64, what: &str) {
        let scale = b.to_cols_array().iter().fold(0.0f64, |m, v| m.max(v.abs()));
        let diff = (a - b)
            .to_cols_array()
            .iter()
            .fold(0.0f64, |m, v| m.max(v.abs()));
        assert!(diff <= rel * scale, "{what}: {a:?}\nvs\n{b:?}");
    }

    fn ply() -> Material {
        Material::bulk("plywood", 630.0)
    }

    fn rectangle(chord: f64, span: f64) -> FinPlanform {
        FinPlanform::Trapezoidal {
            root_chord_m: chord,
            tip_chord_m: chord,
            span_m: span,
            sweep_m: 0.0,
        }
    }

    fn set(count: u32, planform: FinPlanform, thickness: f64) -> FinSet {
        FinSet {
            count,
            planform,
            thickness_m: thickness,
            cross_section: FinCrossSection::Square,
            tab: None,
            fillet: None,
            cant_rad: 0.0,
            base_angle_rad: 0.0,
            material: ply(),
        }
    }

    #[test]
    fn airfoil_constants_are_the_integrals() {
        // With ξ = u², dξ = 2u du, and the integrands are polynomials in u.
        let tol = Tolerance::default();
        let p = naca_polynomial;
        for (k, want) in AIRFOIL_PHI.iter().enumerate() {
            let power = i32::try_from(2 * k).unwrap();
            let got = 10.0
                * integrate_scalar(|u| u.powi(power) * p(u * u) * 2.0 * u, 0.0, 1.0, tol).unwrap();
            close(got, *want, 1e-14, "phi");
        }
        let psi = 1000.0 * integrate_scalar(|u| p(u * u).powi(3) * 2.0 * u, 0.0, 1.0, tol).unwrap();
        close(psi, AIRFOIL_PSI, 1e-14, "psi");
    }

    #[test]
    fn chord_moments_match_quadrature_of_each_section() {
        // The references ask for the tightest tolerance the quadrature allows (50 ε relative).
        let tol = Tolerance {
            relative: 0.0,
            absolute: 0.0,
            max_intervals: 4000,
        };
        let t: f64 = 0.004;
        // (a, b) wider than the thickness, and a sliver narrower than it (a circle of diameter c).
        for (a, b) in [(0.02, 0.17), (0.05, 0.053)] {
            let c: f64 = b - a;
            let rounded = |x: f64| {
                let mid = t.min(c);
                let r = 0.5 * mid;
                let u = (x - a).min(b - x);
                if u >= r {
                    mid
                } else {
                    2.0 * (r * r - (r - u).powi(2)).max(0.0).sqrt()
                }
            };
            let airfoil = |x: f64| 10.0 * t * naca_polynomial((x - a) / c);
            let square = |_x: f64| t;
            let sections: [(FinCrossSection, &dyn Fn(f64) -> f64); 3] = [
                (FinCrossSection::Square, &square),
                (FinCrossSection::Rounded, &rounded),
                (FinCrossSection::Airfoil, &airfoil),
            ];
            // Exact references: split where the edges meet the flat, and substitute x = a + v²
            // (and x = b − v²) at the ends, where the thickness grows like a square root.
            let edge = 0.5 * t.min(c);
            let reference = |f: &dyn Fn(f64) -> f64| -> f64 {
                let fore = integrate_scalar(|v| f(a + v * v) * 2.0 * v, 0.0, edge.sqrt(), tol);
                let aft = integrate_scalar(|v| f(b - v * v) * 2.0 * v, 0.0, edge.sqrt(), tol);
                let middle = integrate_scalar(f, a + edge, b - edge, tol);
                fore.unwrap() + aft.unwrap() + middle.unwrap()
            };
            for (section, thickness) in sections {
                let got = chord_moments(section, a, b, t);
                for (k, &moment) in got.iter().take(3).enumerate() {
                    let power = i32::try_from(k).unwrap();
                    let want = reference(&|x: f64| x.powi(power) * thickness(x));
                    close(
                        moment,
                        want,
                        1e-13,
                        &format!("{section:?} M{k} on [{a}, {b}]"),
                    );
                }
                let want = reference(&|x: f64| thickness(x).powi(3) / 12.0);
                close(got[3], want, 1e-13, &format!("{section:?} T on [{a}, {b}]"));
            }
        }
        // The rounded section of a wide chord loses (1 − π/4) t² of area.
        let m = chord_moments(FinCrossSection::Rounded, 0.0, 0.1, t);
        close(
            m[0],
            0.1 * t - (1.0 - PI / 4.0) * t * t,
            1e-14,
            "rounded area",
        );
    }

    #[test]
    fn planform_areas_and_centroids_are_the_closed_forms() {
        let (cr, ct, s, xt) = (0.15, 0.06, 0.1, 0.07);
        let trapezoid = FinPlanform::Trapezoidal {
            root_chord_m: cr,
            tip_chord_m: ct,
            span_m: s,
            sweep_m: xt,
        };
        let g = trapezoid.geometry().unwrap();
        close(g.area_m2, 0.5 * s * (cr + ct), 1e-13, "trapezoid area");
        close(
            g.centroid_x_m,
            (xt * (cr + 2.0 * ct) + cr * cr + cr * ct + ct * ct) / (3.0 * (cr + ct)),
            1e-13,
            "trapezoid centroid x",
        );
        close(
            g.centroid_span_m,
            s * (cr + 2.0 * ct) / (3.0 * (cr + ct)),
            1e-13,
            "trapezoid centroid h",
        );
        // The same outline as a freeform polygon.
        let freeform = FinPlanform::Freeform {
            points_m: vec![[0.0, 0.0], [xt, s], [xt + ct, s], [cr, 0.0]],
            root_m: Vec::new(),
        };
        let f = freeform.geometry().unwrap();
        close(f.area_m2, g.area_m2, 1e-13, "freeform area");
        close(f.centroid_x_m, g.centroid_x_m, 1e-13, "freeform centroid x");
        close(
            f.centroid_span_m,
            g.centroid_span_m,
            1e-13,
            "freeform centroid h",
        );
        // Half ellipse: area π c_r s / 4, centroid at mid-chord and 4s/3π out.
        let e = FinPlanform::Elliptical {
            root_chord_m: cr,
            span_m: s,
        }
        .geometry()
        .unwrap();
        close(e.area_m2, PI * cr * s / 4.0, 1e-12, "ellipse area");
        close(e.centroid_x_m, cr / 2.0, 1e-12, "ellipse centroid x");
        close(
            e.centroid_span_m,
            4.0 * s / (3.0 * PI),
            1e-12,
            "ellipse centroid h",
        );
    }

    #[test]
    fn a_rectangular_fin_set_matches_box_inertia_by_hand() {
        let (c, s, t, rb, rho) = (0.12, 0.08, 0.003, 0.04, 630.0);
        let m = rho * c * s * t;
        let d = rb + s / 2.0;
        // One fin: a box with s along x, t along y, c along z, centered at (d, 0, −c/2).
        let one = set(1, rectangle(c, s), t).single_fin(rb).unwrap();
        close(one.mass_kg, m, 1e-13, "mass");
        assert!((one.cg_m - DVec3::new(d, 0.0, -c / 2.0)).length() < 1e-15);
        let box_inertia = DMat3::from_diagonal(DVec3::new(
            m * (t * t + c * c) / 12.0,
            m * (s * s + c * c) / 12.0,
            m * (s * s + t * t) / 12.0,
        ));
        mat_close(one.inertia_kg_m2, box_inertia, 1e-12, "one fin");
        // Four fins: two along ±x and two along ±y.
        let four = set(4, rectangle(c, s), t).mass_properties(rb).unwrap();
        close(four.mass_kg, 4.0 * m, 1e-13, "four fins");
        assert!(four.cg_m.truncate().length() < 1e-15);
        let transverse =
            2.0 * m * (t * t + c * c) / 12.0 + 2.0 * (m * (s * s + c * c) / 12.0 + m * d * d);
        let axial = 4.0 * (m * (s * s + t * t) / 12.0 + m * d * d);
        mat_close(
            four.inertia_kg_m2,
            DMat3::from_diagonal(DVec3::new(transverse, transverse, axial)),
            1e-12,
            "four fins",
        );
        // Three fins are isotropic across the axis too; two are not.
        let three = set(3, rectangle(c, s), t)
            .mass_properties(rb)
            .unwrap()
            .inertia_kg_m2;
        close(
            three.x_axis.x,
            three.y_axis.y,
            1e-12,
            "three fins isotropic",
        );
        assert!(three.y_axis.x.abs() < 1e-12 * three.x_axis.x);
        let two = set(2, rectangle(c, s), t)
            .mass_properties(rb)
            .unwrap()
            .inertia_kg_m2;
        assert!((two.x_axis.x - two.y_axis.y).abs() > 1e-4 * two.y_axis.y);

        // Roots on the axis, as on a pod of no length (M1.13b2): the same boxes, `d = s/2` out.
        let d = s / 2.0;
        let on_axis = set(4, rectangle(c, s), t).mass_properties(0.0).unwrap();
        let transverse =
            2.0 * m * (t * t + c * c) / 12.0 + 2.0 * (m * (s * s + c * c) / 12.0 + m * d * d);
        let axial = 4.0 * (m * (s * s + t * t) / 12.0 + m * d * d);
        mat_close(
            on_axis.inertia_kg_m2,
            DMat3::from_diagonal(DVec3::new(transverse, transverse, axial)),
            1e-12,
            "four fins from the axis",
        );
    }

    #[test]
    fn cant_turns_the_fin_about_its_span_axis() {
        let (c, s, t, rb, rho) = (0.12, 0.08, 0.003, 0.04, 630.0);
        let m = rho * c * s * t;
        // Canted 90°, the chord lies along y: a box with s along x, c along y, t along z, whose
        // center stays at the pivot's station because the pivot is the root mid-chord.
        let mut fins = set(1, rectangle(c, s), t);
        fins.cant_rad = PI / 2.0;
        let turned = fins.single_fin(rb).unwrap();
        assert!((turned.cg_m - DVec3::new(rb + s / 2.0, 0.0, -c / 2.0)).length() < 1e-15);
        let expected = DMat3::from_diagonal(DVec3::new(
            m * (c * c + t * t) / 12.0,
            m * (s * s + t * t) / 12.0,
            m * (s * s + c * c) / 12.0,
        ));
        mat_close(turned.inertia_kg_m2, expected, 1e-12, "canted 90°");
        // A small cant only couples x and... keeps mass and center and the trace.
        fins.cant_rad = 0.05;
        let small = fins.single_fin(rb).unwrap();
        let flat = set(1, rectangle(c, s), t).single_fin(rb).unwrap();
        close(small.mass_kg, flat.mass_kg, 1e-15, "mass");
        let trace = |i: DMat3| i.x_axis.x + i.y_axis.y + i.z_axis.z;
        close(
            trace(small.inertia_kg_m2),
            trace(flat.inertia_kg_m2),
            1e-12,
            "trace",
        );
    }

    #[test]
    fn positive_cant_turns_the_leading_edge_toward_negative_y() {
        // A swept fin's centroid lies aft of the root mid-chord, so turning it by δ about the span
        // axis through the mid-chord moves the centroid to y = (x̄ − c_r/2) sin δ and
        // z = −c_r/2 − (x̄ − c_r/2) cos δ; the leading edge (x = 0) goes to y = −(c_r/2) sin δ.
        let (cr, ct, s, xt) = (0.15, 0.05, 0.1, 0.09);
        let planform = FinPlanform::Trapezoidal {
            root_chord_m: cr,
            tip_chord_m: ct,
            span_m: s,
            sweep_m: xt,
        };
        let centroid = planform.geometry().unwrap().centroid_x_m;
        assert!(centroid > cr / 2.0);
        let mut fins = set(1, planform, 0.003);
        let delta: f64 = 0.2;
        fins.cant_rad = delta;
        let fin = fins.single_fin(0.04).unwrap();
        close(
            fin.cg_m.y,
            (centroid - cr / 2.0) * delta.sin(),
            1e-12,
            "centroid y",
        );
        close(
            fin.cg_m.z,
            -cr / 2.0 - (centroid - cr / 2.0) * delta.cos(),
            1e-12,
            "centroid z",
        );
        assert!(fin.cg_m.y > 0.0);
        let leading_edge = DQuat::from_rotation_x(delta) * DVec3::new(0.04, 0.0, cr / 2.0);
        assert!(leading_edge.y < 0.0);
    }

    #[test]
    fn tabs_must_fit_the_root_and_the_body() {
        let mut fins = set(3, rectangle(0.1, 0.05), 0.003);
        // Along the root is checked without a body; the depth needs one.
        fins.tab = Some(FinTab {
            height_m: 0.01,
            length_m: 0.08,
            offset_m: 0.03,
        });
        assert!(matches!(fins.validate(), Err(DesignError::Geometry(_))));
        for tab in [
            FinTab {
                height_m: 0.05,
                length_m: 0.02,
                offset_m: 0.0,
            },
            FinTab {
                height_m: 0.01,
                length_m: 0.02,
                offset_m: -0.01,
            },
            FinTab {
                height_m: 0.01,
                length_m: 0.08,
                offset_m: 0.03,
            },
        ] {
            fins.tab = Some(tab);
            assert!(
                matches!(fins.mass_properties(0.03), Err(DesignError::Geometry(_))),
                "{tab:?}"
            );
        }
        fins.tab = Some(FinTab {
            height_m: 0.03,
            length_m: 0.1,
            offset_m: 0.0,
        });
        fins.mass_properties(0.03).unwrap();
    }

    #[test]
    fn a_swept_fin_has_the_parallel_axis_product_of_inertia() {
        // A single trapezoidal fin: I_xz about the origin is ∫ r x dm, so about the center it is
        // ∫ r x dm − m r̄ x̄ (with z = −x, I_xz = −∫ x_B z_B dm). Check against direct quadrature of
        // the planform.
        let (cr, ct, s, xt, t, rb, rho) = (0.15, 0.05, 0.1, 0.09, 0.004, 0.05, 1800.0);
        let planform = FinPlanform::Trapezoidal {
            root_chord_m: cr,
            tip_chord_m: ct,
            span_m: s,
            sweep_m: xt,
        };
        let mut fins = set(1, planform, t);
        fins.material = Material::bulk("G10", rho);
        let fin = fins.single_fin(rb).unwrap();
        let tol = Tolerance::default();
        let lead = |h: f64| xt * h / s;
        let chord = |h: f64| cr + (ct - cr) * h / s;
        let area = integrate_scalar(chord, 0.0, s, tol).unwrap();
        let mx = integrate_scalar(|h| chord(h) * (lead(h) + chord(h) / 2.0), 0.0, s, tol).unwrap();
        let mr = integrate_scalar(|h| chord(h) * (rb + h), 0.0, s, tol).unwrap();
        let mrx = integrate_scalar(
            |h| (rb + h) * chord(h) * (lead(h) + chord(h) / 2.0),
            0.0,
            s,
            tol,
        )
        .unwrap();
        let m = rho * t * area;
        close(fin.mass_kg, m, 1e-13, "mass");
        let product = rho * t * mrx - m * (mr / area) * (mx / area);
        close(fin.inertia_kg_m2.z_axis.x, product, 1e-11, "I_xz");
        close(fin.inertia_kg_m2.x_axis.z, product, 1e-11, "I_zx");
    }

    #[test]
    fn a_closing_ring_of_tubes_touches_the_body_and_its_neighbours() {
        let body = 0.05;
        for count in 3..=12 {
            let r = TubeFinSet::closing_radius_m(body, count);
            assert!(r > 0.0 && r.is_finite(), "{count}: {r}");
            // Neighbouring axes, on the circle of radius R + r, are 2r apart: the tubes touch.
            let step = 2.0 * PI / f64::from(count);
            let (a, b) = (
                DVec3::new(body + r, 0.0, 0.0),
                DVec3::new((body + r) * step.cos(), (body + r) * step.sin(), 0.0),
            );
            close(a.distance(b), 2.0 * r, 1e-14, "neighbours touch");
        }
        // Six tubes are as wide as the body; four much wider, eight narrower.
        close(TubeFinSet::closing_radius_m(body, 6), body, 1e-15, "six");
        close(
            TubeFinSet::closing_radius_m(body, 4),
            body * (1.0 + 2f64.sqrt()),
            1e-15,
            "four",
        );
        // No finite ring closes with one or two: the body's radius, as OpenRocket 24.12 gives.
        assert_eq!(TubeFinSet::closing_radius_m(body, 1), body);
        assert_eq!(TubeFinSet::closing_radius_m(body, 2), body);
    }

    #[test]
    fn tube_fins_are_hollow_cylinders_around_the_body() {
        let tubes = TubeFinSet {
            count: 6,
            length_m: 0.1,
            outer_radius_m: 0.012,
            thickness_m: 0.001,
            base_angle_rad: 0.0,
            material: Material::bulk("cardboard", 790.0),
        };
        let rb = 0.02;
        let g = tubes.mass_properties(rb).unwrap();
        let (ro, ri) = (0.012, 0.011);
        let m1 = 790.0 * PI * (ro * ro - ri * ri) * 0.1;
        close(g.mass_kg, 6.0 * m1, 1e-13, "mass");
        let d = rb + ro;
        let axial = 6.0 * (0.5 * m1 * (ro * ro + ri * ri) + m1 * d * d);
        close(g.inertia_kg_m2.z_axis.z, axial, 1e-13, "axial");
        let transverse = 6.0 * m1 * ((ro * ro + ri * ri) / 4.0 + 0.01 / 12.0) + 3.0 * m1 * d * d;
        close(g.inertia_kg_m2.x_axis.x, transverse, 1e-12, "transverse");
        close(g.cg_m.z, -0.05, 1e-15, "center");
    }

    /// OpenRocket's *Pods--airframes and winglets* example's cockpit: one balsa fin, 6.985 mm
    /// thick, a triangle on its tangent ogive's last 50 mm, its root through `pieces` straight
    /// pieces on the surface. Returns the set and its root leading edge's station on the nose.
    fn cockpit(pieces: u32) -> (FinSet, Profile, f64) {
        let nose = Profile::nose(
            crate::NoseShape::Ogive { radius_ratio: 1.0 },
            0.136525,
            0.0168275,
        )
        .unwrap();
        let (chord, fore) = (0.05, 0.136525 - 0.05);
        let base = nose.radius_m(fore);
        let rise = |x: f64| nose.radius_m(fore + x) - base;
        let root_m = (1..pieces)
            .rev()
            .map(|i| {
                let x = chord * f64::from(i) / f64::from(pieces);
                [x, rise(x)]
            })
            .collect();
        let planform = FinPlanform::Freeform {
            points_m: vec![
                [0.0, 0.0],
                [0.009347826086956524, 0.006956521739130436],
                [chord, 0.0022276567072510058],
            ],
            root_m,
        };
        let mut set = set(1, planform, 0.006985);
        set.material = Material::bulk("balsa", 170.0);
        (set, nose, fore)
    }

    /// M4.5g4 (ADR-166): a fin's root along a nose cone is read as OpenRocket 24.12 reads it.
    /// OpenRocket draws the cockpit's root through 21 points, 2.5 mm apart: through the same
    /// points, hpr's area, mass and center of mass are OpenRocket's (its `getPlanformArea`,
    /// `getMass` and `getComponentCG`, to the six digits the oracle printed). hpr's own 64
    /// pieces sit closer to the curve, which bulges out of each chord: 2.9e-4 less area than
    /// OpenRocket's 20, and 3.1e-5 more than the curve itself.
    #[test]
    fn a_root_along_a_nose_cone_is_openrockets() {
        let (set, nose, fore) = cockpit(20);
        let radius = set.root_radius_on(&nose, fore).unwrap();
        close(
            radius,
            0.0145998,
            1e-5,
            "body radius at the root leading edge",
        );
        let g = set.planform.geometry().unwrap();
        close(g.area_m2, 1.449544e-4, 1e-6, "OpenRocket's planform area");
        let m = set.mass_properties(radius).unwrap();
        close(m.mass_kg, 1.72126e-4, 1e-5, "OpenRocket's mass");
        // Fin 0 at roll 0: its points at (r, τ, −x) about the root leading edge on the axis.
        close(
            -m.cg_m.z,
            0.0191163,
            1e-5,
            "OpenRocket's center of mass, aft",
        );
        close(m.cg_m.x, 0.017882, 1e-5, "OpenRocket's center of mass, out");
        let fine = |pieces| cockpit(pieces).0.planform.geometry().unwrap().area_m2;
        let (ours, curve) = (fine(64), fine(1024));
        assert!(
            ours < g.area_m2 && curve < ours,
            "{curve} < {ours} < {}",
            g.area_m2
        );
        close(ours, g.area_m2, 3.0e-4, "64 pieces against OpenRocket's 20");
        close(ours, curve, 3.5e-5, "64 pieces against the curve");
        // The slivers ROOT_SLIVER_SHARE's doc quotes, and the chord refused.
        let share = |pieces| {
            let (set, nose, fore) = cockpit(pieces);
            set.root_on_surface(&nose, fore).1 / set.planform.geometry().unwrap().area_m2
        };
        for (pieces, want) in [(20, 3.2e-4), (64, 3.1e-5), (1, 0.114)] {
            close(share(pieces), want, 0.1, "sliver share");
        }
        let (chord, nose, fore) = cockpit(1);
        let err = chord.root_radius_on(&nose, fore).unwrap_err();
        assert!(
            matches!(&err, DesignError::Geometry(m) if m.contains("cuts across")),
            "{err}"
        );
        // The span is the outline's highest point above the root leading edge.
        close(set.planform.span_m(), 0.006956521739130436, 1e-15, "span");
    }

    /// M4.5g4 (ADR-166): a fin set sits on a nose cone or a transition only with its root on the
    /// surface, within a micron, along the body and with no tab or fillet; each refusal is named.
    #[test]
    fn a_root_off_the_surface_is_refused_by_name() {
        // A straight cone's surface is straight: an outline that ends on it needs no root points.
        let cone =
            Profile::transition(crate::NoseShape::Conical {}, 0.2, 0.02, 0.04, false).unwrap();
        let rise = 0.02 / 0.2 * 0.05;
        let on_cone = |end_h: f64| {
            set(
                3,
                FinPlanform::Freeform {
                    points_m: vec![[0.0, 0.0], [0.03, 0.04], [0.05, end_h]],
                    root_m: Vec::new(),
                },
                0.003,
            )
        };
        close(
            on_cone(rise).root_radius_on(&cone, 0.1).unwrap(),
            0.03,
            1e-14,
            "radius at x = 0.1 m",
        );
        let says = |set: FinSet, fore: f64, what: &str| {
            let err = set.root_radius_on(&cone, fore).unwrap_err();
            assert!(
                matches!(&err, DesignError::Geometry(m) if m.contains(what)),
                "{what}: {err}"
            );
        };
        // Just inside and just outside the micron.
        assert!(on_cone(rise + 0.9e-6).root_radius_on(&cone, 0.1).is_ok());
        says(on_cone(rise + 1.1e-6), 0.1, "stands");
        says(on_cone(rise - 1.1e-6), 0.1, "stands");
        // A level root on a cone stands off it at its trailing edge.
        says(set(3, rectangle(0.05, 0.04), 0.003), 0.1, "stands");
        // Along the body, a micron's slack at either end.
        says(on_cone(rise), -2e-6, "must stay on it");
        says(on_cone(rise), 0.15 + 2e-6, "must stay on it");
        assert!(on_cone(rise).root_radius_on(&cone, 0.15).is_ok());
        let mut tabbed = on_cone(rise);
        tabbed.tab = Some(FinTab {
            height_m: 0.005,
            length_m: 0.02,
            offset_m: 0.01,
        });
        says(tabbed, 0.1, "tab or a fillet");
        let mut filleted = on_cone(rise);
        filleted.fillet = Some(FinFillet {
            radius_m: 0.005,
            material: ply(),
        });
        says(filleted, 0.1, "tab or a fillet");
    }

    #[test]
    fn bad_outlines_and_dimensions_are_rejected() {
        let bowtie = FinPlanform::Freeform {
            points_m: vec![[0.0, 0.0], [0.1, 0.1], [0.0, 0.1], [0.1, 0.0]],
            root_m: Vec::new(),
        };
        assert!(matches!(bowtie.validate(), Err(DesignError::Geometry(_))));
        let below = FinPlanform::Freeform {
            points_m: vec![[0.0, 0.0], [0.05, -0.01], [0.1, 0.0]],
            root_m: Vec::new(),
        };
        assert!(below.validate().is_err());
        // An outline that ends above the root leading edge has a root that rises: a planform, but
        // not one a body tube holds (ADR-166); the tree refuses it there (`tree.rs`'s tests).
        let open = FinPlanform::Freeform {
            points_m: vec![[0.0, 0.0], [0.05, 0.05], [0.1, 0.02]],
            root_m: Vec::new(),
        };
        assert!(open.validate().is_ok());
        let err = set(1, open, 0.003).check_level_root().unwrap_err();
        assert!(
            matches!(&err, DesignError::Geometry(m) if m.contains("needs a level root")),
            "{err}"
        );
        // Root points run fore, strictly between the root's ends, and stay at or above its
        // leading edge.
        let outline = vec![[0.0, 0.0], [0.05, 0.05], [0.1, 0.02]];
        let with_root = |root_m: Vec<[f64; 2]>| FinPlanform::Freeform {
            points_m: outline.clone(),
            root_m,
        };
        assert!(
            with_root(vec![[0.07, 0.012], [0.03, 0.004]])
                .validate()
                .is_ok()
        );
        for (root, says) in [
            (vec![[0.03, 0.004], [0.07, 0.012]], "must run fore"),
            (vec![[0.1, 0.02]], "must run fore"),
            (vec![[0.12, 0.02]], "must run fore"),
            (vec![[0.0, 0.0]], "must run fore"),
            (vec![[0.05, 0.06]], "runs below its root"),
            (vec![[0.07, 0.04], [0.03, 0.004]], "crosses itself"),
        ] {
            let err = with_root(root.clone()).validate().unwrap_err();
            assert!(
                matches!(&err, DesignError::Geometry(m) if m.contains(says)),
                "{root:?}: {err}"
            );
        }
        assert!(matches!(
            with_root(vec![[0.05, -0.001]]).validate(),
            Err(DesignError::Domain {
                what: "freeform fin span",
                ..
            })
        ));
        // The root leading edge is the origin, and the outline runs forward to aft.
        let shifted = FinPlanform::Freeform {
            points_m: vec![[0.03, 0.0], [0.08, 0.05], [0.13, 0.0]],
            root_m: Vec::new(),
        };
        assert!(shifted.validate().is_err());
        let backwards = FinPlanform::Freeform {
            points_m: vec![[0.0, 0.0], [-0.05, 0.05], [-0.1, 0.0]],
            root_m: Vec::new(),
        };
        assert!(backwards.validate().is_err());
        assert!(
            set(0, rectangle(0.1, 0.1), 0.003)
                .mass_properties(0.03)
                .is_err()
        );
        assert!(
            set(3, rectangle(0.1, 0.1), 0.0)
                .mass_properties(0.03)
                .is_err()
        );
        let mut fabric = set(3, rectangle(0.1, 0.1), 0.003);
        fabric.material = Material::surface("ripstop", 0.04);
        assert!(matches!(
            fabric.mass_properties(0.03),
            Err(DesignError::MaterialKind { .. })
        ));
    }

    /// `[A, S_x, S_xx, S_yy]` of one fillet's section by quadrature, column by column: from the
    /// fillet circle's tangent point on the body out to the body's radius, the section runs from
    /// the body's circle up to the fillet's; beyond, from the fin's plane up to the fillet's. The
    /// body's circle is integrated in `x = R_b − v²`, which takes out its vertical tangent.
    fn fillet_section_by_quadrature(rb: f64, r: f64) -> [f64; 4] {
        let tol = Tolerance {
            relative: 1e-12,
            absolute: 1e-24,
            ..Tolerance::default()
        };
        let c = (rb * rb + 2.0 * rb * r).sqrt();
        let joint = |x: f64| r - (r * r - (x - c) * (x - c)).max(0.0).sqrt();
        let body = |x: f64| (rb * rb - x * x).max(0.0).sqrt();
        let moments = |bottom: f64, top: f64, x: f64| {
            [
                top - bottom,
                x * (top - bottom),
                x * x * (top - bottom),
                (top.powi(3) - bottom.powi(3)) / 3.0,
            ]
        };
        let near = integrate(
            |v| {
                let x = rb - v * v;
                moments(body(x), joint(x), x).map(|m| m * 2.0 * v)
            },
            0.0,
            (rb - c * rb / (rb + r)).sqrt(),
            tol,
        )
        .unwrap();
        let far = integrate(|x| moments(0.0, joint(x), x), rb, c, tol).unwrap();
        std::array::from_fn(|k| near.value[k] + far.value[k])
    }

    #[test]
    fn fillet_section_is_its_region_by_quadrature() {
        for (rb, r) in [
            (0.05, 0.005),
            (0.05, 0.01),
            (0.05, 0.03),
            (0.1, 0.001),
            (0.02, 0.05),
        ] {
            let got = fillet_section(rb, r);
            let want = fillet_section_by_quadrature(rb, r);
            for (k, what) in ["A", "S_x", "S_xx", "S_yy"].into_iter().enumerate() {
                close(got[k], want[k], 1e-11, what);
            }
        }
        // No body, no fillet.
        assert_eq!(fillet_section(0.0, 0.01)[0], 0.0);
        // The series meets the difference where it hands over.
        let x: f64 = 0.25;
        close(
            span_less_sine(x.next_down()),
            x - x.sin(),
            1e-13,
            "x − sin x",
        );
        close(
            span_less_sine(1e-3),
            1e-9 / 6.0 - 1e-15 / 120.0 + 1e-21 / 5040.0,
            1e-15,
            "x − sin x",
        );
    }

    #[test]
    fn a_fillet_on_a_flat_body_is_a_square_less_a_quarter_circle() {
        // The curvature's share is of order `r/R_b`, 1e-8 here.
        let r = 0.01;
        close(
            fillet_section(1e6, r)[0],
            r * r * (1.0 - PI / 4.0),
            1e-6,
            "flat",
        );
    }

    /// The worked example in `docs/physics/mass.md`: a 5 mm fillet on a 30 mm tube is 4.253 mm²
    /// (the triangle 86.603 mm² less sectors of 64.506 and 17.843 mm²), and six of them along a
    /// 100 mm root in cardboard, 680 kg/m³, weigh 1.735 g.
    #[test]
    fn the_worked_fillet_example_is_the_docs() {
        let [a, ..] = fillet_section(0.030, 0.005);
        assert!((a * 1e6 - 4.253).abs() < 5e-4, "{}", a * 1e6);
        let planform = FinPlanform::Trapezoidal {
            root_chord_m: 0.1,
            tip_chord_m: 0.05,
            span_m: 0.05,
            sweep_m: 0.05,
        };
        let bare = set(3, planform.clone(), 0.003);
        let mut filleted = set(3, planform, 0.003);
        filleted.fillet = Some(FinFillet {
            radius_m: 0.005,
            material: Material::bulk("cardboard", 680.0),
        });
        let mass_g = (filleted.mass_properties(0.030).unwrap().mass_kg
            - bare.mass_properties(0.030).unwrap().mass_kg)
            * 1e3;
        assert!((mass_g - 1.735).abs() < 5e-4, "{mass_g}");
    }

    /// Far from real fillets the section's closed form cancels away: a fillet more than
    /// [`FILLET_RATIO_MAX`] times the body radius is refused, not weighed on noise; one under a
    /// millionth of it, or on a body of no radius, weighs nothing.
    #[test]
    fn a_fillet_far_from_the_body_s_size_is_refused_or_none() {
        let planform = FinPlanform::Trapezoidal {
            root_chord_m: 0.1,
            tip_chord_m: 0.05,
            span_m: 0.05,
            sweep_m: 0.05,
        };
        let filleted = |radius_m: f64| {
            let mut fins = set(3, planform.clone(), 0.003);
            fins.fillet = Some(FinFillet {
                radius_m,
                material: Material::bulk("epoxy", 1200.0),
            });
            fins
        };
        let bare = set(3, planform.clone(), 0.003);
        let rb = 0.05;
        let widest = FILLET_RATIO_MAX * rb;
        for radius_m in [widest * (1.0 + 1e-12), 1e10, 1e12, 1e100] {
            let error = filleted(radius_m).single_fin(rb).unwrap_err();
            assert!(
                matches!(error, DesignError::Domain { what, value }
                    if what.starts_with("fin fillet radius") && value == radius_m / rb),
                "{radius_m}: {error:?}"
            );
        }
        let mass = |fins: &FinSet, rb: f64| fins.single_fin(rb).unwrap().mass_kg;
        assert!(mass(&filleted(widest), rb) > mass(&bare, rb));
        assert!(mass(&filleted(1e-6 * rb), rb) > mass(&bare, rb));
        for (radius_m, body_m) in [(1e-6 * rb * (1.0 - 1e-12), rb), (1e-20, rb), (0.005, 0.0)] {
            assert_eq!(
                mass(&filleted(radius_m), body_m),
                mass(&bare, body_m),
                "{radius_m} on {body_m}"
            );
        }
        // A fillet of a weightless material weighs nothing, and is no error.
        let mut fins = set(3, planform, 0.003);
        fins.fillet = Some(FinFillet {
            radius_m: 0.005,
            material: Material::bulk("air", 0.0),
        });
        assert_eq!(
            fins.single_fin(0.05).unwrap().mass_kg,
            set(3, fins.planform.clone(), 0.003)
                .single_fin(0.05)
                .unwrap()
                .mass_kg
        );
    }

    #[test]
    fn fillets_are_a_prism_of_their_section_along_the_root() {
        // Moments by quadrature, assembled by the equations in `FinSet::fillets`' docs.
        let (rb, r, rho) = (0.05, 0.005, 1200.0);
        let planform = FinPlanform::Trapezoidal {
            root_chord_m: 0.1,
            tip_chord_m: 0.05,
            span_m: 0.05,
            sweep_m: 0.05,
        };
        let mut fins = set(3, planform, 0.003);
        let bare = fins.single_fin(rb).unwrap();
        fins.fillet = Some(FinFillet {
            radius_m: r,
            material: Material::bulk("epoxy", rho),
        });
        let with = fins.single_fin(rb).unwrap();
        let pair = with.without_part(&bare);
        let [a, sx, sxx, syy] = fillet_section_by_quadrature(rb, r);
        let l = 0.1;
        close(pair.mass_kg, 2.0 * rho * a * l, 1e-10, "mass");
        assert!((pair.cg_m - DVec3::new(sx / a, 0.0, -0.5 * l)).length() < 1e-12);
        let about_origin = with.inertia_about(DVec3::ZERO) - bare.inertia_about(DVec3::ZERO);
        let want = DMat3::from_cols(
            DVec3::new(
                2.0 * rho * (syy * l + a * l.powi(3) / 3.0),
                0.0,
                rho * sx * l * l,
            ),
            DVec3::new(0.0, 2.0 * rho * (sxx * l + a * l.powi(3) / 3.0), 0.0),
            DVec3::new(rho * sx * l * l, 0.0, 2.0 * rho * (sxx + syy) * l),
        );
        mat_close(about_origin, want, 1e-9, "fillets about the origin");
        // A fillet of no radius is none, and a negative one is refused.
        fins.fillet = Some(FinFillet {
            radius_m: 0.0,
            material: Material::bulk("epoxy", rho),
        });
        assert_eq!(fins.single_fin(rb).unwrap(), bare);
        fins.fillet = Some(FinFillet {
            radius_m: -0.001,
            material: Material::bulk("epoxy", rho),
        });
        assert!(matches!(
            fins.single_fin(rb),
            Err(DesignError::Domain {
                what: "fin fillet radius",
                ..
            })
        ));
    }

    /// A fin set's count is 1 to [`FinSet::MAX_COUNT`]: none and one over are refused by the
    /// count check itself (its `what` and the count), so four billion fins are an error rather
    /// than a loop that never returns; the bound itself is weighed, as 64 fins of one fin each.
    #[test]
    fn a_fin_count_over_the_bound_is_refused_and_the_bound_is_weighed() {
        let rb = 0.028;
        let planform = FinPlanform::Trapezoidal {
            root_chord_m: 0.05,
            tip_chord_m: 0.02,
            span_m: 0.03,
            sweep_m: 0.02,
        };
        for count in [0, FinSet::MAX_COUNT + 1, 4_000_000_000, u32::MAX] {
            let fins = set(count, planform.clone(), 0.003);
            for result in [fins.validate(), fins.mass_properties(rb).map(|_| ())] {
                match result {
                    Err(DesignError::Domain { what, value }) => {
                        assert_eq!(what, "fin count (1 to 64)", "{count}");
                        assert_eq!(value, f64::from(count), "{count}");
                    }
                    other => panic!("{count} fins: {other:?}"),
                }
            }
        }
        let most = set(FinSet::MAX_COUNT, planform.clone(), 0.003);
        most.validate().unwrap();
        let one = set(1, planform, 0.003).mass_properties(rb).unwrap();
        let all = most.mass_properties(rb).unwrap();
        close(all.mass_kg, 64.0 * one.mass_kg, 1e-13, "64 fins' mass");
        // Evenly spaced, the fins' center lies on the axis.
        assert!(
            all.cg_m.x.abs() < 1e-15 && all.cg_m.y.abs() < 1e-15,
            "{:?}",
            all.cg_m
        );
    }

    /// A tube fin set's count is 1 to [`TubeFinSet::MAX_COUNT`], refused by the count check
    /// itself above it, and weighed at it.
    #[test]
    fn a_tube_fin_count_over_the_bound_is_refused_and_the_bound_is_weighed() {
        let tubes = |count| TubeFinSet {
            count,
            length_m: 0.1,
            outer_radius_m: 0.012,
            thickness_m: 0.001,
            base_angle_rad: 0.0,
            material: Material::bulk("cardboard", 790.0),
        };
        for count in [0, TubeFinSet::MAX_COUNT + 1, 4_000_000_000, u32::MAX] {
            match tubes(count).mass_properties(0.02) {
                Err(DesignError::Domain { what, value }) => {
                    assert_eq!(what, "tube fin count (1 to 64)", "{count}");
                    assert_eq!(value, f64::from(count), "{count}");
                }
                other => panic!("{count} tubes: {other:?}"),
            }
        }
        let one = tubes(1).mass_properties(0.02).unwrap();
        let all = tubes(TubeFinSet::MAX_COUNT).mass_properties(0.02).unwrap();
        close(all.mass_kg, 64.0 * one.mass_kg, 1e-13, "64 tubes' mass");
    }
}