squitter 0.1.1

//! Compact Position Reporting (CPR) latitude/longitude decoding, shared by
//! surface and airborne position messages.
//!
//! Algorithm, the `NL(lat)` boundary table, and the real-position fixture
//! below are cross-checked against pyModeS's `position/_cpr.py` (itself
//! derived from RTCA DO-260B §A.1.7).
//!
//! `core` has no `f64::floor`/`rem_euclid` without a `libm` dependency (only
//! basic arithmetic operators and `abs`, which is bit-manipulation rather
//! than a libm call); the private `floor`/`rem_euclid` helpers below
//! hand-roll exactly what this module needs, matching this crate's (and
//! `aivdm`'s) "hand-roll the small pieces standard library conveniences
//! don't cover in `no_std`" approach rather than pulling in a full math
//! library for two functions.

const CPR_DENOM: f64 = 131_072.0; // 2^17

/// Which CPR encoding (even/odd latitude-zone size) a position report used.
/// A single frame's position is ambiguous on its own; pairing one even and
/// one odd frame resolves it (see [`super::AirbornePosition`]'s
/// `PositionPair` and this module's [`global_decode`]/[`surface_global_decode`]).
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum CprFormat {
    /// Even-format CPR encoding.
    Even,
    /// Odd-format CPR encoding.
    Odd,
}

impl CprFormat {
    pub(crate) const fn from_raw(v: u8) -> Self {
        if v == 0 { Self::Even } else { Self::Odd }
    }

    pub(crate) const fn to_raw(self) -> u8 {
        match self {
            Self::Even => 0,
            Self::Odd => 1,
        }
    }
}

/// `core`-only floor for f64. Correct for any `x` within `i64` range, which
/// every value this module computes always is (CPR math never produces
/// magnitudes anywhere near that large).
fn floor(x: f64) -> f64 {
    #[allow(
        clippy::cast_possible_truncation,
        clippy::cast_precision_loss,
        reason = "CPR intermediate values are always tiny (well under i64::MAX \
                  or f64's 52-bit mantissa range)"
    )]
    let truncated = x as i64 as f64;
    if x < truncated {
        truncated - 1.0
    } else {
        truncated
    }
}

/// `core`-only euclidean remainder for f64 (`y` must be positive), matching
/// Python's `%` operator semantics used throughout the reference algorithm.
fn rem_euclid(x: f64, y: f64) -> f64 {
    let r = x % y;
    if r < 0.0 { r + y } else { r }
}

/// Latitude boundaries (decimal degrees) where `NL(lat)` steps down by one,
/// in ascending order. Entry `i` is where `NL` transitions from `59 - i` to
/// `59 - i - 1`.
// Values copied verbatim (no manual reformatting) from pyModeS's
// `position/_cpr.py` `_NL_BOUNDARIES` to avoid transcription errors -- these
// are cross-referenced against the closed-form trig formula in pyModeS's own
// test suite, so treat any edit here as needing the same re-verification.
#[rustfmt::skip]
#[allow(clippy::unreadable_literal, reason = "copied verbatim from a verified source, not hand-written")]
const NL_BOUNDARIES: [f64; 58] = [
    10.47047129996848,  // NL=59->58
    14.828174368686794, // NL=58->57
    18.186263570713354, // NL=57->56
    21.029394926028463, // NL=56->55
    23.545044865570706, // NL=55->54
    25.829247070587755, // NL=54->53
    27.938987101219045, // NL=53->52
    29.911356857318083, // NL=52->51
    31.77209707681077,  // NL=51->50
    33.53993436298484,  // NL=50->49
    35.22899597796385,  // NL=49->48
    36.85025107593526,  // NL=48->47
    38.41241892412256,  // NL=47->46
    39.922566843338615, // NL=46->45
    41.38651832260239,  // NL=45->44
    42.80914012243555,  // NL=44->43
    44.194549514192744, // NL=43->42
    45.546267226602346, // NL=42->41
    46.867332524987454, // NL=41->40
    48.160391280966216, // NL=40->39
    49.42776439255687,  // NL=39->38
    50.67150165553835,  // NL=38->37
    51.893424691687684, // NL=37->36
    53.09516152796003,  // NL=36->35
    54.278174722729,    // NL=35->34
    55.44378444495043,  // NL=34->33
    56.59318756205918,  // NL=33->32
    57.72747353866114,  // NL=32->31
    58.84763776148457,  // NL=31->30
    59.954592766940294, // NL=30->29
    61.04917774246351,  // NL=29->28
    62.13216659210329,  // NL=28->27
    63.20427479381928,  // NL=27->26
    64.2661652256744,   // NL=26->25
    65.31845309682089,  // NL=25->24
    66.36171008382617,  // NL=24->23
    67.39646774084667,  // NL=23->22
    68.4232202208333,   // NL=22->21
    69.44242631144024,  // NL=21->20
    70.454510749876,    // NL=20->19
    71.45986473028982,  // NL=19->18
    72.45884544728945,  // NL=18->17
    73.45177441667865,  // NL=17->16
    74.43893415725137,  // NL=16->15
    75.42056256653356,  // NL=15->14
    76.39684390794469,  // NL=14->13
    77.36789461328188,  // NL=13->12
    78.33374082922747,  // NL=12->11
    79.29428225456925,  // NL=11->10
    80.24923213280512,  // NL=10->9
    81.19801349271948,  // NL=9->8
    82.13956980510606,  // NL=8->7
    83.07199444719814,  // NL=7->6
    83.99173562980565,  // NL=6->5
    84.89166190702085,  // NL=5->4
    85.75541620944418,  // NL=4->3
    86.535369975121,    // NL=3->2
    87.0,               // NL=2->1
];

/// Number of CPR longitude zones (`1..=59`) for `lat` (decimal degrees).
///
/// Per RTCA DO-260B §A.1.7.2: at the equator `NL=59`; at `±87°` `NL=2`;
/// beyond `±87°` `NL=1`. Uses a hardcoded boundary table for a binary-search
/// lookup rather than the closed-form `arccos`-based formula, since `core`
/// has no `f64::acos`/`cos` without `libm` -- this is also what real-world
/// decoders (`dump1090`, `pyModeS`) do, for the same reason plus speed.
#[must_use]
pub fn nl(lat: f64) -> u8 {
    let abs_lat = lat.abs();
    if abs_lat > 87.0 {
        return 1;
    }
    // exactly 87 deg is a documented special case (the boundary table's
    // last entry is also 87.0, which a plain partition_point lookup would
    // otherwise resolve to NL=1, not the correct NL=2). An exact
    // (non-epsilon) comparison is intentional: this mirrors the reference
    // algorithm's own exact check, not a tolerance-based float comparison.
    #[allow(clippy::float_cmp, reason = "intentional exact match, see above")]
    if abs_lat == 87.0 {
        return 2;
    }
    #[allow(
        clippy::cast_possible_truncation,
        reason = "idx is at most 58 (NL_BOUNDARIES.len()), fits u8"
    )]
    let idx = NL_BOUNDARIES.partition_point(|&boundary| boundary <= abs_lat) as u8;
    59 - idx
}

/// Resolves an absolute latitude/longitude from a paired even-format and
/// odd-format CPR frame (the raw 17-bit `LAT-CPR`/`LON-CPR` fields).
///
/// Per RTCA DO-260B §A.1.7.3. `even_is_newer` selects which frame's
/// latitude zone defines the reported position (the more recently received
/// frame). The two frames must have been broadcast close enough in time
/// (in practice, within about 10 seconds) that the aircraft stayed in the
/// same latitude zone; this function performs an `NL`-equality check as a
/// lightweight sanity guard, returning `None` if it fails, but cannot
/// detect a same-zone false positive caused by a wide time gap -- that's
/// the caller's responsibility.
///
/// # Errors
/// Returns `None` if the frames fall in different CPR latitude zones, or if
/// the resolved position is outside the physically valid range (typically
/// indicating the frames shouldn't have been paired).
#[must_use]
pub fn global_decode(
    lat_cpr_even_raw: u32,
    lon_cpr_even_raw: u32,
    lat_cpr_odd_raw: u32,
    lon_cpr_odd_raw: u32,
    even_is_newer: bool,
) -> Option<(f64, f64)> {
    let cpr_lat_even = f64::from(lat_cpr_even_raw) / CPR_DENOM;
    let cpr_lon_even = f64::from(lon_cpr_even_raw) / CPR_DENOM;
    let cpr_lat_odd = f64::from(lat_cpr_odd_raw) / CPR_DENOM;
    let cpr_lon_odd = f64::from(lon_cpr_odd_raw) / CPR_DENOM;

    let j = floor(59.0 * cpr_lat_even - 60.0 * cpr_lat_odd + 0.5);

    let mut lat_even = (360.0 / 60.0) * (rem_euclid(j, 60.0) + cpr_lat_even);
    let mut lat_odd = (360.0 / 59.0) * (rem_euclid(j, 59.0) + cpr_lat_odd);
    if lat_even >= 270.0 {
        lat_even -= 360.0;
    }
    if lat_odd >= 270.0 {
        lat_odd -= 360.0;
    }

    if nl(lat_even) != nl(lat_odd) {
        return None;
    }

    let lat = if even_is_newer { lat_even } else { lat_odd };
    let nl_here = f64::from(nl(lat));
    let m = floor(cpr_lon_even * (nl_here - 1.0) - cpr_lon_odd * nl_here + 0.5);

    let (ni, cpr_lon_selected) = if even_is_newer {
        (nl_here.max(1.0), cpr_lon_even)
    } else {
        ((nl_here - 1.0).max(1.0), cpr_lon_odd)
    };

    let mut lon = (360.0 / ni) * (rem_euclid(m, ni) + cpr_lon_selected);
    if lon > 180.0 {
        lon -= 360.0;
    }

    if lat.abs() > 90.0 || lon.abs() > 180.0 {
        return None;
    }

    Some((lat, lon))
}

/// Resolves an absolute latitude/longitude from a paired even-format and
/// odd-format *surface* CPR frame (the raw 17-bit `LAT-CPR`/`LON-CPR`
/// fields from a [`super::SurfacePosition`] message pair).
///
/// Per RTCA DO-260B §A.1.7.4. Surface CPR spans only a 90° latitude zone
/// (versus 360° for airborne) and four longitude quadrants, so -- unlike
/// [`global_decode`] -- a nearby reference position (typically the
/// receiving station's own location) is required to disambiguate which
/// hemisphere and quadrant the aircraft is actually in; the reference must
/// be within about 45 nautical miles of the true position.
///
/// # Errors
/// Returns `None` if the frames fall in different CPR latitude zones.
#[must_use]
#[allow(
    clippy::similar_names,
    reason = "north/south suffixes are the clearest naming here"
)]
pub fn surface_global_decode(
    lat_cpr_even_raw: u32,
    lon_cpr_even_raw: u32,
    lat_cpr_odd_raw: u32,
    lon_cpr_odd_raw: u32,
    lat_ref: f64,
    lon_ref: f64,
    even_is_newer: bool,
) -> Option<(f64, f64)> {
    let cpr_lat_even = f64::from(lat_cpr_even_raw) / CPR_DENOM;
    let cpr_lon_even = f64::from(lon_cpr_even_raw) / CPR_DENOM;
    let cpr_lat_odd = f64::from(lat_cpr_odd_raw) / CPR_DENOM;
    let cpr_lon_odd = f64::from(lon_cpr_odd_raw) / CPR_DENOM;

    let j = floor(59.0 * cpr_lat_even - 60.0 * cpr_lat_odd + 0.5);

    // surface CPR's 90-degree zone is ambiguous between two hemispheres;
    // resolve using whichever is actually closer to the reference position,
    // not the reference's own hemisphere (a receiver and aircraft can
    // straddle the equator while still within the 45 NM limit).
    let lat_even_n = (90.0 / 60.0) * (rem_euclid(j, 60.0) + cpr_lat_even);
    let lat_odd_n = (90.0 / 59.0) * (rem_euclid(j, 59.0) + cpr_lat_odd);
    let lat_even_s = lat_even_n - 90.0;
    let lat_odd_s = lat_odd_n - 90.0;

    let newer_n = if even_is_newer { lat_even_n } else { lat_odd_n };
    let newer_s = if even_is_newer { lat_even_s } else { lat_odd_s };
    let (lat_even, lat_odd) = if (lat_ref - newer_n).abs() <= (lat_ref - newer_s).abs() {
        (lat_even_n, lat_odd_n)
    } else {
        (lat_even_s, lat_odd_s)
    };

    if nl(lat_even) != nl(lat_odd) {
        return None;
    }

    let lat = if even_is_newer { lat_even } else { lat_odd };
    let nl_here = f64::from(nl(lat));
    let m = floor(cpr_lon_even * (nl_here - 1.0) - cpr_lon_odd * nl_here + 0.5);

    let (ni, cpr_lon_selected) = if even_is_newer {
        (nl_here.max(1.0), cpr_lon_even)
    } else {
        ((nl_here - 1.0).max(1.0), cpr_lon_odd)
    };
    let lon_base = (90.0 / ni) * (rem_euclid(m, ni) + cpr_lon_selected);

    // four candidate longitudes (one per 90-degree quadrant), each wrapped
    // to [-180, 180]; pick whichever is circularly closest to lon_ref.
    let mut best_lon = lon_base;
    let mut best_dist = f64::INFINITY;
    for q in [0.0, 90.0, 180.0, 270.0] {
        let candidate = rem_euclid(lon_base + q + 180.0, 360.0) - 180.0;
        let dist = (rem_euclid(candidate - lon_ref + 180.0, 360.0) - 180.0).abs();
        if dist < best_dist {
            best_dist = dist;
            best_lon = candidate;
        }
    }

    Some((lat, best_lon))
}

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

    #[test]
    fn nl_matches_known_reference_values() {
        assert_eq!(nl(0.0), 59);
        assert_eq!(nl(87.0), 2);
        assert_eq!(nl(87.1), 1);
        assert_eq!(nl(-87.0), 2);
        assert_eq!(nl(10.0), 59);
        assert_eq!(nl(10.5), 58);
    }

    #[test]
    fn floor_matches_std_for_representative_values() {
        assert!((floor(2.3) - 2.0).abs() < 1e-9);
        assert!((floor(-2.3) - -3.0).abs() < 1e-9);
        assert!((floor(-2.0) - -2.0).abs() < 1e-9);
        assert!((floor(0.0) - 0.0).abs() < 1e-9);
    }

    #[test]
    fn rem_euclid_is_always_nonnegative() {
        assert!((rem_euclid(-1.0, 60.0) - 59.0).abs() < 1e-9);
        assert!((rem_euclid(59.0, 60.0) - 59.0).abs() < 1e-9);
        assert!((rem_euclid(61.0, 60.0) - 1.0).abs() < 1e-9);
    }

    // Real even/odd airborne position pair (DF17, ICAO 40058B), CPR fields
    // extracted from hex messages 8D40058B58C901375147EFD09357 (even) and
    // 8D40058B58C904A87F402D3B8C59 (odd); expected results independently
    // verified by pyModeS's own `airborne_position_pair` test suite.
    const EVEN_LAT_CPR: u32 = 39848;
    const EVEN_LON_CPR: u32 = 83951;
    const ODD_LAT_CPR: u32 = 21567;
    const ODD_LON_CPR: u32 = 81965;

    #[test]
    fn decodes_a_real_captured_position_pair_odd_newer() {
        let (lat, lon) =
            global_decode(EVEN_LAT_CPR, EVEN_LON_CPR, ODD_LAT_CPR, ODD_LON_CPR, false).unwrap();
        assert!((lat - 49.817_55).abs() < 1e-3);
        assert!((lon - 6.084_42).abs() < 1e-3);
    }

    #[test]
    fn decodes_a_real_captured_position_pair_even_newer() {
        let (lat, lon) =
            global_decode(EVEN_LAT_CPR, EVEN_LON_CPR, ODD_LAT_CPR, ODD_LON_CPR, true).unwrap();
        assert!((lat - 49.824_10).abs() < 1e-3);
        assert!((lon - 6.067_85).abs() < 1e-3);
    }

    #[test]
    fn rejects_mismatched_zones() {
        // real vector from pyModeS's test suite: cpr_lat_even=8192 (lat
        // ~66.38 deg, NL=23) paired with cpr_lat_odd=114688 (lat ~66.36
        // deg, NL=24) -- different latitude zones, must be rejected.
        assert_eq!(global_decode(8192, 0, 114_688, 0, true), None);
    }

    #[test]
    fn rejects_impossible_latitude() {
        // real DF17 pair from OpenSky (ICAO 485A33) that slips past the
        // NL-equality check but resolves to |lat| > 90 -- the final sanity
        // guard must reject it. Regression fixture for a real bug pyModeS
        // hit: prior to the guard this resolved to lat ~113 deg.
        assert_eq!(global_decode(113_939, 119_217, 72_761, 59_498, false), None);
    }

    #[test]
    fn surface_decode_resolves_a_real_lfbo_taxiway_pair() {
        // real even/odd surface pair (DF17, ICAO 3A23FF taxiing at LFBO
        // Toulouse-Blagnac), CPR fields extracted from hex messages
        // 903a23ff426a38565950432ebf95 (even) and 903a23ff426a4e65f7487a775d17
        // (odd); expected result independently verified by pyModeS's own
        // `surface_position_pair` test suite, resolved against the airport
        // reference position.
        let (lat, lon) =
            surface_global_decode(11052, 86083, 78587, 84090, 43.63, 1.37, false).unwrap();
        assert!((lat - 43.626_46).abs() < 1e-3);
        assert!((lon - 1.374_76).abs() < 1e-3);
    }
}