geo.distance

Imports name this capability’s declared dependencies, which fune builds next to it in your project; each one links to its page.

use super::funejson::Value;  ← the fune runtime: the JSON value the test vectors use; fune build keeps it only where a signature takes one

// Why this does not call f64::sin, cos or atan2: those come from the platform
// maths library, which may differ in the last bit between Rust, JavaScript and
// Python, and a value next to a rounding boundary can then round differently.
// The trigonometry below uses only operations IEEE 754 requires to be
// correctly rounded (+ - * / sqrt floor), in the same order as the other two
// implementations, so the unrounded result is the same double everywhere.
// Rust never contracts a * b + c into a fused multiply-add on its own.

/// IUGG mean radius R1 = (2a + b) / 3 of the GRS 80 ellipsoid (Moritz,
/// "Geodetic Reference System 1980", Journal of Geodesy 74 (2000) 128-133).
const EARTH_RADIUS_METRES: f64 = 6371008.8;

const PI: f64 = 3.141592653589793;
const HALF_PI: f64 = PI / 2.0;
const DEGREES: f64 = PI / 180.0;

const S3: f64 = -1.0 / 6.0;
const S5: f64 = 1.0 / 120.0;
const S7: f64 = -1.0 / 5040.0;
const S9: f64 = 1.0 / 362880.0;
const S11: f64 = -1.0 / 39916800.0;
const S13: f64 = 1.0 / 6227020800.0;
const S15: f64 = -1.0 / 1307674368000.0;
const S17: f64 = 1.0 / 355687428096000.0;
const C2: f64 = -1.0 / 2.0;
const C4: f64 = 1.0 / 24.0;
const C6: f64 = -1.0 / 720.0;
const C8: f64 = 1.0 / 40320.0;
const C10: f64 = -1.0 / 3628800.0;
const C12: f64 = 1.0 / 479001600.0;
const C14: f64 = -1.0 / 87178291200.0;
const C16: f64 = 1.0 / 20922789888000.0;
const C18: f64 = -1.0 / 6402373705728000.0;
const A3: f64 = -1.0 / 3.0;
const A5: f64 = 1.0 / 5.0;
const A7: f64 = -1.0 / 7.0;
const A9: f64 = 1.0 / 9.0;
const A11: f64 = -1.0 / 11.0;
const A13: f64 = 1.0 / 13.0;
const A15: f64 = -1.0 / 15.0;
const A17: f64 = 1.0 / 17.0;
const A19: f64 = -1.0 / 19.0;
const A21: f64 = 1.0 / 21.0;

fn sin_small(r: f64) -> f64 {
    let s = r * r;
    let mut p = S17;
    p = S15 + s * p;
    p = S13 + s * p;
    p = S11 + s * p;
    p = S9 + s * p;
    p = S7 + s * p;
    p = S5 + s * p;
    p = S3 + s * p;
    r + r * s * p
}

fn cos_small(r: f64) -> f64 {
    let s = r * r;
    let mut p = C18;
    p = C16 + s * p;
    p = C14 + s * p;
    p = C12 + s * p;
    p = C10 + s * p;
    p = C8 + s * p;
    p = C6 + s * p;
    p = C4 + s * p;
    p = C2 + s * p;
    1.0 + s * p
}

fn sin_cos(x: f64) -> (f64, f64) {
    let k = (x / HALF_PI + 0.5).floor();
    let r = x - k * HALF_PI;
    let sr = sin_small(r);
    let cr = cos_small(r);
    match (k as i64).rem_euclid(4) {
        0 => (sr, cr),
        1 => (cr, -sr),
        2 => (-sr, -cr),
        _ => (-cr, sr),
    }
}

fn atan_unit(u: f64) -> f64 {
    let mut v = u;
    v = v / (1.0 + (1.0 + v * v).sqrt());
    v = v / (1.0 + (1.0 + v * v).sqrt());
    v = v / (1.0 + (1.0 + v * v).sqrt());
    let s = v * v;
    let mut p = A21;
    p = A19 + s * p;
    p = A17 + s * p;
    p = A15 + s * p;
    p = A13 + s * p;
    p = A11 + s * p;
    p = A9 + s * p;
    p = A7 + s * p;
    p = A5 + s * p;
    p = A3 + s * p;
    8.0 * (v + v * s * p)
}

fn atan2_positive(y: f64, x: f64) -> f64 {
    if x == 0.0 {
        return if y == 0.0 { 0.0 } else { HALF_PI };
    }
    let t = y / x;
    if t > 1.0 {
        HALF_PI - atan_unit(x / y)
    } else {
        atan_unit(t)
    }
}

/// Half away from zero on the binary64 value; + 0.0 turns -0.0 into 0.0.
fn round_to(x: f64, scale: f64) -> f64 {
    let y = x.abs() * scale;
    let mut r = y.floor();
    if y - r >= 0.5 {
        r += 1.0;
    }
    let out = r / scale;
    (if x < 0.0 { -out } else { out }) + 0.0
}

fn check_latitude(value: f64) {
    if !value.is_finite() {
        panic!("latitude must be a finite number of degrees");
    }
    if value < -90.0 || value > 90.0 {
        panic!("latitude must be between -90 and 90 degrees");
    }
}

fn check_longitude(value: f64) {
    if !value.is_finite() {
        panic!("longitude must be a finite number of degrees");
    }
    if value < -180.0 || value > 180.0 {
        panic!("longitude must be between -180 and 180 degrees");
    }
}

/// Great-circle distance in metres (haversine, IUGG mean radius), to the millimetre.
///
/// # Panics
/// Panics if a coordinate is not finite or out of range.
pub fn distance(from_lat: f64, from_lng: f64, to_lat: f64, to_lng: f64) -> f64 {
    check_latitude(from_lat);
    check_longitude(from_lng);
    check_latitude(to_lat);
    check_longitude(to_lng);

    let phi1 = from_lat * DEGREES;
    let phi2 = to_lat * DEGREES;
    // A longitude difference of 359 degrees is a 1 degree step across the
    // antimeridian; sin^2 of the half-angle is the same either way.
    let sin_half_lat = sin_cos(((to_lat - from_lat) * DEGREES) / 2.0).0;
    let sin_half_lng = sin_cos(((to_lng - from_lng) * DEGREES) / 2.0).0;
    let cos1 = sin_cos(phi1).1;
    let cos2 = sin_cos(phi2).1;

    let mut a = sin_half_lat * sin_half_lat + cos1 * cos2 * sin_half_lng * sin_half_lng;
    if a < 0.0 {
        a = 0.0;
    }
    if a > 1.0 {
        a = 1.0;
    }
    // atan2 rather than asin(sqrt(a)): stays accurate for antipodal points too.
    let c = 2.0 * atan2_positive(a.sqrt(), (1.0 - a).sqrt());
    round_to(EARTH_RADIUS_METRES * c, 1000.0)
}

pub fn fune_vector(args: &[Value]) -> Value {
    // Refuse what the typed signature cannot hold, with the wording TypeScript
    // and Python use, rather than let the conversion below quietly change it.
    if !matches!(args[0], Value::Int(_) | Value::Float(_)) {
        panic!("latitude must be a finite number of degrees");
    }
    if !matches!(args[1], Value::Int(_) | Value::Float(_)) {
        panic!("longitude must be a finite number of degrees");
    }
    if !matches!(args[2], Value::Int(_) | Value::Float(_)) {
        panic!("latitude must be a finite number of degrees");
    }
    if !matches!(args[3], Value::Int(_) | Value::Float(_)) {
        panic!("longitude must be a finite number of degrees");
    }
    Value::Float(distance(
        args[0].as_f64(),
        args[1].as_f64(),
        args[2].as_f64(),
        args[3].as_f64(),
    ))
}