astrodynamics-gnss 0.9.5

//! Ionospheric delay models.
//!
//! This module exposes the single-frequency ionospheric group-delay models used
//! to correct GNSS pseudoranges. The GPS broadcast Klobuchar model and the IONEX
//! vertical-TEC grid path are both reached through the same [`ionosphere_delay`]
//! entry, and the IONEX grid parser is exposed directly as [`Ionex`].
//!
//! All delays returned are group delays and are positive: they increase the
//! measured pseudorange (the carrier-phase advance is the negation of this
//! value). The ionosphere is dispersive, so a delay reported on a carrier other
//! than the model's native L1 is the L1 delay scaled by `(f_l1 / f)^2`.

mod grid;
mod klobuchar;
mod slant;

#[cfg(test)]
mod tests;

use core::f64::consts::PI;

use astrodynamics::time::model::{Instant, InstantRepr};

use crate::frame::Wgs84Geodetic;

pub use grid::Ionex;

pub(crate) use klobuchar::klobuchar_l1_components;
use klobuchar::F_L1;

/// Broadcast Klobuchar alpha/beta coefficients.
///
/// `alpha` are the four coefficients of the cosine-amplitude polynomial (in
/// seconds and seconds-per-semicircle powers); `beta` are the four coefficients
/// of the period polynomial (in seconds and seconds-per-semicircle powers).
/// These are the eight values transmitted in the GPS navigation message.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct KlobucharParams {
    /// Cosine-amplitude polynomial coefficients (a0..a3).
    pub alpha: [f64; 4],
    /// Period polynomial coefficients (b0..b3).
    pub beta: [f64; 4],
}

/// Selects which ionospheric model produces the delay.
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum IonoModel {
    /// GPS broadcast Klobuchar model with its eight alpha/beta coefficients.
    Klobuchar(KlobucharParams),
}

/// Single-frequency ionospheric group delay (code, positive meters).
///
/// Dispatches on `model`. `frequency_hz` is the carrier on which the delay is
/// reported; the model is dispersive, so the delay scales as `1 / f^2`. The
/// returned value is positive meters that increase the pseudorange.
pub fn ionosphere_delay(
    receiver: Wgs84Geodetic,
    elevation_rad: f64,
    azimuth_rad: f64,
    epoch: Instant,
    frequency_hz: f64,
    model: &IonoModel,
) -> f64 {
    match model {
        IonoModel::Klobuchar(params) => klobuchar(
            params,
            receiver,
            elevation_rad,
            azimuth_rad,
            epoch,
            frequency_hz,
        ),
    }
}

/// GPS broadcast Klobuchar ionospheric group delay (positive meters).
///
/// Lower-level entry the Klobuchar arm of [`ionosphere_delay`] calls; exposed so
/// the broadcast model can be used directly. The receiver geodetic
/// latitude/longitude and the satellite azimuth/elevation are converted from
/// radians to the model's published degree boundary, and the GPS second-of-day
/// is taken from `epoch`. The model evaluates the L1 group delay; the result is
/// then scaled to `frequency_hz` by the dispersive `(f_l1 / f)^2` factor.
///
/// Note on bit-exactness: this wrapper converts both angle (radians -> degrees)
/// and time (`epoch` -> GPS second-of-day) at its boundary; both are
/// representation-bound, so the wrapper is NOT bit-exact to a golden expressed
/// in the kernel's native units (degrees and an exact second-of-day) — the
/// difference is at the nanometre level. The 0-ULP parity contract is on the
/// model kernel in those native units; this convenience entry agrees with it to
/// within that conversion bound.
pub fn klobuchar(
    params: &KlobucharParams,
    receiver: Wgs84Geodetic,
    elevation_rad: f64,
    azimuth_rad: f64,
    epoch: Instant,
    frequency_hz: f64,
) -> f64 {
    klobuchar_native(
        params,
        receiver.lat_rad * 180.0 / PI,
        receiver.lon_rad * 180.0 / PI,
        azimuth_rad * 180.0 / PI,
        elevation_rad * 180.0 / PI,
        gps_second_of_day(epoch),
        frequency_hz,
    )
}

/// GPS broadcast Klobuchar group delay in the model's native input units
/// (positive meters).
///
/// Latitude/longitude and azimuth/elevation are in **degrees** (the model's
/// published boundary) and `t_gps_s` is the GPS **second-of-day** in
/// `[0, 86400)`. This is the bit-exact (0-ULP) entry: it feeds the model kernel
/// directly with no angle or time conversion, so a caller holding native inputs
/// (for example the Elixir wrapper, which already has degrees and an integer
/// time of day) gets exactly the reference result. The L1 delay is scaled to
/// `frequency_hz` by the dispersive `(f_l1 / f)^2` factor.
pub fn klobuchar_native(
    params: &KlobucharParams,
    lat_deg: f64,
    lon_deg: f64,
    az_deg: f64,
    el_deg: f64,
    t_gps_s: f64,
    frequency_hz: f64,
) -> f64 {
    let c = klobuchar_l1_components(
        lat_deg,
        lon_deg,
        az_deg,
        el_deg,
        t_gps_s,
        params.alpha,
        params.beta,
    );

    let ratio = F_L1 / frequency_hz;
    c.delay_l1_m * (ratio * ratio)
}

/// IONEX vertical-TEC-grid slant ionospheric group delay (positive meters).
///
/// Maps the parsed [`Ionex`] vertical-TEC grid to the line of sight in the
/// single-layer-model convention: a single-layer pierce point at the product's
/// shell height, an explicit four-term bilinear VTEC per map, a linear-in-time
/// blend between the two maps bracketing `epoch_j2000_s` (holding the endpoint
/// map outside coverage), the `1/sqrt(1 - s^2)` obliquity factor, and the
/// dispersive `40.3e16 / f^2` frequency scaling.
///
/// The receiver geodetic latitude/longitude come from `receiver` (height is
/// unused: the pierce point rides on the IONEX shell, not the antenna height).
/// The epoch is taken as integer J2000 seconds so it lands exactly on the
/// product's own epoch axis, with no float-rounded time entering the temporal
/// bracket. `frequency_hz` is the carrier on which the delay is reported. The
/// returned value is positive meters that increase the pseudorange.
pub fn ionex_slant_delay(
    ionex: &Ionex,
    receiver: Wgs84Geodetic,
    elevation_rad: f64,
    azimuth_rad: f64,
    epoch_j2000_s: i64,
    frequency_hz: f64,
) -> f64 {
    slant::slant_delay_components(
        receiver.lat_rad,
        receiver.lon_rad,
        azimuth_rad,
        elevation_rad,
        frequency_hz,
        ionex.base_radius_km(),
        ionex.shell_height_km(),
        epoch_j2000_s,
        ionex.map_epochs_s(),
        ionex.tec_maps(),
        ionex.lat_nodes_deg(),
        ionex.lon_nodes_deg(),
        ionex.dlat_deg(),
        ionex.dlon_deg(),
    )
    .delay_m
}

/// GPS second-of-day in `[0, 86400)` carried by an instant.
///
/// The Klobuchar diurnal term needs the local-solar-time argument, built from
/// the GPS second-of-day. A Julian date's civil day begins at noon, so the
/// midnight day fraction is `(jd + 0.5)` modulo one.
///
/// Precision: for a split-Julian-date instant the second-of-day is
/// `day_fraction * 86400`, and `day_fraction` is itself a rounded binary
/// fraction of a day, so this recovers the second-of-day only to within the
/// float granularity of a day fraction (a few microseconds) — a
/// sub-nanometre-to-nanometre perturbation in the delay. The bit-exact (0-ULP)
/// contract is on the model kernel evaluated at an exact second-of-day, not on
/// this convenience conversion. An integer-nanosecond instant is exact (it
/// reduces by the seconds-per-day modulus).
fn gps_second_of_day(epoch: Instant) -> f64 {
    match epoch.repr {
        InstantRepr::JulianDate(jd) => {
            // jd_whole is on a *.0 or *.5 boundary; (jd_whole + 0.5) shifts the
            // day origin from noon to midnight. Keep only the part within a day.
            let from_midnight = jd.jd_whole + 0.5 + jd.fraction;
            let day_fraction = from_midnight - from_midnight.floor();
            day_fraction * 86400.0
        }
        InstantRepr::Nanos(nanos) => {
            let ns_per_day: i128 = 86_400 * 1_000_000_000;
            let mut rem = nanos % ns_per_day;
            if rem < 0 {
                rem += ns_per_day;
            }
            rem as f64 / 1.0e9
        }
    }
}