ballistics_engine/

hold_curve.rs

//! `HoldCurve`: a solved, sampled drop-vs-range curve read by interpolation, and the
//! sampled-trajectory helpers it is built on (MBA-1361/MBA-1362).
//!
//! Promoted out of the CLI binary (0.33.0 decision-support Task 8) so a library consumer can
//! reuse it without the `cli` feature -- starting with the constant-drop range-card work
//! [`HoldCurve`]'s own doc comment already nominates as its next consumer. No feature gate:
//! must compile for wasm32 (pure math over already-resolved inputs; no fs, no clap -- CLI
//! argument resolution, `InverseSolverLoadArgs::resolve`, stays in `main.rs` and constructs
//! [`HoldCurveLoad`] there).

use crate::atmosphere::StationAtmospherePresence;
use crate::constants::GRAINS_TO_KG;
use crate::trajectory_observation::{bracket_param, Bracket};
use crate::truing::fallback_bullet_length_m;
use crate::truing_dsf::DsfTable;
use crate::{
    trajectory_sampling, AtmosphericConditions, BCSegmentData, BallisticInputs, DragModel,
    TrajectorySolver, WindConditions,
};
use std::error::Error;

/// Shared: build BallisticInputs + atmosphere + wind from common parameters (all in metric)
#[allow(
    clippy::too_many_arguments,
    reason = "flat arguments preserve the shared CLI trajectory compatibility helper"
)]
pub fn build_trajectory_components(
    velocity: f64,
    bc: f64,
    mass: f64,
    diameter: f64,
    drag_model: DragModel,
    sight_height: f64,
    temperature: f64,
    pressure: f64,
    humidity: f64,
    altitude: f64,
    wind_speed: f64,
    wind_direction: f64,
    max_range: f64,
    sample_interval: f64,
    // MBA-1323 Phase 2: saved-profile velocity-BC segments / Mach-Cd drag curve, already
    // resolved to engine shapes by the caller (bc_segments_from_profile /
    // drag_table_from_profile). `None` for every caller that does not (yet) source these from
    // a saved profile — see handle_come_ups/handle_lead for the callers that do.
    bc_segments_data: Option<Vec<BCSegmentData>>,
    custom_drag_table: Option<crate::drag::DragTable>,
    // MBA-1359: deliberate POI offset at the zero range (meters), carried on the inputs so
    // the zero solves that reuse these components inherit it. 0.0 = no offset.
    zero_poi_vertical_m: f64,
    zero_poi_horizontal_m: f64,
    // MBA-1396: lateral sight-to-bore mount offset (meters, positive = sight right of
    // bore). 0.0 = sight directly above the bore.
    sight_offset_lateral_m: f64,
) -> (BallisticInputs, WindConditions, AtmosphericConditions) {
    let drag_model_enum = drag_model;
    let wind_direction_rad = wind_direction.to_radians();
    let use_bc_segments = bc_segments_data.is_some();

    let inputs = BallisticInputs {
        bc_value: bc,
        bc_type: drag_model_enum,
        bullet_mass: mass,
        muzzle_velocity: velocity,
        bullet_diameter: diameter,
        bullet_length: fallback_bullet_length_m(diameter, mass),
        muzzle_angle: 0.0,
        target_distance: max_range,
        sight_height,
        sight_offset_lateral_m,
        zero_poi_vertical_m,
        zero_poi_horizontal_m,
        altitude,
        temperature,
        pressure,
        humidity,
        wind_speed,
        wind_angle: wind_direction_rad,
        use_rk4: true,           // Required for non-Euler solver
        use_adaptive_rk45: true, // Use RK45 adaptive (default solver)
        enable_trajectory_sampling: true,
        sample_interval,
        caliber_inches: diameter / 0.0254,
        weight_grains: mass / GRAINS_TO_KG,
        twist_rate: 12.0,
        is_twist_right: true,
        use_bc_segments,
        bc_segments_data,
        custom_drag_table,
        ..Default::default()
    };

    // wind_direction enters from the CLI in degrees; both engine structures use radians.
    let wind = WindConditions {
        speed: wind_speed,
        direction: wind_direction_rad,
        ..Default::default()
    };

    let atmosphere = AtmosphericConditions {
        temperature,
        pressure,
        humidity,
        altitude,
    };

    (inputs, wind, atmosphere)
}

/// A solved flight: the sampled grid, plus how far the bullet ACTUALLY got (MBA-1476).
///
/// `reach_m` is the solved trajectory's own terminal downrange distance — where it hit the
/// ground, ran out of flight time, or was stopped by the caller's `max_range` — not the last
/// point on the sampling grid. The two differ by up to one grid cell, and only `reach_m` is
/// independent of the grid the caller chose: it is the figure a card quotes when it has to
/// tell a shooter the bullet does not get out to the row they asked for. Quoting the last
/// SAMPLE there instead made that number a function of the request (`--end` sets the sampled
/// span), so the same load's "reaches only …" figure moved when nothing about the load had.
pub struct SampledFlight {
    pub samples: Vec<trajectory_sampling::TrajectorySample>,
    /// Terminal downrange distance of the solved flight, meters.
    pub reach_m: f64,
}

/// Run a trajectory and return sampled points at the given zero angle.
///
/// Discards the flight's own reach; [`run_sampled_flight`] keeps it, and every card surface
/// uses that one because a truncated card has to state how far the load actually flew.
///
/// Presence-blind: `temperature` and `pressure` go through the legacy default-sentinel rule
/// (see [`TrajectorySolver::new`]). A caller that knows which were supplied uses
/// [`run_sampled_flight_with_station_presence`] and reads its `.samples`.
#[allow(
    clippy::too_many_arguments,
    reason = "flat arguments preserve the shared sampled-trajectory compatibility helper"
)]
pub fn run_sampled_trajectory(
    velocity: f64,
    bc: f64,
    mass: f64,
    diameter: f64,
    drag_model: DragModel,
    sight_height: f64,
    temperature: f64,
    pressure: f64,
    humidity: f64,
    altitude: f64,
    wind_speed: f64,
    wind_direction: f64,
    max_range: f64,
    sample_interval: f64,
    zero_angle_rad: f64,
    bc_segments_data: Option<Vec<BCSegmentData>>,
    custom_drag_table: Option<crate::drag::DragTable>,
    dsf_table: Option<&DsfTable>,
    zero_poi_vertical_m: f64,
    zero_poi_horizontal_m: f64,
    sight_offset_lateral_m: f64,
    zero_solve_distance_m: Option<f64>,
) -> Result<Vec<trajectory_sampling::TrajectorySample>, Box<dyn Error>> {
    Ok(run_sampled_flight(
        velocity,
        bc,
        mass,
        diameter,
        drag_model,
        sight_height,
        temperature,
        pressure,
        humidity,
        altitude,
        wind_speed,
        wind_direction,
        max_range,
        sample_interval,
        zero_angle_rad,
        bc_segments_data,
        custom_drag_table,
        dsf_table,
        zero_poi_vertical_m,
        zero_poi_horizontal_m,
        sight_offset_lateral_m,
        zero_solve_distance_m,
    )?
    .samples)
}

/// Run a trajectory and return its sampled points AND its terminal reach at the given zero
/// angle. See [`SampledFlight`] for why the reach is carried out alongside the samples.
///
/// Presence-blind: `temperature` and `pressure` go through the legacy default-sentinel rule.
/// See [`run_sampled_flight_with_station_presence`].
#[allow(
    clippy::too_many_arguments,
    reason = "flat arguments preserve the shared sampled-trajectory compatibility helper"
)]
pub fn run_sampled_flight(
    velocity: f64,
    bc: f64,
    mass: f64,
    diameter: f64,
    drag_model: DragModel,
    sight_height: f64,
    temperature: f64,
    pressure: f64,
    humidity: f64,
    altitude: f64,
    wind_speed: f64,
    wind_direction: f64,
    max_range: f64,
    sample_interval: f64,
    zero_angle_rad: f64,
    bc_segments_data: Option<Vec<BCSegmentData>>,
    custom_drag_table: Option<crate::drag::DragTable>,
    dsf_table: Option<&DsfTable>,
    zero_poi_vertical_m: f64,
    zero_poi_horizontal_m: f64,
    sight_offset_lateral_m: f64,
    zero_solve_distance_m: Option<f64>,
) -> Result<SampledFlight, Box<dyn Error>> {
    run_sampled_flight_with_station_presence(
        velocity,
        bc,
        mass,
        diameter,
        drag_model,
        sight_height,
        temperature,
        pressure,
        humidity,
        altitude,
        wind_speed,
        wind_direction,
        max_range,
        sample_interval,
        zero_angle_rad,
        bc_segments_data,
        custom_drag_table,
        dsf_table,
        zero_poi_vertical_m,
        zero_poi_horizontal_m,
        sight_offset_lateral_m,
        zero_solve_distance_m,
        StationAtmospherePresence::NEITHER,
    )
}

/// [`run_sampled_flight`] knowing which of `temperature` and `pressure` the caller supplied
/// (0.49.0): a supplied value flies as given, an omitted one keeps the legacy default-sentinel
/// rule (see [`TrajectorySolver::new_with_station_presence`]).
/// [`StationAtmospherePresence::NEITHER`] is bit-identical to [`run_sampled_flight`].
#[allow(
    clippy::too_many_arguments,
    reason = "flat arguments preserve the shared sampled-trajectory compatibility helper"
)]
pub fn run_sampled_flight_with_station_presence(
    velocity: f64,
    bc: f64,
    mass: f64,
    diameter: f64,
    drag_model: DragModel,
    sight_height: f64,
    temperature: f64,
    pressure: f64,
    humidity: f64,
    altitude: f64,
    wind_speed: f64,
    wind_direction: f64,
    max_range: f64,
    sample_interval: f64,
    zero_angle_rad: f64,
    // MBA-1323 Phase 2: see build_trajectory_components's doc comment on these two.
    bc_segments_data: Option<Vec<BCSegmentData>>,
    custom_drag_table: Option<crate::drag::DragTable>,
    // MBA-1357: saved profile's DSF table, already validated by the caller. `None` for
    // every call site except come-ups (the only sampled-trajectory command with a DSF
    // auto-apply story so far).
    dsf_table: Option<&DsfTable>,
    // MBA-1359: deliberate POI offset at the zero range (meters) plus the distance the
    // caller's zero solve used (Some iff a zero was solved). The vertical offset already
    // rides inside `zero_angle_rad` (the zero solve returns it biased); the horizontal
    // offset becomes an azimuth bias here, exactly as calculate_and_set_zero_angle applies
    // it to its own solver. (0.0, 0.0, _) and (_, _, None) are exact no-ops.
    zero_poi_vertical_m: f64,
    zero_poi_horizontal_m: f64,
    // MBA-1396: lateral sight-mount offset (meters). Physically displaces the initial
    // lateral position, and (via windage_zero_bias_rad) joins the azimuth convergence
    // below when a zero was solved.
    sight_offset_lateral_m: f64,
    zero_solve_distance_m: Option<f64>,
    station_presence: StationAtmospherePresence,
) -> Result<SampledFlight, Box<dyn Error>> {
    let (mut inputs, wind, atmosphere) = build_trajectory_components(
        velocity,
        bc,
        mass,
        diameter,
        drag_model,
        sight_height,
        temperature,
        pressure,
        humidity,
        altitude,
        wind_speed,
        wind_direction,
        max_range,
        sample_interval,
        bc_segments_data,
        custom_drag_table,
        zero_poi_vertical_m,
        zero_poi_horizontal_m,
        sight_offset_lateral_m,
    );
    inputs.muzzle_angle = zero_angle_rad;
    if let Some(zero_distance_m) = zero_solve_distance_m {
        inputs.azimuth_angle += inputs.windage_zero_bias_rad(zero_distance_m);
    }

    let mut solver =
        TrajectorySolver::new_with_station_presence(inputs, wind, atmosphere, station_presence);
    solver.set_max_range(max_range);
    solver.set_time_step(0.001);
    let result = solver.solve()?;

    let mut samples = result.sampled_points.unwrap_or_default();
    // MBA-1357: drop-only correction, same convention as apply_dsf — TrajectorySample's
    // `drop_m` uses the identical `LOS - actual` sign (see sample_trajectory's own doc
    // comment in trajectory_sampling.rs), so scaling it directly by the DSF factor is
    // the exact per-point transform apply_dsf performs on a full TrajectoryResult.
    // Per-point Mach uses the SAME frozen station speed of sound divisor apply_dsf uses
    // (velocity_mps / station_speed_of_sound_mps), not a per-altitude recompute.
    if let Some(table) = dsf_table {
        let station_sos = result.station_speed_of_sound_mps;
        for s in samples.iter_mut() {
            let mach = if station_sos > 0.0 {
                s.velocity_mps / station_sos
            } else {
                0.0
            };
            s.drop_m *= table.factor_at(mach);
        }
    }

    Ok(SampledFlight {
        samples,
        // `TrajectoryResult::max_range` is the terminal point's downrange distance, not the
        // ceiling that was requested — the flight's own reach.
        reach_m: result.max_range,
    })
}

/// Distance interval every DOPE-card surface samples its trajectory on, meters (~1 yard).
///
/// **Deliberately independent of the card's row spacing** (MBA-1476). Sampling on a grid
/// whose spacing was the card's own `--step` made the value printed against a range depend
/// on the step, the start and the end of the card that asked for it: the grid is anchored at
/// zero, so a card starting anywhere that is not a multiple of its step put every requested
/// row between two samples, and [`sample_at_range`]'s predecessor -- a nearest-sample search
/// with a tolerance of one and a half whole steps -- then substituted a neighbouring range's
/// data without saying so. A shooter dialing `300 yd` off a `--start 300 --step 200` card was
/// reading the 200-yard line. The grid spacing is now a constant of the engine, and the value
/// at a range is read off it by interpolation, so a range's row is the same row on every card
/// that contains it.
///
/// Same grid as [`HoldCurve::SAMPLE_INTERVAL_M`], and for the same reason: fine enough that
/// linear interpolation between neighbours is orders of magnitude below what any card
/// resolves, coarse enough that a 1500 m card is a couple of thousand points. Sampling is
/// pure post-integration interpolation over the solver's knots -- it does not re-integrate --
/// so refining it is measurably free next to the solve that produced those knots.
pub const CARD_SAMPLE_INTERVAL_M: f64 = HoldCurve::SAMPLE_INTERVAL_M;

/// Read a sampled trajectory at an EXACT range by linear interpolation between the two
/// samples that bracket it (MBA-1476).
///
/// `None` when `range_m` lies outside the sampled span (or the span is degenerate) — that is
/// a row the solved flight genuinely cannot supply, and every caller reports it as an error.
/// It must never be answered with the closest sample lying around: the substitution is
/// invisible in the output, and a card row is a number someone dials.
///
/// The returned sample carries `distance_m == range_m` — the range asked for, which is also
/// the range the row is labelled with — so the angular conversion downstream divides by the
/// same distance the shooter reads. `flags` are dropped: they mark grid points (zero
/// crossing, apex, Mach transition), not arbitrary interpolated ranges.
pub fn sample_at_range(
    samples: &[trajectory_sampling::TrajectorySample],
    range_m: f64,
) -> Option<trajectory_sampling::TrajectorySample> {
    if !range_m.is_finite() {
        return None;
    }
    let Bracket::Inside { lo, t } =
        bracket_param(samples.len(), |i| samples[i].distance_m, range_m)
    else {
        return None;
    };
    let hi = lo + 1;
    let lerp = |a: f64, b: f64| a + (b - a) * t;
    Some(trajectory_sampling::TrajectorySample {
        distance_m: range_m,
        drop_m: lerp(samples[lo].drop_m, samples[hi].drop_m),
        wind_drift_m: lerp(samples[lo].wind_drift_m, samples[hi].wind_drift_m),
        velocity_mps: lerp(samples[lo].velocity_mps, samples[hi].velocity_mps),
        energy_j: lerp(samples[lo].energy_j, samples[hi].energy_j),
        time_s: lerp(samples[lo].time_s, samples[hi].time_s),
        flags: Vec::new(),
    })
}

/// How far out a card's trajectory must be solved and sampled, given the furthest row it
/// asked for (`end_m`). Meters in, meters out.
///
/// **The headroom past the last row is load-bearing, not a leftover.** The sampling grid is
/// anchored at zero with a fixed [`CARD_SAMPLE_INTERVAL_M`] spacing, so the last grid point
/// inside a span is the largest multiple of that interval the span contains: a span ending
/// only 10 % past the last requested row leaves less than one whole cell of headroom as soon
/// as `end_m` drops below ten intervals (~10 yd), and the last grid point then lands SHORT of
/// the last row. That row is genuinely unsampled, and the card correctly declines to invent
/// it — which is how a request as ordinary as a 5.1 yd card ended up with no rows at all.
///
/// Trimming the 10 % — reading it as vestigial because the row-spacing grid it once fed is
/// gone — widens that hole rather than closing it. The floor is what closes it: at least one
/// whole grid cell beyond the last row, so a flight that gets there is always sampled there.
pub fn card_sample_max_range_m(end_m: f64) -> f64 {
    (end_m * 1.1).max(end_m + CARD_SAMPLE_INTERVAL_M)
}

/// The one wording every card surface uses when its solved flight supplies NONE of the rows
/// asked for (MBA-1476) — stated in the caller's display unit, because that is the unit the
/// request named the range in.
///
/// A flight that supplies some but not all of them is not an error: the card prints the rows
/// it has and says what it left out, via [`CardTruncation`].
pub fn range_not_sampled_message(requested_display: f64, reach_display: f64, unit: &str) -> String {
    format!(
        "no trajectory sample at {requested_display:.0} {unit}: this load's solved flight \
         reaches only {reach_display:.0} {unit}"
    )
}

/// A card that stopped short of the rows it was asked for, in the caller's DISPLAY unit
/// (MBA-1476 follow-up).
///
/// Refusing to FABRICATE a row the flight never reached is right, and stays. Refusing to
/// print the rows it did reach is not: a `.22 LR` asked for the default 1200 yd card should
/// see its hundred-yard lines and learn where the bullet stops, not a bare exit code. So the
/// card truncates, and this is the record of that — carried to stderr for a CLI reader and as
/// a structured field on every JSON/bridge surface, because a truncation a caller cannot
/// detect is its own silent failure.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct CardTruncation {
    /// The furthest row the card asked for (its `--end`).
    pub requested_end: f64,
    /// The last row it could actually print.
    pub last_row: f64,
    /// How far the solved flight itself got — a property of the load, not of the request.
    pub reach: f64,
}

impl CardTruncation {
    /// The one wording every truncated card surface uses. `unit` labels the display distance
    /// unit ("yd" / "m").
    pub fn message(&self, unit: &str) -> String {
        let Self {
            requested_end,
            last_row,
            reach,
        } = *self;
        format!(
            "card truncated at {last_row:.0} {unit}: {requested_end:.0} {unit} was requested, \
             but this load's solved flight reaches only {reach:.0} {unit}; the rows past that \
             are left out rather than filled in from a range the bullet does reach"
        )
    }

    /// The same notice sized for PAPER: the one line a printed card carries in its footer
    /// (MBA-1477).
    ///
    /// [`Self::message`] is written for a terminal, where a 180-character sentence wraps and
    /// costs nothing. A dope card footer is a single centred line about eighty characters
    /// wide at the fixed footer size, drawn on every page beside the engine/table provenance,
    /// so the printed form states the same three facts — where the card stops, what was
    /// asked for, how far the load flies — and drops the sentence explaining that the missing
    /// rows were not fabricated. Both live on this type so the paper and the terminal cannot
    /// come to say different things about the same card.
    pub fn printed_note(&self, unit: &str) -> String {
        let Self {
            requested_end,
            last_row,
            reach,
        } = *self;
        format!(
            "TRUNCATED at {last_row:.0} {unit}: {requested_end:.0} {unit} requested, this load \
             reaches only {reach:.0} {unit}"
        )
    }
}

/// Everything one sampled hold curve needs, already in METRIC (MBA-1361/MBA-1362).
///
/// Flat and small on purpose: it is built once at a CLI boundary and handed to
/// [`HoldCurve::solve`], which owns the zero solve so every consumer zeroes identically.
#[derive(Debug, Clone)]
pub struct HoldCurveLoad {
    pub velocity_mps: f64,
    pub bc: f64,
    pub mass_kg: f64,
    pub diameter_m: f64,
    pub drag_model: DragModel,
    pub sight_height_m: f64,
    pub zero_distance_m: f64,
    pub temperature_c: f64,
    pub pressure_hpa: f64,
    pub humidity: f64,
    pub altitude_m: f64,
    pub wind_speed_mps: f64,
    pub wind_direction_deg: f64,
    /// Which of `temperature_c` / `pressure_hpa` were supplied (0.49.0). A supplied value flies
    /// as given; [`StationAtmospherePresence::NEITHER`] keeps the legacy default-sentinel rule,
    /// under which a sea-level default at a real altitude flies ICAO standard air.
    pub station_presence: StationAtmospherePresence,
}

/// One point read off a [`HoldCurve`], in the ANGULAR units a hold is expressed in.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct HoldPoint {
    pub range_m: f64,
    /// Milliradians below the line of sight (positive = the shooter holds over).
    pub drop_mil: f64,
    /// Milliradians of lateral deflection, positive = the bullet goes RIGHT.
    pub wind_mil: f64,
    pub velocity_mps: f64,
    pub energy_j: f64,
    pub time_s: f64,
}

/// A solved, sampled drop-vs-range curve expressed in ANGULAR units (MBA-1361/MBA-1362).
///
/// One solve, sampled finely, then read by interpolation — the same shape `come-ups` uses,
/// but keyed on angle rather than linear drop because that is what a reticle mark is.
///
/// **This is THE drop-vs-range helper**, forward ([`Self::at_range`]) and inverse
/// ([`Self::range_for_angular_drop_mil`]), and all four consumers go through it:
/// `reticle hold --range` (MBA-1361) plus `mark-to-range`, `bdc-match` and `optimal-zero`
/// (MBA-1362). Sharing it is the point — three separately-written root finds would be
/// three chances for "the drop at 500 yards" to mean three different things. The future
/// constant-drop range-card ticket is mark-to-range with uniformly spaced marks and should
/// reuse it too rather than growing a fourth.
pub struct HoldCurve {
    samples: Vec<trajectory_sampling::TrajectorySample>,
    /// The solved flight's own terminal distance, meters — see [`Self::reach_m`].
    reach_m: f64,
}

/// Result of inverting the hold curve for one mark subtension (MBA-1362).
///
/// The two failure arms are REPORTED, never silently dropped: a mark the load cannot reach
/// is information a shooter needs (it is the reticle's usable range limit), and dropping it
/// would quietly shorten the answer table.
#[derive(Debug, Clone, PartialEq)]
pub enum MarkToRangeOutcome {
    /// The mark's subtension is matched at this range.
    Reached(Box<HoldPoint>),
    /// The subtension corresponds to no range past the far zero. Angular drop is exactly
    /// zero AT the far zero and only grows past it, so a mark at or ABOVE the optical
    /// center (a non-positive subtension) is never matched — those are under-holds for
    /// ranges inside the zero, where drop-vs-range is not monotone and "the range for this
    /// mark" is not well defined.
    InsideZero { far_zero_range_m: f64 },
    /// Angular drop never grows to the subtension within the searched trajectory.
    BeyondSearch { max_range_m: f64, max_drop_mil: f64 },
}

impl HoldCurve {
    /// Sample interval used by every hold curve, meters (~1 yard).
    ///
    /// Fine enough that linear interpolation between neighbours is well below the
    /// resolution any reticle can be read to, and coarse enough that a 1500 m curve is a
    /// few thousand points.
    pub const SAMPLE_INTERVAL_M: f64 = 0.9144;

    /// Solve once and sample out to `max_range_m`.
    pub fn solve(load: &HoldCurveLoad, max_range_m: f64) -> Result<Self, Box<dyn Error>> {
        if !max_range_m.is_finite() || max_range_m <= 0.0 {
            return Err("hold curve max range must be finite and greater than zero".into());
        }
        let zero_inputs = BallisticInputs {
            bc_value: load.bc,
            bc_type: load.drag_model,
            bullet_mass: load.mass_kg,
            muzzle_velocity: load.velocity_mps,
            bullet_diameter: load.diameter_m,
            bullet_length: fallback_bullet_length_m(load.diameter_m, load.mass_kg),
            sight_height: load.sight_height_m,
            use_rk4: true,
            ..Default::default()
        };
        let atmosphere = AtmosphericConditions {
            temperature: load.temperature_c,
            pressure: load.pressure_hpa,
            humidity: load.humidity,
            altitude: load.altitude_m,
        };
        let zero_angle = crate::calculate_zero_angle_with_station_presence(
            zero_inputs,
            load.zero_distance_m,
            load.sight_height_m,
            WindConditions::default(),
            atmosphere,
            load.station_presence,
        )?;

        // MBA-1478: the FLIGHT, not just its samples. A card built on this curve has to be
        // able to quote how far the load actually flew (see `SampledFlight::reach_m`), and
        // the last sample is not that figure — it moves with the span the caller asked for.
        let flight = run_sampled_flight_with_station_presence(
            load.velocity_mps,
            load.bc,
            load.mass_kg,
            load.diameter_m,
            load.drag_model,
            load.sight_height_m,
            load.temperature_c,
            load.pressure_hpa,
            load.humidity,
            load.altitude_m,
            load.wind_speed_mps,
            load.wind_direction_deg,
            max_range_m,
            Self::SAMPLE_INTERVAL_M,
            zero_angle,
            None,
            None,
            None,
            0.0,
            0.0,
            0.0,
            Some(load.zero_distance_m),
            load.station_presence,
        )?;
        let SampledFlight { samples, reach_m } = flight;
        if samples.len() < 2 {
            return Err(
                "the trajectory produced too few sampled points to read a hold from".into(),
            );
        }
        Ok(Self { samples, reach_m })
    }

    /// The furthest range this curve reaches, meters.
    pub fn max_sampled_range_m(&self) -> f64 {
        self.samples.last().map_or(0.0, |s| s.distance_m)
    }

    /// The solved flight's own terminal downrange distance, meters (MBA-1478).
    ///
    /// Distinct from [`Self::max_sampled_range_m`], which is the last point on the sampling
    /// grid and therefore a function of the span the CALLER asked to sample. This one is a
    /// property of the load: it is the figure a truncated card quotes when it tells a shooter
    /// the bullet does not get out to the row they asked for. See [`SampledFlight`], whose
    /// doc comment records why quoting the last sample there was wrong.
    pub fn reach_m(&self) -> f64 {
        self.reach_m
    }

    /// This curve's own sample ranges, in order, meters.
    ///
    /// Exact multiples of [`Self::SAMPLE_INTERVAL_M`] (`i as f64 * SAMPLE_INTERVAL_M` for
    /// `i = 0..N`), the same arithmetic sequence a caller would otherwise have to reproduce by
    /// hand to reason about where this curve was actually verified -- an additive accessor so
    /// no consumer needs read access to the private `samples` field just to answer "which
    /// ranges did this curve solve at."
    pub fn sample_ranges_m(&self) -> Vec<f64> {
        self.samples.iter().map(|s| s.distance_m).collect()
    }

    /// Linearly interpolate the angular hold at `range_m`.
    ///
    /// `None` when the range is outside the sampled span or non-positive (an angular drop
    /// is undefined at the muzzle — it divides by the range).
    pub fn at_range(&self, range_m: f64) -> Option<HoldPoint> {
        if !range_m.is_finite() || range_m <= 0.0 {
            return None;
        }
        let samples = &self.samples;
        let Bracket::Inside { lo, t } =
            bracket_param(samples.len(), |i| samples[i].distance_m, range_m)
        else {
            return None;
        };
        let hi = lo + 1;
        let lerp = |a: f64, b: f64| a + (b - a) * t;
        let drop_m = lerp(samples[lo].drop_m, samples[hi].drop_m);
        let drift_m = lerp(samples[lo].wind_drift_m, samples[hi].wind_drift_m);
        Some(HoldPoint {
            range_m,
            // Milliradian small-angle definition: 1 mil subtends 1/1000 of the range.
            drop_mil: drop_m / range_m * 1000.0,
            wind_mil: drift_m / range_m * 1000.0,
            velocity_mps: lerp(samples[lo].velocity_mps, samples[hi].velocity_mps),
            energy_j: lerp(samples[lo].energy_j, samples[hi].energy_j),
            time_s: lerp(samples[lo].time_s, samples[hi].time_s),
        })
    }

    /// The downrange distance of the FAR zero crossing, meters — the point past which
    /// angular drop grows monotonically with range.
    ///
    /// Angular drop is not monotone over the whole flight: it starts large and positive at
    /// the muzzle (the bullet is a sight height below the line of sight, divided by a tiny
    /// range), falls through zero at the near zero, goes negative while the bullet rides
    /// above the line of sight, and returns through zero at the far zero. Only past that
    /// second crossing is the inverse below single-valued, so the search domain starts
    /// there rather than at the muzzle.
    pub fn far_zero_range_m(&self) -> f64 {
        let mut far = self.samples.first().map_or(0.0, |s| s.distance_m);
        for sample in &self.samples {
            if sample.distance_m > 0.0 && sample.drop_m <= 0.0 {
                far = sample.distance_m;
            }
        }
        far
    }

    /// Bisection cap for the inverse below. The curve is monotone on the searched interval,
    /// so bisection halves the bracket every step; 80 iterations takes any realistic
    /// bracket far below the tolerance and exists only as a runaway guard.
    const INVERSE_MAX_ITERATIONS: u32 = 80;

    /// Bracket width (meters) at which the inverse stops regardless of residual — a
    /// hundredth of a millimeter, four orders below anything a range card resolves.
    const INVERSE_TOLERANCE_M: f64 = 1.0e-5;

    /// Invert the curve: the range at which the angular drop equals `target_mil`.
    ///
    /// Bisection over the interpolated curve on `[far zero, furthest sample]`, where drop
    /// is monotone increasing in range. Both out-of-domain cases come back as their own
    /// outcome rather than as a clamped range.
    pub fn range_for_angular_drop_mil(&self, target_mil: f64) -> MarkToRangeOutcome {
        let far_zero_range_m = self.far_zero_range_m();
        let max_range_m = self.max_sampled_range_m();
        let drop_at = |range_m: f64| self.at_range(range_m).map(|p| p.drop_mil);
        let max_drop_mil = drop_at(max_range_m).unwrap_or(f64::NEG_INFINITY);

        if !target_mil.is_finite() || target_mil <= 0.0 {
            return MarkToRangeOutcome::InsideZero { far_zero_range_m };
        }
        if target_mil > max_drop_mil {
            return MarkToRangeOutcome::BeyondSearch {
                max_range_m,
                max_drop_mil,
            };
        }
        // Start the bracket just past the far zero: exactly at it the drop is 0, and any
        // positive target is therefore bracketed.
        let mut lo = far_zero_range_m.max(Self::SAMPLE_INTERVAL_M);
        if drop_at(lo).is_none_or(|drop| drop >= target_mil) {
            return MarkToRangeOutcome::InsideZero { far_zero_range_m };
        }
        let mut hi = max_range_m;
        for _ in 0..Self::INVERSE_MAX_ITERATIONS {
            if hi - lo <= Self::INVERSE_TOLERANCE_M {
                break;
            }
            let mid = 0.5 * (lo + hi);
            match drop_at(mid) {
                Some(drop) if drop < target_mil => lo = mid,
                Some(_) => hi = mid,
                None => break,
            }
        }
        match self.at_range(0.5 * (lo + hi)) {
            Some(point) => MarkToRangeOutcome::Reached(Box::new(point)),
            None => MarkToRangeOutcome::BeyondSearch {
                max_range_m,
                max_drop_mil,
            },
        }
    }
}

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

    /// A representative metric load (a .308-class 168 gr match bullet at a 100 m zero),
    /// solved with a nonzero crosswind so both `drop_mil` and `wind_mil` are exercised.
    fn oracle_test_load() -> HoldCurveLoad {
        HoldCurveLoad {
            velocity_mps: 800.0,
            bc: 0.223,
            mass_kg: 0.0109,
            diameter_m: 0.00782,
            drag_model: DragModel::G7,
            sight_height_m: 0.045,
            zero_distance_m: 100.0,
            temperature_c: 15.0,
            pressure_hpa: 1013.25,
            humidity: 50.0,
            altitude_m: 0.0,
            wind_speed_mps: 3.0,
            wind_direction_deg: 90.0,
            station_presence: StationAtmospherePresence::NEITHER,
        }
    }

    /// Reproduces exactly what [`HoldCurve::solve`] does internally, as a SEPARATE, from-
    /// scratch call -- it never reads any `HoldCurve`'s private `samples` field. This gives
    /// `at_range` an independent direct-solve oracle to be checked against, rather than a
    /// self-consistency check of `HoldCurve` against its own stored state.
    fn independent_direct_solve(
        load: &HoldCurveLoad,
        max_range_m: f64,
    ) -> Vec<trajectory_sampling::TrajectorySample> {
        let zero_inputs = BallisticInputs {
            bc_value: load.bc,
            bc_type: load.drag_model,
            bullet_mass: load.mass_kg,
            muzzle_velocity: load.velocity_mps,
            bullet_diameter: load.diameter_m,
            bullet_length: fallback_bullet_length_m(load.diameter_m, load.mass_kg),
            sight_height: load.sight_height_m,
            use_rk4: true,
            ..Default::default()
        };
        let atmosphere = AtmosphericConditions {
            temperature: load.temperature_c,
            pressure: load.pressure_hpa,
            humidity: load.humidity,
            altitude: load.altitude_m,
        };
        let zero_angle = crate::calculate_zero_angle_with_conditions(
            zero_inputs,
            load.zero_distance_m,
            load.sight_height_m,
            WindConditions::default(),
            atmosphere,
        )
        .expect("zero solve should succeed for a realistic load");

        run_sampled_trajectory(
            load.velocity_mps,
            load.bc,
            load.mass_kg,
            load.diameter_m,
            load.drag_model,
            load.sight_height_m,
            load.temperature_c,
            load.pressure_hpa,
            load.humidity,
            load.altitude_m,
            load.wind_speed_mps,
            load.wind_direction_deg,
            max_range_m,
            HoldCurve::SAMPLE_INTERVAL_M,
            zero_angle,
            None,
            None,
            None,
            0.0,
            0.0,
            0.0,
            Some(load.zero_distance_m),
        )
        .expect("sampled trajectory should succeed for a realistic load")
    }

    #[test]
    fn at_range_matches_an_independent_direct_solves_sampled_observation() {
        let load = oracle_test_load();
        let max_range_m = 1500.0;

        let hold_curve = HoldCurve::solve(&load, max_range_m).expect("hold curve should solve");
        let oracle_samples = independent_direct_solve(&load, max_range_m);

        // A probe range read off the oracle's own grid rather than a hand-picked constant --
        // it only needs to land exactly on a sample node so `at_range`'s bracket search
        // resolves with t == 1.0 (the exact-endpoint-hit convention Task 7 disclosed: an exact
        // key hit lands on the upper endpoint of its bracket, not the lower one) and its own
        // linear interpolation is not itself a source of disagreement with the oracle.
        let probe = &oracle_samples[oracle_samples.len() / 2];
        let probe_range_m = probe.distance_m;
        assert!(probe_range_m > 0.0 && probe_range_m.is_finite());

        let point = hold_curve
            .at_range(probe_range_m)
            .expect("probe range must be inside the sampled span");

        let expected_drop_mil = probe.drop_m / probe_range_m * 1000.0;
        let expected_wind_mil = probe.wind_drift_m / probe_range_m * 1000.0;

        assert!((point.range_m - probe_range_m).abs() < 1e-9);
        assert!((point.drop_mil - expected_drop_mil).abs() < 1e-9);
        assert!((point.wind_mil - expected_wind_mil).abs() < 1e-9);
        assert!((point.velocity_mps - probe.velocity_mps).abs() < 1e-9);
        assert!((point.energy_j - probe.energy_j).abs() < 1e-9);
        assert!((point.time_s - probe.time_s).abs() < 1e-9);
    }

    #[test]
    fn at_range_outside_the_sampled_span_returns_none() {
        let load = oracle_test_load();
        let max_range_m = 1500.0;
        let hold_curve = HoldCurve::solve(&load, max_range_m).expect("hold curve should solve");

        let beyond = hold_curve.max_sampled_range_m() + 10_000.0;
        assert_eq!(hold_curve.at_range(beyond), None);
    }

    /// `sample_ranges_m` must reproduce the exact arithmetic sequence `HoldCurve::solve` built
    /// internally: index 0 at the muzzle (`0.0`), every later index an exact multiple of
    /// `SAMPLE_INTERVAL_M`, and the last entry equal to `max_sampled_range_m()`.
    #[test]
    fn sample_ranges_m_is_the_exact_multiples_of_the_sample_interval() {
        let load = oracle_test_load();
        let hold_curve = HoldCurve::solve(&load, 1500.0).expect("hold curve should solve");

        let ranges = hold_curve.sample_ranges_m();
        assert!(!ranges.is_empty());
        assert_eq!(ranges[0], 0.0);
        assert_eq!(*ranges.last().unwrap(), hold_curve.max_sampled_range_m());
        for (i, &r) in ranges.iter().enumerate() {
            let expected = i as f64 * HoldCurve::SAMPLE_INTERVAL_M;
            assert!(
                (r - expected).abs() < 1e-9,
                "index {i}: got {r}, expected {expected} ({} * SAMPLE_INTERVAL_M)",
                i
            );
        }
    }
}