diff --git a/docs/src/api.md b/docs/src/api.md index de5dd2c..5051cd1 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -6,6 +6,7 @@ SatelliteState PVTSolution SatInfo +InterFrequencyBias ``` ## PVT Computation @@ -62,3 +63,10 @@ receiver-clock grouping (one bias per time system); `GNSSSignals.get_band_id` (e `GNSSSignals.get_signal_id` (e.g. `:GPSL1CA`) is the per-signal identity used in the `sats` key of [`PVTSolution`](@ref). This package forwards each of the latter two to a [`SatelliteState`](@ref). + +The receiver inter-frequency biases and their reference bands (reported per +[`InterFrequencyBias`](@ref)) are laid out from the constellation × band coverage graph: + +```@docs +PositionVelocityTime.band_ifb_layout +``` diff --git a/src/PositionVelocityTime.jl b/src/PositionVelocityTime.jl index 66b7756..c0f1ef9 100755 --- a/src/PositionVelocityTime.jl +++ b/src/PositionVelocityTime.jl @@ -11,12 +11,13 @@ using CoordinateTransformations, Unitful, Dates -using Unitful: s, Hz +using Unitful: s, Hz, m, °, ustrip using Dictionaries: Dictionary export calc_pvt, PVTSolution, SatInfo, + InterFrequencyBias, SatelliteState, get_LLA, get_sat_info, @@ -129,14 +130,34 @@ satellite used in the fix). # Fields - `position::ECEF`: Satellite ECEF position at transmit time (metres). - `time::Float64`: Satellite transmit time (system time of week, seconds). -- `residual::Float64`: Post-fit least-squares pseudorange residual (metres) — the +- `residual::typeof(1.0m)`: Post-fit least-squares pseudorange residual (metres) — the modeled minus the (atmosphere-corrected) measured pseudorange. A per-satellite fit-quality / outlier indicator. """ struct SatInfo position::ECEF time::Float64 - residual::Float64 + residual::typeof(1.0m) +end + +""" + InterFrequencyBias + +A receiver inter-frequency bias attached to a [`PVTSolution`](@ref): the differential +hardware delay of one frequency band's RF chain, together with the band it is measured +against. + +# Fields +- `value::typeof(1.0m)`: The bias (metres), relative to `reference` — how much longer + this band's receiver chain delay is than the reference band's. +- `reference::Symbol`: The reference band (e.g. `:L1`; see `get_band_id`) whose delay is + folded into the clock biases and against which `value` is measured. Chosen per + coverage component (see [`band_ifb_layout`](@ref)), so different biases in a + disconnected solution can in principle carry different references. +""" +struct InterFrequencyBias + value::typeof(1.0m) + reference::Symbol end """ @@ -147,7 +168,14 @@ Complete Position, Velocity, and Time solution from GNSS measurements. # Fields - `position::ECEF`: User position in ECEF coordinates (meters) - `velocity::ECEF`: User velocity in ECEF coordinates (m/s) -- `time_correction::Float64`: Estimated receiver clock bias of the reference GNSS +- `course_over_ground::typeof(1.0°)`: Horizontal direction of travel (degrees), the + azimuth of the velocity vector in the local East-North-Up frame at `position` — + measured clockwise from true North and wrapped to `[0, 360)°`, following the GNSS + course-over-ground (COG) convention. Derived from the velocity alone, so this is the + direction of motion, not vehicle heading (which a single-antenna receiver cannot + observe). The vertical (Up) velocity component is ignored. `0°` when the horizontal + velocity is zero (stationary or purely vertical), where course is undefined. +- `time_correction::typeof(1.0m)`: Estimated receiver clock bias of the reference GNSS (meters). For a multi-GNSS solution this is the bias of `reference_system`; other systems' biases are `time_correction + inter_system_biases[system]`. - `time::Union{TAIEpoch{Float64}, Nothing}`: Estimated time as a TAI epoch @@ -162,7 +190,7 @@ Complete Position, Velocity, and Time solution from GNSS measurements. see `reference_system`. - `reference_system::Union{GNSSSignals.TimeSystem, Nothing}`: GNSS time system (e.g. `GNSSSignals.GPST()`, `GST()`) that `time` and `time_correction` are referenced to. -- `inter_system_biases::Dict{GNSSSignals.TimeSystem, Float64}`: For each GNSS time system +- `inter_system_biases::Dict{GNSSSignals.TimeSystem, typeof(1.0m)}`: For each GNSS time system other than `reference_system`, the offset of that system's time scale relative to the reference system's (meters) — the inter-system bias. This is a **system / space-segment** effect, not a receiver one (this receiver has no inter-system hardware bias): it is the GNSS @@ -172,26 +200,33 @@ Complete Position, Velocity, and Time solution from GNSS measurements. GGTO-aided collapse), as that broadcast value may be erroneous. Reference-independent (the difference of two entries is the offset between those two systems); empty for a single-system solution. -- `inter_frequency_biases::Dict{Symbol, Float64}`: For each frequency band other than - the reference band, the receiver inter-frequency bias relative to it (meters) — the +- `inter_frequency_biases::Dict{Symbol, InterFrequencyBias}`: For each frequency band + other than the reference band, the receiver inter-frequency bias relative to it — the differential hardware delay of that band's RF chain, estimated as an extra unknown when satellites are processed on more than one band. The key is the band - (e.g. `:L1`, `:L5`; see `get_band_id`) and the value is shared across - all constellations on that band. Empty for a single-band solution. A solution then - needs `n ≥ 3 + M + B` satellites for `M` time systems and `B` extra bands. + (e.g. `:L1`, `:L5`; see `get_band_id`) and the value is shared across all + constellations on that band. Each [`InterFrequencyBias`](@ref) carries both the bias + (in metres) and the reference band it is measured against — the anchor whose delay is + folded into the clock biases — so each reported bias's anchor is explicit rather than + implicit. The reference is chosen per coverage component (see + [`band_ifb_layout`](@ref)), so a disconnected constellation can in principle yield + several references. Empty for a single-band solution. A solution then needs + `n ≥ 3 + M + B` satellites for `M` time systems and `B` extra bands. """ @kwdef struct PVTSolution position::ECEF = ECEF(0, 0, 0) velocity::ECEF = ECEF(0, 0, 0) - time_correction::Float64 = 0 + course_over_ground::typeof(1.0°) = 0.0° + time_correction::typeof(1.0m) = 0.0m time::Union{TAIEpoch{Float64},Nothing} = nothing relative_clock_drift::Float64 = 0 dop::Union{DOP,Nothing} = nothing sats::Dictionary{Tuple{Symbol,Int},SatInfo} = Dictionary{Tuple{Symbol,Int},SatInfo}() reference_system::Union{GNSSSignals.TimeSystem,Nothing} = nothing - inter_system_biases::Dict{GNSSSignals.TimeSystem,Float64} = - Dict{GNSSSignals.TimeSystem,Float64}() - inter_frequency_biases::Dict{Symbol,Float64} = Dict{Symbol,Float64}() + inter_system_biases::Dict{GNSSSignals.TimeSystem,typeof(1.0m)} = + Dict{GNSSSignals.TimeSystem,typeof(1.0m)}() + inter_frequency_biases::Dict{Symbol,InterFrequencyBias} = + Dict{Symbol,InterFrequencyBias}() end """ @@ -233,7 +268,26 @@ get_sat_enu(user_pos_ecef::ECEF, sat_pos_ecef::ECEF) = get_sat_enu(ENUfromECEF(user_pos_ecef, wgs84), sat_pos_ecef) """ - band_ifb_layout(systems, bands) -> (ifb_indices, extra_bands, num_components) + calc_course_over_ground(position::ECEF, velocity::ECEF) -> typeof(1.0°) + +Course over ground (degrees): the azimuth of `velocity` in the local East-North-Up +frame at `position`, measured clockwise from true North and wrapped to `[0, 360)°`, +following the GNSS COG convention. This is the direction of motion derived from the +velocity, not vehicle heading. The vertical (Up) component is ignored. Returns `0°` +when the horizontal velocity is zero (stationary or purely vertical), where course is +undefined. + +`ENUfromECEF` is affine with `position` as its origin (which it maps to zero), so +applying it to `position + velocity` rotates the velocity vector into the ENU frame +without introducing any translation. +""" +function calc_course_over_ground(position::ECEF, velocity::ECEF) + enu_velocity = ENUfromECEF(position, wgs84)(ECEF(position + velocity)) + rad2deg(mod2pi(atan(enu_velocity[1], enu_velocity[2]))) * ° +end + +""" + band_ifb_layout(systems, bands) -> (ifb_indices, extra_bands, reference_bands, num_components) Lay out the receiver inter-frequency biases from the (constellation × band) coverage graph. Two bands share a *coverage component* iff some constellation is tracked on @@ -247,8 +301,9 @@ matrix is always full rank. `ifb_indices[j]` is satellite `j`'s IFB column (1…`length(extra_bands)`), or `0` for a per-component reference band; `extra_bands[i]` is the band of IFB column `i` (ordered -deterministically); `num_components` is the number of coverage components (`1` ⇔ the -graph is connected). +deterministically) and `reference_bands[i]` is the reference band of that column's +coverage component (the anchor its IFB is measured against); `num_components` is the +number of coverage components (`1` ⇔ the graph is connected). """ function band_ifb_layout(systems, bands) unique_bands = unique(bands) @@ -278,14 +333,15 @@ function band_ifb_layout(systems, bands) end end extra_bands = filter(b -> reference_of[root(b)] != b, unique_bands) + reference_bands = [reference_of[root(b)] for b in extra_bands] band_column = Dict(b => i for (i, b) in enumerate(extra_bands)) ifb_indices = [get(band_column, b, 0) for b in bands] - return ifb_indices, extra_bands, length(reference_of) + return ifb_indices, extra_bands, reference_bands, length(reference_of) end """ decide_bias_layout(states, systems, bands, times) - -> (bias_columns::BiasColumns, extra_bands, ggto_offsets) + -> (bias_columns::BiasColumns, extra_bands, reference_bands, ggto_offsets) Decide the full least-squares bias layout — one clock column per GNSS time system plus the per-band inter-frequency-bias columns from [`band_ifb_layout`](@ref) — and the @@ -323,15 +379,17 @@ function decide_bias_layout(states, systems, bands, times) unique_effective = unique(effective_systems) index = Dict(sys => i for (i, sys) in enumerate(unique_effective)) clock_bias_indices = [index[sys] for sys in effective_systems] - ifb_indices, extra_bands, num_components = band_ifb_layout(effective_systems, bands) + ifb_indices, extra_bands, reference_bands, num_components = + band_ifb_layout(effective_systems, bands) (; clock_bias_indices, num_clock_biases = length(unique_effective), ifb_indices, - extra_bands, num_components) + extra_bands, reference_bands, num_components) end as_result(layout, ggto_offsets) = ( BiasColumns(layout.clock_bias_indices, layout.num_clock_biases, layout.ifb_indices, length(layout.extra_bands)), layout.extra_bands, + layout.reference_bands, ggto_offsets, ) @@ -514,7 +572,7 @@ function calc_pvt( bias_layout = decide_bias_layout(healthy_states, systems, bands, times) isnothing(bias_layout) && return prev_pvt - bias_columns, extra_bands, ggto_offsets = bias_layout + bias_columns, extra_bands, reference_bands, ggto_offsets = bias_layout (; clock_bias_indices, num_clock_biases) = bias_columns # Propagating the ephemerides. @@ -551,14 +609,22 @@ function calc_pvt( # absolute bias from the reference bias plus its stored inter-system bias. prev_abs_bias(sys) = sys == prev_pvt.reference_system ? prev_pvt.time_correction : - prev_pvt.time_correction + get(prev_pvt.inter_system_biases, sys, 0.0) + prev_pvt.time_correction + get(prev_pvt.inter_system_biases, sys, 0.0m) prev_ξ = zeros(num_lsq_params(bias_columns)) prev_ξ[1], prev_ξ[2], prev_ξ[3] = prev_pvt.position for j in 1:num_sats - prev_ξ[3+clock_bias_indices[j]] = prev_abs_bias(systems[j]) + prev_ξ[3+clock_bias_indices[j]] = ustrip(m, prev_abs_bias(systems[j])) end + # Only warm-start an IFB column from the previous solution when that band's bias was + # measured against the same reference band; a reference change (the anchor differs + # across epochs) makes the stored value refer to a different quantity, so seed at 0 + # instead. The IFB enters the design linearly, so a stale seed only costs iterations, + # but this keeps the starting point meaningful. for (i, band) in enumerate(extra_bands) - prev_ξ[3+num_clock_biases+i] = get(prev_pvt.inter_frequency_biases, band, 0.0) + prev_ifb = get(prev_pvt.inter_frequency_biases, band, nothing) + prev_ξ[3+num_clock_biases+i] = + !isnothing(prev_ifb) && prev_ifb.reference == reference_bands[i] ? + ustrip(m, prev_ifb.value) : 0.0 end # Atmospheric corrections, summed per satellite and subtracted from the @@ -611,6 +677,7 @@ function calc_pvt( user_velocity_and_clock_drift[3], ) relative_clock_drift = user_velocity_and_clock_drift[4] / SPEEDOFLIGHT + course_over_ground = calc_course_over_ground(position, velocity) time_correction = ξ[3+primary_clock_index] # The estimated time correction is negative # See https://github.com/JuliaGNSS/PositionVelocityTime.jl/issues/8 @@ -625,7 +692,7 @@ function calc_pvt( corrected_reference_time - floor(Int, corrected_reference_time), ) - sat_infos = SatInfo.(sat_positions, times, residuals) + sat_infos = SatInfo.(sat_positions, times, residuals .* m) dop = calc_DOP(H, position, primary_clock_index) if dop.GDOP < 0 @@ -637,19 +704,20 @@ function calc_pvt( # GGTO-collapsed system this is the broadcast offset −c·Δt_systems, read from the # system's first satellite (the per-satellite offsets differ only by the GGTO # drift A_1G over the transmit-time spread — sub-millimetre). - inter_system_biases = Dict{GNSSSignals.TimeSystem,Float64}() + inter_system_biases = Dict{GNSSSignals.TimeSystem,typeof(1.0m)}() for sys in unique_systems sys == primary_system && continue j = findfirst(==(sys), systems) inter_system_biases[sys] = - ξ[3+clock_bias_indices[j]] + ggto_offsets[j] - time_correction + (ξ[3+clock_bias_indices[j]] + ggto_offsets[j] - time_correction) * m end # Receiver inter-frequency biases relative to the reference band, in meters # (the reference band is omitted; its bias is folded into the clock biases). - inter_frequency_biases = Dict{Symbol,Float64}() + inter_frequency_biases = Dict{Symbol,InterFrequencyBias}() for (i, band) in enumerate(extra_bands) - inter_frequency_biases[band] = ξ[3+num_clock_biases+i] + inter_frequency_biases[band] = + InterFrequencyBias(ξ[3+num_clock_biases+i] * m, reference_bands[i]) end # Per-satellite `sats` key: (signal id, PRN) — signal-level (not time system), so a @@ -661,7 +729,8 @@ function calc_pvt( PVTSolution( position, velocity, - time_correction, + course_over_ground, + time_correction * m, time, relative_clock_drift, dop, diff --git a/src/sat_time.jl b/src/sat_time.jl index c2d4078..a08e0f1 100644 --- a/src/sat_time.jl +++ b/src/sat_time.jl @@ -24,7 +24,15 @@ function calc_uncorrected_time(state::SatelliteState) t_bits = state.decoder.num_bits_after_valid_syncro_sequence / get_data_frequency(state.decoder) * Hz - t_code_phase = state.code_phase / get_code_frequency(system) * Hz + # The code phase is shared across the same-band signals of a jointly tracked + # satellite and can run up to a whole data symbol of chips. That whole-symbol part + # is already counted by `t_bits` (which uses the decoder's data-symbol rate), so + # reduce the code phase to the residual within one data symbol to avoid + # double-counting it. A no-op for a single-signal sat, whose code phase is already + # below one data symbol of chips. + chips_per_symbol = + ustrip(Hz, get_code_frequency(system)) / ustrip(Hz, get_data_frequency(state.decoder)) + t_code_phase = mod(state.code_phase, chips_per_symbol) / get_code_frequency(system) * Hz t_carrier_phase = state.carrier_phase / get_center_frequency(system) * Hz t_tow + t_bits + t_code_phase + t_carrier_phase diff --git a/test/inter_frequency_bias.jl b/test/inter_frequency_bias.jl index d20b687..3cbb2cb 100644 --- a/test/inter_frequency_bias.jl +++ b/test/inter_frequency_bias.jl @@ -6,6 +6,21 @@ kw = (; approximate_year = 2021, enable_ionospheric_correction = false, enable_tropospheric_correction = false) + # Decompose an uncorrected transmit time `target` (s) into a realistic measurement on + # `system` with time-of-week `tow`: whole data symbols go into the bit count and only + # the residual (< one data symbol) into the code phase — matching + # `calc_uncorrected_time`, which reduces the shared code phase modulo one data symbol. + # Packing the whole sub-TOW interval into `code_phase` (as a naive copy would) is + # reduced away by that mod, so the copy must spread it across the bit count. + function bits_and_code_phase(system, tow, target) + datafreq = Float64(GNSSSignals.get_data_frequency(system) / Hz) + codefreq = Float64(GNSSSignals.get_code_frequency(system) / Hz) + elapsed = target - tow + num_bits = floor(Int, elapsed * datafreq) + code_phase = (elapsed - num_bits / datafreq) * codefreq + (num_bits, code_phase) + end + # Build an L5 (E5a) measurement for a Galileo E1B satellite, transmit-time-consistent # with it (same ephemeris/clock, BGD matched, observables solved to the same epoch). # `ifb_shift_s` injects a uniform receiver L5 delay (seconds); `ggto` (A_0G, seconds) @@ -21,11 +36,11 @@ broadcast_group_delay_e1_e5a = d.broadcast_group_delay_e1_e5b, signal_health_e5a = GNSSDecoder.signal_ok, data_validity_status_e5a = GNSSDecoder.navigation_data_valid, ggto_fields...) - dec = GNSSDecoder.GNSSDecoderState(GNSSDecoder.GalileoE5aDecoderState(state.decoder.prn); - data = e5a_data, raw_data = e5a_data, num_bits_after_valid_syncro_sequence = 0) target = PositionVelocityTime.calc_uncorrected_time(state) + ifb_shift_s - codefreq = Float64(GNSSSignals.get_code_frequency(GalileoE5aI()) / Hz) - code_phase = (target - PositionVelocityTime.get_tow(dec)) * codefreq + num_bits, code_phase = bits_and_code_phase(GalileoE5aI(), e5a_data.TOW, target) + dec = GNSSDecoder.GNSSDecoderState(GNSSDecoder.GalileoE5aDecoderState(state.decoder.prn); + data = e5a_data, raw_data = e5a_data, + num_bits_after_valid_syncro_sequence = num_bits) SatelliteState(; decoder = dec, system = GalileoE5aI(), code_phase = code_phase, carrier_doppler = 0.0Hz, carrier_phase = 0.0) end @@ -46,11 +61,11 @@ C_rc = d.C_rc, C_rs = d.C_rs, C_ic = d.C_ic, C_is = d.C_is, t_0c = d.t_0c, a_f0 = d.a_f0, a_f1 = d.a_f1, a_f2 = d.a_f2, T_GD = d.T_GD, ISC_L2C = 0.0, l2_health = false) - dec = GNSSDecoder.GNSSDecoderState(GNSSDecoder.GPSL2CMDecoderState(state.decoder.prn); - data = l2c_data, raw_data = l2c_data, num_bits_after_valid_syncro_sequence = 0) target = PositionVelocityTime.calc_uncorrected_time(state) + ifb_shift_s - codefreq = Float64(GNSSSignals.get_code_frequency(GPSL2CM()) / Hz) - code_phase = (target - PositionVelocityTime.get_tow(dec)) * codefreq + num_bits, code_phase = bits_and_code_phase(GPSL2CM(), l2c_data.TOW, target) + dec = GNSSDecoder.GNSSDecoderState(GNSSDecoder.GPSL2CMDecoderState(state.decoder.prn); + data = l2c_data, raw_data = l2c_data, + num_bits_after_valid_syncro_sequence = num_bits) SatelliteState(; decoder = dec, system = GPSL2CM(), code_phase = code_phase, carrier_doppler = 0.0Hz, carrier_phase = 0.0) end @@ -94,18 +109,32 @@ bcl = PositionVelocityTime.band_ifb_layout # Disjoint: GPS only on L5, Galileo only on L1 ⇒ two components, no IFB column # (each band is the sole band of its component, so its bias folds into a clock). - ifb, extra, ncomp = bcl([GPST(), GPST(), GST(), GST()], [:L5, :L5, :L1, :L1]) + ifb, extra, refs, ncomp = bcl([GPST(), GPST(), GST(), GST()], [:L5, :L5, :L1, :L1]) @test ncomp == 2 @test isempty(extra) + @test isempty(refs) @test all(==(0), ifb) # Connected via a constellation spanning both bands ⇒ one component, one IFB. - _, extra2, ncomp2 = bcl([GPST(), GPST(), GST(), GST()], [:L1, :L5, :L1, :L5]) + # L1 and L5 are equally populated ⇒ first-seen (L1) is the reference. + _, extra2, refs2, ncomp2 = bcl([GPST(), GPST(), GST(), GST()], [:L1, :L5, :L1, :L5]) @test ncomp2 == 1 - @test length(extra2) == 1 + @test extra2 == [:L5] + @test refs2 == [:L1] # Single constellation on two bands ⇒ connected by its shared clock ⇒ one IFB. - _, extra3, ncomp3 = bcl([GPST(), GPST()], [:L1, :L5]) + _, extra3, refs3, ncomp3 = bcl([GPST(), GPST()], [:L1, :L5]) @test ncomp3 == 1 @test length(extra3) == 1 + @test refs3 == [:L1] + # Disconnected yet one component carries an IFB: GPS on L1+L2 (connected) with + # Galileo stranded on L5. The L2 bias is anchored to L1; L5 folds into its clock. + ifb4, extra4, refs4, ncomp4 = bcl( + [GPST(), GPST(), GPST(), GPST(), GST(), GST()], + [:L1, :L1, :L2, :L2, :L5, :L5], + ) + @test ncomp4 == 2 + @test extra4 == [:L2] + @test refs4 == [:L1] + @test ifb4 == [0, 0, 1, 1, 0, 0] end @testset "layout count gate accounts for the extra band" begin @@ -145,10 +174,11 @@ pvt0 = calc_pvt([e1b; map(as_e5a, e1b)]; kw...) @test pvt0.reference_system == GST() @test haskey(pvt0.inter_frequency_biases, :L5) + @test pvt0.inter_frequency_biases[:L5].reference == :L1 @test length(pvt0.sats) == 2 * length(ref.sats) @test norm(pvt0.position - ref.position) < 1e-2 - @test abs(pvt0.inter_frequency_biases[:L5]) < 1e-2 - @test maximum(abs, [info.residual for info in values(pvt0.sats)]) ≈ base atol = 1e-2 + @test abs(pvt0.inter_frequency_biases[:L5].value) < 1e-2m + @test maximum(abs, [info.residual for info in values(pvt0.sats)]) ≈ base atol = 1e-2m # A uniform 12 m receiver L5 delay (signals appear farther ⇒ earlier transmit # time) is absorbed by the IFB, leaving the position and the residuals at the @@ -156,9 +186,9 @@ # residuals instead. δ = 12.0 pvtδ = calc_pvt([e1b; map(s -> as_e5a(s; ifb_shift_s = -δ / C), e1b)]; kw...) - @test pvtδ.inter_frequency_biases[:L5] ≈ δ atol = 0.05 + @test pvtδ.inter_frequency_biases[:L5].value ≈ δ * m atol = 0.05m @test norm(pvtδ.position - ref.position) < 1e-2 - @test maximum(abs, [info.residual for info in values(pvtδ.sats)]) ≈ base atol = 0.05 + @test maximum(abs, [info.residual for info in values(pvtδ.sats)]) ≈ base atol = 0.05m end # GPS two-band fix: L1 C/A + L2C. The L2C copies reproduce each L1 satellite's @@ -175,18 +205,19 @@ pvt0 = calc_pvt([gps; map(as_l2c, gps)]; kw...) @test pvt0.reference_system == GPST() @test haskey(pvt0.inter_frequency_biases, :L2) + @test pvt0.inter_frequency_biases[:L2].reference == :L1 @test length(pvt0.sats) == 2 * length(ref.sats) @test norm(pvt0.position - ref.position) < 1e-2 - @test abs(pvt0.inter_frequency_biases[:L2]) < 1e-2 - @test maximum(abs, [info.residual for info in values(pvt0.sats)]) ≈ base atol = 1e-2 + @test abs(pvt0.inter_frequency_biases[:L2].value) < 1e-2m + @test maximum(abs, [info.residual for info in values(pvt0.sats)]) ≈ base atol = 1e-2m # A uniform 12 m receiver L2 delay is absorbed by the IFB, leaving the position # and residuals at the baseline. δ = 12.0 pvtδ = calc_pvt([gps; map(s -> as_l2c(s; ifb_shift_s = -δ / C), gps)]; kw...) - @test pvtδ.inter_frequency_biases[:L2] ≈ δ atol = 0.05 + @test pvtδ.inter_frequency_biases[:L2].value ≈ δ * m atol = 0.05m @test norm(pvtδ.position - ref.position) < 1e-2 - @test maximum(abs, [info.residual for info in values(pvtδ.sats)]) ≈ base atol = 0.05 + @test maximum(abs, [info.residual for info in values(pvtδ.sats)]) ≈ base atol = 0.05m end # Regression for the disjoint-band bug: GPS on L1 only + Galileo on L5 only makes a @@ -217,12 +248,14 @@ # With a (correct) broadcast GGTO, the collapse reconnects the bands: a clean # inter-frequency bias is recovered and the inter-system bias comes from the GGTO. true_isb = connected.inter_system_biases[GST()] - gal_l5_ggto = map(s -> as_e5a(s; ggto = -true_isb / C), galileo_e1b_states(0.0Hz)) + gal_l5_ggto = + map(s -> as_e5a(s; ggto = -ustrip(m, true_isb) / C), galileo_e1b_states(0.0Hz)) pvt_ggto = calc_pvt([gps_l1; gal_l5_ggto]; kw...) @test pvt_ggto.reference_system == GPST() # Galileo collapsed onto GPS @test haskey(pvt_ggto.inter_frequency_biases, :L5) # reconnected ⇒ IFB observable - @test pvt_ggto.inter_system_biases[GST()] ≈ true_isb atol = 5 - @test abs(pvt_ggto.inter_frequency_biases[:L5]) < 5 + @test pvt_ggto.inter_frequency_biases[:L5].reference == :L1 + @test pvt_ggto.inter_system_biases[GST()] ≈ true_isb atol = 5m + @test abs(pvt_ggto.inter_frequency_biases[:L5].value) < 5m @test isfinite(pvt_ggto.dop.GDOP) && 0 < pvt_ggto.dop.GDOP < 1e4 @test norm(pvt_ggto.position - connected.position) < 10 end diff --git a/test/pvt.jl b/test/pvt.jl index a892e77..5dc825a 100644 --- a/test/pvt.jl +++ b/test/pvt.jl @@ -144,7 +144,7 @@ end resids = [info.residual for info in values(pvt.sats)] @test length(resids) == length(pvt.sats) @test all(isfinite, resids) - @test maximum(abs, resids) < 10 + @test maximum(abs, resids) < 10m 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) @@ -164,7 +164,7 @@ end # reference_system is GPS (most satellites), so the Galileo inter-system bias # is Galileo − GPS = −c·(GST − GPST). The physical broadcast GGTO that would # reproduce it is its negation, in seconds. - Δ = -reference.inter_system_biases[GST()] / C + Δ = -ustrip(m, reference.inter_system_biases[GST()]) / C # 3 GPS + 1 Galileo: the independent solve needs 3 + 2 = 5 satellites, so # without GGTO the constellation is under-determined and calc_pvt returns the @@ -197,7 +197,7 @@ end # Galileo measurement is moved into the GPS frame by subtracting GST − GPST. A_big = 1.0e-6 big = calc_pvt([gps[1:3]; [with_ggto(gal[1]; A_0G = A_big)]]; approximate_year = 2021, enable_ionospheric_correction = false, enable_tropospheric_correction = false) - @test big.inter_system_biases[GST()] ≈ -C * A_big rtol = 1e-6 + @test big.inter_system_biases[GST()] ≈ -C * A_big * m rtol = 1e-6 # calc_ggto_offset evaluates the OS SIS ICD word-type-10 polynomial, taking # the reference week difference modulo 64. @@ -253,3 +253,26 @@ end get_sat_info(pvt, :GalileoE1B, shared_prn) @test length(pvt.sats) == length(baseline.sats) end + +# Course over ground: azimuth of the velocity vector in the local ENU frame, degrees +# clockwise from North. The recorded fixtures are stationary, so exercise the +# computation directly with known ENU velocities rotated into ECEF at a fixed position. +@testset "calc_course_over_ground is the ENU velocity azimuth clockwise from North" begin + user = ECEF(ECEFfromLLA(wgs84)(LLA(50.1, 8.7, 120.0))) + ecef_from_enu = ECEFfromENU(user, wgs84) + # ENU velocity → ECEF velocity vector: rotate the tip and drop the origin translation. + ecef_velocity(e, n, u) = ECEF(ecef_from_enu(ENU(e, n, u))) - user + cog(e, n, u) = PositionVelocityTime.calc_course_over_ground(user, ecef_velocity(e, n, u)) + + @test cog(0.0, 10.0, 0.0) ≈ 0° atol = 1e-6° # North + @test cog(10.0, 0.0, 0.0) ≈ 90° atol = 1e-6° # East + @test cog(0.0, -10.0, 0.0) ≈ 180° atol = 1e-6° # South + @test cog(-10.0, 0.0, 0.0) ≈ 270° atol = 1e-6° # West + @test cog(10.0, 10.0, 0.0) ≈ 45° atol = 1e-6° # North-East + # A stationary receiver has no horizontal velocity ⇒ course defined as 0°. + @test PositionVelocityTime.calc_course_over_ground(user, ECEF(0.0, 0.0, 0.0)) == 0.0° + # The vertical component does not affect the horizontal course. + @test cog(10.0, 0.0, 7.0) ≈ cog(10.0, 0.0, 0.0) + # Always wrapped to [0, 360)°. + @test 0° ≤ cog(-3.0, -4.0, 0.0) < 360° +end diff --git a/test/runtests.jl b/test/runtests.jl index c7ea0da..cb14d18 100755 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,6 @@ using Test, PositionVelocityTime, GNSSDecoder, AstroTime, GNSSSignals, Geodesy, Dates, LinearAlgebra -using Unitful: Hz +using Unitful: Hz, m, °, ustrip include("aqua.jl") include("fixtures.jl")