navcore_math/

rhumb.rs

//! Rhumb line (loxodrome) sailing -- the path of constant bearing.
//!
//! Longer than the great circle, but the constant-bearing path actually
//! steered. A rhumb line is also the path a leg drawn on a Mercator
//! chart represents: a straight line on that projection is a constant
//! bearing, the property Mercator was built for. Measuring a leg drawn
//! on a Mercator map as a great circle gives a distance that disagrees
//! with the line under it, so route legs should use this module
//! throughout.

use crate::{EARTH_RADIUS_NM, Position, angle};

/// Guards the `Δφ/Δψ` quotient near the equator, where both go to zero
/// together and the ratio becomes 0/0. Below this the limit `cos φ` is used
/// instead, which is what the quotient converges to.
const MERIDIONAL_EPS: f64 = 1e-12;

/// The inverse Gudermannian-ish term: the Mercator projected latitude,
/// which is what makes a constant bearing a straight line.
fn projected_lat(lat_rad: f64) -> f64 {
    (std::f64::consts::FRAC_PI_4 + lat_rad / 2.0).tan().ln()
}

/// Distance along the rhumb line between two positions, in nautical miles.
#[must_use]
pub fn distance_nm(from: Position, to: Position) -> f64 {
    let (lat1, lat2) = (from.lat_rad(), to.lat_rad());
    let d_lat = lat2 - lat1;
    let d_lon = shortest_d_lon(from.lon_rad(), to.lon_rad());
    let d_proj = projected_lat(lat2) - projected_lat(lat1);

    // The ratio is the local scale between projected and real latitude; at
    // the equator it degenerates and the limit cos(lat) takes over.
    let q = if d_proj.abs() > MERIDIONAL_EPS {
        d_lat / d_proj
    } else {
        lat1.cos()
    };

    (d_lat * d_lat + q * q * d_lon * d_lon).sqrt() * EARTH_RADIUS_NM
}

/// The constant bearing that runs from `from` to `to`, in degrees true.
///
/// Unlike [`crate::great_circle::initial_bearing_deg`] this bearing holds
/// for the whole leg -- steer it and you arrive.
#[must_use]
pub fn bearing_deg(from: Position, to: Position) -> f64 {
    let d_lon = shortest_d_lon(from.lon_rad(), to.lon_rad());
    let d_proj = projected_lat(to.lat_rad()) - projected_lat(from.lat_rad());
    angle::norm_360(d_lon.atan2(d_proj).to_degrees())
}

/// The position reached by holding `bearing_deg` for `distance_nm`.
#[must_use]
pub fn destination(from: Position, bearing_deg: f64, distance_nm: f64) -> Position {
    let angular = distance_nm / EARTH_RADIUS_NM;
    let brg = bearing_deg.to_radians();
    let (lat1, lon1) = (from.lat_rad(), from.lon_rad());

    let d_lat = angular * brg.cos();
    let mut lat2 = lat1 + d_lat;

    // Sailing over a pole on a constant bearing is possible only due north
    // or due south; anything else spirals in without reaching it. Reflecting
    // keeps a nonsensical input from producing a position off the sphere.
    if lat2.abs() > std::f64::consts::FRAC_PI_2 {
        lat2 = if lat2 > 0.0 {
            std::f64::consts::PI - lat2
        } else {
            -std::f64::consts::PI - lat2
        };
    }

    let d_proj = projected_lat(lat2) - projected_lat(lat1);
    let q = if d_proj.abs() > MERIDIONAL_EPS {
        d_lat / d_proj
    } else {
        lat1.cos()
    };
    let lon2 = lon1 + angular * brg.sin() / q;

    Position::new(lat2.to_degrees(), angle::norm_180(lon2.to_degrees()))
}

/// Longitude difference taken the short way round, in radians.
fn shortest_d_lon(lon1_rad: f64, lon2_rad: f64) -> f64 {
    let d = lon2_rad - lon1_rad;
    if d.abs() > std::f64::consts::PI {
        if d > 0.0 {
            d - std::f64::consts::TAU
        } else {
            d + std::f64::consts::TAU
        }
    } else {
        d
    }
}

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

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

    #[test]
    fn along_a_meridian_it_agrees_with_the_great_circle() {
        // A meridian *is* a great circle, so the two sailings must agree
        // exactly -- a good check that the projected-latitude term is right.
        let a = Position::new(45.0, 13.0);
        let b = Position::new(46.0, 13.0);
        assert!(close(distance_nm(a, b), great_circle::distance_nm(a, b), 1e-6));
        assert!(close(bearing_deg(a, b), 0.0, 1e-9));
    }

    #[test]
    fn along_the_equator_it_agrees_too() {
        let a = Position::new(0.0, 10.0);
        let b = Position::new(0.0, 11.0);
        assert!(close(distance_nm(a, b), great_circle::distance_nm(a, b), 1e-6));
        assert!(close(bearing_deg(a, b), 90.0, 1e-9));
    }

    #[test]
    fn along_a_parallel_it_is_longer_than_the_great_circle() {
        // The point of the whole module: holding 090 at 60N costs distance
        // against the great circle that bends north of it.
        let a = Position::new(60.0, -10.0);
        let b = Position::new(60.0, 10.0);
        assert!(distance_nm(a, b) > great_circle::distance_nm(a, b));
        assert!(close(bearing_deg(a, b), 90.0, 1e-9));
    }

    #[test]
    fn destination_inverts_bearing_and_distance() {
        let start = Position::new(45.55, 13.73);
        for bearing in [0.0, 37.5, 90.0, 180.0, 271.3, 359.0] {
            for dist in [0.5, 12.0, 300.0] {
                let end = destination(start, bearing, dist);
                assert!(
                    close(distance_nm(start, end), dist, 1e-6),
                    "distance differs for {bearing}/{dist}"
                );
                assert!(
                    close(bearing_deg(start, end), bearing, 1e-6),
                    "bearing differs for {bearing}/{dist}"
                );
            }
        }
    }

    #[test]
    fn bearing_holds_all_the_way_along_the_leg() {
        // The defining property: sample the leg and every remaining stretch
        // has the same bearing to the end.
        let start = Position::new(43.0, 9.0);
        let end = destination(start, 62.0, 400.0);
        for frac in [0.1, 0.25, 0.5, 0.9] {
            let mid = destination(start, 62.0, 400.0 * frac);
            assert!(
                close(bearing_deg(mid, end), 62.0, 1e-6),
                "bearing drifted at {frac}"
            );
        }
    }

    #[test]
    fn crossing_the_antimeridian_takes_the_short_way() {
        let a = Position::new(10.0, 179.0);
        let b = Position::new(10.0, -179.0);
        let d = distance_nm(a, b);
        assert!(d < 130.0, "took the long way round: {d}");
        assert!(close(bearing_deg(a, b), 90.0, 1e-6));
    }
}