astrodynamics-gnss 0.8.0

//! Tropospheric delay model.
//!
//! Computes the neutral-atmosphere (tropospheric) delay on a GNSS signal as a
//! Saastamoinen (1972) zenith hydrostatic and wet delay, driven by supplied
//! surface meteorology, mapped to the line of sight by the Niell (1996)
//! continued-fraction mapping functions (NMF). The zenith primitives and the
//! mapping primitives are exposed separately so a caller can apply the distinct
//! hydrostatic and wet mappings, and a convenience entry composes the full
//! slant delay.
//!
//! The tropospheric delay is non-dispersive: it has the same sign and magnitude
//! for code and carrier phase, and it is a positive additive range error. The
//! returned delays are positive meters that increase the measured pseudorange;
//! `delay_m > 0` means the signal arrived later and the pseudorange is too long
//! by `delay_m`.
//!
//! Angles are radians internally (`_rad`); height is the WGS84 ellipsoidal
//! height in meters carried by [`Wgs84Geodetic`]. Pressure is hectopascals,
//! temperature is kelvin, and relative humidity is a unit fraction in `[0, 1]`.

mod saastamoinen;

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

use crate::frame::Wgs84Geodetic;

pub(crate) use saastamoinen::slant_components;

/// Surface meteorological conditions at the receiver.
///
/// These are the inputs to the Saastamoinen zenith delays. Pressure is in
/// hectopascals (millibars) because that is the unit the troposphere formulas
/// and meteorological products use; temperature is absolute (kelvin) to avoid a
/// Celsius zero-point slip; relative humidity is a unit fraction in `[0, 1]`,
/// not a percentage.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Met {
    /// Total atmospheric pressure in hectopascals (= millibars).
    pub pressure_hpa: f64,
    /// Ambient temperature in kelvin.
    pub temperature_k: f64,
    /// Relative humidity as a unit fraction in `[0, 1]` (0.5 == 50%).
    pub relative_humidity: f64,
}

impl Met {
    /// Construct surface meteorology from explicit values.
    pub const fn new(pressure_hpa: f64, temperature_k: f64, relative_humidity: f64) -> Self {
        Self {
            pressure_hpa,
            temperature_k,
            relative_humidity,
        }
    }

    /// Standard-atmosphere pressure and temperature for an ellipsoidal height,
    /// carrying the supplied relative humidity through unchanged.
    ///
    /// A convenience for callers without live meteorology. The height is clamped
    /// to be non-negative before the standard pressure and temperature lapses
    /// are applied, so a below-sea-level height yields the sea-level values.
    pub fn standard(height_m: f64, relative_humidity: f64) -> Self {
        let s = saastamoinen::standard_atmosphere(height_m, relative_humidity);
        Self {
            pressure_hpa: s.pressure_hpa,
            temperature_k: s.temperature_k,
            relative_humidity: s.relative_humidity,
        }
    }
}

/// Tropospheric zenith-delay model selector.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum TropoModel {
    /// Saastamoinen (1972) zenith hydrostatic + wet delays.
    Saastamoinen,
}

/// Tropospheric mapping-function selector.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum MappingModel {
    /// Niell (1996) mapping functions (NMF).
    Niell,
}

/// Zenith tropospheric delay split into its hydrostatic and wet parts.
///
/// The two components are returned separately because the Niell mapping applies
/// a distinct hydrostatic and wet mapping factor; the total slant delay is
/// `dry_m * mapping.dry + wet_m * mapping.wet`.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct ZenithDelay {
    /// Zenith hydrostatic (dry) delay in positive meters.
    pub dry_m: f64,
    /// Zenith wet delay in positive meters.
    pub wet_m: f64,
}

/// Dimensionless mapping factors at a given elevation.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct MappingFactors {
    /// Hydrostatic mapping factor (includes the height correction).
    pub dry: f64,
    /// Wet mapping factor.
    pub wet: f64,
}

/// Zenith tropospheric delay (hydrostatic + wet) from supplied meteorology.
///
/// The receiver geodetic latitude and ellipsoidal height come from `receiver`;
/// the surface pressure, temperature, and humidity come from `met`. Both
/// components are returned as positive meters. The possibly-negative ellipsoidal
/// height is used with its sign.
pub fn tropo_zenith(model: TropoModel, receiver: Wgs84Geodetic, met: Met) -> ZenithDelay {
    match model {
        TropoModel::Saastamoinen => {
            let z = saastamoinen::zenith_delays(
                met.pressure_hpa,
                met.temperature_k,
                met.relative_humidity,
                receiver.lat_rad,
                receiver.height_m,
            );
            ZenithDelay {
                dry_m: z.zhd_m,
                wet_m: z.zwd_m,
            }
        }
    }
}

/// Niell hydrostatic and wet mapping factors at the given elevation.
///
/// The mapping depends on the receiver geodetic latitude and ellipsoidal height
/// (via `receiver`) and on the fractional day-of-year (derived from `epoch`),
/// hence both are arguments. The hydrostatic mapping carries the height
/// correction; the wet mapping does not.
pub fn tropo_mapping(
    model: MappingModel,
    elevation_rad: f64,
    receiver: Wgs84Geodetic,
    epoch: Instant,
) -> MappingFactors {
    match model {
        MappingModel::Niell => {
            let doy = fractional_day_of_year(epoch);
            let m = saastamoinen::niell_mapping(
                elevation_rad,
                receiver.lat_rad,
                receiver.height_m,
                doy,
            );
            MappingFactors {
                dry: m.mh,
                wet: m.mw,
            }
        }
    }
}

/// Full slant tropospheric delay in positive meters.
///
/// Composes the Saastamoinen zenith delays and the Niell mapping into the total
/// line-of-sight delay `dry_m * mapping.dry + wet_m * mapping.wet`. The delay is
/// zero for an elevation at or below the horizon and for a height outside the
/// model's validity range; inside that range the result is positive.
///
/// Note on bit-exactness: the fractional day-of-year is derived from `epoch` at
/// this boundary. Its integer day is exact (from the calendar date); the
/// within-day fraction carries the float granularity of a split Julian date,
/// which enters the Niell seasonal term so weakly that it is below the last bit
/// in practice. The 0-ULP parity contract is on the kernel evaluated at an
/// exact day-of-year; this wrapper agrees with it to within that bound.
pub fn tropo_slant(elevation_rad: f64, receiver: Wgs84Geodetic, met: Met, epoch: Instant) -> f64 {
    let doy = fractional_day_of_year(epoch);
    slant_components(
        elevation_rad,
        receiver.lat_rad,
        receiver.lon_rad,
        receiver.height_m,
        met.pressure_hpa,
        met.temperature_k,
        met.relative_humidity,
        doy,
    )
    .slant_m
}

/// Fractional day-of-year carried by an instant, Jan 1 00:00:00 = 1.0.
///
/// The Niell seasonal term needs the continuous day-of-year (it carries the
/// fractional time of day). This converts the instant to a calendar date and
/// returns the day number plus the within-day fraction. The conversion is the
/// standard civil-date algorithm from the Julian day number; the within-day
/// fraction comes from the day part shifted from the noon Julian-date origin to
/// midnight.
fn fractional_day_of_year(epoch: Instant) -> f64 {
    let jd = match epoch.repr {
        InstantRepr::JulianDate(split) => split.jd_whole + split.fraction,
        InstantRepr::Nanos(nanos) => {
            // Nanoseconds since the J2000 Julian-date origin 2451545.0.
            (nanos as f64) / 1.0e9 / 86_400.0 + 2_451_545.0
        }
    };

    // Shift the noon origin to midnight, then split into an integer day and the
    // within-day fraction.
    let jd_midnight = jd + 0.5;
    let day_floor = jd_midnight.floor();
    let day_fraction = jd_midnight - day_floor;

    let (year, month, day) = civil_from_jd_floor(day_floor);
    let doy_integer = day_of_year(year, month, day);
    doy_integer as f64 + day_fraction
}

/// Calendar (year, month, day) from the floor of a midnight-shifted Julian date.
///
/// This is the standard Fliegel-Van Flandern conversion operating on the integer
/// Julian day number.
fn civil_from_jd_floor(day_floor: f64) -> (i64, i64, i64) {
    let jdn = day_floor as i64;
    let l = jdn + 68_569;
    let n = (4 * l) / 146_097;
    let l = l - (146_097 * n + 3) / 4;
    let i = (4_000 * (l + 1)) / 1_461_001;
    let l = l - (1_461 * i) / 4 + 31;
    let j = (80 * l) / 2_447;
    let day = l - (2_447 * j) / 80;
    let l = j / 11;
    let month = j + 2 - 12 * l;
    let year = 100 * (n - 49) + i + l;
    (year, month, day)
}

/// Day-of-year (Jan 1 = 1) for a civil date, Gregorian leap-year rule.
fn day_of_year(year: i64, month: i64, day: i64) -> i64 {
    const CUM: [i64; 12] = [0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334];
    let leap = (year % 4 == 0 && year % 100 != 0) || (year % 400 == 0);
    let mut doy = CUM[(month - 1) as usize] + day;
    if leap && month > 2 {
        doy += 1;
    }
    doy
}

#[cfg(test)]
mod tests;