diff --git a/docs/src/api.md b/docs/src/api.md index 53b4f08..23ff545 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -22,6 +22,20 @@ calc_satellite_position_and_velocity get_sat_enu ``` +## Atmospheric Corrections + +These corrections are applied automatically by [`calc_pvt`](@ref); they are +documented here for reference and for diagnostic use. + +```@docs +PositionVelocityTime.select_ionospheric_correction +PositionVelocityTime.ionospheric_delay +PositionVelocityTime.tropospheric_delay +PositionVelocityTime.KlobucharParams +PositionVelocityTime.NTCMGParams +PositionVelocityTime._elevation_azimuth +``` + ## Dilution of Precision ```@docs diff --git a/src/PositionVelocityTime.jl b/src/PositionVelocityTime.jl index 39075d4..b4e50a6 100755 --- a/src/PositionVelocityTime.jl +++ b/src/PositionVelocityTime.jl @@ -186,7 +186,9 @@ end """ calc_pvt(states::AbstractVector{<:SatelliteState}, prev_pvt::PVTSolution = PVTSolution(); - approximate_year::Integer = year(now(UTC))) -> PVTSolution + approximate_year::Integer = year(now(UTC)), + enable_ionospheric_correction::Bool = true, + enable_tropospheric_correction::Bool = true) -> PVTSolution Calculate Position, Velocity, and Time (PVT) from GNSS satellite measurements. @@ -194,6 +196,18 @@ Requires at least 4 healthy satellites from the same GNSS system. Uses least-squ estimation for position and time, and solves for velocity and clock drift from carrier Doppler measurements. +Unless disabled via `enable_ionospheric_correction`, the ionospheric delay is +corrected automatically using only the coefficients decoded from the navigation +messages. A single model is chosen for the whole solve and applied to every +satellite: NTCM-G if Galileo Effective Ionisation Level coefficients have been +decoded (the more accurate model), otherwise the Klobuchar model if GPS L1 α/β +have been decoded, otherwise no correction. See +[`select_ionospheric_correction`](@ref) and [`ionospheric_delay`](@ref). + +Unless disabled via `enable_tropospheric_correction`, the tropospheric delay is +corrected with the blind Saastamoinen model (no broadcast coefficients needed). +See [`tropospheric_delay`](@ref). + # Arguments - `states`: Vector of [`SatelliteState`](@ref) for observed satellites - `prev_pvt`: Previous PVT solution used as initial guess (default: origin) @@ -206,6 +220,12 @@ carrier Doppler measurements. ±9 years of the actual observation date works. Defaults to the current UTC year, which is correct for live signals; for processing archived recordings, pass the rough year of the recording. +- `enable_ionospheric_correction`: when `true` (default), apply the automatic + ionospheric correction described above. Set to `false` to skip it entirely + and solve from the raw pseudoranges (e.g. for diagnostics or when an external + correction is applied elsewhere). +- `enable_tropospheric_correction`: when `true` (default), apply the Saastamoinen + tropospheric correction. Set to `false` to skip it. # Returns A [`PVTSolution`](@ref) containing position, velocity, time, DOP values, and @@ -219,6 +239,8 @@ function calc_pvt( states::AbstractVector{<:SatelliteState}, prev_pvt::PVTSolution = PVTSolution(); approximate_year::Integer = year(now(UTC)), + enable_ionospheric_correction::Bool = true, + enable_tropospheric_correction::Bool = true, ) length(states) < 4 && throw(ArgumentError("You'll need at least 4 satellites to calculate PVT")) @@ -238,7 +260,62 @@ function calc_pvt( ) sat_positions = map(get_sat_position, sat_positions_and_velocities) pseudo_ranges, reference_time = calc_pseudo_ranges(times) - ξ, rmse = user_position(sat_positions, pseudo_ranges, prev_ξ) + # Atmospheric corrections, summed per satellite and subtracted from the + # pseudoranges. The ionospheric model is chosen for the whole solve from the + # coefficients decoded across the constellation (NTCM-G if Galileo coefficients + # are available, else Klobuchar if GPS α/β are available, else none; see + # `select_ionospheric_correction`); the troposphere uses the blind Saastamoinen + # model (see `tropospheric_delay`). A single corrected solve is enough: the + # delays depend on position only through the satellite elevation/azimuth (and, + # for the troposphere, the user height), and ∂delay/∂position is negligible over + # the metre-level position uncertainty (a 15 m shift moves the elevation by + # ~1e-5°), so delays predicted at a nearby position are accurate to well under a + # millimetre — no iterate-to-convergence needed. + correction = + enable_ionospheric_correction ? select_ionospheric_correction(healthy_states) : + nothing + function predict_atmospheric_delays(position) + user_pos = ECEF(position[1], position[2], position[3]) + # The user position is the same for every satellite, so the geodetic + # coordinates and the ENU transform are computed once per epoch and reused. + user_lla = LLAfromECEF(wgs84)(user_pos) + enu_from_ecef = ENUfromECEF(user_pos, wgs84) + map(healthy_states, sat_positions) do state, sat_pos + elevation, azimuth = _elevation_azimuth(enu_from_ecef, sat_pos) + iono = ionospheric_delay( + correction, + state.system, + elevation, + azimuth, + user_lla, + reference_time, + ) + tropo = + enable_tropospheric_correction ? + tropospheric_delay(elevation, user_lla) : 0.0 + iono + tropo + end + end + ξ, rmse = if iszero(prev_ξ) + # Cold start: no prior position, and the Klobuchar model is undefined near + # the geocenter, so first obtain an approximate fix from an uncorrected + # solve, then re-solve once with the delay-corrected pseudoranges (only if + # there is anything to correct, so the uncorrected case stays a single solve). + ξ_uncorrected, rmse_uncorrected = + user_position(sat_positions, pseudo_ranges, prev_ξ) + atmospheric_delays = predict_atmospheric_delays(ξ_uncorrected) + any(!iszero, atmospheric_delays) ? + user_position(sat_positions, pseudo_ranges .- atmospheric_delays, ξ_uncorrected) : + (ξ_uncorrected, rmse_uncorrected) + else + # Warm start: predict the delays from the previous (already metre-accurate) + # position before solving, so ξ never needs a post-solve correction. + user_position( + sat_positions, + pseudo_ranges .- predict_atmospheric_delays(prev_ξ), + prev_ξ, + ) + end sat_positions_mat = reduce(hcat, sat_positions) H = calc_H(sat_positions_mat, ξ) user_velocity_and_clock_drift = calc_user_velocity_and_clock_drift( @@ -360,4 +437,6 @@ end include("user_position.jl") include("sat_time.jl") include("sat_position.jl") +include("ionosphere.jl") +include("troposphere.jl") end diff --git a/src/ionosphere.jl b/src/ionosphere.jl new file mode 100644 index 0000000..d828f9a --- /dev/null +++ b/src/ionosphere.jl @@ -0,0 +1,333 @@ +# =========================================================================== +# Ionospheric correction — constellation-wide model selection +# +# The broadcast ionospheric coefficients are global to a GNSS (every GPS +# satellite broadcasts the same Klobuchar α/β; every Galileo satellite the same +# NTCM-G Effective Ionisation Level a_i0..a_i2). So rather than choosing a model +# per satellite, a single model is chosen for the whole solve from whatever the +# healthy decoders have delivered, and applied to *all* satellites: +# +# - No coefficients decoded → no correction (0 m for every sat). +# - Only Klobuchar (GPS) decoded → Klobuchar for every sat. +# - Only NTCM-G (Galileo) decoded → NTCM-G for every sat. +# - Both decoded → NTCM-G for every sat (more accurate). +# +# Only data actually decoded from the navigation message is used; there are no +# user-supplied fallback coefficients. The ionospheric delay scales as 1/f², so a +# model is evaluated for the constellation it came from and then rescaled to each +# satellite's actual carrier frequency. This makes a single set of coefficients +# applicable across constellations and frequency bands (L1/E1, L2, L5/E5a, E6, …), +# not just the band it was broadcast on. +# =========================================================================== + +# The frequency the Klobuchar group delay refers to (GPS L1), taken from GNSSSignals +# so it stays consistent with the per-satellite carrier frequencies used below. Kept +# as a Hz quantity: only the (f_L1/f) ratio is used, which is dimensionless. +const _GPS_L1_FREQUENCY = GNSSSignals.get_center_frequency(GPSL1CA()) + +""" + KlobucharParams(α_0, α_1, α_2, α_3, β_0, β_1, β_2, β_3) + +The eight Klobuchar ionospheric coefficients decoded from a GPS L1 navigation +message, in IS-GPS-200 SI units (seconds and seconds·semicircle⁻ⁿ). The field +names mirror `GNSSDecoder.GPSL1Data` (`α_0…α_3`, `β_0…β_3`). +""" +struct KlobucharParams + α_0::Float64 + α_1::Float64 + α_2::Float64 + α_3::Float64 + β_0::Float64 + β_1::Float64 + β_2::Float64 + β_3::Float64 +end + +""" + NTCMGParams(a_i0, a_i1, a_i2, week_number::Integer) + +The broadcast Galileo Effective Ionisation Level coefficients `a_i0`/`a_i1`/`a_i2` +(sfu, sfu/deg, sfu/deg²) decoded from an E1B navigation message, together with the +GST week number needed to derive the day of year and universal time for NTCM-G. +""" +struct NTCMGParams + a_i0::Float64 + a_i1::Float64 + a_i2::Float64 + week_number::Int +end + +""" + klobuchar_params(decoder) -> Union{KlobucharParams,Nothing} + +Klobuchar α/β decoded from a GPS L1 navigation message, or `nothing` if they have +not been broadcast yet (subframe 4, page 18) or the decoder is not GPS L1. +""" +klobuchar_params(decoder) = nothing +function klobuchar_params(decoder::GNSSDecoder.GNSSDecoderState{<:GNSSDecoder.GPSL1Data}) + d = decoder.data + (isnothing(d.α_0) || isnothing(d.β_0)) && return nothing + return KlobucharParams(d.α_0, d.α_1, d.α_2, d.α_3, d.β_0, d.β_1, d.β_2, d.β_3) +end + +""" + ntcm_g_params(decoder) -> Union{NTCMGParams,Nothing} + +NTCM-G Effective Ionisation Level coefficients decoded from a Galileo E1B +navigation message, or `nothing` if they (or the week number) have not been +decoded yet or the decoder is not Galileo E1B. +""" +ntcm_g_params(decoder) = nothing +function ntcm_g_params(decoder::GNSSDecoder.GNSSDecoderState{<:GNSSDecoder.GalileoE1BData}) + d = decoder.data + (isnothing(d.a_i0) || isnothing(d.a_i1) || isnothing(d.a_i2) || isnothing(d.WN)) && + return nothing + return NTCMGParams(d.a_i0, d.a_i1, d.a_i2, d.WN) +end + +""" + select_ionospheric_correction(states) -> Union{KlobucharParams,NTCMGParams,Nothing} + +Scan all (healthy) satellite decoders and pick the single ionospheric correction +to apply to the whole solve. NTCM-G is preferred whenever Galileo coefficients are +available (it is the more accurate model); otherwise Klobuchar is used if GPS +coefficients are available; otherwise `nothing` (no correction). The coefficients +are global to a constellation, so the first decoder that carries each set is used. +""" +function select_ionospheric_correction(states) + klobuchar = nothing + ntcm_g = nothing + for state in states + klobuchar === nothing && (klobuchar = klobuchar_params(state.decoder)) + ntcm_g === nothing && (ntcm_g = ntcm_g_params(state.decoder)) + end + ntcm_g !== nothing && return ntcm_g # most accurate model when available + return klobuchar # KlobucharParams, or nothing if neither +end + +""" + ionospheric_delay(correction, system, elevation, azimuth, lla, gps_time) -> Float64 + +Slant ionospheric group delay in metres for one satellite (`system` is the +satellite's GNSS, used for its carrier frequency), using the constellation-wide +`correction` returned by [`select_ionospheric_correction`](@ref): + +- `::Nothing` → `0.0` (no coefficients were decoded). +- [`KlobucharParams`](@ref) → Klobuchar model (IS-GPS-200, Fig. 20-4). +- [`NTCMGParams`](@ref) → NTCM-G model. + +The line of sight is given by the satellite `elevation`/`azimuth` (radians) and +the user geodetic position `lla` (a `Geodesy.LLA`); `gps_time` is the system time +of week of the measurement in seconds. The geometry is taken precomputed — and +shared across satellites and with [`tropospheric_delay`](@ref) — so a whole-epoch +correction does the user geodetic conversion only once. Derive the geometry from +ECEF with `LLAfromECEF(wgs84)(user)` and +[`_elevation_azimuth`](@ref)`(ENUfromECEF(user, wgs84), sat)`. +""" +ionospheric_delay(::Nothing, system, elevation, azimuth, lla, gps_time) = 0.0 + +function ionospheric_delay(p::KlobucharParams, system, elevation, azimuth, lla, gps_time) + # IS-GPS-200 works in semicircles: lat/lon in deg/180, elevation/azimuth in rad/π. + l1_seconds = klobuchar_group_delay( + lla.lat / 180, + lla.lon / 180, + elevation / π, + azimuth / π, + gps_time, + (p.α_0, p.α_1, p.α_2, p.α_3), + (p.β_0, p.β_1, p.β_2, p.β_3), + ) + # The Klobuchar broadcast coefficients define the group delay at the GPS L1 + # frequency (IS-GPS-200). The ionospheric delay scales as 1/f², so rescale it + # to this satellite's actual carrier frequency — yielding the correct (larger) + # delay on lower bands such as L2, L5/E5a and E6. For an L1/E1 signal the + # factor is exactly 1. The ratio of two Hz quantities is dimensionless. + f = GNSSSignals.get_center_frequency(system) + return SPEEDOFLIGHT * l1_seconds * (_GPS_L1_FREQUENCY / f)^2 +end + +function ionospheric_delay(p::NTCMGParams, system, elevation, azimuth, lla, gps_time) + doy, ut = _galileo_doy_and_ut(p.week_number, gps_time) + stec = ntcm_g_stec(elevation, azimuth, lla, doy, ut, p.a_i0, p.a_i1, p.a_i2) # TECU + f = ustrip(Hz, GNSSSignals.get_center_frequency(system)) + # Eq. 1: group delay [m] = 40.3 / f² · STEC, with STEC in electrons/m² (1 TECU = 1e16). + return 40.3 / f^2 * stec * 1.0e16 +end + +""" + klobuchar_group_delay(φ_u, λ_u, E, A, gps_time, α, β) + +Klobuchar single-frequency ionospheric group delay for GPS L1 (IS-GPS-200N, +Fig. 20-4), returned in **seconds**. All angles are in **semicircles**: +`φ_u`/`λ_u` are the user geodetic latitude/longitude, `E`/`A` the satellite +elevation/azimuth. `gps_time` is GPS system time in seconds. `α`/`β` are the +4-element Klobuchar coefficient tuples (SI units). +""" +function klobuchar_group_delay(φ_u, λ_u, E, A, gps_time, α, β) + # Earth-centred angle between user and ionospheric pierce point (semicircles) + ψ = 0.0137 / (E + 0.11) - 0.022 + # Geodetic latitude of the ionospheric pierce point (IPP), clamped per ICD + φ_i = clamp(φ_u + ψ * cos(A * π), -0.416, 0.416) + # Geodetic longitude of the IPP + λ_i = λ_u + ψ * sin(A * π) / cos(φ_i * π) + # Geomagnetic latitude of the IPP + φ_m = φ_i + 0.064 * cos((λ_i - 1.617) * π) + # Local time at the IPP (seconds), wrapped to [0, 86400) + t = mod(4.32e4 * λ_i + gps_time, 86400.0) + # Obliquity / slant factor + F = 1.0 + 16.0 * (0.53 - E)^3 + # Amplitude (s) and period (s) of the cosine model, with ICD floors + AMP = max(α[1] + φ_m * (α[2] + φ_m * (α[3] + φ_m * α[4])), 0.0) + PER = max(β[1] + φ_m * (β[2] + φ_m * (β[3] + φ_m * β[4])), 72000.0) + x = 2π * (t - 50400.0) / PER + return abs(x) < 1.57 ? F * (5.0e-9 + AMP * (1 - x^2 / 2 + x^4 / 24)) : F * 5.0e-9 +end + +""" + _elevation_azimuth(enu_from_ecef::ENUfromECEF, sat_position) -> (elevation, azimuth) + +Elevation and azimuth (radians) of `sat_position` (ECEF) in the local East-North-Up +frame defined by `enu_from_ecef = ENUfromECEF(user_position, wgs84)`. Azimuth is +measured clockwise from North. The transform is taken precomputed so it can be built +once per user position and reused across satellites. +""" +function _elevation_azimuth(enu_from_ecef::ENUfromECEF, sat_position) + enu = enu_from_ecef(ECEF(sat_position)) + elevation = atan(enu.u, hypot(enu.e, enu.n)) + azimuth = atan(enu.e, enu.n) + return elevation, azimuth +end + +# =========================================================================== +# NTCM-G — Galileo single-frequency ionospheric model +# (European GNSS (Galileo) NTCM-G Ionospheric Model Description, Issue 1.0, +# May 2022). Driven by the broadcast Effective Ionisation Level coefficients +# a_i0, a_i1, a_i2. Returns slant TEC (TECU); equation numbers below refer to +# that document. +# =========================================================================== + +# NTCM-G model coefficients k1..k12 (Table 3) +const _NTCM_K = ( + 0.92519, + 0.16951, + 0.00443, + 0.06626, + 0.00899, + 0.21289, + -0.15414, + -0.38439, + 1.14023, + 1.20556, + 1.41808, + 0.13985, +) +const _NTCM_RE = 6371.0 # Earth mean radius [km] (Table 2) +const _NTCM_HI = 450.0 # ionospheric pierce point height [km] +const _NTCM_GNP_LAT = deg2rad(79.74) # geomagnetic North pole latitude +const _NTCM_GNP_LON = deg2rad(-71.78) # geomagnetic North pole longitude + +# Effective Ionisation Level Azpar [sfu] from the broadcast coefficients (Eq. 2). +function _azpar(a_i0, a_i1, a_i2) + radicand = a_i0^2 + 1633.33 * a_i1^2 + 4802000.0 * a_i2^2 + 3266.67 * a_i0 * a_i2 + return sqrt(max(radicand, 0.0)) +end + +# Ionospheric pierce point geographic latitude/longitude [rad] (Eq. 24-26), +# given user geodetic lat/lon and satellite elevation/azimuth [rad]. +function _pierce_point(φ_u, λ_u, elevation, azimuth) + ψ = π / 2 - elevation - asin(_NTCM_RE / (_NTCM_RE + _NTCM_HI) * cos(elevation)) + φ_pp = asin(sin(φ_u) * cos(ψ) + cos(φ_u) * sin(ψ) * cos(azimuth)) + λ_pp = λ_u + asin(sin(ψ) * sin(azimuth) / cos(φ_pp)) + return φ_pp, λ_pp +end + +# Sun's declination [rad] for the day of year (Eq. 28). +_sun_declination(doy) = deg2rad(23.44) * sin(deg2rad(0.9856 * (doy - 80.7))) + +# Modified Single Layer Model mapping function (Eq. 32-33). +function _mslm_mapping_function(elevation) + sinz = _NTCM_RE / (_NTCM_RE + _NTCM_HI) * sin(0.9782 * (π / 2 - elevation)) + return 1.0 / sqrt(1.0 - sinz^2) +end + +""" + ntcm_g_vtec(φ_pp, λ_pp, doy, ut, azpar) -> Float64 + +Vertical TEC in TECU at the ionospheric pierce point (geographic latitude/longitude +`φ_pp`/`λ_pp` in radians) for day of year `doy`, universal time `ut` in hours, and +Effective Ionisation Level `azpar` in solar flux units. Implements the NTCM-G model +`VTEC = F1·F2·F3·F4·F5` (Eq. 3-15). +""" +function ntcm_g_vtec(φ_pp, λ_pp, doy, ut, azpar) + k = _NTCM_K + lt = ut + rad2deg(λ_pp) / 15 # local time [h] (Eq. 27) + δ = _sun_declination(doy) + + # Solar zenith angle dependence (Eq. 29-30) + cosχ3 = cos(φ_pp - δ) + 0.4 # cosχ*** (PF1 = 0.4) + cosχ2 = cos(φ_pp - δ) - (2 / π) * φ_pp * sin(δ) # cosχ** + + # F1 — local-time dependency (Eq. 4-7) + V_D = 2π * (lt - 14) / 24 + V_SD = 2π * lt / 12 + V_TD = 2π * lt / 8 + F1 = + cosχ3 + + cosχ2 * ( + k[1] * cos(V_D) + + k[2] * cos(V_SD) + + k[3] * sin(V_SD) + + k[4] * cos(V_TD) + + k[5] * sin(V_TD) + ) + + # F2 — seasonal dependency (Eq. 8-10) + V_A = 2π * (doy - 18) / 365.25 + V_SA = 4π * (doy - 6) / 365.25 + F2 = 1 + k[6] * cos(V_A) + k[7] * cos(V_SA) + + # Geomagnetic latitude of the pierce point (Eq. 31) + φ_m = asin( + sin(φ_pp) * sin(_NTCM_GNP_LAT) + + cos(φ_pp) * cos(_NTCM_GNP_LAT) * cos(λ_pp - _NTCM_GNP_LON), + ) + + # F3 — geomagnetic field dependency (Eq. 11), φ_m in radians + F3 = 1 + k[8] * cos(φ_m) + + # F4 — equatorial (Appleton) anomaly dependency (Eq. 12-14), φ_m in degrees + φ_m_deg = rad2deg(φ_m) + EC1 = -(φ_m_deg - 16.0)^2 / (2 * 12.0^2) + EC2 = -(φ_m_deg + 10.0)^2 / (2 * 13.0^2) + F4 = 1 + k[9] * exp(EC1) + k[10] * exp(EC2) + + # F5 — solar activity dependency (Eq. 15) + F5 = k[11] + k[12] * azpar + + return max(F1 * F2 * F3 * F4 * F5, 0.0) +end + +""" + ntcm_g_stec(elevation, azimuth, lla, doy, ut, a_i0, a_i1, a_i2) -> Float64 + +Slant TEC in TECU along the user→satellite line of sight using NTCM-G, for the +satellite `elevation`/`azimuth` (radians) seen from the user geodetic position +`lla` (a `Geodesy.LLA`), day of year `doy`, universal time `ut` (hours), and the +broadcast Galileo Effective Ionisation Level coefficients `a_i0`/`a_i1`/`a_i2`. +""" +function ntcm_g_stec(elevation, azimuth, lla, doy, ut, a_i0, a_i1, a_i2) + φ_pp, λ_pp = _pierce_point(deg2rad(lla.lat), deg2rad(lla.lon), elevation, azimuth) + vtec = ntcm_g_vtec(φ_pp, λ_pp, doy, ut, _azpar(a_i0, a_i1, a_i2)) + return _mslm_mapping_function(elevation) * vtec +end + +# Day of year and universal time (hours) from a Galileo System Time (GST) week +# number and time of week [s]. GST epoch is 1999-08-22 00:00:00 UTC; GST is +# continuous (offset from UTC by leap seconds, ~18 s), which is negligible for +# the day-of-year / UT inputs of NTCM-G. +function _galileo_doy_and_ut(week_number, time_of_week) + total_seconds = week_number * 604800 + time_of_week + days = floor(Int, total_seconds / 86400) + ut = (total_seconds - days * 86400) / 3600 + return dayofyear(Date(1999, 8, 22) + Day(days)), ut +end diff --git a/src/troposphere.jl b/src/troposphere.jl new file mode 100644 index 0000000..5987d29 --- /dev/null +++ b/src/troposphere.jl @@ -0,0 +1,56 @@ +# =========================================================================== +# Tropospheric correction — Saastamoinen model +# +# The tropospheric delay is non-dispersive — the same on every GNSS frequency. +# Unlike the ionosphere there are no broadcast coefficients: the delay is a +# "blind" function of the user height and a standard atmosphere, mapped to the +# line of sight by the satellite elevation. This is the Saastamoinen model as +# implemented in RTKLIB's `tropmodel` (and used by GNSS-SDR / PocketSDR for +# single-point positioning), with the relative humidity fixed at 70 % and a +# simple 1/cos(z) obliquity mapping. +# =========================================================================== + +const _DEFAULT_RELATIVE_HUMIDITY = 0.7 # GNSS-SDR / RTKLIB default + +""" + tropospheric_delay(elevation, lla; humidity = 0.7) -> Float64 + +Slant tropospheric group delay in metres for a satellite at `elevation` (radians) +seen from the user geodetic position `lla` (a `Geodesy.LLA`), using the +Saastamoinen model driven by a standard atmosphere. The delay is non-dispersive +(frequency independent). `humidity` is the relative humidity (0…1) used for the +wet component. Returns `0.0` when the satellite is at or below the horizon or the +user height is outside the model's valid range (−100 m … 10 km). + +The geometry is taken precomputed so a whole-epoch correction shares the user +geodetic conversion and the elevation with [`ionospheric_delay`](@ref). +""" +function tropospheric_delay(elevation, lla; humidity = _DEFAULT_RELATIVE_HUMIDITY) + return saastamoinen_delay(deg2rad(lla.lat), lla.alt, elevation, humidity) +end + +""" + saastamoinen_delay(latitude, height, elevation, humidity) -> Float64 + +Saastamoinen slant tropospheric delay in metres. `latitude`/`elevation` are in +radians and `height` is the geodetic (ellipsoidal) height in metres. The standard +atmosphere (pressure, temperature, water-vapour partial pressure) is derived from +the height with a 15 °C sea-level temperature, and the hydrostatic and wet zenith +delays are mapped to the slant direction by `1/cos(z)` with `z = π/2 − elevation`. +Mirrors RTKLIB's `tropmodel`. +""" +function saastamoinen_delay(latitude, height, elevation, humidity) + (height < -100.0 || height > 1.0e4 || elevation <= 0.0) && return 0.0 + hgt = height < 0.0 ? 0.0 : height + # Standard atmosphere + pressure = 1013.25 * (1.0 - 2.2557e-5 * hgt)^5.2568 # hPa + temperature = 15.0 - 6.5e-3 * hgt + 273.16 # K (15 °C at sea level) + e = 6.108 * humidity * exp((17.15 * temperature - 4684.0) / (temperature - 38.45)) # hPa + # Saastamoinen hydrostatic and wet slant delays (1/cos(z) obliquity mapping) + z = π / 2 - elevation + trph = + 0.0022768 * pressure / + (1.0 - 0.00266 * cos(2.0 * latitude) - 0.00028 * hgt / 1.0e3) / cos(z) + trpw = 0.002277 * (1255.0 / temperature + 0.05) * e / cos(z) + return trph + trpw +end diff --git a/test/ionosphere.jl b/test/ionosphere.jl new file mode 100644 index 0000000..2c751d2 --- /dev/null +++ b/test/ionosphere.jl @@ -0,0 +1,265 @@ +# Independent RTKLIB-style reference implementation (inputs in radians) used to +# cross-check the package's semicircle-based `klobuchar_group_delay`. +function klobuchar_reference(lat, lon, az, el, gps_time, α, β) + ψ = 0.0137 / (el / π + 0.11) - 0.022 + φ = clamp(lat / π + ψ * cos(az), -0.416, 0.416) + λ = lon / π + ψ * sin(az) / cos(φ * π) + φ += 0.064 * cos((λ - 1.617) * π) + tt = 43200.0 * λ + gps_time + tt -= floor(tt / 86400.0) * 86400.0 + f = 1.0 + 16.0 * (0.53 - el / π)^3 + amp = max(α[1] + φ * (α[2] + φ * (α[3] + φ * α[4])), 0.0) + per = max(β[1] + φ * (β[2] + φ * (β[3] + φ * β[4])), 72000.0) + x = 2π * (tt - 50400.0) / per + return f * (abs(x) < 1.57 ? 5e-9 + amp * (1 - x^2 / 2 + x^4 / 24) : 5e-9) +end + +@testset "Klobuchar ionospheric model" begin + # Exemplary Klobuchar coefficients (IS-GPS-200 SI units) + α = (4.6566129e-09, 1.4901161e-08, -5.96046e-08, -5.96046e-08) + β = (79872.0, 65536.0, -65536.0, -393216.0) + + @testset "klobuchar_group_delay matches independent reference" begin + for lat in deg2rad.((10.0, 48.0, -33.0)), + lon in deg2rad.((-120.0, 11.0, 150.0)), + az in deg2rad.((0.0, 95.0, 270.0)), + el in deg2rad.((5.0, 25.0, 89.0)), + t in (0.0, 50400.0, 80000.0) + + got = PositionVelocityTime.klobuchar_group_delay( + lat / π, + lon / π, + el / π, + az / π, + t, + α, + β, + ) + ref = klobuchar_reference(lat, lon, az, el, t, α, β) + @test got ≈ ref rtol = 1e-12 + end + end + + @testset "physical behaviour" begin + latsc, lonsc, t = 48 / 180, 11 / 180, 50400.0 # local ~14:00 → near diurnal peak + d90 = + PositionVelocityTime.klobuchar_group_delay(latsc, lonsc, 90 / 180, 0.0, t, α, β) + d30 = + PositionVelocityTime.klobuchar_group_delay(latsc, lonsc, 30 / 180, 0.0, t, α, β) + d10 = + PositionVelocityTime.klobuchar_group_delay(latsc, lonsc, 10 / 180, 0.0, t, α, β) + @test d90 > 0 + @test d30 > d90 # obliquity increases slant delay at low elevation + @test d10 > d30 + @test 299792458.0 * d90 < 15 # zenith L1 delay is a few metres, not tens + end + + # --- helpers: build decoders carrying broadcast coefficients --- + function gps_decoder_with(α, β) + dec = GNSSDecoderState(GPSL1CA(), 1) + GNSSDecoder.GNSSDecoderState( + dec; + data = GNSSDecoder.GPSL1Data( + dec.data; + α_0 = α[1], + α_1 = α[2], + α_2 = α[3], + α_3 = α[4], + β_0 = β[1], + β_1 = β[2], + β_2 = β[3], + β_3 = β[4], + ), + ) + end + function galileo_decoder_with(a_i0, a_i1, a_i2, WN) + dec = GNSSDecoderState(GalileoE1B(), 1) + GNSSDecoder.GNSSDecoderState( + dec; + data = GNSSDecoder.GalileoE1BData(dec.data; a_i0, a_i1, a_i2, WN), + ) + end + mkstate(dec, sys) = SatelliteState(; + decoder = dec, + system = sys, + code_phase = 0.0, + carrier_doppler = 0.0Hz, + ) + + @testset "parameter extraction from decoders" begin + # Fresh decoders carry no coefficients + @test PositionVelocityTime.klobuchar_params(GNSSDecoderState(GPSL1CA(), 1)) === + nothing + @test PositionVelocityTime.ntcm_g_params(GNSSDecoderState(GalileoE1B(), 1)) === + nothing + # Broadcast α/β are extracted into KlobucharParams + kp = PositionVelocityTime.klobuchar_params(gps_decoder_with(α, β)) + @test kp isa PositionVelocityTime.KlobucharParams + @test (kp.α_0, kp.α_1, kp.α_2, kp.α_3) == α + @test (kp.β_0, kp.β_1, kp.β_2, kp.β_3) == β + # Broadcast a_i0..a_i2 + WN are extracted into NTCMGParams + np = PositionVelocityTime.ntcm_g_params( + galileo_decoder_with(236.831641, -0.39362878, 0.00402826613, 1100), + ) + @test np isa PositionVelocityTime.NTCMGParams + @test np.a_i0 == 236.831641 + @test np.week_number == 1100 + # Wrong-system decoders never yield the other model's params + @test PositionVelocityTime.ntcm_g_params(gps_decoder_with(α, β)) === nothing + @test PositionVelocityTime.klobuchar_params( + galileo_decoder_with(100.0, 0.0, 0.0, 1100), + ) === nothing + end + + @testset "constellation-wide model selection" begin + gps_bare = mkstate(GNSSDecoderState(GPSL1CA(), 1), GPSL1CA()) + gps_klob = mkstate(gps_decoder_with(α, β), GPSL1CA()) + gal_bare = mkstate(GNSSDecoderState(GalileoE1B(), 1), GalileoE1B()) + gal_ntcm = mkstate(galileo_decoder_with(121.13, 0.35, 0.013, 1100), GalileoE1B()) + + # Nothing decoded → no correction + @test PositionVelocityTime.select_ionospheric_correction([gps_bare, gal_bare]) === + nothing + # Only Klobuchar → Klobuchar + @test PositionVelocityTime.select_ionospheric_correction([gps_klob, gps_bare]) isa + PositionVelocityTime.KlobucharParams + # Only Galileo → NTCM-G + @test PositionVelocityTime.select_ionospheric_correction([gal_ntcm, gal_bare]) isa + PositionVelocityTime.NTCMGParams + # Both available → NTCM-G wins (more accurate) + @test PositionVelocityTime.select_ionospheric_correction([gps_klob, gal_ntcm]) isa + PositionVelocityTime.NTCMGParams + end + + @testset "selected model applied to every satellite" begin + user = ECEFfromLLA(wgs84)(LLA(48.0, 11.0, 550.0)) + sat = ECEF(user[1] + 1.0e7, user[2] + 1.0e7, user[3] + 2.0e7) + klob = PositionVelocityTime.KlobucharParams(α..., β...) + # Line-of-sight geometry (same user/sat for every call below) + el, az = PositionVelocityTime._elevation_azimuth(ENUfromECEF(user, wgs84), sat) + lla = LLAfromECEF(wgs84)(user) + + # No correction → exactly zero for any system + @test PositionVelocityTime.ionospheric_delay( + nothing, + GPSL1CA(), + el, + az, + lla, + 50400.0, + ) == 0.0 + @test PositionVelocityTime.ionospheric_delay( + nothing, + GalileoE1B(), + el, + az, + lla, + 50400.0, + ) == 0.0 + # Klobuchar applied to GPS *and* Galileo (E1 shares the L1 frequency, so + # the delay is identical for both systems) + d_gps = + PositionVelocityTime.ionospheric_delay(klob, GPSL1CA(), el, az, lla, 50400.0) + d_gal = + PositionVelocityTime.ionospheric_delay(klob, GalileoE1B(), el, az, lla, 50400.0) + @test d_gps > 0.0 + @test d_gps ≈ d_gal rtol = 1e-12 + + # The delay scales as 1/f², so the same coefficients applied on a lower band + # (here GPS L5, 1176.45 MHz) give the correct larger delay, not the L1 value. + d_l5 = PositionVelocityTime.ionospheric_delay(klob, GPSL5I(), el, az, lla, 50400.0) + # ratio of two Hz quantities is dimensionless + f_ratio = + GNSSSignals.get_center_frequency(GPSL1CA()) / + GNSSSignals.get_center_frequency(GPSL5I()) + @test d_l5 ≈ d_gps * f_ratio^2 rtol = 1e-12 + @test d_l5 > d_gps + end +end + +@testset "NTCM-G ionospheric model (Galileo)" begin + toecef(lat, lon, h) = ECEFfromLLA(wgs84)(LLA(lat, lon, h)) + HI = (236.831641, -0.39362878, 0.00402826613) # high solar activity + MED = (121.129893, 0.351254133, 0.0134635348) # medium + LOW = (2.580271, 0.127628236, 0.0252748384) # low + # Official Input/Output verification data (NTCM-G description, Annex D): + # (a0, a1, a2, doy, ut, stn_lon, stn_lat, stn_h, sat_lon, sat_lat, sat_h, STEC) + vectors = [ + (HI..., 105, 0, -62.34, 82.49, 78.11, 8.23, 54.29, 20281546.18, 33.7567), + (HI..., 105, 12, -62.34, 82.49, 78.11, 81.09, 35.20, 20278071.03, 65.0500), + (HI..., 105, 20, -52.81, 5.25, -25.76, 10.94, 44.72, 20450566.19, 252.0204), + (HI..., 105, 16, -52.81, 5.25, -25.76, -70.26, 50.63, 20043030.72, 216.2278), + (MED..., 105, 4, 40.19, -3.00, -23.32, 107.19, -10.65, 19943686.06, 108.8940), + (MED..., 105, 20, 115.89, -31.80, 12.78, 131.65, -31.56, 20066111.12, 7.5508), + (LOW..., 105, 0, 141.13, 39.14, 117.00, 165.14, -13.93, 20181976.50, 51.5270), + (LOW..., 105, 20, -155.46, 19.80, 3754.69, -82.52, 20.64, 19937791.48, 67.4750), + ] + @testset "official Annex D test vectors" begin + for v in vectors + a0, a1, a2, doy, ut, slon, slat, sh, blon, blat, bh, expected = v + u = toecef(slat, slon, sh) + s = toecef(blat, blon, bh) + el, az = PositionVelocityTime._elevation_azimuth(ENUfromECEF(u, wgs84), s) + stec = PositionVelocityTime.ntcm_g_stec( + el, + az, + LLAfromECEF(wgs84)(u), + doy, + ut, + a0, + a1, + a2, + ) + @test stec ≈ expected atol = 1e-3 # published values rounded to 4 dp + end + end + + @testset "Azpar (Eq. 2)" begin + # |√(a0² + 1633.33·a1² + 4802000·a2² + 3266.67·a0·a2)|, in sfu + @test PositionVelocityTime._azpar(HI...) ≈ 244.007 atol = 1e-2 + @test PositionVelocityTime._azpar(0.0, 0.0, 0.0) == 0.0 + end + + @testset "GST week/TOW → day-of-year and UT" begin + # 1 week + 12 h past the GST epoch (1999-08-22 00:00 UTC) → 1999-08-29, 12:00 + doy, ut = PositionVelocityTime._galileo_doy_and_ut(1, 12 * 3600) + @test doy == dayofyear(Date(1999, 8, 29)) + @test ut ≈ 12.0 + end + + @testset "selection and delay from a Galileo decoder" begin + user = ECEFfromLLA(wgs84)(LLA(48.0, 11.0, 550.0)) + sat = ECEF(user[1] + 1.0e7, user[2] + 1.0e7, user[3] + 2.0e7) + # Galileo decoder carrying broadcast a_i0..a_i2 and a week number + gal = GNSSDecoderState(GalileoE1B(), 1) + gal = GNSSDecoder.GNSSDecoderState( + gal; + data = GNSSDecoder.GalileoE1BData( + gal.data; + a_i0 = 236.831641, + a_i1 = -0.39362878, + a_i2 = 0.00402826613, + WN = 1100, + ), + ) + state = SatelliteState(; + decoder = gal, + system = GalileoE1B(), + code_phase = 0.0, + carrier_doppler = 0.0Hz, + ) + correction = PositionVelocityTime.select_ionospheric_correction([state]) + @test correction isa PositionVelocityTime.NTCMGParams + el, az = PositionVelocityTime._elevation_azimuth(ENUfromECEF(user, wgs84), sat) + delay = PositionVelocityTime.ionospheric_delay( + correction, + GalileoE1B(), + el, + az, + LLAfromECEF(wgs84)(user), + 200000.0, + ) + @test delay > 0.0 + @test delay < 100.0 # a sane L1/E1 ionospheric delay magnitude (metres) + end +end diff --git a/test/pvt.jl b/test/pvt.jl index f66f90f..99e9f2f 100644 --- a/test/pvt.jl +++ b/test/pvt.jl @@ -494,14 +494,25 @@ BitIntegers.@define_integers 320 # this test stays valid as the wall clock drifts past 2030 (when the # default `now()` anchor would start picking the wrong rollover cycle # for these archived fixtures). - pvt = calc_pvt(states; approximate_year = 2021) + pvt = calc_pvt( + states; + approximate_year = 2021, + enable_ionospheric_correction = false, + enable_tropospheric_correction = false, + ) expected_lla = LLA(; lat = 50.778851672464015, lon = 6.065568885758519, alt = 289.4069805158367) @test pvt.position ≈ ECEFfromLLA(wgs84)(expected_lla) rtol = 1e-8 @test pvt.time ≈ TAIEpoch(2021, 5, 31, 12, 53, 14.1183385390904732) @test pvt.velocity ≈ ECEF(0.0, 0.0, 0.0) atol = 9 @test get_frequency_offset(pvt, get_center_frequency(galileo_e1b)) ≈ -(1675.63Hz + freq_offset) atol = 0.01Hz - warm_pvt = calc_pvt(states, pvt; approximate_year = 2021) + warm_pvt = calc_pvt( + states, + pvt; + approximate_year = 2021, + enable_ionospheric_correction = false, + enable_tropospheric_correction = false, + ) @test get_LLA(warm_pvt) ≈ get_LLA(pvt) @test warm_pvt.time ≈ pvt.time @test warm_pvt.velocity ≈ pvt.velocity atol = 1e-6 @@ -1432,14 +1443,25 @@ end # Fixture data was recorded on 2021-05-31. See the note on the # Galileo testset above for why we pin `approximate_year`. - pvt = calc_pvt(states; approximate_year = 2021) + pvt = calc_pvt( + states; + approximate_year = 2021, + enable_ionospheric_correction = false, + enable_tropospheric_correction = false, + ) expected_lla = LLA(; lat = 50.77885249310784, lon = 6.0656199911189175, alt = 291.95658091689086) @test pvt.position ≈ ECEFfromLLA(wgs84)(expected_lla) rtol = 1e-8 @test pvt.time ≈ TAIEpoch(2021, 5, 31, 12, 53, 14.1491024351271335) @test pvt.velocity ≈ ECEF(0.0, 0.0, 0.0) atol = 2.5 @test get_frequency_offset(pvt, get_center_frequency(gpsl1)) ≈ -(1632.59Hz + freq_offset) atol = 0.01Hz - warm_pvt = calc_pvt(states, pvt; approximate_year = 2021) + warm_pvt = calc_pvt( + states, + pvt; + approximate_year = 2021, + enable_ionospheric_correction = false, + enable_tropospheric_correction = false, + ) @test get_LLA(warm_pvt) ≈ get_LLA(pvt) @test warm_pvt.time ≈ pvt.time @test warm_pvt.velocity ≈ pvt.velocity atol = 1e-6 diff --git a/test/pvt_iono_tropo.jl b/test/pvt_iono_tropo.jl new file mode 100644 index 0000000..237fb55 --- /dev/null +++ b/test/pvt_iono_tropo.jl @@ -0,0 +1,945 @@ +using LinearAlgebra: norm + +# Exact-width unsigned types backing the captured decoder bit buffers +# (288-bit for Galileo E1B, 320-bit for GPS L1 C/A). Guarded so the file is +# self-contained yet does not clash with the same definitions in `pvt.jl`. +isdefined(@__MODULE__, :UInt288) || BitIntegers.@define_integers 288 +isdefined(@__MODULE__, :UInt320) || BitIntegers.@define_integers 320 + +# Ground truth position +const IONO_TROPO_GROUND_TRUTH = ECEFfromLLA(wgs84)(LLA(48.0, 11.0, 550.0)) + +@testset "calc_pvt with iono + tropo corrections" begin + @testset "GPS L1 C/A — Klobuchar" begin + gpsl1 = GPSL1CA() + states = [ + SatelliteState(; + decoder = GNSSDecoderState( + 9, + uint320"0xfffffac0000014ffffaf741b36875df7ac0111755517d5cdf8e80000ed4d11000f800004f20e38e3", + uint320"0xf208baaa8beae71d80000000000000f8d15980b47fffffd6000000a7fffd7ba0d9b43aefbd60088b", + GNSSDecoder.GPSL1Data( + last_subframe_id = 2, + integrity_status_flag = false, + TOW = 259242, + alert_flag = false, + anti_spoof_flag = false, + trans_week = 58, + codeonl2 = 1, + ura = 2.0, + svhealth = "000000", + IODC = "0000000000", + l2pcode = false, + T_GD = 0.0, + t_0c = 266400, + a_f2 = 0.0, + a_f1 = 0.0, + a_f0 = 0.0, + IODE_Sub_2 = "00000000", + C_rs = 0.0, + Δn = 0.0, + M_0 = 1.5441823417813125, + C_uc = 0.0, + e = 0.0, + C_us = 0.0, + sqrt_A = 5153.700811386108, + t_0e = 266400, + fit_interval = false, + AODO = 31, + C_ic = 0.0, + Ω_0 = 2.905155877506184, + C_is = 0.0, + i_0 = 0.9599310884343369, + C_rc = 0.0, + ω = 1.4594111442521003, + Ω_dot = 0.0, + IODE_Sub_3 = "00000000", + i_dot = 0.0, + α_0 = nothing, + α_1 = nothing, + α_2 = nothing, + α_3 = nothing, + β_0 = nothing, + β_1 = nothing, + β_2 = nothing, + β_3 = nothing, + ), + GNSSDecoder.GPSL1Data( + last_subframe_id = 2, + integrity_status_flag = false, + TOW = 259242, + alert_flag = false, + anti_spoof_flag = false, + trans_week = 58, + codeonl2 = 1, + ura = 2.0, + svhealth = "000000", + IODC = "0000000000", + l2pcode = false, + T_GD = 0.0, + t_0c = 266400, + a_f2 = 0.0, + a_f1 = 0.0, + a_f0 = 0.0, + IODE_Sub_2 = "00000000", + C_rs = 0.0, + Δn = 0.0, + M_0 = 1.5441823417813125, + C_uc = 0.0, + e = 0.0, + C_us = 0.0, + sqrt_A = 5153.700811386108, + t_0e = 266400, + fit_interval = false, + AODO = 31, + C_ic = 0.0, + Ω_0 = 2.905155877506184, + C_is = 0.0, + i_0 = 0.9599310884343369, + C_rc = 0.0, + ω = 1.4594111442521003, + Ω_dot = 0.0, + IODE_Sub_3 = "00000000", + i_dot = 0.0, + α_0 = 4.6566129e-9, + α_1 = 1.4901161e-8, + α_2 = -5.96046e-8, + α_3 = -5.96046e-8, + β_0 = 79872.0, + β_1 = 65536.0, + β_2 = -65536.0, + β_3 = -393216.0, + ), + GNSSDecoder.GPSL1Constants( + 300, + 0x8b, + 8, + 30, + 3.1415926535898, + 7.2921151467e-5, + 2.99792458e8, + 3.986005e14, + -4.442807633e-10, + ), + GNSSDecoder.GPSL1Cache(), + 173, + 173, + false, + ), + system = gpsl1, + code_phase = 7457.6350011159975, + carrier_doppler = 945.8683607586818Hz, + carrier_phase = -0.22095823550737081, + ), + SatelliteState(; + decoder = GNSSDecoderState( + 11, + uint320"0x00000000000000000050d81b3686120853ef6e8aaae82a320717fffe20cd0c5197800004f20e38e3", + uint320"0xf208baaa8beae71d800000000000017f5f5201ecfffffffffffffffffffd793f264bcf6fbd60848b", + GNSSDecoder.GPSL1Data( + last_subframe_id = 2, + integrity_status_flag = false, + TOW = 259242, + alert_flag = false, + anti_spoof_flag = false, + trans_week = 58, + codeonl2 = 1, + ura = 2.0, + svhealth = "000000", + IODC = "0000000000", + l2pcode = false, + T_GD = 0.0, + t_0c = 266400, + a_f2 = 0.0, + a_f1 = 0.0, + a_f0 = 0.0, + IODE_Sub_2 = "00000000", + C_rs = 0.0, + Δn = 0.0, + M_0 = 2.3326933798024854, + C_uc = 0.0, + e = 0.0, + C_us = 0.0, + sqrt_A = 5153.700811386108, + t_0e = 266400, + fit_interval = false, + AODO = 31, + C_ic = 0.0, + Ω_0 = -0.4081865199555507, + C_is = 0.0, + i_0 = 0.9599310884343369, + C_rc = 0.0, + ω = 0.6840685603956145, + Ω_dot = 0.0, + IODE_Sub_3 = "00000000", + i_dot = 0.0, + α_0 = nothing, + α_1 = nothing, + α_2 = nothing, + α_3 = nothing, + β_0 = nothing, + β_1 = nothing, + β_2 = nothing, + β_3 = nothing, + ), + GNSSDecoder.GPSL1Data( + last_subframe_id = 2, + integrity_status_flag = false, + TOW = 259242, + alert_flag = false, + anti_spoof_flag = false, + trans_week = 58, + codeonl2 = 1, + ura = 2.0, + svhealth = "000000", + IODC = "0000000000", + l2pcode = false, + T_GD = 0.0, + t_0c = 266400, + a_f2 = 0.0, + a_f1 = 0.0, + a_f0 = 0.0, + IODE_Sub_2 = "00000000", + C_rs = 0.0, + Δn = 0.0, + M_0 = 2.3326933798024854, + C_uc = 0.0, + e = 0.0, + C_us = 0.0, + sqrt_A = 5153.700811386108, + t_0e = 266400, + fit_interval = false, + AODO = 31, + C_ic = 0.0, + Ω_0 = -0.4081865199555507, + C_is = 0.0, + i_0 = 0.9599310884343369, + C_rc = 0.0, + ω = 0.6840685603956145, + Ω_dot = 0.0, + IODE_Sub_3 = "00000000", + i_dot = 0.0, + α_0 = 4.6566129e-9, + α_1 = 1.4901161e-8, + α_2 = -5.96046e-8, + α_3 = -5.96046e-8, + β_0 = 79872.0, + β_1 = 65536.0, + β_2 = -65536.0, + β_3 = -393216.0, + ), + GNSSDecoder.GPSL1Constants( + 300, + 0x8b, + 8, + 30, + 3.1415926535898, + 7.2921151467e-5, + 2.99792458e8, + 3.986005e14, + -4.442807633e-10, + ), + GNSSDecoder.GPSL1Cache(), + 173, + 173, + true, + ), + system = gpsl1, + code_phase = 6078.218551898637, + carrier_doppler = -3276.264988960822Hz, + carrier_phase = 0.7638002718956874, + ), + SatelliteState(; + decoder = GNSSDecoderState( + 17, + uint320"0xffffffffffffffffffaf27e4c979edf7ac1091755517d5cdf8e800019ad6016e98bffffb190e38e3", + uint320"0xf208baaa8beae71d8000000000000367d702b006fffffffffffffffffffd793f264bcf6fbd60848b", + GNSSDecoder.GPSL1Data( + last_subframe_id = 2, + integrity_status_flag = false, + TOW = 259242, + alert_flag = false, + anti_spoof_flag = false, + trans_week = 58, + codeonl2 = 1, + ura = 2.0, + svhealth = "000000", + IODC = "0000000000", + l2pcode = false, + T_GD = 0.0, + t_0c = 266400, + a_f2 = 0.0, + a_f1 = 0.0, + a_f0 = 0.0, + IODE_Sub_2 = "00000000", + C_rs = 0.0, + Δn = 0.0, + M_0 = -0.9434142898148061, + C_uc = 0.0, + e = 0.0, + C_us = 0.0, + sqrt_A = 5153.700811386108, + t_0e = 266400, + fit_interval = false, + AODO = 31, + C_ic = 0.0, + Ω_0 = -1.2333033919680705, + C_is = 0.0, + i_0 = 0.9599310884343369, + C_rc = 0.0, + ω = -2.750026464429168, + Ω_dot = 0.0, + IODE_Sub_3 = "00000000", + i_dot = 0.0, + α_0 = nothing, + α_1 = nothing, + α_2 = nothing, + α_3 = nothing, + β_0 = nothing, + β_1 = nothing, + β_2 = nothing, + β_3 = nothing, + ), + GNSSDecoder.GPSL1Data( + last_subframe_id = 2, + integrity_status_flag = false, + TOW = 259242, + alert_flag = false, + anti_spoof_flag = false, + trans_week = 58, + codeonl2 = 1, + ura = 2.0, + svhealth = "000000", + IODC = "0000000000", + l2pcode = false, + T_GD = 0.0, + t_0c = 266400, + a_f2 = 0.0, + a_f1 = 0.0, + a_f0 = 0.0, + IODE_Sub_2 = "00000000", + C_rs = 0.0, + Δn = 0.0, + M_0 = -0.9434142898148061, + C_uc = 0.0, + e = 0.0, + C_us = 0.0, + sqrt_A = 5153.700811386108, + t_0e = 266400, + fit_interval = false, + AODO = 31, + C_ic = 0.0, + Ω_0 = -1.2333033919680705, + C_is = 0.0, + i_0 = 0.9599310884343369, + C_rc = 0.0, + ω = -2.750026464429168, + Ω_dot = 0.0, + IODE_Sub_3 = "00000000", + i_dot = 0.0, + α_0 = 4.6566129e-9, + α_1 = 1.4901161e-8, + α_2 = -5.96046e-8, + α_3 = -5.96046e-8, + β_0 = 79872.0, + β_1 = 65536.0, + β_2 = -65536.0, + β_3 = -393216.0, + ), + GNSSDecoder.GPSL1Constants( + 300, + 0x8b, + 8, + 30, + 3.1415926535898, + 7.2921151467e-5, + 2.99792458e8, + 3.986005e14, + -4.442807633e-10, + ), + GNSSDecoder.GPSL1Cache(), + 173, + 173, + false, + ), + system = gpsl1, + code_phase = 17299.53480240955, + carrier_doppler = -189.2252575679783Hz, + carrier_phase = 2.7456602492891107, + ), + SatelliteState(; + decoder = GNSSDecoderState( + 27, + uint320"0xfffffac0000014ffffaf741b36875df7ac0111755517d5cdf8e80000eafce6d996400004e6f1c71c", + uint320"0xf208baaa8beae71d80000000000003d8fb25beea7fffffd6000000a7fffd7ba0d9b43aefbd60088b", + GNSSDecoder.GPSL1Data( + last_subframe_id = 2, + integrity_status_flag = false, + TOW = 259242, + alert_flag = false, + anti_spoof_flag = false, + trans_week = 58, + codeonl2 = 1, + ura = 2.0, + svhealth = "000000", + IODC = "0000000000", + l2pcode = false, + T_GD = 0.0, + t_0c = 266400, + a_f2 = 0.0, + a_f1 = 0.0, + a_f0 = 0.0, + IODE_Sub_2 = "00000000", + C_rs = 0.0, + Δn = 0.0, + M_0 = -0.23799319949775619, + C_uc = 0.0, + e = 0.0, + C_us = 0.0, + sqrt_A = 5153.700811386108, + t_0e = 266400, + fit_interval = false, + AODO = 31, + C_ic = 0.0, + Ω_0 = 2.881117323198194, + C_is = 0.0, + i_0 = 0.9599310884343369, + C_rc = 0.0, + ω = -1.651172995948652, + Ω_dot = 0.0, + IODE_Sub_3 = "00000000", + i_dot = 0.0, + α_0 = nothing, + α_1 = nothing, + α_2 = nothing, + α_3 = nothing, + β_0 = nothing, + β_1 = nothing, + β_2 = nothing, + β_3 = nothing, + ), + GNSSDecoder.GPSL1Data( + last_subframe_id = 2, + integrity_status_flag = false, + TOW = 259242, + alert_flag = false, + anti_spoof_flag = false, + trans_week = 58, + codeonl2 = 1, + ura = 2.0, + svhealth = "000000", + IODC = "0000000000", + l2pcode = false, + T_GD = 0.0, + t_0c = 266400, + a_f2 = 0.0, + a_f1 = 0.0, + a_f0 = 0.0, + IODE_Sub_2 = "00000000", + C_rs = 0.0, + Δn = 0.0, + M_0 = -0.23799319949775619, + C_uc = 0.0, + e = 0.0, + C_us = 0.0, + sqrt_A = 5153.700811386108, + t_0e = 266400, + fit_interval = false, + AODO = 31, + C_ic = 0.0, + Ω_0 = 2.881117323198194, + C_is = 0.0, + i_0 = 0.9599310884343369, + C_rc = 0.0, + ω = -1.651172995948652, + Ω_dot = 0.0, + IODE_Sub_3 = "00000000", + i_dot = 0.0, + α_0 = 4.6566129e-9, + α_1 = 1.4901161e-8, + α_2 = -5.96046e-8, + α_3 = -5.96046e-8, + β_0 = 79872.0, + β_1 = 65536.0, + β_2 = -65536.0, + β_3 = -393216.0, + ), + GNSSDecoder.GPSL1Constants( + 300, + 0x8b, + 8, + 30, + 3.1415926535898, + 7.2921151467e-5, + 2.99792458e8, + 3.986005e14, + -4.442807633e-10, + ), + GNSSDecoder.GPSL1Cache(), + 173, + 173, + false, + ), + system = gpsl1, + code_phase = 4175.34652009604, + carrier_doppler = -3554.3272562888696Hz, + carrier_phase = 1.7547749028909547, + ), + ] + # Constellation-wide model selection picks Klobuchar: the GPS decoders + # carry the broadcast α/β, no Galileo coefficients are present. + @test PositionVelocityTime.select_ionospheric_correction(states) isa + PositionVelocityTime.KlobucharParams + + # Corrections are on by default — the fix lands within 0.5 m of the ground truth. + pvt = calc_pvt(states; approximate_year = 2020) + @test norm(pvt.position - IONO_TROPO_GROUND_TRUTH) < 0.5 + + # Disabling the corrections moves the fix several metres off truth, so + # the corrections demonstrably improve the solution. + uncorrected = calc_pvt( + states; + approximate_year = 2020, + enable_ionospheric_correction = false, + enable_tropospheric_correction = false, + ) + @test norm(uncorrected.position - IONO_TROPO_GROUND_TRUTH) > 5.0 + @test norm(pvt.position - IONO_TROPO_GROUND_TRUTH) < + norm(uncorrected.position - IONO_TROPO_GROUND_TRUTH) + end + + @testset "Galileo E1B — NTCM-G (decoded coefficients)" begin + galileo_e1b = GalileoE1B() + states = [ + SatelliteState(; + decoder = GNSSDecoderState( + 5, + uint288"0xbfffffed00000047fffffe10000004bfffffeb04200f832fffe1e3c9ff809f5803ef0640", + uint288"0xc00bdee5802000000c7fffffc4000000f7fffffda0000008ffffffc200000097fffffd60", + GNSSDecoder.GalileoE1BData( + WN = 1082, + TOW = 259279, + t_0e = 259200.0, + M_0 = 1.9477946226528733, + e = 0.0, + sqrt_A = 5440.588203430176, + Ω_0 = -1.3732156716246846, + i_0 = 0.9773843813769011, + ω = 0.0, + i_dot = 0.0, + Ω_dot = 0.0, + Δn = 0.0, + C_uc = 0.0, + C_us = 0.0, + C_rc = 0.0, + C_rs = 0.0, + C_ic = 0.0, + C_is = 0.0, + t_0c = 259200.0, + a_f0 = 0.0, + a_f1 = 0.0, + a_f2 = 0.0, + IOD_nav1 = 0x0000000000000030, + IOD_nav2 = 0x0000000000000030, + IOD_nav3 = 0x0000000000000030, + IOD_nav4 = 0x0000000000000030, + num_pages_after_last_TOW = 4, + num_bits_after_valid_syncro_sequence_after_last_TOW = 1010, + signal_health_e1b = GNSSDecoder.signal_ok, + signal_health_e5b = GNSSDecoder.signal_ok, + data_validity_status_e1b = GNSSDecoder.navigation_data_valid, + data_validity_status_e5b = GNSSDecoder.navigation_data_valid, + broadcast_group_delay_e1_e5a = 0.0, + broadcast_group_delay_e1_e5b = 0.0, + a_i0 = 100.0, + a_i1 = 0.0, + a_i2 = 0.0, + ), + GNSSDecoder.GalileoE1BData( + WN = 1082, + TOW = 259279, + t_0e = 259200.0, + M_0 = 1.9477946226528733, + e = 0.0, + sqrt_A = 5440.588203430176, + Ω_0 = -1.3732156716246846, + i_0 = 0.9773843813769011, + ω = 0.0, + i_dot = 0.0, + Ω_dot = 0.0, + Δn = 0.0, + C_uc = 0.0, + C_us = 0.0, + C_rc = 0.0, + C_rs = 0.0, + C_ic = 0.0, + C_is = 0.0, + t_0c = 259200.0, + a_f0 = 0.0, + a_f1 = 0.0, + a_f2 = 0.0, + IOD_nav1 = 0x0000000000000030, + IOD_nav2 = 0x0000000000000030, + IOD_nav3 = 0x0000000000000030, + IOD_nav4 = 0x0000000000000030, + num_pages_after_last_TOW = 4, + num_bits_after_valid_syncro_sequence_after_last_TOW = 1010, + signal_health_e1b = GNSSDecoder.signal_ok, + signal_health_e5b = GNSSDecoder.signal_ok, + data_validity_status_e1b = GNSSDecoder.navigation_data_valid, + data_validity_status_e5b = GNSSDecoder.navigation_data_valid, + broadcast_group_delay_e1_e5a = 0.0, + broadcast_group_delay_e1_e5b = 0.0, + a_i0 = 100.0, + a_i1 = 0.0, + a_i2 = 0.0, + ), + GNSSDecoder.GalileoE1BConstants( + 250, + 0x0160, + 10, + 3.1415926535898, + 7.2921151467e-5, + 2.99792458e8, + 3.986004418e14, + -4.442807309e-10, + ), + GNSSDecoder.GalileoE1BCache(), + 141, + 1391, + false, + ), + system = galileo_e1b, + code_phase = 1134.9227691611316, + carrier_doppler = -1048.4103530464388Hz, + carrier_phase = -1.1559529083133402, + ), + SatelliteState(; + decoder = GNSSDecoderState( + 6, + uint288"0xf7fffffda0000008ffffffc200000097fffffd608401f065fffc3c793ff013eb007de0c8", + uint288"0xc00be565802000000c7fffffc4000000f7fffffda0000008ffffffc200000097fffffd60", + GNSSDecoder.GalileoE1BData( + WN = 1082, + TOW = 259279, + t_0e = 259200.0, + M_0 = 2.7331927860503233, + e = 0.0, + sqrt_A = 5440.588203430176, + Ω_0 = -1.3732156716246846, + i_0 = 0.9773843813769011, + ω = 0.0, + i_dot = 0.0, + Ω_dot = 0.0, + Δn = 0.0, + C_uc = 0.0, + C_us = 0.0, + C_rc = 0.0, + C_rs = 0.0, + C_ic = 0.0, + C_is = 0.0, + t_0c = 259200.0, + a_f0 = 0.0, + a_f1 = 0.0, + a_f2 = 0.0, + IOD_nav1 = 0x0000000000000030, + IOD_nav2 = 0x0000000000000030, + IOD_nav3 = 0x0000000000000030, + IOD_nav4 = 0x0000000000000030, + num_pages_after_last_TOW = 4, + num_bits_after_valid_syncro_sequence_after_last_TOW = 1010, + signal_health_e1b = GNSSDecoder.signal_ok, + signal_health_e5b = GNSSDecoder.signal_ok, + data_validity_status_e1b = GNSSDecoder.navigation_data_valid, + data_validity_status_e5b = GNSSDecoder.navigation_data_valid, + broadcast_group_delay_e1_e5a = 0.0, + broadcast_group_delay_e1_e5b = 0.0, + a_i0 = 100.0, + a_i1 = 0.0, + a_i2 = 0.0, + ), + GNSSDecoder.GalileoE1BData( + WN = 1082, + TOW = 259279, + t_0e = 259200.0, + M_0 = 2.7331927860503233, + e = 0.0, + sqrt_A = 5440.588203430176, + Ω_0 = -1.3732156716246846, + i_0 = 0.9773843813769011, + ω = 0.0, + i_dot = 0.0, + Ω_dot = 0.0, + Δn = 0.0, + C_uc = 0.0, + C_us = 0.0, + C_rc = 0.0, + C_rs = 0.0, + C_ic = 0.0, + C_is = 0.0, + t_0c = 259200.0, + a_f0 = 0.0, + a_f1 = 0.0, + a_f2 = 0.0, + IOD_nav1 = 0x0000000000000030, + IOD_nav2 = 0x0000000000000030, + IOD_nav3 = 0x0000000000000030, + IOD_nav4 = 0x0000000000000030, + num_pages_after_last_TOW = 4, + num_bits_after_valid_syncro_sequence_after_last_TOW = 1010, + signal_health_e1b = GNSSDecoder.signal_ok, + signal_health_e5b = GNSSDecoder.signal_ok, + data_validity_status_e1b = GNSSDecoder.navigation_data_valid, + data_validity_status_e5b = GNSSDecoder.navigation_data_valid, + broadcast_group_delay_e1_e5a = 0.0, + broadcast_group_delay_e1_e5b = 0.0, + a_i0 = 100.0, + a_i1 = 0.0, + a_i2 = 0.0, + ), + GNSSDecoder.GalileoE1BConstants( + 250, + 0x0160, + 10, + 3.1415926535898, + 7.2921151467e-5, + 2.99792458e8, + 3.986004418e14, + -4.442807309e-10, + ), + GNSSDecoder.GalileoE1BCache(), + 138, + 1388, + false, + ), + system = galileo_e1b, + code_phase = 552.8323205662532, + carrier_doppler = -2401.7199003375345Hz, + carrier_phase = -0.5865229251135221, + ), + SatelliteState(; + decoder = GNSSDecoderState( + 20, + uint288"0x080000025ffffff70000003dffffff680000029f7bfe0f9a0003c386c00fec14ff821f37", + uint288"0xc00bc465802000000c7fffffc4000000f7fffffda0000008ffffffc200000097fffffd60", + GNSSDecoder.GalileoE1BData( + WN = 1082, + TOW = 259279, + t_0e = 259200.0, + M_0 = 1.685995235829002, + e = 0.0, + sqrt_A = 5440.588203430176, + Ω_0 = 2.815574532186437, + i_0 = 0.9773843813769011, + ω = 0.0, + i_dot = 0.0, + Ω_dot = 0.0, + Δn = 0.0, + C_uc = 0.0, + C_us = 0.0, + C_rc = 0.0, + C_rs = 0.0, + C_ic = 0.0, + C_is = 0.0, + t_0c = 259200.0, + a_f0 = 0.0, + a_f1 = 0.0, + a_f2 = 0.0, + IOD_nav1 = 0x0000000000000030, + IOD_nav2 = 0x0000000000000030, + IOD_nav3 = 0x0000000000000030, + IOD_nav4 = 0x0000000000000030, + num_pages_after_last_TOW = 4, + num_bits_after_valid_syncro_sequence_after_last_TOW = 1010, + signal_health_e1b = GNSSDecoder.signal_ok, + signal_health_e5b = GNSSDecoder.signal_ok, + data_validity_status_e1b = GNSSDecoder.navigation_data_valid, + data_validity_status_e5b = GNSSDecoder.navigation_data_valid, + broadcast_group_delay_e1_e5a = 0.0, + broadcast_group_delay_e1_e5b = 0.0, + a_i0 = 100.0, + a_i1 = 0.0, + a_i2 = 0.0, + ), + GNSSDecoder.GalileoE1BData( + WN = 1082, + TOW = 259279, + t_0e = 259200.0, + M_0 = 1.685995235829002, + e = 0.0, + sqrt_A = 5440.588203430176, + Ω_0 = 2.815574532186437, + i_0 = 0.9773843813769011, + ω = 0.0, + i_dot = 0.0, + Ω_dot = 0.0, + Δn = 0.0, + C_uc = 0.0, + C_us = 0.0, + C_rc = 0.0, + C_rs = 0.0, + C_ic = 0.0, + C_is = 0.0, + t_0c = 259200.0, + a_f0 = 0.0, + a_f1 = 0.0, + a_f2 = 0.0, + IOD_nav1 = 0x0000000000000030, + IOD_nav2 = 0x0000000000000030, + IOD_nav3 = 0x0000000000000030, + IOD_nav4 = 0x0000000000000030, + num_pages_after_last_TOW = 4, + num_bits_after_valid_syncro_sequence_after_last_TOW = 1010, + signal_health_e1b = GNSSDecoder.signal_ok, + signal_health_e5b = GNSSDecoder.signal_ok, + data_validity_status_e1b = GNSSDecoder.navigation_data_valid, + data_validity_status_e5b = GNSSDecoder.navigation_data_valid, + broadcast_group_delay_e1_e5a = 0.0, + broadcast_group_delay_e1_e5b = 0.0, + a_i0 = 100.0, + a_i1 = 0.0, + a_i2 = 0.0, + ), + GNSSDecoder.GalileoE1BConstants( + 250, + 0x0160, + 10, + 3.1415926535898, + 7.2921151467e-5, + 2.99792458e8, + 3.986004418e14, + -4.442807309e-10, + ), + GNSSDecoder.GalileoE1BCache(), + 138, + 1388, + true, + ), + system = galileo_e1b, + code_phase = 1998.2750592825137, + carrier_doppler = 1342.86940202154Hz, + carrier_phase = -1.9500644306669892, + ), + SatelliteState(; + decoder = GNSSDecoderState( + 22, + uint288"0x080000025ffffff70000003dffffff680000029f7bfe0f9a0003c386c00fec14ff821f37", + uint288"0xc00bf825802000000c7fffffc4000000f7fffffda0000008ffffffc200000097fffffd60", + GNSSDecoder.GalileoE1BData( + WN = 1082, + TOW = 259279, + t_0e = 259200.0, + M_0 = -3.026393744555698, + e = 0.0, + sqrt_A = 5440.588203430176, + Ω_0 = 2.815574532186437, + i_0 = 0.9773843813769011, + ω = 0.0, + i_dot = 0.0, + Ω_dot = 0.0, + Δn = 0.0, + C_uc = 0.0, + C_us = 0.0, + C_rc = 0.0, + C_rs = 0.0, + C_ic = 0.0, + C_is = 0.0, + t_0c = 259200.0, + a_f0 = 0.0, + a_f1 = 0.0, + a_f2 = 0.0, + IOD_nav1 = 0x0000000000000030, + IOD_nav2 = 0x0000000000000030, + IOD_nav3 = 0x0000000000000030, + IOD_nav4 = 0x0000000000000030, + num_pages_after_last_TOW = 4, + num_bits_after_valid_syncro_sequence_after_last_TOW = 1010, + signal_health_e1b = GNSSDecoder.signal_ok, + signal_health_e5b = GNSSDecoder.signal_ok, + data_validity_status_e1b = GNSSDecoder.navigation_data_valid, + data_validity_status_e5b = GNSSDecoder.navigation_data_valid, + broadcast_group_delay_e1_e5a = 0.0, + broadcast_group_delay_e1_e5b = 0.0, + a_i0 = 100.0, + a_i1 = 0.0, + a_i2 = 0.0, + ), + GNSSDecoder.GalileoE1BData( + WN = 1082, + TOW = 259279, + t_0e = 259200.0, + M_0 = -3.026393744555698, + e = 0.0, + sqrt_A = 5440.588203430176, + Ω_0 = 2.815574532186437, + i_0 = 0.9773843813769011, + ω = 0.0, + i_dot = 0.0, + Ω_dot = 0.0, + Δn = 0.0, + C_uc = 0.0, + C_us = 0.0, + C_rc = 0.0, + C_rs = 0.0, + C_ic = 0.0, + C_is = 0.0, + t_0c = 259200.0, + a_f0 = 0.0, + a_f1 = 0.0, + a_f2 = 0.0, + IOD_nav1 = 0x0000000000000030, + IOD_nav2 = 0x0000000000000030, + IOD_nav3 = 0x0000000000000030, + IOD_nav4 = 0x0000000000000030, + num_pages_after_last_TOW = 4, + num_bits_after_valid_syncro_sequence_after_last_TOW = 1010, + signal_health_e1b = GNSSDecoder.signal_ok, + signal_health_e5b = GNSSDecoder.signal_ok, + data_validity_status_e1b = GNSSDecoder.navigation_data_valid, + data_validity_status_e5b = GNSSDecoder.navigation_data_valid, + broadcast_group_delay_e1_e5a = 0.0, + broadcast_group_delay_e1_e5b = 0.0, + a_i0 = 100.0, + a_i1 = 0.0, + a_i2 = 0.0, + ), + GNSSDecoder.GalileoE1BConstants( + 250, + 0x0160, + 10, + 3.1415926535898, + 7.2921151467e-5, + 2.99792458e8, + 3.986004418e14, + -4.442807309e-10, + ), + GNSSDecoder.GalileoE1BCache(), + 138, + 1388, + true, + ), + system = galileo_e1b, + code_phase = 2721.0058329387184, + carrier_doppler = -3042.4604297211754Hz, + carrier_phase = -0.5604621503191517, + ), + ] + # NTCM-G is selected from the Effective Ionisation Level coefficients + # decoded from the E1B navigation message. + @test PositionVelocityTime.select_ionospheric_correction(states) isa + PositionVelocityTime.NTCMGParams + + pvt = calc_pvt(states; approximate_year = 2020) + @test norm(pvt.position - IONO_TROPO_GROUND_TRUTH) < 0.5 + + uncorrected = calc_pvt( + states; + approximate_year = 2020, + enable_ionospheric_correction = false, + enable_tropospheric_correction = false, + ) + @test norm(uncorrected.position - IONO_TROPO_GROUND_TRUTH) > 5.0 + @test norm(pvt.position - IONO_TROPO_GROUND_TRUTH) < + norm(uncorrected.position - IONO_TROPO_GROUND_TRUTH) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 88b89e4..fb329a2 100755 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,5 +5,8 @@ using Unitful: Hz include("aqua.jl") include("sat_time.jl") include("pvt.jl") +include("pvt_iono_tropo.jl") include("get_week.jl") -include("tracking_ext.jl") \ No newline at end of file +include("tracking_ext.jl") +include("ionosphere.jl") +include("troposphere.jl") diff --git a/test/troposphere.jl b/test/troposphere.jl new file mode 100644 index 0000000..edab603 --- /dev/null +++ b/test/troposphere.jl @@ -0,0 +1,73 @@ +# Independent reimplementation of the Saastamoinen model (RTKLIB `tropmodel`), +# written with sin(elevation) instead of cos(z) and explicit intermediate terms, +# to cross-check the package's `saastamoinen_delay`. +function saastamoinen_reference(lat, height, elevation, humidity) + (height < -100.0 || height > 1.0e4 || elevation <= 0.0) && return 0.0 + h = max(height, 0.0) + P = 1013.25 * (1 - 2.2557e-5 * h)^5.2568 + T = 288.16 - 6.5e-3 * h + e = 6.108 * humidity * exp((17.15 * T - 4684.0) / (T - 38.45)) + mapping = 1 / sin(elevation) # = 1/cos(z), z = π/2 - el + zhd = 0.0022768 * P / (1 - 0.00266 * cos(2lat) - 0.00028 * h / 1000) + zwd = 0.002277 * (1255.0 / T + 0.05) * e + return (zhd + zwd) * mapping +end + +@testset "Saastamoinen tropospheric model" begin + @testset "matches independent reference" begin + for lat in deg2rad.((0.0, 48.0, -67.0)), + h in (0.0, 550.0, 3000.0, 9000.0), + el in deg2rad.((5.0, 15.0, 45.0, 90.0)), + humi in (0.0, 0.5, 0.7, 1.0) + + got = PositionVelocityTime.saastamoinen_delay(lat, h, el, humi) + ref = saastamoinen_reference(lat, h, el, humi) + @test got ≈ ref rtol = 1e-12 + end + end + + @testset "physical behaviour" begin + lat, humi = deg2rad(48.0), 0.7 + zenith = PositionVelocityTime.saastamoinen_delay(lat, 0.0, deg2rad(90.0), humi) + # Zenith total delay at sea level is a couple of metres + @test 2.2 < zenith < 2.6 + # Obliquity: lower elevation → larger slant delay (≈ 1/sin(el)) + d30 = PositionVelocityTime.saastamoinen_delay(lat, 0.0, deg2rad(30.0), humi) + d10 = PositionVelocityTime.saastamoinen_delay(lat, 0.0, deg2rad(10.0), humi) + @test d30 > zenith + @test d10 > d30 + @test d30 ≈ zenith / sin(deg2rad(30.0)) rtol = 1e-12 + # Delay decreases with user height (thinner atmosphere above) + @test PositionVelocityTime.saastamoinen_delay(lat, 3000.0, deg2rad(90.0), humi) < + zenith + # Wet component grows with humidity; dry-only is the floor + @test PositionVelocityTime.saastamoinen_delay(lat, 0.0, deg2rad(90.0), 1.0) > + PositionVelocityTime.saastamoinen_delay(lat, 0.0, deg2rad(90.0), 0.0) + end + + @testset "edge cases return zero" begin + lat, humi = deg2rad(48.0), 0.7 + # Satellite at/below the horizon, or user height outside the model's valid + # range (−100 m … 10 km), returns exactly zero. The four cases below are, in + # order: at horizon, below horizon, too high, too low. + @test PositionVelocityTime.saastamoinen_delay(lat, 0.0, 0.0, humi) == 0.0 + @test PositionVelocityTime.saastamoinen_delay(lat, 0.0, deg2rad(-5.0), humi) == 0.0 + @test PositionVelocityTime.saastamoinen_delay(lat, 1.5e4, deg2rad(45.0), humi) == 0.0 + @test PositionVelocityTime.saastamoinen_delay(lat, -200.0, deg2rad(45.0), humi) == 0.0 + end + + @testset "tropospheric_delay from line-of-sight geometry" begin + user = ECEFfromLLA(wgs84)(LLA(48.0, 11.0, 550.0)) + lla = LLAfromECEF(wgs84)(user) + # Satellite roughly overhead → near-zenith, small slant delay (~2 m) + sat_up = ECEFfromLLA(wgs84)(LLA(48.0, 11.0, 2.0e7)) + el_up, _ = PositionVelocityTime._elevation_azimuth(ENUfromECEF(user, wgs84), sat_up) + d_up = PositionVelocityTime.tropospheric_delay(el_up, lla) + @test 2.0 < d_up < 2.6 + # A low-elevation satellite has a larger slant delay than the near-zenith one + sat_low = ECEF(user[1] + 2.0e7, user[2] + 2.0e6, user[3] + 2.0e6) + el_low, _ = PositionVelocityTime._elevation_azimuth(ENUfromECEF(user, wgs84), sat_low) + @test PositionVelocityTime.tropospheric_delay(el_low, lla) > d_up + # Frequency independence is implicit: tropospheric_delay takes no system/freq. + end +end