BREP_kernel 0.3.0

use crate::analytic_surface::{circumcenter, AnalyticSurface};
use crate::topology::{BrepSolid, CoedgeRecord, EdgeRecord, FaceRecord, VertexRecord};
use crate::{make_arc, KernelTolerances, NurbsCurve, NurbsSurface, Vec3};
use rustc_hash::FxHashMap as HashMap;

#[path = "step/pcurve.rs"]
mod pcurve;
use pcurve::{build_pcurve, EmittedCurve, EmittedFrame, EmittedSurface, Pcurve2d};

fn step_string(value: &str) -> String {
    value.replace('\'', "''")
}

fn real(value: f64) -> Result<String, String> {
    if !value.is_finite() {
        return Err(format!("export_step: non-finite number {value}"));
    }
    // Analytic frame axes come from cross products that can round to -0.0;
    // normalize so directions never print a negative zero component.
    if value == 0.0 {
        return Ok("0.".into());
    }
    if value.fract() == 0.0 && value.abs() < 1e15 {
        return Ok(format!("{value:.0}."));
    }
    let mut output = format!("{value:.15}");
    while output.ends_with('0') {
        output.pop();
    }
    if output.ends_with('.') {
        output.push('0');
    }
    if output == "-0.0" {
        output = "0.0".into();
    }
    Ok(output)
}

fn knot_runs(knots: &[f64]) -> (Vec<f64>, Vec<usize>) {
    let mut values = Vec::new();
    let mut multiplicities = Vec::new();
    for &knot in knots {
        if values
            .last()
            .is_some_and(|previous: &f64| (*previous - knot).abs() <= 1e-12)
        {
            *multiplicities.last_mut().unwrap() += 1;
        } else {
            values.push(knot);
            multiplicities.push(1);
        }
    }
    (values, multiplicities)
}

pub(crate) fn edge_subcurve(edge: &EdgeRecord) -> Result<NurbsCurve, String> {
    let [start, end] = edge.curve.domain()?;
    let epsilon = (1e-9 * (end - start)).max(2e-9);
    let mut curve = edge.curve.clone();
    if edge.t0 > start + epsilon && edge.t0 < end - epsilon {
        curve = curve.split(edge.t0)?.1;
    }
    let domain = curve.domain()?;
    if edge.t1 < domain[1] - epsilon && edge.t1 > domain[0] + epsilon {
        curve = curve.split(edge.t1)?.0;
    }
    Ok(curve)
}

#[derive(Default)]
struct StepWriter {
    lines: Vec<String>,
}

impl StepWriter {
    fn add(&mut self, body: impl Into<String>) -> usize {
        let id = self.lines.len() + 1;
        self.lines.push(format!("#{id}={};", body.into()));
        id
    }

    fn data(&self) -> String {
        self.lines.join("\n")
    }
}

fn write_point(writer: &mut StepWriter, point: Vec3) -> Result<usize, String> {
    Ok(writer.add(format!(
        "CARTESIAN_POINT('',({},{},{}))",
        real(point.x)?,
        real(point.y)?,
        real(point.z)?
    )))
}

fn id_list(ids: &[usize]) -> String {
    format!(
        "({})",
        ids.iter()
            .map(|id| format!("#{id}"))
            .collect::<Vec<_>>()
            .join(",")
    )
}

fn write_direction(writer: &mut StepWriter, direction: Vec3) -> Result<usize, String> {
    Ok(writer.add(format!(
        "DIRECTION('',({},{},{}))",
        real(direction.x)?,
        real(direction.y)?,
        real(direction.z)?
    )))
}

fn write_placement(
    writer: &mut StepWriter,
    origin: Vec3,
    axis: Vec3,
    ref_direction: Vec3,
) -> Result<usize, String> {
    let origin = write_point(writer, origin)?;
    let axis = write_direction(writer, axis)?;
    let ref_direction = write_direction(writer, ref_direction)?;
    Ok(writer.add(format!(
        "AXIS2_PLACEMENT_3D('',#{origin},#{axis},#{ref_direction})"
    )))
}

/// Emit the analytic AP214 surface entity (PLANE / CYLINDRICAL_SURFACE /
/// CONICAL_SURFACE / SPHERICAL_SURFACE / TOROIDAL_SURFACE) for a recognized
/// carrier, or `None` when the surface must stay a B-spline. The second value
/// reports whether the STEP-standard orientation of the emitted entity (plane
/// normal along the placement axis; revolution normal outward) is the REVERSE
/// of the stored NURBS orientation, so ADVANCED_FACE can invert `same_sense`
/// and the face normal survives the round trip. The third describes the
/// emitted entity's own (u,v) parameterization for the pcurve writer.
fn write_analytic_surface(
    writer: &mut StepWriter,
    surface: &NurbsSurface,
) -> Result<Option<(usize, bool, EmittedSurface)>, String> {
    let Some(analytic) = surface.analytic() else {
        return Ok(None);
    };
    match analytic {
        AnalyticSurface::Plane {
            origin,
            u_dir,
            v_dir,
            ..
        } => {
            // STEP planes are unbounded; the importer re-sizes the patch from
            // the face's edges, so only origin/normal/ref matter.
            let (Ok(normal), Ok(x_axis)) = (u_dir.cross(*v_dir).normalized(), u_dir.normalized())
            else {
                return Ok(None);
            };
            let placement = write_placement(writer, *origin, normal, x_axis)?;
            // The reader derives the plane's second parameter axis as
            // normal × ref_direction, so S(x,y) = origin + x·x̂ + y·ŷ with
            // those exact vectors — an orthonormal frame regardless of how the
            // stored patch scaled or sheared its own u_dir/v_dir.
            Ok(Some((
                writer.add(format!("PLANE('',#{placement})")),
                false,
                EmittedSurface::Plane {
                    origin: *origin,
                    x_axis,
                    y_axis: normal.cross(x_axis),
                },
            )))
        }
        AnalyticSurface::RuledRevolution {
            frame,
            rho0,
            rho1,
            height,
        } => {
            // The recognizer allows height < 0 (descending generatrix), whose
            // normal is the reverse of the standard outward convention the
            // importer reconstructs; report that so the face sense compensates.
            let flipped = *height < 0.0;
            let radius_scale = 1.0 + rho0.abs().max(rho1.abs());
            if (rho1 - rho0).abs() <= 1e-9 * radius_scale {
                if *rho0 <= 0.0 {
                    return Ok(None);
                }
                let placement = write_placement(writer, frame.origin, frame.axis, frame.x_axis)?;
                return Ok(Some((
                    writer.add(format!(
                        "CYLINDRICAL_SURFACE('',#{placement},{})",
                        real(*rho0)?
                    )),
                    flipped,
                    EmittedSurface::Cylinder {
                        frame: EmittedFrame {
                            origin: frame.origin,
                            x_axis: frame.x_axis,
                            y_axis: frame.y_axis,
                            axis: frame.axis,
                            azimuth_sign: 1.0,
                        },
                        radius: *rho0,
                    },
                )));
            }
            // Cone. The importer re-sizes the carrier from edge samples with a
            // 1e-4-scaled axial margin and clamps a negative extended-end
            // radius to zero — which BENDS the rebuilt slope when the apex sits
            // at an end of the face's axial range. Keep apex-touching cones as
            // exact NURBS instead of exporting a distorted carrier.
            let slope = (rho1 - rho0) / height;
            let apex_margin = 2.0 * (1e-4 * height.abs().max(1.0) + 1e-9) * slope.abs();
            if rho0.min(*rho1) <= apex_margin {
                return Ok(None);
            }
            // Orient the placement axis so the radius grows along +axis: STEP
            // semi-angles are positive. Radius at the placement origin stays
            // rho0 either way because the origin is on-axis at the v = 0 base.
            let axis = if slope >= 0.0 {
                frame.axis
            } else {
                frame.axis.scale(-1.0)
            };
            let placement = write_placement(writer, frame.origin, axis, frame.x_axis)?;
            Ok(Some((
                writer.add(format!(
                    "CONICAL_SURFACE('',#{placement},{},{})",
                    real(*rho0)?,
                    real(slope.abs().atan())?
                )),
                flipped,
                EmittedSurface::Cone {
                    frame: EmittedFrame {
                        origin: frame.origin,
                        x_axis: frame.x_axis,
                        // Reversing the placement axis reverses the derived
                        // second axis with it, so the emitted azimuth runs
                        // opposite the stored carrier's u.
                        y_axis: axis.cross(frame.x_axis),
                        axis,
                        azimuth_sign: if slope >= 0.0 { 1.0 } else { -1.0 },
                    },
                    radius: *rho0,
                    semi_angle: slope.abs().atan(),
                },
            )))
        }
        AnalyticSurface::Sphere { frame, radius } => {
            // Recognition template and importer reconstruction share the same
            // south-to-north meridian construction, so the rebuild is exact.
            let placement = write_placement(writer, frame.origin, frame.axis, frame.x_axis)?;
            Ok(Some((
                writer.add(format!(
                    "SPHERICAL_SURFACE('',#{placement},{})",
                    real(*radius)?
                )),
                false,
                EmittedSurface::Sphere {
                    frame: EmittedFrame {
                        origin: frame.origin,
                        x_axis: frame.x_axis,
                        y_axis: frame.y_axis,
                        axis: frame.axis,
                        azimuth_sign: 1.0,
                    },
                    radius: *radius,
                },
            )))
        }
        AnalyticSurface::Torus {
            frame,
            major_radius,
            minor_radius,
        } => {
            let placement = write_placement(writer, frame.origin, frame.axis, frame.x_axis)?;
            Ok(Some((
                writer.add(format!(
                    "TOROIDAL_SURFACE('',#{placement},{},{})",
                    real(*major_radius)?,
                    real(*minor_radius)?
                )),
                false,
                EmittedSurface::Torus {
                    frame: EmittedFrame {
                        origin: frame.origin,
                        x_axis: frame.x_axis,
                        y_axis: frame.y_axis,
                        axis: frame.axis,
                        azimuth_sign: 1.0,
                    },
                    major_radius: *major_radius,
                    minor_radius: *minor_radius,
                },
            )))
        }
        // No SURFACE_OF_REVOLUTION reader exists yet; general revolutions keep
        // their exact NURBS form.
        AnalyticSurface::Revolution { .. } => Ok(None),
    }
}

/// A curve recognized as an exact `make_arc` product: a circular arc of
/// `radius` about `axis`, starting at angle 0 on `x_axis` and travelling
/// counterclockwise through `sweep` — exactly the CIRCLE parameterization the
/// importer trims between the edge's vertices.
struct CircularArc {
    center: Vec3,
    axis: Vec3,
    x_axis: Vec3,
    /// axis × x_axis — the second placement axis a STEP reader derives, so the
    /// emitted parameterization is C(a) = center + r(cos a·x̂ + sin a·ŷ).
    y_axis: Vec3,
    radius: f64,
    /// Total swept angle: the emitted entity's parameter range is [0, sweep].
    sweep: f64,
    /// Rational-quadratic span count of the kernel arc this was recognized
    /// from — the bridge from the emitted ANGLE back to the kernel's own
    /// parameter, which a fitted pcurve on a B-spline carrier needs.
    spans: usize,
}

fn curve_scale(curve: &NurbsCurve) -> f64 {
    curve
        .control_points
        .iter()
        .map(|p| (p.x.abs() / p.w).max(p.y.abs() / p.w).max(p.z.abs() / p.w))
        .fold(0.0, f64::max)
}

/// Homogeneous-net equality to a scale-relative tolerance (the same
/// reconstruction contract analytic_surface.rs uses): matching nets mean the
/// curves are the SAME exact rational arc, not merely close.
fn curves_match(a: &NurbsCurve, b: &NurbsCurve, scale: f64) -> bool {
    if a.degree != b.degree
        || a.knots.len() != b.knots.len()
        || a.control_points.len() != b.control_points.len()
    {
        return false;
    }
    if a.knots
        .iter()
        .zip(&b.knots)
        .any(|(x, y)| (x - y).abs() > 1e-12)
    {
        return false;
    }
    let tolerance = 1e-9 * scale.max(1.0);
    a.control_points
        .iter()
        .zip(&b.control_points)
        .all(|(p, q)| {
            (p.x - q.x).abs() <= tolerance
                && (p.y - q.y).abs() <= tolerance
                && (p.z - q.z).abs() <= tolerance
                && (p.w - q.w).abs() <= 1e-9
        })
}

/// Recognition by exact reconstruction: extract a candidate circle from three
/// curve points, rebuild it with `make_arc`, and demand the identical net.
/// Split subranges of a circle (whose knots are no longer the pristine
/// make_arc pattern) are rejected and honestly stay NURBS.
fn recognize_circular_arc(curve: &NurbsCurve) -> Option<CircularArc> {
    if curve.degree != 2
        || curve.control_points.len() < 3
        || curve.control_points.len() % 2 == 0
        || (curve.control_points.len() - 1) / 2 > 4
    {
        return None;
    }
    let [t0, t1] = curve.domain().ok()?;
    let at = |fraction: f64| curve.evaluate(t0 + (t1 - t0) * fraction);
    // Three points at < 74% of the sweep apart, so consecutive pairs subtend
    // less than pi and the cross product below gives the travel direction.
    let p0 = at(0.0).ok()?;
    let pa = at(0.35).ok()?;
    let pb = at(0.7).ok()?;
    let center = circumcenter(p0, pa, pb)?;
    let radial = p0.sub(center);
    let radius = radial.length();
    let scale = curve_scale(curve);
    if radius <= 1e-9 * scale.max(1.0) {
        return None;
    }
    let x_axis = radial.scale(1.0 / radius);
    let axis = radial.cross(pa.sub(center)).normalized().ok()?;
    let y_axis = axis.cross(x_axis);
    let p_end = at(1.0).ok()?;
    let sweep = if p_end.sub(p0).length() <= 1e-9 * (1.0 + radius) {
        std::f64::consts::TAU
    } else {
        let closing = p_end.sub(center);
        let mut angle = closing.dot(y_axis).atan2(closing.dot(x_axis));
        if angle < 0.0 {
            angle += std::f64::consts::TAU;
        }
        angle
    };
    let rebuilt = make_arc(center, x_axis, y_axis, radius, 0.0, sweep).ok()?;
    curves_match(curve, &rebuilt, scale).then_some(CircularArc {
        center,
        axis,
        x_axis,
        y_axis,
        radius,
        sweep,
        spans: (curve.control_points.len() - 1) / 2,
    })
}

/// Emit LINE or CIRCLE for a recognized analytic edge curve (already oriented
/// start-to-end by `edge_subcurve`), or `None` for the B-spline fallback.
///
/// The second value describes the parameterization ISO 10303-42 gives the
/// entity that was written, because a pcurve on this edge has to share THAT
/// parameter — not the kernel's knot values (see `step/pcurve.rs`).
fn write_analytic_curve(
    writer: &mut StepWriter,
    curve: &NurbsCurve,
) -> Result<Option<(usize, EmittedCurve)>, String> {
    if curve.degree == 1
        && curve.control_points.len() == 2
        && curve
            .control_points
            .iter()
            .all(|control| (control.w - 1.0).abs() <= 1e-12)
    {
        let start = curve.control_points[0].point()?;
        let end = curve.control_points[1].point()?;
        let Ok(direction) = end.sub(start).normalized() else {
            return Ok(None);
        };
        let point = write_point(writer, start)?;
        let step_direction = write_direction(writer, direction)?;
        let vector = writer.add(format!(
            "VECTOR('',#{step_direction},{})",
            real(end.sub(start).length())?
        ));
        return Ok(Some((
            writer.add(format!("LINE('',#{point},#{vector})")),
            // The VECTOR carries the full chord length, so the STEP parameter
            // of this line is the fraction along the chord: C(s) = start +
            // s·(end − start) over [0, 1].
            EmittedCurve::Line { start, end },
        )));
    }
    if let Some(arc) = recognize_circular_arc(curve) {
        // ref_direction points at the edge's start vertex and the arc runs
        // counterclockwise about the axis, so the importer's vertex-trimmed
        // CCW rebuild reproduces the same directed curve with sense .T.
        let placement = write_placement(writer, arc.center, arc.axis, arc.x_axis)?;
        return Ok(Some((
            writer.add(format!("CIRCLE('',#{placement},{})", real(arc.radius)?)),
            EmittedCurve::Circle {
                center: arc.center,
                x_axis: arc.x_axis,
                y_axis: arc.y_axis,
                radius: arc.radius,
                sweep: arc.sweep,
                spans: arc.spans,
            },
        )));
    }
    Ok(None)
}

fn write_curve(writer: &mut StepWriter, curve: &NurbsCurve) -> Result<usize, String> {
    let points = curve
        .control_points
        .iter()
        .map(|control| write_point(writer, control.point()?))
        .collect::<Result<Vec<_>, _>>()?;
    write_bspline_curve(writer, curve, &points)
}

/// The B_SPLINE_CURVE_WITH_KNOTS entity (or its rational complex form) over
/// control points that have ALREADY been written — shared by the 3D edge
/// curves and by the 2D pcurves, whose CARTESIAN_POINTs carry two coordinates
/// instead of three but whose degree/knots/weights are written identically.
fn write_bspline_curve(
    writer: &mut StepWriter,
    curve: &NurbsCurve,
    points: &[usize],
) -> Result<usize, String> {
    let (knot_values, multiplicities) = knot_runs(&curve.knots);
    let multiplicities = format!(
        "({})",
        multiplicities
            .iter()
            .map(usize::to_string)
            .collect::<Vec<_>>()
            .join(",")
    );
    let knots = format!(
        "({})",
        knot_values
            .iter()
            .map(|value| real(*value))
            .collect::<Result<Vec<_>, _>>()?
            .join(",")
    );
    let rational = curve
        .control_points
        .iter()
        .any(|control| (control.w - 1.0).abs() > 1e-12);
    if !rational {
        return Ok(writer.add(format!(
            "B_SPLINE_CURVE_WITH_KNOTS('',{},{},.UNSPECIFIED.,.F.,.F.,{multiplicities},{knots},.UNSPECIFIED.)",
            curve.degree,
            id_list(points),
        )));
    }
    let weights = format!(
        "({})",
        curve
            .control_points
            .iter()
            .map(|control| real(control.w))
            .collect::<Result<Vec<_>, _>>()?
            .join(",")
    );
    Ok(writer.add(format!(
        "(BOUNDED_CURVE()B_SPLINE_CURVE({},{},.UNSPECIFIED.,.F.,.F.)\
         B_SPLINE_CURVE_WITH_KNOTS({multiplicities},{knots},.UNSPECIFIED.)\
         CURVE()GEOMETRIC_REPRESENTATION_ITEM()RATIONAL_B_SPLINE_CURVE({weights})\
         REPRESENTATION_ITEM(''))",
        curve.degree,
        id_list(points),
    )))
}

fn write_point_2d(writer: &mut StepWriter, point: [f64; 2]) -> Result<usize, String> {
    Ok(writer.add(format!(
        "CARTESIAN_POINT('',({},{}))",
        real(point[0])?,
        real(point[1])?
    )))
}

fn write_direction_2d(writer: &mut StepWriter, direction: [f64; 2]) -> Result<usize, String> {
    Ok(writer.add(format!(
        "DIRECTION('',({},{}))",
        real(direction[0])?,
        real(direction[1])?
    )))
}

/// Write the 2D geometry of one pcurve.
///
/// The `VECTOR` of a 2D LINE carries its TRUE magnitude rather than being
/// normalized to 1: ISO 10303-42 parameterizes a line as pnt + u·(magnitude ×
/// orientation), so the magnitude is what makes the 2D parameter equal the 3D
/// entity's parameter exactly.  (OpenCASCADE normalizes every VECTOR it writes
/// and then loses that agreement wherever the two speeds differ — its own cone
/// ruling pcurves are off by cos(semi-angle) — so this is strictly the more
/// faithful of the two conventions, and identical wherever the speeds agree.)
fn write_pcurve_geometry(writer: &mut StepWriter, curve: &Pcurve2d) -> Result<usize, String> {
    match curve {
        Pcurve2d::Line { point, vector } => {
            let magnitude = vector[0].hypot(vector[1]);
            if magnitude <= 0.0 {
                return Err("export_step: degenerate 2D line pcurve".into());
            }
            let point_id = write_point_2d(writer, *point)?;
            let direction =
                write_direction_2d(writer, [vector[0] / magnitude, vector[1] / magnitude])?;
            let vector_id = writer.add(format!(
                "VECTOR('',#{direction},{})",
                real(magnitude)?
            ));
            Ok(writer.add(format!("LINE('',#{point_id},#{vector_id})")))
        }
        Pcurve2d::Circle {
            center,
            ref_direction,
            radius,
        } => {
            let center_id = write_point_2d(writer, *center)?;
            let direction = write_direction_2d(writer, *ref_direction)?;
            let placement = writer.add(format!(
                "AXIS2_PLACEMENT_2D('',#{center_id},#{direction})"
            ));
            Ok(writer.add(format!("CIRCLE('',#{placement},{})", real(*radius)?)))
        }
        Pcurve2d::Spline(spline) => {
            let points = spline
                .control_points
                .iter()
                .map(|control| {
                    let point = control.point()?;
                    write_point_2d(writer, [point.x, point.y])
                })
                .collect::<Result<Vec<_>, String>>()?;
            write_bspline_curve(writer, spline, &points)
        }
    }
}

/// `PCURVE('', surface, DEFINITIONAL_REPRESENTATION('', (2d curve), ctx))` —
/// the association of one 2D curve with the surface it parameterizes.
fn write_pcurve_entity(
    writer: &mut StepWriter,
    surface_id: usize,
    context_2d: usize,
    curve: &Pcurve2d,
) -> Result<usize, String> {
    let geometry = write_pcurve_geometry(writer, curve)?;
    let representation = writer.add(format!(
        "DEFINITIONAL_REPRESENTATION('',(#{geometry}),#{context_2d})"
    ));
    Ok(writer.add(format!("PCURVE('',#{surface_id},#{representation})")))
}

fn write_surface(writer: &mut StepWriter, surface: &NurbsSurface) -> Result<usize, String> {
    let rows = surface
        .control_points
        .iter()
        .map(|row| {
            row.iter()
                .map(|control| write_point(writer, control.point()?))
                .collect::<Result<Vec<_>, _>>()
                .map(|ids| id_list(&ids))
        })
        .collect::<Result<Vec<_>, _>>()?;
    let grid = format!("({})", rows.join(","));
    let (u_values, u_multiplicities) = knot_runs(&surface.knots_u);
    let (v_values, v_multiplicities) = knot_runs(&surface.knots_v);
    let multiplicities = |values: &[usize]| {
        format!(
            "({})",
            values
                .iter()
                .map(usize::to_string)
                .collect::<Vec<_>>()
                .join(",")
        )
    };
    let knots = |values: &[f64]| -> Result<String, String> {
        Ok(format!(
            "({})",
            values
                .iter()
                .map(|value| real(*value))
                .collect::<Result<Vec<_>, _>>()?
                .join(",")
        ))
    };
    let u_mults = multiplicities(&u_multiplicities);
    let v_mults = multiplicities(&v_multiplicities);
    let u_knots = knots(&u_values)?;
    let v_knots = knots(&v_values)?;
    let rational = surface
        .control_points
        .iter()
        .flatten()
        .any(|control| (control.w - 1.0).abs() > 1e-12);
    if !rational {
        return Ok(writer.add(format!(
            "B_SPLINE_SURFACE_WITH_KNOTS('',{},{},{grid},.UNSPECIFIED.,.F.,.F.,.F.,\
             {u_mults},{v_mults},{u_knots},{v_knots},.UNSPECIFIED.)",
            surface.degree_u, surface.degree_v,
        )));
    }
    let weights = format!(
        "({})",
        surface
            .control_points
            .iter()
            .map(|row| {
                row.iter()
                    .map(|control| real(control.w))
                    .collect::<Result<Vec<_>, _>>()
                    .map(|values| format!("({})", values.join(",")))
            })
            .collect::<Result<Vec<_>, _>>()?
            .join(",")
    );
    Ok(writer.add(format!(
        "(BOUNDED_SURFACE()B_SPLINE_SURFACE({},{},{grid},.UNSPECIFIED.,.F.,.F.,.F.)\
         B_SPLINE_SURFACE_WITH_KNOTS({u_mults},{v_mults},{u_knots},{v_knots},.UNSPECIFIED.)\
         GEOMETRIC_REPRESENTATION_ITEM()RATIONAL_B_SPLINE_SURFACE({weights})\
         REPRESENTATION_ITEM('')SURFACE())",
        surface.degree_u, surface.degree_v,
    )))
}

fn write_length_unit(writer: &mut StepWriter, unit: &str) -> Result<usize, String> {
    let normalized = unit.to_lowercase();
    if normalized == "meter" || normalized == "metre" {
        return Ok(writer.add("(LENGTH_UNIT()NAMED_UNIT(*)SI_UNIT($,.METRE.))"));
    }
    if normalized == "centimeter" || normalized == "centimetre" {
        return Ok(writer.add("(LENGTH_UNIT()NAMED_UNIT(*)SI_UNIT(.CENTI.,.METRE.))"));
    }
    if matches!(normalized.as_str(), "micron" | "micrometer" | "micrometre") {
        return Ok(writer.add("(LENGTH_UNIT()NAMED_UNIT(*)SI_UNIT(.MICRO.,.METRE.))"));
    }
    if normalized == "inch" || normalized == "foot" {
        let metre = writer.add("(LENGTH_UNIT()NAMED_UNIT(*)SI_UNIT($,.METRE.))");
        let (factor, name) = if normalized == "inch" {
            (0.0254, "INCH")
        } else {
            (0.3048, "FOOT")
        };
        let measure = writer.add(format!(
            "LENGTH_MEASURE_WITH_UNIT(LENGTH_MEASURE({}),#{metre})",
            real(factor)?
        ));
        return Ok(writer.add(format!(
            "(CONVERSION_BASED_UNIT('{name}',#{measure})LENGTH_UNIT()NAMED_UNIT(*))"
        )));
    }
    Ok(writer.add("(LENGTH_UNIT()NAMED_UNIT(*)SI_UNIT(.MILLI.,.METRE.))"))
}

fn vertex_for(solid: &BrepSolid, id: u64) -> Result<&VertexRecord, String> {
    solid
        .vertices
        .iter()
        .find(|vertex| vertex.id == id)
        .ok_or_else(|| format!("export_step: missing vertex {id}"))
}

fn edge_for(solid: &BrepSolid, id: u64) -> Result<&EdgeRecord, String> {
    solid
        .edges
        .iter()
        .find(|edge| edge.id == id)
        .ok_or_else(|| format!("export_step: missing edge {id}"))
}

fn surface_key(face: &FaceRecord) -> usize {
    face as *const FaceRecord as usize
}

/// One use of an edge by a face's loop.  Curve-on-surface geometry is written
/// per USE — a pcurve names the surface it lives on — so the writer needs the
/// full adjacency of an edge before it can emit that edge.
struct CoedgeUse<'a> {
    surface_key: usize,
    face: &'a FaceRecord,
    coedge: &'a CoedgeRecord,
}

/// The AP214 document plus what the writer measured while producing it.
///
/// The pcurve counters are the honest half of the verify-or-omit contract in
/// `step/pcurve.rs`: a pcurve is written only when the emitted 2D geometry was
/// proved to reproduce the emitted 3D curve, and every one that could not be
/// proved is COUNTED here rather than guessed into the file.  A reader
/// reprojects an omitted pcurve exactly as it must reproject every pcurve in
/// the files this exporter wrote before curve-on-surface geometry existed.
#[derive(Clone, Debug, Default)]
pub struct StepExportReport {
    /// The Part 21 text.
    pub text: String,
    /// `PCURVE` entities written — at most one per coedge use of a
    /// non-degenerate edge.
    pub pcurves_written: usize,
    /// Coedge uses left without a pcurve because no candidate 2D geometry
    /// verified inside the export band.
    pub pcurves_omitted: usize,
    /// Edges written as `SURFACE_CURVE` (their two uses sit on two surfaces).
    pub surface_curves: usize,
    /// Edges written as `SEAM_CURVE` (both uses on ONE surface — a periodic
    /// seam, the case an importer otherwise has to re-detect geometrically).
    pub seam_curves: usize,
    /// Edges that kept a bare 3D curve for `edge_geometry`, because no pcurve
    /// survived (or a seam lost one of the pair it is required to carry).
    pub bare_curves: usize,
    /// Collapsed face boundaries written as `VERTEX_LOOP` — cone apexes and
    /// sphere poles, which this writer used to drop entirely.
    pub vertex_loops: usize,
    /// Largest verified ‖S(c₂d(s)) − c₃d(s)‖ among the pcurves actually
    /// written, in model units.
    pub max_pcurve_deviation: f64,
    /// Largest deviation among the BEST candidate for each OMITTED pcurve —
    /// how far the closest miss was, so an omission can be diagnosed as "a
    /// candidate applied and missed the band by this much" rather than only
    /// counted. Zero when nothing was omitted.
    ///
    /// It is infinite when ANY omission had no candidate proposed at all, which
    /// then hides the finite misses behind it — the two omission classes share
    /// one counter. Splitting them is a follow-up; the max is the useful half,
    /// because it is the finite value that says whether widening the band would
    /// have helped.
    pub worst_omitted_deviation: f64,
}

/// Serialize exact NURBS BREP topology as an AP214 STEP Part 21 document.
pub fn export_step(
    solids: &[BrepSolid],
    name: &str,
    unit: &str,
    timestamp: &str,
) -> Result<String, String> {
    export_step_report(solids, name, unit, timestamp).map(|report| report.text)
}

/// [`export_step`] plus the pcurve-coverage measurements behind the file.
pub fn export_step_report(
    solids: &[BrepSolid],
    name: &str,
    unit: &str,
    timestamp: &str,
) -> Result<StepExportReport, String> {
    if solids.is_empty() {
        return Err("export_step: at least one solid is required".into());
    }
    for solid in solids {
        let policy = KernelTolerances::for_solid(solid, 1e-7);
        let issues = solid.validate_with_tolerances(&KernelTolerances {
            pcurve_consistency: policy.export_knit,
            ..policy
        });
        if !issues.is_empty() {
            return Err(format!("export_step: invalid solid: {issues:?}"));
        }
    }
    let mut report = StepExportReport::default();
    let mut writer = StepWriter::default();
    let safe_name = step_string(name);
    let application = writer.add("APPLICATION_CONTEXT('automotive design')");
    writer.add(format!(
        "APPLICATION_PROTOCOL_DEFINITION('','automotive_design',2010,#{application})"
    ));
    let product_context = writer.add(format!("PRODUCT_CONTEXT('',#{application},'mechanical')"));
    let product = writer.add(format!(
        "PRODUCT('{safe_name}','{safe_name}','',(#{product_context}))"
    ));
    let formation = writer.add(format!("PRODUCT_DEFINITION_FORMATION('','',#{product})"));
    let definition_context = writer.add(format!(
        "PRODUCT_DEFINITION_CONTEXT('part definition',#{application},'design')"
    ));
    let definition = writer.add(format!(
        "PRODUCT_DEFINITION('design','',#{formation},#{definition_context})"
    ));
    let product_shape = writer.add(format!("PRODUCT_DEFINITION_SHAPE('','',#{definition})"));
    let length_unit = write_length_unit(&mut writer, unit)?;
    let angle_unit = writer.add("(NAMED_UNIT(*)PLANE_ANGLE_UNIT()SI_UNIT($,.RADIAN.))");
    let solid_angle_unit = writer.add("(NAMED_UNIT(*)SI_UNIT($,.STERADIAN.)SOLID_ANGLE_UNIT())");
    let uncertainty = writer.add(format!(
        "UNCERTAINTY_MEASURE_WITH_UNIT(LENGTH_MEASURE(1.E-6),#{length_unit},'distance_accuracy_value','')"
    ));
    let geometry_context = writer.add(format!(
        "(GEOMETRIC_REPRESENTATION_CONTEXT(3)\
         GLOBAL_UNCERTAINTY_ASSIGNED_CONTEXT((#{uncertainty}))\
         GLOBAL_UNIT_ASSIGNED_CONTEXT((#{length_unit},#{angle_unit},#{solid_angle_unit}))\
         REPRESENTATION_CONTEXT('',''))"
    ));
    // ONE parametric context for every DEFINITIONAL_REPRESENTATION in the file.
    // It carries no unit assignment, so the importer's file-scale search
    // (`derive_length_scale_mm`, which keys on GLOBAL_UNIT_ASSIGNED_CONTEXT /
    // LENGTH_UNIT) cannot mistake it for the geometric context.
    let parametric_context = writer.add(
        "(GEOMETRIC_REPRESENTATION_CONTEXT(2)\
         PARAMETRIC_REPRESENTATION_CONTEXT()\
         REPRESENTATION_CONTEXT('2D SPACE',''))",
    );
    let origin = write_point(&mut writer, Vec3::default())?;
    let direction_z = writer.add("DIRECTION('',(0.,0.,1.))");
    let direction_x = writer.add("DIRECTION('',(1.,0.,0.))");
    let axis = writer.add(format!(
        "AXIS2_PLACEMENT_3D('',#{origin},#{direction_z},#{direction_x})"
    ));

    let mut solid_ids = Vec::new();
    for solid in solids {
        // The same band the export gate above already held every stored pcurve
        // to, so a pcurve that survives here is no looser than the data it was
        // derived from.
        let band = KernelTolerances::for_solid(solid, 1e-7).export_knit;
        let mut vertex_ids = HashMap::<u64, usize>::default();
        let mut edge_ids = HashMap::<u64, usize>::default();
        let mut surfaces = HashMap::<usize, (usize, bool, EmittedSurface)>::default();
        for shell in &solid.shells {
            // Pass 1 — SURFACES. A pcurve references the surface entity it
            // parameterizes, so every surface id has to exist before the first
            // edge is written. (Before curve-on-surface geometry the writer
            // could emit surfaces after the loops; it no longer can.)
            for face in &shell.faces {
                let key = surface_key(face);
                if surfaces.contains_key(&key) {
                    continue;
                }
                let entry = match write_analytic_surface(&mut writer, &face.surface)? {
                    Some(triple) => triple,
                    None => (
                        write_surface(&mut writer, &face.surface)?,
                        false,
                        EmittedSurface::Spline,
                    ),
                };
                surfaces.insert(key, entry);
            }

            // Pass 2 — ADJACENCY. Which faces use each edge, in first-encounter
            // order so the emitted entity numbering stays deterministic.
            let mut edge_uses = HashMap::<u64, Vec<CoedgeUse>>::default();
            let mut edge_order: Vec<u64> = Vec::new();
            for face in &shell.faces {
                for loop_record in &face.loops {
                    for coedge in &loop_record.coedges {
                        let edge = edge_for(solid, coedge.edge_id)?;
                        if edge.degenerate {
                            continue;
                        }
                        let uses = edge_uses.entry(edge.id).or_default();
                        if uses.is_empty() {
                            edge_order.push(edge.id);
                        }
                        uses.push(CoedgeUse {
                            surface_key: surface_key(face),
                            face,
                            coedge,
                        });
                    }
                }
            }

            // Pass 3 — EDGES with their curve-on-surface geometry.
            for edge_id in &edge_order {
                if edge_ids.contains_key(edge_id) {
                    continue;
                }
                let edge = edge_for(solid, *edge_id)?;
                let subcurve = edge_subcurve(edge)?;
                let (curve, emitted_curve) = match write_analytic_curve(&mut writer, &subcurve)? {
                    Some(pair) => pair,
                    None => (
                        write_curve(&mut writer, &subcurve)?,
                        EmittedCurve::Spline { curve: subcurve },
                    ),
                };
                let uses = &edge_uses[edge_id];
                // Two uses on ONE surface is a periodic seam: the two halves of
                // the edge sit on opposite domain boundaries and the reader is
                // told so explicitly instead of having to rediscover it.
                let seam = uses.len() == 2 && uses[0].surface_key == uses[1].surface_key;
                let mut pcurves: Vec<(usize, Pcurve2d)> = Vec::new();
                let mut omitted = 0usize;
                if uses.len() == 2 {
                    // Slot order for a seam: the FORWARD-oriented coedge first,
                    // the reversed one second — the pairing OpenCASCADE writes
                    // and reads (BRep_CurveOnClosedSurface's PCurve1/PCurve2).
                    let mut ordered: Vec<&CoedgeUse> = uses.iter().collect();
                    if seam && !ordered[0].coedge.forward {
                        ordered.swap(0, 1);
                    }
                    for coedge_use in ordered {
                        // The stored pcurve runs in LOOP direction; the edge's
                        // geometry runs in EDGE direction. Fraction-synchronised
                        // reversal is exact, so a reversed coedge's pcurve is
                        // simply flipped before it is expressed.
                        let oriented = if coedge_use.coedge.forward {
                            coedge_use.coedge.pcurve.clone()
                        } else {
                            coedge_use.coedge.pcurve.reversed()?
                        };
                        let (surface_id, _, emitted_surface) = &surfaces[&coedge_use.surface_key];
                        let outcome = build_pcurve(
                            &coedge_use.face.surface,
                            emitted_surface,
                            &emitted_curve,
                            &oriented,
                            band,
                        )?;
                        match outcome.curve {
                            Some(curve_2d) => {
                                report.max_pcurve_deviation =
                                    report.max_pcurve_deviation.max(outcome.deviation);
                                pcurves.push((*surface_id, curve_2d));
                            }
                            None => {
                                omitted += 1;
                                if outcome.deviation.is_finite() {
                                    report.worst_omitted_deviation =
                                        report.worst_omitted_deviation.max(outcome.deviation);
                                } else {
                                    report.worst_omitted_deviation = f64::INFINITY;
                                }
                            }
                        }
                    }
                }
                // A SEAM_CURVE is required to carry BOTH pcurves on the one
                // surface, so a seam that lost one falls all the way back to a
                // bare 3D curve rather than to a half-described seam.
                if seam && pcurves.len() != 2 {
                    omitted += pcurves.len();
                    pcurves.clear();
                }
                report.pcurves_omitted += omitted;
                report.pcurves_written += pcurves.len();
                let geometry = if pcurves.is_empty() {
                    report.bare_curves += 1;
                    curve
                } else {
                    let ids = pcurves
                        .iter()
                        .map(|(surface_id, curve_2d)| {
                            write_pcurve_entity(
                                &mut writer,
                                *surface_id,
                                parametric_context,
                                curve_2d,
                            )
                        })
                        .collect::<Result<Vec<_>, String>>()?;
                    let keyword = if seam {
                        report.seam_curves += 1;
                        "SEAM_CURVE"
                    } else {
                        report.surface_curves += 1;
                        "SURFACE_CURVE"
                    };
                    // master_representation is .CURVE_3D.: our 3D curves are
                    // the exact authority the whole kernel and the export gate
                    // treat them as, and the pcurves above are verified AGAINST
                    // them. Promoting an approximation to master would be the
                    // one claim this writer must not make.
                    writer.add(format!(
                        "{keyword}('',#{curve},{},.CURVE_3D.)",
                        id_list(&ids)
                    ))
                };
                let start = vertex_step_id(&mut writer, &mut vertex_ids, solid, edge.start_vertex_id)?;
                let end = vertex_step_id(&mut writer, &mut vertex_ids, solid, edge.end_vertex_id)?;
                let step_id =
                    writer.add(format!("EDGE_CURVE('',#{start},#{end},#{geometry},.T.)"));
                edge_ids.insert(edge.id, step_id);
            }

            // Pass 4 — LOOPS and FACES over the ids the passes above fixed.
            let mut face_ids = Vec::new();
            for face in &shell.faces {
                let mut bound_ids = Vec::new();
                for (loop_index, loop_record) in face.loops.iter().enumerate() {
                    let mut oriented_edges = Vec::new();
                    for coedge in &loop_record.coedges {
                        let edge = edge_for(solid, coedge.edge_id)?;
                        if edge.degenerate {
                            continue;
                        }
                        let edge_id = *edge_ids
                            .get(&edge.id)
                            .ok_or_else(|| format!("export_step: unwritten edge {}", edge.id))?;
                        let orientation = if coedge.forward { ".T." } else { ".F." };
                        oriented_edges.push(
                            writer.add(format!("ORIENTED_EDGE('',*,*,#{edge_id},{orientation})")),
                        );
                    }
                    let kind = if loop_index == 0 {
                        "FACE_OUTER_BOUND"
                    } else {
                        "FACE_BOUND"
                    };
                    if oriented_edges.is_empty() {
                        // A bound left empty by the degenerate-edge skip is a
                        // COLLAPSED boundary — a cone apex, a sphere pole —
                        // and AP214 spells that `VERTEX_LOOP`, not nothing.
                        // Dropping it (what this writer did before) deletes the
                        // apex from the face's parameter domain, so an imported
                        // pointed cone lost its apex bound on re-export and
                        // every reader had to re-synthesise it from the
                        // surface's own degeneracy. A MIXED loop keeps the
                        // skip: that is OCCT's own write convention, both our
                        // importer and OCC's ShapeFix rebuild those, and the
                        // advanced-face conformance class permits vertex loops
                        // but not zero-length edge curves.
                        let Some(coedge) = loop_record.coedges.first() else {
                            continue;
                        };
                        let collapsed = edge_for(solid, coedge.edge_id)?;
                        let vertex = vertex_step_id(
                            &mut writer,
                            &mut vertex_ids,
                            solid,
                            collapsed.start_vertex_id,
                        )?;
                        let vertex_loop = writer.add(format!("VERTEX_LOOP('',#{vertex})"));
                        report.vertex_loops += 1;
                        bound_ids.push(writer.add(format!("{kind}('',#{vertex_loop},.T.)")));
                        continue;
                    }
                    let edge_loop =
                        writer.add(format!("EDGE_LOOP('',{})", id_list(&oriented_edges)));
                    bound_ids.push(writer.add(format!("{kind}('',#{edge_loop},.T.)")));
                }
                let (surface, flipped, _) = surfaces[&surface_key(face)];
                // `flipped` analytic entities are written with the reverse of
                // the stored NURBS orientation, so invert the flag to keep the
                // face normal identical through the round trip.
                let sense = if face.same_sense != flipped {
                    ".T."
                } else {
                    ".F."
                };
                face_ids.push(writer.add(format!(
                    "ADVANCED_FACE('',{},#{surface},{sense})",
                    id_list(&bound_ids)
                )));
            }
            let closed_shell = writer.add(format!("CLOSED_SHELL('',{})", id_list(&face_ids)));
            solid_ids.push(writer.add(format!(
                "MANIFOLD_SOLID_BREP('{safe_name}',#{closed_shell})"
            )));
        }
    }
    let mut items = vec![axis];
    items.extend(&solid_ids);
    let representation = writer.add(format!(
        "ADVANCED_BREP_SHAPE_REPRESENTATION('',{},#{geometry_context})",
        id_list(&items)
    ));
    writer.add(format!(
        "SHAPE_DEFINITION_REPRESENTATION(#{product_shape},#{representation})"
    ));
    let safe_timestamp = step_string(timestamp);
    let output = [
        "ISO-10303-21;".to_string(),
        "HEADER;".to_string(),
        "FILE_DESCRIPTION((''),'2;1');".to_string(),
        format!(
            "FILE_NAME('{safe_name}.step','{safe_timestamp}',(''),(''),'brep-kernel-rs','brep-kernel-rs','');"
        ),
        "FILE_SCHEMA(('AUTOMOTIVE_DESIGN { 1 0 10303 214 1 1 1 1 }'));".to_string(),
        "ENDSEC;".to_string(),
        "DATA;".to_string(),
        writer.data(),
        "ENDSEC;".to_string(),
        "END-ISO-10303-21;".to_string(),
        String::new(),
    ]
    .join("\n");
    let manifold_issues = audit_step_manifold(&output);
    if !manifold_issues.is_empty() {
        return Err(format!(
            "export_step: emitted AP214 manifold audit failed: {}",
            manifold_issues.join("; ")
        ));
    }
    let pcurve_issues = audit_step_pcurves(&output);
    if !pcurve_issues.is_empty() {
        return Err(format!(
            "export_step: emitted AP214 pcurve audit failed: {}",
            pcurve_issues.join("; ")
        ));
    }
    report.text = output;
    Ok(report)
}

/// VERTEX_POINT for a kernel vertex, written once per solid.
fn vertex_step_id(
    writer: &mut StepWriter,
    vertex_ids: &mut HashMap<u64, usize>,
    solid: &BrepSolid,
    id: u64,
) -> Result<usize, String> {
    if let Some(step_id) = vertex_ids.get(&id) {
        return Ok(*step_id);
    }
    let point = write_point(writer, vertex_for(solid, id)?.point)?;
    let step_id = writer.add(format!("VERTEX_POINT('',#{point})"));
    vertex_ids.insert(id, step_id);
    Ok(step_id)
}

/// Entity id -> body, for the serialized-text audits.  One entity per line is
/// this writer's own invariant, so the "parse" is a split.
fn step_entity_bodies(step: &str) -> HashMap<u64, &str> {
    step.lines()
        .filter_map(|line| {
            let rest = line.strip_prefix('#')?;
            let (digits, body) = rest.split_once('=')?;
            Some((
                digits.parse::<u64>().ok()?,
                body.trim_end().trim_end_matches(';'),
            ))
        })
        .collect()
}

/// Every `#N` reference in an entity body, in order.  For the entities audited
/// below that order IS the attribute order (`SURFACE_CURVE` puts its 3D curve
/// first and its pcurves after; `PCURVE` puts its surface first).
fn step_entity_refs(body: &str) -> Vec<u64> {
    let mut refs = Vec::new();
    let bytes = body.as_bytes();
    let mut index = 0;
    while index < bytes.len() {
        if bytes[index] == b'#' {
            let start = index + 1;
            let mut end = start;
            while end < bytes.len() && bytes[end].is_ascii_digit() {
                end += 1;
            }
            if end > start {
                if let Ok(id) = body[start..end].parse::<u64>() {
                    refs.push(id);
                }
            }
            index = end;
        } else {
            index += 1;
        }
    }
    refs
}

/// Audit the curve-on-surface half of the emitted graph, on the serialized
/// text, for the same reason `audit_step_manifold` exists: valid in-memory
/// intent is not proof of a correctly written file.
///
/// An `EDGE_CURVE` whose geometry is a bare 3D curve is NOT an issue — that is
/// the counted, deliberate outcome of the verify-or-omit gate. What is an
/// issue is a malformed bundle: a `SURFACE_CURVE` with no or more than two
/// associated geometries, a `SEAM_CURVE` that does not carry exactly two, a
/// `SEAM_CURVE` whose two pcurves name DIFFERENT surfaces (it is by definition
/// one surface's seam), a `SURFACE_CURVE` whose two pcurves name the SAME
/// surface (that is a seam, and must say so), or a `PCURVE` that does not
/// point at both a surface and a `DEFINITIONAL_REPRESENTATION`.
pub fn audit_step_pcurves(step: &str) -> Vec<String> {
    let bodies = step_entity_bodies(step);
    let mut issues = Vec::new();
    for (id, body) in &bodies {
        if !body.starts_with("EDGE_CURVE(") {
            continue;
        }
        let refs = step_entity_refs(body);
        let Some(geometry) = refs.get(2) else {
            issues.push(format!("EDGE_CURVE #{id} has no edge_geometry"));
            continue;
        };
        let Some(wrapper) = bodies.get(geometry) else {
            issues.push(format!("EDGE_CURVE #{id} references missing #{geometry}"));
            continue;
        };
        let seam = wrapper.starts_with("SEAM_CURVE(");
        if !seam && !wrapper.starts_with("SURFACE_CURVE(") {
            continue;
        }
        let wrapper_refs = step_entity_refs(wrapper);
        let pcurves = &wrapper_refs[wrapper_refs.len().min(1)..];
        if pcurves.is_empty() || pcurves.len() > 2 || (seam && pcurves.len() != 2) {
            issues.push(format!(
                "#{geometry} carries {} associated geometries",
                pcurves.len()
            ));
            continue;
        }
        let mut surfaces = Vec::new();
        for pcurve in pcurves {
            let Some(pcurve_body) = bodies.get(pcurve) else {
                issues.push(format!("#{geometry} references missing #{pcurve}"));
                continue;
            };
            if !pcurve_body.starts_with("PCURVE(") {
                issues.push(format!("#{geometry} associate #{pcurve} is not a PCURVE"));
                continue;
            }
            let pcurve_refs = step_entity_refs(pcurve_body);
            let representation = pcurve_refs.get(1).and_then(|id| bodies.get(id));
            if !representation.is_some_and(|body| body.starts_with("DEFINITIONAL_REPRESENTATION("))
            {
                issues.push(format!(
                    "PCURVE #{pcurve} has no DEFINITIONAL_REPRESENTATION"
                ));
            }
            if let Some(surface) = pcurve_refs.first() {
                surfaces.push(*surface);
            }
        }
        if surfaces.len() == 2 && (surfaces[0] == surfaces[1]) != seam {
            issues.push(format!(
                "#{geometry} pcurves name {} surface(s) but it is a {}",
                if surfaces[0] == surfaces[1] { 1 } else { 2 },
                if seam { "SEAM_CURVE" } else { "SURFACE_CURVE" }
            ));
        }
    }
    issues.sort();
    issues
}

/// Audit the serialized entity graph rather than assuming that valid
/// in-memory topology was necessarily written correctly.  Every EDGE_CURVE
/// in a closed shell must have exactly two ORIENTED_EDGE users with opposite
/// senses.
pub fn audit_step_manifold(step: &str) -> Vec<String> {
    let marker = "ORIENTED_EDGE('',*,*,#";
    let mut uses = HashMap::<u64, Vec<bool>>::default();
    for line in step.lines() {
        let Some(offset) = line.find(marker) else {
            continue;
        };
        let rest = &line[offset + marker.len()..];
        let digits = rest
            .chars()
            .take_while(|character| character.is_ascii_digit())
            .collect::<String>();
        let Ok(edge_id) = digits.parse::<u64>() else {
            continue;
        };
        let suffix = &rest[digits.len()..];
        let sense = suffix.starts_with(",.T.");
        uses.entry(edge_id).or_default().push(sense);
    }
    let mut issues = uses
        .into_iter()
        .filter_map(|(edge, senses)| {
            (senses.len() != 2 || senses[0] == senses[1]).then(|| {
                format!(
                    "EDGE_CURVE #{edge} has {} uses with senses {:?}",
                    senses.len(),
                    senses
                )
            })
        })
        .collect::<Vec<_>>();
    issues.sort();
    issues
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::{
        boolean_operation, import_step, make_box_brep, make_cone_brep, make_cylinder_brep,
        make_cylinder_surface, make_sphere_brep, make_torus_brep, solid_mass_properties,
        BooleanOperation, BooleanOptions,
    };

    #[test]
    fn box_step_contains_exact_manifold_topology() {
        let box_solid = make_box_brep(Vec3::default(), 2.0, 3.0, 4.0).unwrap();
        let step = export_step(&[box_solid], "box", "millimeter", "2026-07-27T00:00:00").unwrap();
        assert!(step.starts_with("ISO-10303-21;\nHEADER;"));
        assert!(audit_step_manifold(&step).is_empty());
        assert!(step.contains("MANIFOLD_SOLID_BREP('box'"));
        assert_eq!(step.matches("ADVANCED_FACE(").count(), 6);
        assert_eq!(step.matches("EDGE_CURVE(").count(), 12);
        assert!(step.ends_with("END-ISO-10303-21;\n"));
    }

    #[test]
    fn step_manifold_audit_rejects_single_and_same_sense_uses() {
        let single = "#1=ORIENTED_EDGE('',*,*,#9,.T.);";
        assert_eq!(audit_step_manifold(single).len(), 1);
        let same = "#1=ORIENTED_EDGE('',*,*,#9,.T.);\n\
                    #2=ORIENTED_EDGE('',*,*,#9,.T.);";
        assert_eq!(audit_step_manifold(same).len(), 1);
        let good = "#1=ORIENTED_EDGE('',*,*,#9,.T.);\n\
                    #2=ORIENTED_EDGE('',*,*,#9,.F.);";
        assert!(audit_step_manifold(good).is_empty());
    }

    #[test]
    fn unrecognized_surfaces_and_curves_still_write_rational_complex_entities() {
        // The all-NURBS fallback writers must stay intact for carriers no
        // analytic entity covers (general revolutions, split arc subranges).
        let mut writer = StepWriter::default();
        let cylinder =
            make_cylinder_surface(Vec3::default(), Vec3::new(0.0, 0.0, 1.0), 2.0, 5.0).unwrap();
        write_surface(&mut writer, &cylinder).unwrap();
        let split_arc = make_arc(
            Vec3::default(),
            Vec3::new(1.0, 0.0, 0.0),
            Vec3::new(0.0, 1.0, 0.0),
            2.0,
            0.0,
            std::f64::consts::TAU,
        )
        .unwrap()
        .split(0.37)
        .unwrap()
        .1;
        assert!(
            recognize_circular_arc(&split_arc).is_none(),
            "a split subrange is not the pristine make_arc net"
        );
        assert!(write_analytic_curve(&mut writer, &split_arc)
            .unwrap()
            .is_none());
        write_curve(&mut writer, &split_arc).unwrap();
        let data = writer.data();
        assert!(data.contains("RATIONAL_B_SPLINE_SURFACE"));
        assert!(data.contains("RATIONAL_B_SPLINE_CURVE"));
    }

    /// Export → assert the analytic entities appear (and, when the solid is
    /// fully analytic, that NO B-spline entity remains) → import → validate,
    /// match volume to 1e-6, and re-recognize every face's carrier.
    fn assert_analytic_round_trip(
        label: &str,
        original: &BrepSolid,
        expected_markers: &[&str],
        forbid_nurbs: bool,
    ) -> BrepSolid {
        let step = export_step(std::slice::from_ref(original), label, "millimeter", "fixed")
            .expect("export");
        for marker in expected_markers {
            assert!(step.contains(marker), "{label}: missing {marker}");
        }
        if forbid_nurbs {
            assert!(
                !step.contains("B_SPLINE"),
                "{label}: expected a fully analytic export"
            );
        }
        assert!(audit_step_manifold(&step).is_empty(), "{label}: audit");
        let imported = import_step(&step).expect("import");
        assert_eq!(imported.len(), 1, "{label}: one solid");
        let solid = imported.into_iter().next().unwrap();
        assert!(
            solid.validate().is_empty(),
            "{label}: imported solid invalid: {:?}",
            solid.validate()
        );
        let original_volume = solid_mass_properties(original).unwrap().volume;
        let volume = solid_mass_properties(&solid).unwrap().volume;
        let relative = ((volume - original_volume) / original_volume).abs();
        assert!(
            relative < 1e-6,
            "{label}: volume {volume} vs {original_volume} (rel {relative:.3e})"
        );
        for shell in &solid.shells {
            for face in &shell.faces {
                assert!(
                    face.surface.analytic().is_some(),
                    "{label}: imported face {} did not re-recognize as analytic",
                    face.id
                );
            }
        }
        solid
    }

    #[test]
    fn box_round_trips_through_plane_and_line_entities() {
        let solid = make_box_brep(Vec3::new(-1.0, 0.5, 2.0), 2.0, 3.0, 4.0).unwrap();
        let step = export_step(std::slice::from_ref(&solid), "box", "millimeter", "fixed").unwrap();
        assert_eq!(step.matches("PLANE(").count(), 6);
        // 12 three-dimensional edge lines plus the 24 two-dimensional pcurve
        // lines that now carry each edge on both of the planes it bounds.
        assert_eq!(step.matches("LINE(").count(), 36);
        assert_analytic_round_trip("box", &solid, &["PLANE(", "LINE("], true);
    }

    /// Curve-on-surface coverage on the case where every pcurve is exact: a
    /// box's edges are straight lines on planes, so each of the 12 edges goes
    /// out as a SURFACE_CURVE carrying the 2D LINE it traces on each of its
    /// two planes — no B-spline anywhere, nothing omitted.
    #[test]
    fn box_pcurves_cover_every_edge() {
        let solid = make_box_brep(Vec3::new(-1.0, 0.5, 2.0), 2.0, 3.0, 4.0).unwrap();
        let report =
            export_step_report(std::slice::from_ref(&solid), "box", "millimeter", "fixed").unwrap();
        assert_eq!(report.pcurves_written, 24);
        assert_eq!(report.pcurves_omitted, 0);
        assert_eq!(report.surface_curves, 12);
        assert_eq!(report.seam_curves, 0);
        assert_eq!(report.bare_curves, 0);
        assert_eq!(report.text.matches("SURFACE_CURVE(").count(), 12);
        assert_eq!(report.text.matches("PCURVE(").count(), 24);
        assert_eq!(
            report.text.matches("DEFINITIONAL_REPRESENTATION(").count(),
            24
        );
        // ONE parametric context shared by all 24 definitional representations.
        assert_eq!(
            report
                .text
                .matches("PARAMETRIC_REPRESENTATION_CONTEXT()")
                .count(),
            1
        );
        assert!(
            !report.text.contains("B_SPLINE"),
            "an all-planar solid must stay B-spline-free on both sides"
        );
        assert!(audit_step_pcurves(&report.text).is_empty());
        // Every plane frame here is axis aligned, so reconstructing a sampled
        // point through (p−o)·x̂, (p−o)·ŷ and back is exact in binary: the
        // measured worst deviation is 0.
        assert_eq!(report.max_pcurve_deviation, 0.0);
    }

    /// The same coverage on a solid with NO axis-aligned plane frame, so the
    /// verify gate reports a real floating-point residual rather than an
    /// exact-by-luck zero. Observed worst deviation 3.1e-15 mm; the band this
    /// is checked against is `export_knit` = 4e-3 mm, twelve orders of
    /// magnitude larger, so the assertion below is deliberately tight at 1e-13
    /// (32x the observation) instead of at the shipping band.
    #[test]
    fn rotated_box_pcurves_stay_exact_off_axis() {
        let (sin, cos) = 0.7_f64.sin_cos();
        let axis = Vec3::new(0.3, 0.7, 0.2).normalized().unwrap();
        let (ax, ay, az) = (axis.x, axis.y, axis.z);
        let one = 1.0 - cos;
        let rotation = crate::AffineTransform::new([
            cos + ax * ax * one,
            ax * ay * one - az * sin,
            ax * az * one + ay * sin,
            0.0,
            ay * ax * one + az * sin,
            cos + ay * ay * one,
            ay * az * one - ax * sin,
            0.0,
            az * ax * one - ay * sin,
            az * ay * one + ax * sin,
            cos + az * az * one,
            0.0,
            0.0,
            0.0,
            0.0,
            1.0,
        ])
        .unwrap();
        let solid = crate::transform_brep(
            &make_box_brep(Vec3::new(-1.0, 0.5, 2.0), 2.0, 3.0, 4.0).unwrap(),
            rotation,
            false,
        )
        .unwrap();
        let report =
            export_step_report(std::slice::from_ref(&solid), "tilted", "millimeter", "fixed")
                .unwrap();
        assert_eq!(report.pcurves_written, 24);
        assert_eq!(report.pcurves_omitted, 0);
        assert!(
            report.max_pcurve_deviation < 1e-13,
            "off-axis pcurve deviation {:.3e} exceeded 1e-13",
            report.max_pcurve_deviation
        );
        assert!(audit_step_pcurves(&report.text).is_empty());
        assert!(import_step(&report.text).is_ok());
    }

    /// Every edge of a cylinder and a frustum now carries curve-on-surface
    /// geometry, and every piece of it is EXACT: rim circles become 2D lines
    /// on the revolution (the STEP angle is the surface's own u) and 2D
    /// circles on the cap planes, the seam ruling becomes a pair of 2D lines
    /// on one surface. Nothing is omitted and no B-spline appears, so the
    /// analytic exports stay analytic on both sides.
    ///
    /// `frustum` shrinks along +axis and `frustum_growing` grows: the shrinking
    /// one is written with a REVERSED placement axis (STEP semi-angles are
    /// positive), which reverses the emitted azimuth against the stored
    /// carrier's u. Both are here because only the reversed one exercises that.
    #[test]
    fn revolution_carriers_cover_every_edge_exactly() {
        let axis = Vec3::new(0.0, 0.0, 1.0);
        for (label, solid) in [
            (
                "cylinder",
                make_cylinder_brep(Vec3::new(1.0, -2.0, 0.5), axis, 2.0, 5.0).unwrap(),
            ),
            (
                "cylinder_reversed",
                make_cylinder_brep(Vec3::new(0.0, 0.0, 5.0), Vec3::new(0.0, 0.0, -1.0), 2.0, 5.0)
                    .unwrap(),
            ),
            (
                "frustum",
                make_cone_brep(Vec3::new(0.5, 0.5, -1.0), axis, 3.0, 1.5, 5.0).unwrap(),
            ),
            (
                "frustum_growing",
                make_cone_brep(Vec3::new(0.5, 0.5, -1.0), axis, 1.5, 3.0, 5.0).unwrap(),
            ),
        ] {
            let report =
                export_step_report(std::slice::from_ref(&solid), label, "millimeter", "fixed")
                    .unwrap();
            // Two rim circles (two uses each) plus the seam ruling (two uses).
            assert_eq!(report.pcurves_written, 6, "{label}: written");
            assert_eq!(report.pcurves_omitted, 0, "{label}: omitted");
            assert_eq!(report.surface_curves, 2, "{label}: surface curves");
            assert_eq!(report.seam_curves, 1, "{label}: seam curves");
            assert_eq!(report.bare_curves, 0, "{label}: bare curves");
            // Both cap rims: a 2D CIRCLE apiece on the cap planes.
            assert_eq!(
                report.text.matches("AXIS2_PLACEMENT_2D(").count(),
                2,
                "{label}: 2D circles"
            );
            assert!(
                !report.text.contains("B_SPLINE"),
                "{label}: an analytic solid must export analytic on both sides"
            );
            assert!(audit_step_pcurves(&report.text).is_empty(), "{label}: audit");
            // Observed worst deviations: cylinder 8.9e-16, reversed 6.3e-16,
            // frustum 2.1e-15, growing frustum 1.9e-15 mm. The shipping band is
            // export_knit = 4e-3 mm; assert at 1e-12 (about 500x the worst
            // observation) so this stays a claim about exactness.
            assert!(
                report.max_pcurve_deviation < 1e-12,
                "{label}: pcurve deviation {:.3e} exceeded 1e-12",
                report.max_pcurve_deviation
            );
        }
    }

    /// A sphere's only non-degenerate edge is its pole-to-pole seam, used twice
    /// by the one spherical face. Both ends of that edge sit ON the axis, where
    /// the azimuth is undefined and the surface collapses — the sampler fills
    /// those from a neighbouring sample, and the result is still a pair of
    /// exact 2D lines at u = 0 and u = 2π.
    #[test]
    fn sphere_pole_and_seam() {
        let solid =
            make_sphere_brep(Vec3::new(2.0, 1.0, -1.0), 3.0, Vec3::new(0.0, 0.0, 1.0)).unwrap();
        let report =
            export_step_report(std::slice::from_ref(&solid), "sphere", "millimeter", "fixed")
                .unwrap();
        assert_eq!(report.seam_curves, 1);
        assert_eq!(report.pcurves_written, 2);
        assert_eq!(report.pcurves_omitted, 0);
        assert_eq!(report.text.matches("SEAM_CURVE(").count(), 1);
        assert!(!report.text.contains("B_SPLINE"));
        assert!(audit_step_pcurves(&report.text).is_empty());
        // Observed 4.768e-15 mm on a radius-3 sphere; band 1e-12 is 210x that.
        assert!(
            report.max_pcurve_deviation < 1e-12,
            "sphere pcurve deviation {:.3e}",
            report.max_pcurve_deviation
        );
    }

    /// A torus is closed in BOTH directions, so it has two seam edges and no
    /// ordinary surface curve at all — the case where a reader with no pcurves
    /// has the most to re-derive.
    #[test]
    fn torus_seams() {
        let solid =
            make_torus_brep(Vec3::new(0.0, 0.0, 1.0), Vec3::new(0.0, 0.0, 1.0), 5.0, 1.5).unwrap();
        let report =
            export_step_report(std::slice::from_ref(&solid), "torus", "millimeter", "fixed")
                .unwrap();
        assert_eq!(report.seam_curves, 2);
        assert_eq!(report.surface_curves, 0);
        assert_eq!(report.pcurves_written, 4);
        assert_eq!(report.pcurves_omitted, 0);
        assert!(!report.text.contains("B_SPLINE"));
        assert!(audit_step_pcurves(&report.text).is_empty());
        // Observed 1.592e-15 mm; band 1e-12 is 630x that.
        assert!(
            report.max_pcurve_deviation < 1e-12,
            "torus pcurve deviation {:.3e}",
            report.max_pcurve_deviation
        );
    }

    /// Every number in a body of a written entity, ignoring `#N` references.
    /// Only tokens containing a decimal point count, which is exactly what
    /// `real` emits and never an entity's dimension count or spline degree.
    fn step_numbers(body: &str) -> Vec<f64> {
        let characters: Vec<char> = body.chars().collect();
        let mut values = Vec::new();
        let mut index = 0;
        while index < characters.len() {
            if characters[index] == '#' {
                index += 1;
                while index < characters.len() && characters[index].is_ascii_digit() {
                    index += 1;
                }
                continue;
            }
            let signed = characters[index] == '-'
                && index + 1 < characters.len()
                && characters[index + 1].is_ascii_digit();
            if !characters[index].is_ascii_digit() && !signed {
                index += 1;
                continue;
            }
            let start = index;
            if signed {
                index += 1;
            }
            while index < characters.len() && characters[index].is_ascii_digit() {
                index += 1;
            }
            if index < characters.len() && characters[index] == '.' {
                index += 1;
                while index < characters.len() && characters[index].is_ascii_digit() {
                    index += 1;
                }
                values.push(
                    characters[start..index]
                        .iter()
                        .collect::<String>()
                        .parse()
                        .expect("number"),
                );
            }
        }
        values
    }

    /// Read the emitted file back and re-evaluate it against ISO 10303-42's own
    /// formulas — written HERE, from the standard, not shared with the writer —
    /// to prove the 2D and 3D geometry of an edge agree at EQUAL parameters.
    ///
    /// This is the check the internal verify-or-omit gate cannot make about
    /// itself: if the writer's evaluators were self-consistently wrong, the gate
    /// would pass every candidate and this test would fail.
    #[test]
    fn pcurve_parameter_shared_with_analytic_curve() {
        let solid =
            make_cylinder_brep(Vec3::default(), Vec3::new(0.0, 0.0, 1.0), 2.0, 5.0).unwrap();
        let step =
            export_step(std::slice::from_ref(&solid), "cylinder", "millimeter", "fixed").unwrap();
        let bodies = step_entity_bodies(&step);
        let body = |id: u64| -> &str { bodies[&id] };
        let refs = |id: u64| step_entity_refs(body(id));
        let numbers = |id: u64| step_numbers(body(id));
        let vector3 = |id: u64| {
            let values = numbers(id);
            Vec3::new(values[0], values[1], values[2])
        };

        // The cylindrical surface: placement (origin, axis, ref) + radius.
        let surface_id = bodies
            .iter()
            .find(|(_, text)| text.starts_with("CYLINDRICAL_SURFACE("))
            .map(|(id, _)| *id)
            .expect("cylindrical surface");
        let radius = numbers(surface_id)[0];
        let placement = refs(refs(surface_id)[0]);
        let origin = vector3(placement[0]);
        let axis = vector3(placement[1]);
        let x_axis = vector3(placement[2]);
        let y_axis = axis.cross(x_axis);
        // ISO 10303-42 cylindrical_surface.
        let evaluate = |u: f64, v: f64| {
            origin
                .add(x_axis.scale(radius * u.cos()))
                .add(y_axis.scale(radius * u.sin()))
                .add(axis.scale(v))
        };

        let mut checked = 0;
        let mut worst: f64 = 0.0;
        for (id, text) in &bodies {
            if !text.starts_with("SURFACE_CURVE(") && !text.starts_with("SEAM_CURVE(") {
                continue;
            }
            let bundle = refs(*id);
            // The emitted 3D curve: LINE(pnt, VECTOR(dir, magnitude)) with
            // C(s) = pnt + s·magnitude·dir over [0,1], or CIRCLE(placement, r)
            // with C(a) = c + r(cos a·x̂ + sin a·ŷ) over [0, 2π].
            let curve_id = bundle[0];
            let curve_body = body(curve_id);
            let (domain, curve_3d): (f64, Box<dyn Fn(f64) -> Vec3>) =
                if curve_body.starts_with("LINE(") {
                    let parts = refs(curve_id);
                    let start = vector3(parts[0]);
                    let vector = refs(parts[1]);
                    let magnitude = numbers(parts[1])[0];
                    let direction = vector3(vector[0]);
                    (1.0, Box::new(move |s| start.add(direction.scale(magnitude * s))))
                } else {
                    let arc = refs(refs(curve_id)[0]);
                    let radius = numbers(curve_id)[0];
                    let center = vector3(arc[0]);
                    let arc_axis = vector3(arc[1]);
                    let arc_x = vector3(arc[2]);
                    let arc_y = arc_axis.cross(arc_x);
                    (
                        std::f64::consts::TAU,
                        Box::new(move |a: f64| {
                            center
                                .add(arc_x.scale(radius * a.cos()))
                                .add(arc_y.scale(radius * a.sin()))
                        }),
                    )
                };
            for pcurve_id in &bundle[1..] {
                let pcurve = refs(*pcurve_id);
                if pcurve[0] != surface_id {
                    continue; // the cap-plane half; this test grades the cylinder
                }
                let geometry = refs(refs(*pcurve_id)[1])[0];
                assert!(body(geometry).starts_with("LINE("), "iso-lines only");
                let parts = refs(geometry);
                let point = numbers(parts[0]);
                let magnitude = numbers(parts[1])[0];
                let direction = numbers(refs(parts[1])[0]);
                for step_index in 0..=40 {
                    let s = domain * step_index as f64 / 40.0;
                    let u = point[0] + s * magnitude * direction[0];
                    let v = point[1] + s * magnitude * direction[1];
                    worst = worst.max(evaluate(u, v).sub(curve_3d(s)).length());
                }
                checked += 1;
            }
        }
        assert_eq!(checked, 4, "two rim circles and both halves of the seam");
        // Observed 8.882e-16 mm on a radius-2, height-5 cylinder. The plan's
        // bar was 1e-9·scale = 5e-9; this asserts 1e-12, three orders tighter
        // and still ~1000x the observation.
        assert!(
            worst < 1e-12,
            "independent re-evaluation deviated {worst:.3e} mm"
        );
    }

    /// A pointed cone's wall is a full-revolution B-spline patch whose two
    /// generatrix coedges are the SAME edge, so it exercises the seam path:
    /// both uses land on one surface and the edge goes out as a SEAM_CURVE
    /// carrying both pcurves, forward slot first.
    ///
    /// The rim on that B-spline wall is the case a REMAP cannot serve and a
    /// RESAMPLE can: the stored pcurve is synchronised to the kernel's
    /// rational-arc parameter and the emitted circle to the angle, so the
    /// stored net may not be reused, but reading it at
    /// `circle_angle_to_parameter(angle)` — the exact inverse of the `make_arc`
    /// construction the recognizer matched — gives the right point for every
    /// emitted parameter, and the fit through those points carries it.
    #[test]
    fn pointed_cone_seam_exports_seam_curve_with_both_pcurves() {
        let solid =
            make_cone_brep(Vec3::default(), Vec3::new(0.0, 0.0, 1.0), 3.0, 0.0, 5.0).unwrap();
        let report =
            export_step_report(std::slice::from_ref(&solid), "cone", "millimeter", "fixed")
                .unwrap();
        assert_eq!(report.seam_curves, 1);
        assert_eq!(report.surface_curves, 1);
        // Both uses of the seam and both uses of the rim: full coverage.
        assert_eq!(report.pcurves_written, 4);
        assert_eq!(report.pcurves_omitted, 0);
        assert_eq!(report.bare_curves, 0);
        assert_eq!(report.text.matches("SEAM_CURVE(").count(), 1);
        assert!(audit_step_pcurves(&report.text).is_empty());
        // The two seam pcurves are exact (8.6e-16 mm); the fitted rim on the
        // B-spline wall dominates at 3.659e-5 mm after refinement — it started
        // at 2.393e-3 with a single 33-sample fit, which is inside the 4e-3
        // band but only by 1.7x, so the refinement loop earns its keep here.
        // Band 1e-4 is 2.7x the observation; the shipping gate is 40x looser.
        assert!(
            report.max_pcurve_deviation < 1e-4,
            "pcurve deviation {:.3e} exceeded 1e-4",
            report.max_pcurve_deviation
        );
        let imported = import_step(&report.text).expect("import");
        let volume = solid_mass_properties(&imported[0]).unwrap().volume;
        let expected = solid_mass_properties(&solid).unwrap().volume;
        assert!(((volume - expected) / expected).abs() < 1e-6);
    }

    /// Boolean and fillet solids: the edges that are neither iso-lines nor rims.
    ///
    /// `box_minus_cylinder` reaches full coverage only through the fit — one of
    /// its two hole rims runs CLOCKWISE in its cap plane's parameters, and a 2D
    /// `AXIS2_PLACEMENT_2D` is right-handed by definition, so the exact 2D
    /// `CIRCLE` for that one is the MIRROR of the truth and misses by the full
    /// diameter (measured: 3.000e0 mm, refused). The fit carries it instead.
    /// The filleted box has no analytic answer for any of its blend edges.
    ///
    /// "0 omitted" is a coverage claim, not a correctness one — correctness is
    /// the writer's gate, which measured every one of these against the emitted
    /// entities before writing it, `audit_step_pcurves` on the serialized text,
    /// and (for the analytic half) the independent re-evaluation in
    /// `pcurve_parameter_shared_with_analytic_curve`. A fitted 2D B-spline has
    /// no independent evaluator here; `step-validation/export-oracle.mjs` is
    /// where OpenCASCADE grades those.
    #[test]
    fn boolean_and_fillet_solids_reach_full_pcurve_coverage() {
        let axis = Vec3::new(0.0, 0.0, 1.0);
        let block = make_box_brep(Vec3::new(-3.0, -3.0, 0.0), 6.0, 6.0, 4.0).unwrap();
        let drill = make_cylinder_brep(Vec3::new(0.0, 0.0, -1.0), axis, 1.5, 6.0).unwrap();
        let cut = boolean_operation(
            &block,
            &drill,
            BooleanOperation::Subtract,
            &BooleanOptions::default(),
        )
        .unwrap();
        let report =
            export_step_report(std::slice::from_ref(&cut), "cut", "millimeter", "fixed").unwrap();
        assert_eq!(report.pcurves_omitted, 0, "boolean: omitted");
        assert_eq!(report.bare_curves, 0, "boolean: bare");
        assert_eq!(report.pcurves_written, 30, "boolean: written");
        assert_eq!(report.seam_curves, 1, "boolean: the drill's seam");
        assert!(audit_step_pcurves(&report.text).is_empty());
        // Observed 3.885e-6 mm, dominated by the fitted clockwise rim (a single
        // 33-sample fit gave 6.062e-5 before refinement). Band 1e-5 is 2.6x the
        // observation; the shipping gate is 4e-3, 1000x looser.
        assert!(
            report.max_pcurve_deviation < 1e-5,
            "boolean pcurve deviation {:.3e} exceeded 1e-5",
            report.max_pcurve_deviation
        );

        let plain = make_box_brep(Vec3::default(), 6.0, 6.0, 4.0).unwrap();
        let rounded = crate::fillet_edges(&plain, &[Vec3::new(0.0, 0.0, 2.0)], None, 0.8, false, None)
            .expect("fillet");
        let report =
            export_step_report(std::slice::from_ref(&rounded), "fillet", "millimeter", "fixed")
                .unwrap();
        assert_eq!(report.pcurves_omitted, 0, "fillet: omitted");
        assert_eq!(report.pcurves_written, 30, "fillet: written");
        assert!(audit_step_pcurves(&report.text).is_empty());
        // Observed 1.264e-7 mm; band 1e-6 is 8x that.
        assert!(
            report.max_pcurve_deviation < 1e-6,
            "fillet pcurve deviation {:.3e} exceeded 1e-6",
            report.max_pcurve_deviation
        );
        assert!(import_step(&report.text).is_ok());
    }

    /// A boundary collapsed to a single point survives a re-export.
    ///
    /// `abc_00000036` is a vendor file that uses `VERTEX_LOOP`, and our
    /// importer turns one into a face bound holding a single degenerate edge
    /// (`builder/collect.rs`). Its two solids come back carrying two such
    /// bounds between them. Before this the exporter dropped ANY bound whose
    /// coedges were all degenerate, so those collapsed boundaries vanished
    /// from their faces' parameter domains and the next reader had to
    /// re-synthesise them from the surfaces' own degeneracy. They now go back
    /// out as the `VERTEX_LOOP`s they came in as.
    ///
    /// (The plan proposed a pointed-cone fixture for this. There isn't one:
    /// OpenCASCADE writes a pointed cone with the apex as an ordinary vertex
    /// where the seam meets itself, no vertex loop, and our own
    /// `make_cone_brep` puts the apex edge inside a MIXED loop, which keeps the
    /// skip by design. An imported vendor file is the only place the case
    /// arises, so that is what this tests.)
    #[test]
    fn vertex_loops_survive_a_re_export() {
        let text = include_str!(concat!(
            env!("CARGO_MANIFEST_DIR"),
            "/tests/fixtures/step-import/abc_00000036.step"
        ));
        let solids = import_step(text).expect("import");
        assert_eq!(solids.len(), 2);
        let report =
            export_step_report(&solids, "abc_00000036", "millimeter", "fixed").expect("export");
        // Both collapsed bounds are written, and the entity count in the file
        // matches the counter rather than being asserted independently of it.
        assert_eq!(report.vertex_loops, 2);
        assert_eq!(report.text.matches("VERTEX_LOOP(").count(), 2);
        assert!(audit_step_manifold(&report.text).is_empty());
        assert!(audit_step_pcurves(&report.text).is_empty());
        let reimported = import_step(&report.text).expect("re-import");
        assert_eq!(reimported.len(), 2);
        for (index, (before, after)) in solids.iter().zip(&reimported).enumerate() {
            assert!(
                after.validate().is_empty(),
                "solid {index} invalid after re-import: {:?}",
                after.validate()
            );
            let original = solid_mass_properties(before).unwrap().volume;
            let volume = solid_mass_properties(after).unwrap().volume;
            let relative = ((volume - original) / original).abs();
            assert!(
                relative < 1e-6,
                "solid {index} volume {volume} vs {original} (rel {relative:.3e})"
            );
        }
    }

    /// The pcurve audit reads the serialized graph, so it has to catch bundles
    /// this writer would never produce as well as the ones it does.
    #[test]
    fn step_pcurve_audit_rejects_malformed_bundles() {
        // The audit's own contract: one entity per line, unindented, exactly
        // as `StepWriter` emits them.
        let bundle = |surface_curve: &str, pcurve: &str, representation: &str| {
            [
                "#1=PLANE('',#9);",
                representation,
                pcurve,
                surface_curve,
                "#5=EDGE_CURVE('',#10,#11,#4,.T.);",
            ]
            .join("\n")
        };
        let representation = "#2=DEFINITIONAL_REPRESENTATION('',(#8),#7);";
        let pcurve = "#3=PCURVE('',#1,#2);";
        let good = bundle(
            "#4=SURFACE_CURVE('',#6,(#3),.CURVE_3D.);",
            pcurve,
            representation,
        );
        assert!(audit_step_pcurves(&good).is_empty());
        // A SEAM_CURVE must carry exactly two pcurves.
        let short_seam = bundle(
            "#4=SEAM_CURVE('',#6,(#3),.CURVE_3D.);",
            pcurve,
            representation,
        );
        assert_eq!(audit_step_pcurves(&short_seam).len(), 1);
        // Two pcurves naming ONE surface is a seam, and must say so.
        let mislabelled = bundle(
            "#4=SURFACE_CURVE('',#6,(#3,#3),.CURVE_3D.);",
            pcurve,
            representation,
        );
        assert_eq!(audit_step_pcurves(&mislabelled).len(), 1);
        // A PCURVE with no DEFINITIONAL_REPRESENTATION carries no geometry.
        let no_representation = bundle(
            "#4=SURFACE_CURVE('',#6,(#3),.CURVE_3D.);",
            pcurve,
            "#2=REPRESENTATION('',(#8),#7);",
        );
        assert_eq!(audit_step_pcurves(&no_representation).len(), 1);
        // An EDGE_CURVE that kept a bare 3D curve is the counted, deliberate
        // outcome of the omission gate — never an audit issue.
        let bare = "#1=LINE('',#2,#3);\n#5=EDGE_CURVE('',#10,#11,#1,.T.);";
        assert!(audit_step_pcurves(bare).is_empty());
    }


    #[test]
    fn cylinder_round_trips_through_analytic_entities() {
        let solid = make_cylinder_brep(
            Vec3::new(1.0, -2.0, 0.5),
            Vec3::new(0.0, 0.0, 1.0),
            2.0,
            5.0,
        )
        .unwrap();
        assert_analytic_round_trip(
            "cylinder",
            &solid,
            &[
                "CYLINDRICAL_SURFACE(",
                "PLANE(",
                "CIRCLE(",
                "LINE(",
                "VECTOR(",
            ],
            true,
        );
    }

    #[test]
    fn cylinder_export_keeps_unit_conversion_entities() {
        let cylinder =
            make_cylinder_brep(Vec3::default(), Vec3::new(0.0, 0.0, 1.0), 2.0, 5.0).unwrap();
        let step = export_step(&[cylinder], "cylinder", "inch", "fixed").unwrap();
        assert!(step.contains("CYLINDRICAL_SURFACE("));
        assert!(step.contains("CONVERSION_BASED_UNIT('INCH'"));
    }

    #[test]
    fn frustum_round_trips_through_conical_surface() {
        let solid = make_cone_brep(
            Vec3::new(0.5, 0.5, -1.0),
            Vec3::new(0.0, 0.0, 1.0),
            3.0,
            1.5,
            5.0,
        )
        .unwrap();
        assert_analytic_round_trip("frustum", &solid, &["CONICAL_SURFACE(", "PLANE("], true);
    }

    #[test]
    fn pointed_cone_wall_stays_nurbs_but_caps_and_rim_export_analytic() {
        // The importer's axial cover margin clamps a negative apex-end radius
        // to zero, bending the rebuilt slope; an apex-touching CONICAL export
        // would round-trip with ~1e-4 relative volume error, so the wall must
        // honestly stay NURBS while the cap plane and rim circle go analytic.
        let solid =
            make_cone_brep(Vec3::default(), Vec3::new(0.0, 0.0, 1.0), 3.0, 0.0, 5.0).unwrap();
        let step =
            export_step(std::slice::from_ref(&solid), "cone", "millimeter", "fixed").unwrap();
        assert!(!step.contains("CONICAL_SURFACE("));
        assert!(step.contains("B_SPLINE_SURFACE"));
        assert!(step.contains("PLANE("));
        assert!(step.contains("CIRCLE("));
        let imported = import_step(&step).expect("import");
        let volume = solid_mass_properties(&imported[0]).unwrap().volume;
        let expected = solid_mass_properties(&solid).unwrap().volume;
        assert!(((volume - expected) / expected).abs() < 1e-6);
    }

    #[test]
    fn sphere_round_trips_through_spherical_surface() {
        let solid =
            make_sphere_brep(Vec3::new(2.0, 1.0, -1.0), 3.0, Vec3::new(0.0, 0.0, 1.0)).unwrap();
        assert_analytic_round_trip("sphere", &solid, &["SPHERICAL_SURFACE(", "CIRCLE("], true);
    }

    #[test]
    fn torus_round_trips_through_toroidal_surface() {
        let solid =
            make_torus_brep(Vec3::new(0.0, 0.0, 1.0), Vec3::new(0.0, 0.0, 1.0), 5.0, 1.5).unwrap();
        assert_analytic_round_trip("torus", &solid, &["TOROIDAL_SURFACE(", "CIRCLE("], true);
    }

    #[test]
    fn box_minus_cylinder_round_trips_with_analytic_entities() {
        let block = make_box_brep(Vec3::new(-3.0, -3.0, 0.0), 6.0, 6.0, 4.0).unwrap();
        let drill = make_cylinder_brep(
            Vec3::new(0.0, 0.0, -1.0),
            Vec3::new(0.0, 0.0, 1.0),
            1.5,
            6.0,
        )
        .unwrap();
        let cut = boolean_operation(
            &block,
            &drill,
            BooleanOperation::Subtract,
            &BooleanOptions::default(),
        )
        .unwrap();
        // Boolean-produced edges may be split subranges (NURBS fallback), so
        // only the surface entities are required to be analytic here.
        let solid = assert_analytic_round_trip(
            "box_minus_cyl",
            &cut,
            &["CYLINDRICAL_SURFACE(", "PLANE("],
            false,
        );
        assert_eq!(solid.genus, 1, "through-hole genus survives the round trip");
    }
}