navcore_math/

polar.rs

//! Boat performance polars, and the tactical numbers built on one.
//!
//! A polar table states a boat's expected speed for a given wind angle
//! and strength, values that depend on hull and rig and cannot be
//! computed from a formula. [`crate::sailing::tack_legs`] and
//! [`crate::sailing::close_hauled_headings_deg`] take the beat angle as
//! a caller-stated number; [`Polar::best_upwind_angle_deg`] derives it
//! instead, from the polar itself.
//!
//! # What is not modelled
//!
//! Sea state, reefing points, and crew skill are not modelled. A polar
//! is a prediction, usually from a velocity prediction program or a
//! measured racing record, not a live reading. Nothing here compares it
//! against how the boat is actually sailing.

use crate::{angle, sailing};

/// Something wrong with a polar table.
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum PolarError {
    /// The text was not a polar table at all -- too few lines, or a
    /// header or row that did not parse as numbers.
    Malformed(String),
    /// Parsed, but the shape does not make sense as a table -- fewer than
    /// two wind speeds or angles, a row of the wrong length, or an axis
    /// that is not strictly ascending.
    InvalidShape(String),
}

impl std::fmt::Display for PolarError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            Self::Malformed(reason) | Self::InvalidShape(reason) => write!(f, "{reason}"),
        }
    }
}

impl std::error::Error for PolarError {}

/// A boat's expected speed for every true wind angle and speed on a grid.
///
/// Symmetric port/starboard by construction: every lookup folds its own
/// angle to `0..=180` first -- see [`Polar::boat_speed_kn`] -- because a
/// hull does not sail differently on the other tack.
#[derive(Debug, Clone, PartialEq)]
pub struct Polar {
    /// True wind speeds the table is defined at, knots, strictly
    /// ascending.
    tws_kn: Vec<f64>,
    /// True wind angles the table is defined at, degrees `0..=180`,
    /// strictly ascending.
    twa_deg: Vec<f64>,
    /// Boat speed in knots, one row per `twa_deg` entry, one column per
    /// `tws_kn` entry -- the same layout a polar file's own rows and
    /// columns are already in.
    boat_speed_kn: Vec<Vec<f64>>,
}

impl Polar {
    /// Builds a table from its own axes and the speed at every
    /// intersection.
    ///
    /// # Errors
    ///
    /// [`PolarError::InvalidShape`] when either axis has fewer than two
    /// points, either axis is not strictly ascending, `twa_deg` strays
    /// outside `0..=180`, or `boat_speed_kn`'s own shape does not match
    /// `twa_deg.len()` rows of `tws_kn.len()` columns each.
    pub fn new(tws_kn: Vec<f64>, twa_deg: Vec<f64>, boat_speed_kn: Vec<Vec<f64>>) -> Result<Self, PolarError> {
        if tws_kn.len() < 2 {
            return Err(PolarError::InvalidShape("fewer than two wind speeds".to_owned()));
        }
        if twa_deg.len() < 2 {
            return Err(PolarError::InvalidShape("fewer than two wind angles".to_owned()));
        }
        if !tws_kn.windows(2).all(|pair| pair[0] < pair[1]) {
            return Err(PolarError::InvalidShape("wind speeds are not strictly ascending".to_owned()));
        }
        if !twa_deg.windows(2).all(|pair| pair[0] < pair[1]) {
            return Err(PolarError::InvalidShape("wind angles are not strictly ascending".to_owned()));
        }
        if twa_deg[0] < 0.0 || *twa_deg.last().expect("checked len >= 2") > 180.0 {
            return Err(PolarError::InvalidShape("wind angles must lie within 0..=180".to_owned()));
        }
        if boat_speed_kn.len() != twa_deg.len() {
            return Err(PolarError::InvalidShape(format!(
                "expected {} rows of boat speed, one per wind angle, got {}",
                twa_deg.len(),
                boat_speed_kn.len()
            )));
        }
        if let Some(bad_row) = boat_speed_kn.iter().find(|row| row.len() != tws_kn.len()) {
            return Err(PolarError::InvalidShape(format!(
                "every row must have {} boat speeds, one per wind speed, found one with {}",
                tws_kn.len(),
                bad_row.len()
            )));
        }

        Ok(Self { tws_kn, twa_deg, boat_speed_kn })
    }

    /// Parses the common tab/space-separated polar table format: a header
    /// row of wind speeds (its own first cell, conventionally `twa/tws`,
    /// ignored), then one row per wind angle -- the angle itself, then
    /// one boat speed per wind speed column, in the same order as the
    /// header.
    ///
    /// # Errors
    ///
    /// [`PolarError::Malformed`] when the text does not have this shape
    /// at all; [`PolarError::InvalidShape`] when [`Polar::new`]'s own
    /// checks fail on what was read.
    pub fn parse(text: &str) -> Result<Self, PolarError> {
        let mut lines = text.lines().map(str::trim).filter(|line| !line.is_empty());

        let header = lines.next().ok_or_else(|| PolarError::Malformed("empty table".to_owned()))?;
        let mut header_fields = header.split_whitespace();
        // The header's own first cell names the table ("twa/tws"
        // conventionally), not a number -- skipped rather than parsed.
        header_fields.next();
        let tws_kn: Vec<f64> = header_fields
            .map(|field| {
                field.parse::<f64>().map_err(|_| PolarError::Malformed(format!("bad wind speed: {field}")))
            })
            .collect::<Result<_, _>>()?;

        let mut twa_deg = Vec::new();
        let mut boat_speed_kn = Vec::new();
        for line in lines {
            let mut fields = line.split_whitespace();
            let angle_deg = fields
                .next()
                .ok_or_else(|| PolarError::Malformed("row with no wind angle".to_owned()))?
                .parse::<f64>()
                .map_err(|_| PolarError::Malformed(format!("bad wind angle in row: {line}")))?;
            let speeds: Vec<f64> = fields
                .map(|field| {
                    field.parse::<f64>().map_err(|_| PolarError::Malformed(format!("bad boat speed: {field}")))
                })
                .collect::<Result<_, _>>()?;
            twa_deg.push(angle_deg);
            boat_speed_kn.push(speeds);
        }

        Self::new(tws_kn, twa_deg, boat_speed_kn)
    }

    /// The boat's expected speed at `twa_deg`/`tws_kn`, bilinearly
    /// interpolated between the table's own grid points.
    ///
    /// `twa_deg` is folded to `0..=180` first -- see this type's own doc
    /// on why the table itself only ever states one side. A wind speed or
    /// angle outside the table's own range is clamped to its nearest edge
    /// rather than extrapolated: a polar has no honest opinion about a
    /// wind speed nobody sailed it in.
    #[must_use]
    pub fn boat_speed_kn(&self, twa_deg: f64, tws_kn: f64) -> f64 {
        let twa_deg = fold_wind_angle_deg(twa_deg);

        let (twa_lo, twa_hi, twa_t) = bracket(&self.twa_deg, twa_deg);
        let (tws_lo, tws_hi, tws_t) = bracket(&self.tws_kn, tws_kn);

        let at = |twa_index: usize, tws_index: usize| self.boat_speed_kn[twa_index][tws_index];

        // Bilinear: interpolate along wind speed at each of the two
        // bracketing angles, then interpolate that pair along angle.
        let low_angle = lerp(at(twa_lo, tws_lo), at(twa_lo, tws_hi), tws_t);
        let high_angle = lerp(at(twa_hi, tws_lo), at(twa_hi, tws_hi), tws_t);
        lerp(low_angle, high_angle, twa_t)
    }

    /// Velocity made good directly to windward at `twa_deg`/`tws_kn`,
    /// using this table's own predicted boat speed --
    /// [`sailing::vmg_to_wind_kn`] fed from [`Polar::boat_speed_kn`]
    /// rather than a measured one.
    #[must_use]
    pub fn vmg_to_wind_kn(&self, twa_deg: f64, tws_kn: f64) -> f64 {
        sailing::vmg_to_wind_kn(self.boat_speed_kn(twa_deg, tws_kn), twa_deg)
    }

    /// The true wind angle that makes the most progress to windward at
    /// `tws_kn`, searched across the upwind half of the table (`0..90`)
    /// -- the "actual" beat angle real laylines can offer once a polar
    /// exists, in place of a mariner-stated constant.
    ///
    /// Scanned in fine steps rather than solved analytically: a polar's
    /// own boat-speed curve has no closed form, and VMG-to-wind is
    /// unimodal across a beat for every real polar this was checked
    /// against, so a dense scan finds the same peak an analytic solver
    /// would, without needing one.
    #[must_use]
    pub fn best_upwind_angle_deg(&self, tws_kn: f64) -> f64 {
        best_angle_in_range(0.0, 90.0, |twa_deg| self.vmg_to_wind_kn(twa_deg, tws_kn))
    }

    /// The true wind angle that makes the most progress downwind at
    /// `tws_kn`, searched across `90..180` -- the downwind counterpart to
    /// [`Polar::best_upwind_angle_deg`], for a boat whose polar makes
    /// sailing angles downwind faster than a dead run.
    #[must_use]
    pub fn best_downwind_angle_deg(&self, tws_kn: f64) -> f64 {
        // Downwind VMG is negative by vmg_to_wind_kn's own convention
        // (see that function's own doc); "best" downwind means most
        // negative, i.e. minimised rather than maximised.
        best_angle_in_range(90.0, 180.0, |twa_deg| -self.vmg_to_wind_kn(twa_deg, tws_kn))
    }

    /// The wind speeds this table is defined at, knots, ascending.
    ///
    /// For a caller that wants to walk the table's own grid -- writing it
    /// back out, whether as [`Polar::to_text`] or in another shape
    /// entirely, such as Signal K's own `performance.polars`.
    #[must_use]
    pub fn tws_kn(&self) -> &[f64] {
        &self.tws_kn
    }

    /// The wind angles this table is defined at, degrees `0..=180`,
    /// ascending. See [`Polar::tws_kn`] for why this is exposed.
    #[must_use]
    pub fn twa_deg(&self) -> &[f64] {
        &self.twa_deg
    }

    /// Writes this table back out in the same tab-separated format
    /// [`Polar::parse`] reads. `Polar::parse(&polar.to_text())` round
    /// trips, to the precision written -- the format an unrelated
    /// mariner's own sailing software can most likely already open.
    #[must_use]
    pub fn to_text(&self) -> String {
        let mut lines = Vec::with_capacity(self.twa_deg.len() + 1);

        let header: Vec<String> =
            std::iter::once("twa/tws".to_owned()).chain(self.tws_kn.iter().map(|tws| format!("{tws:.2}"))).collect();
        lines.push(header.join("\t"));

        for (row_index, twa_deg) in self.twa_deg.iter().enumerate() {
            let mut fields = vec![format!("{twa_deg:.1}")];
            fields.extend(self.boat_speed_kn[row_index].iter().map(|speed| format!("{speed:.2}")));
            lines.push(fields.join("\t"));
        }

        lines.join("\n")
    }
}

/// Folds any wind angle to the `0..=180` a polar table is stored over --
/// a hull sails identically on either tack, so `-45` and `225` both mean
/// the same 45 degrees off the bow that `45` does.
fn fold_wind_angle_deg(twa_deg: f64) -> f64 {
    let wrapped = angle::norm_360(twa_deg);
    if wrapped > 180.0 { 360.0 - wrapped } else { wrapped }
}

/// The two indices in `axis` bracketing `value`, and where between them
/// it falls, `0.0..=1.0`. Clamped at either end rather than extrapolated
/// -- see [`Polar::boat_speed_kn`]'s own doc.
fn bracket(axis: &[f64], value: f64) -> (usize, usize, f64) {
    if value <= axis[0] {
        return (0, 0, 0.0);
    }
    let last = axis.len() - 1;
    if value >= axis[last] {
        return (last, last, 0.0);
    }
    let hi = axis.iter().position(|&x| x > value).unwrap_or(last);
    let lo = hi - 1;
    let t = (value - axis[lo]) / (axis[hi] - axis[lo]);
    (lo, hi, t)
}

fn lerp(a: f64, b: f64, t: f64) -> f64 {
    a + (b - a) * t
}

/// The angle in `[lo, hi]` (inclusive) that maximises `score`, scanned at
/// a fine, fixed step -- see [`Polar::best_upwind_angle_deg`]'s own doc
/// for why a scan rather than a solver.
fn best_angle_in_range(lo: f64, hi: f64, score: impl Fn(f64) -> f64) -> f64 {
    const STEP_DEG: f64 = 0.1;
    let steps = ((hi - lo) / STEP_DEG).round() as usize;
    (0..=steps)
        .map(|i| lo + f64::from(u32::try_from(i).unwrap_or(u32::MAX)) * STEP_DEG)
        .max_by(|&a, &b| score(a).partial_cmp(&score(b)).expect("polar-derived VMG is never NaN"))
        .unwrap_or(lo)
}

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

    /// A real-shaped cruiser-racer polar, the kind a velocity prediction
    /// program hands back -- zero speed head to wind, rising through the
    /// beat, peaking on a reach, tailing off downwind.
    const SAMPLE: &str = "\
twa/tws\t6\t8\t10\t12\t14\t16\t20
0\t0\t0\t0\t0\t0\t0\t0
52\t5.24\t6.05\t6.41\t6.66\t6.77\t6.81\t6.87
60\t5.55\t6.28\t6.6\t6.83\t6.95\t7.03\t7.06
75\t5.79\t6.48\t6.85\t7.12\t7.28\t7.34\t7.36
90\t5.79\t6.55\t6.98\t7.32\t7.52\t7.6\t7.63
110\t5.59\t6.51\t7.11\t7.55\t7.83\t7.98\t8.09
120\t5.28\t6.19\t6.85\t7.42\t7.83\t8.11\t8.34
135\t4.65\t5.55\t6.19\t6.79\t7.29\t7.67\t8.19
150\t4.03\t4.81\t5.4\t5.94\t6.42\t6.87\t7.62
170\t3.55\t4.24\t4.76\t5.24\t5.66\t6.06\t6.72
180\t3.45\t4.12\t4.62\t5.09\t5.5\t5.89\t6.53";

    fn close(a: f64, b: f64, tol: f64) -> bool {
        (a - b).abs() < tol
    }

    #[test]
    fn the_sample_table_parses_with_its_own_shape() {
        let polar = Polar::parse(SAMPLE).expect("a well-formed table");
        assert_eq!(polar.tws_kn.len(), 7);
        assert_eq!(polar.twa_deg.len(), 11);
        assert_eq!(polar.boat_speed_kn.len(), 11);
        assert!(polar.boat_speed_kn.iter().all(|row| row.len() == 7));
    }

    #[test]
    // A boat speed in the sample table happens to read like an
    // approximation of tau; it is a knots figure from the fixture, not a
    // mistyped constant.
    #[allow(clippy::approx_constant)]
    fn a_grid_point_is_read_back_exactly() {
        let polar = Polar::parse(SAMPLE).expect("a well-formed table");
        assert!(close(polar.boat_speed_kn(90.0, 12.0), 7.32, 1e-9));
        assert!(close(polar.boat_speed_kn(60.0, 8.0), 6.28, 1e-9));
    }

    #[test]
    fn boat_speed_interpolates_between_two_wind_speeds() {
        let polar = Polar::parse(SAMPLE).expect("a well-formed table");
        // Exactly halfway between the 10 and 12 kn columns at TWA 90:
        // 6.98 and 7.32.
        let halfway = polar.boat_speed_kn(90.0, 11.0);
        assert!(close(halfway, (6.98 + 7.32) / 2.0, 1e-9), "{halfway}");
    }

    #[test]
    fn boat_speed_interpolates_between_two_wind_angles() {
        let polar = Polar::parse(SAMPLE).expect("a well-formed table");
        // Exactly halfway between TWA 75 and TWA 90 at 12 kn: 7.12 and
        // 7.32.
        let halfway = polar.boat_speed_kn(82.5, 12.0);
        assert!(close(halfway, (7.12 + 7.32) / 2.0, 1e-9), "{halfway}");
    }

    #[test]
    fn a_wind_speed_below_the_table_is_clamped_not_extrapolated() {
        let polar = Polar::parse(SAMPLE).expect("a well-formed table");
        assert!(close(polar.boat_speed_kn(90.0, 0.0), polar.boat_speed_kn(90.0, 6.0), 1e-9));
    }

    #[test]
    fn a_wind_speed_above_the_table_is_clamped_not_extrapolated() {
        let polar = Polar::parse(SAMPLE).expect("a well-formed table");
        assert!(close(polar.boat_speed_kn(90.0, 100.0), polar.boat_speed_kn(90.0, 20.0), 1e-9));
    }

    #[test]
    fn an_angle_past_180_mirrors_onto_the_stored_side() {
        let polar = Polar::parse(SAMPLE).expect("a well-formed table");
        assert!(close(polar.boat_speed_kn(225.0, 12.0), polar.boat_speed_kn(135.0, 12.0), 1e-9));
    }

    #[test]
    fn a_negative_angle_mirrors_the_same_way() {
        let polar = Polar::parse(SAMPLE).expect("a well-formed table");
        assert!(close(polar.boat_speed_kn(-52.0, 12.0), polar.boat_speed_kn(52.0, 12.0), 1e-9));
    }

    #[test]
    fn vmg_to_wind_matches_the_sailing_module_fed_the_same_speed() {
        let polar = Polar::parse(SAMPLE).expect("a well-formed table");
        let speed = polar.boat_speed_kn(60.0, 12.0);
        assert!(close(polar.vmg_to_wind_kn(60.0, 12.0), sailing::vmg_to_wind_kn(speed, 60.0), 1e-9));
    }

    #[test]
    fn the_best_upwind_angle_actually_beats_its_neighbours() {
        let polar = Polar::parse(SAMPLE).expect("a well-formed table");
        let best = polar.best_upwind_angle_deg(12.0);
        assert!((0.0..90.0).contains(&best), "{best}");

        let best_vmg = polar.vmg_to_wind_kn(best, 12.0);
        assert!(best_vmg >= polar.vmg_to_wind_kn(best - 1.0, 12.0));
        assert!(best_vmg >= polar.vmg_to_wind_kn(best + 1.0, 12.0));
        // And it genuinely beats a coarse sweep across the whole beat,
        // not just its own immediate neighbours.
        for degrees in (0..90).step_by(5) {
            assert!(best_vmg >= polar.vmg_to_wind_kn(f64::from(degrees), 12.0) - 1e-9);
        }
    }

    #[test]
    fn the_best_downwind_angle_actually_beats_its_neighbours() {
        let polar = Polar::parse(SAMPLE).expect("a well-formed table");
        let best = polar.best_downwind_angle_deg(12.0);
        assert!((90.0..180.0).contains(&best), "{best}");

        let best_vmg = polar.vmg_to_wind_kn(best, 12.0);
        for degrees in (90..180).step_by(5) {
            assert!(best_vmg <= polar.vmg_to_wind_kn(f64::from(degrees), 12.0) + 1e-9);
        }
    }

    #[test]
    fn a_stronger_breeze_points_higher_on_this_sample() {
        // The everyday observation this sample's own shape happens to
        // show: more breeze lets the boat point closer to the wind for
        // the same or better VMG. Not a universal law of polars -- just a
        // property of this one, worth a regression test since laylines
        // will lean on it.
        let polar = Polar::parse(SAMPLE).expect("a well-formed table");
        let light = polar.best_upwind_angle_deg(6.0);
        let strong = polar.best_upwind_angle_deg(16.0);
        assert!(strong <= light, "{strong} vs {light}");
    }

    #[test]
    fn the_axes_are_exposed_for_a_caller_writing_the_table_out() {
        let polar = Polar::parse(SAMPLE).expect("a well-formed table");
        assert_eq!(polar.tws_kn(), &[6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 20.0]);
        assert_eq!(polar.twa_deg().first(), Some(&0.0));
        assert_eq!(polar.twa_deg().last(), Some(&180.0));
    }

    #[test]
    fn writing_a_table_out_and_reading_it_back_round_trips() {
        let original = Polar::parse(SAMPLE).expect("a well-formed table");
        let text = original.to_text();
        let read_back = Polar::parse(&text).expect("what this just wrote is itself a well-formed table");

        assert_eq!(read_back.tws_kn(), original.tws_kn());
        assert_eq!(read_back.twa_deg(), original.twa_deg());
        for &twa_deg in original.twa_deg() {
            for &tws_kn in original.tws_kn() {
                // to_text formats to two decimal places, so the round
                // trip is exact to that precision, not bit-exact.
                assert!(
                    close(read_back.boat_speed_kn(twa_deg, tws_kn), original.boat_speed_kn(twa_deg, tws_kn), 1e-2),
                    "twa {twa_deg} tws {tws_kn}"
                );
            }
        }
    }

    #[test]
    fn garbage_text_does_not_parse() {
        assert!(matches!(Polar::parse("not a polar table"), Err(PolarError::Malformed(_))));
        assert!(matches!(Polar::parse(""), Err(PolarError::Malformed(_))));
    }

    #[test]
    fn a_row_with_the_wrong_number_of_speeds_is_rejected() {
        let bad = "twa/tws\t6\t8\n60\t5.5";
        assert!(matches!(Polar::parse(bad), Err(PolarError::InvalidShape(_))));
    }

    #[test]
    fn too_few_wind_speeds_is_rejected() {
        assert!(matches!(
            Polar::new(vec![10.0], vec![0.0, 90.0], vec![vec![0.0], vec![5.0]]),
            Err(PolarError::InvalidShape(_))
        ));
    }

    #[test]
    fn too_few_wind_angles_is_rejected() {
        assert!(matches!(
            Polar::new(vec![10.0, 20.0], vec![90.0], vec![vec![5.0, 6.0]]),
            Err(PolarError::InvalidShape(_))
        ));
    }

    #[test]
    fn a_non_ascending_axis_is_rejected() {
        assert!(matches!(
            Polar::new(vec![10.0, 5.0], vec![0.0, 90.0], vec![vec![0.0, 0.0], vec![5.0, 6.0]]),
            Err(PolarError::InvalidShape(_))
        ));
    }

    #[test]
    fn an_angle_outside_0_to_180_is_rejected() {
        assert!(matches!(
            Polar::new(vec![10.0, 20.0], vec![-10.0, 90.0], vec![vec![0.0, 0.0], vec![5.0, 6.0]]),
            Err(PolarError::InvalidShape(_))
        ));
        assert!(matches!(
            Polar::new(vec![10.0, 20.0], vec![90.0, 190.0], vec![vec![5.0, 6.0], vec![5.0, 6.0]]),
            Err(PolarError::InvalidShape(_))
        ));
    }
}