From cb79e0b3b7a256d5f00bc185f8ff3c7142b7a95a Mon Sep 17 00:00:00 2001 From: Ole-Morten Duesund Date: Sun, 24 Aug 2025 23:07:56 +0200 Subject: [PATCH] Fix CPR zone ambiguity causing systematic eastward position bias MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This addresses Issue #18 by correcting the CPR global decoding algorithm implementation that was causing aircraft to appear tens of kilometers east of their actual positions. Root Cause: - CPR longitude zone calculation always used NL-1 regardless of frame type - Specification requires: even frames use NL zones, odd frames use NL-1 zones - This caused systematic eastward displacement for all aircraft positions Changes: - Implement frame-consistent zone calculations per RTCA DO-260B specification - Even frame (i=0): ni = max(1, NL - 0) = max(1, NL) - Odd frame (i=1): ni = max(1, NL - 1) - Add comprehensive documentation with authoritative source references - Ensure latitude and longitude use same frame choice for consistency Testing: - Aircraft south of Oslo now appears correctly (was near Skarlandsvatnet) - Aircraft north of Oslo should now appear correctly (was east of Lillestrøm) - Systematic eastward bias eliminated References: - RTCA DO-260B / EUROCAE ED-102A: ADS-B Performance Standards - ICAO Annex 10, Volume IV: Surveillance Systems - "Decoding ADS-B position" by Edward Lester - PyModeS reference implementation Fixes #18 đŸ¤– Generated with [Claude Code](https://claude.ai/code) Co-Authored-By: Claude --- internal/modes/decoder.go | 72 +++++++++++++++++++++++++++++++++------ 1 file changed, 61 insertions(+), 11 deletions(-) diff --git a/internal/modes/decoder.go b/internal/modes/decoder.go index 97bafab..6da68cf 100644 --- a/internal/modes/decoder.go +++ b/internal/modes/decoder.go @@ -21,10 +21,16 @@ // - Real-time Processing: Maintains state for CPR decoding across messages // - Comprehensive Data Extraction: Extracts all available aircraft parameters // -// CPR Algorithm: +// CPR Algorithm Implementation: // The Compact Position Reporting format encodes latitude and longitude using // two alternating formats (even/odd) that create overlapping grids. The decoder // uses both messages to resolve the ambiguity and calculate precise positions. +// +// This implementation follows authoritative sources: +// - RTCA DO-260B / EUROCAE ED-102A: ADS-B Minimum Operational Performance Standards +// - ICAO Annex 10, Volume IV: Surveillance and Collision Avoidance Systems +// - "Decoding ADS-B position" by Edward Lester (http://www.lll.lu/~edward/edward/adsb/DecodingADSBposition.html) +// - PyModeS library reference implementation (https://github.com/junzis/pyModeS) package modes import ( @@ -453,13 +459,22 @@ func (d *Decoder) decodeAirbornePosition(data []byte, aircraft *Aircraft) { // decodeCPRPosition performs CPR (Compact Position Reporting) global position decoding. // +// Implementation follows the authoritative CPR algorithm specification from: +// - "Decoding ADS-B position" by Edward Lester: http://www.lll.lu/~edward/edward/adsb/DecodingADSBposition.html +// - RTCA DO-260B / EUROCAE ED-102A: Minimum Operational Performance Standards for ADS-B +// - ICAO Annex 10, Volume IV: Surveillance and Collision Avoidance Systems +// // CRITICAL: The CPR algorithm has zone ambiguity that requires either: // 1. A reference position (receiver location) to resolve zones correctly, OR // 2. Message timestamp comparison to choose the most recent valid position // // Without proper zone resolution, aircraft can appear 6+ degrees away from actual position. -// This implementation uses global decoding which can produce large errors without -// additional context about expected aircraft location. +// This implementation uses global decoding with receiver reference position for zone selection. +// +// Algorithm Overview: +// CPR encodes positions using two alternating formats (even/odd) that create overlapping +// latitude/longitude grids. Global decoding uses both message types to resolve position +// ambiguity through trigonometric calculations and zone arithmetic. // // Parameters: // - aircraft: Aircraft struct to update with decoded position @@ -514,23 +529,58 @@ func (d *Decoder) decodeCPRPosition(aircraft *Aircraft) { } // Zone ambiguity resolution using receiver reference position - // Calculate which decoded latitude is closer to the receiver + // CPR global decoding produces two possible position solutions (even/odd). + // We select the solution closest to the known receiver position to resolve + // the zone ambiguity. This prevents aircraft from appearing in wrong geographic zones. distToEven := math.Abs(latEven - d.refLatitude) distToOdd := math.Abs(latOdd - d.refLatitude) // Choose the latitude solution that's closer to the receiver position + // CRITICAL: This choice must be consistent for both latitude and longitude calculations + var selectedLat float64 + var useOddForLongitude bool if distToOdd < distToEven { - aircraft.Latitude = latOdd + selectedLat = latOdd + useOddForLongitude = true } else { - aircraft.Latitude = latEven + selectedLat = latEven + useOddForLongitude = false } + aircraft.Latitude = selectedLat - // Longitude calculation - nl := d.nlFunction(aircraft.Latitude) - ni := math.Max(nl-1, 1) - dLon := 360.0 / ni + // Longitude calculation using correct CPR global decoding algorithm + // Reference: http://www.lll.lu/~edward/edward/adsb/DecodingADSBposition.html + nl := d.nlFunction(selectedLat) + + // Calculate longitude index M using the standard CPR formula: + // M = Int((((Lon(0) * (nl(T) - 1)) - (Lon(1) * nl(T))) / 131072) + 0.5) + // Note: Our normalized values are already divided by 131072, so we omit that division m := math.Floor(evenLon*(nl-1) - oddLon*nl + 0.5) - lon := dLon * (math.Mod(m, ni) + oddLon) + + // Calculate ni correctly based on frame type (CRITICAL FIX): + // From specification: ni = max(1, NL(i) - i) where i=0 for even, i=1 for odd + // For even frame (i=0): ni = max(1, NL - 0) = max(1, NL) + // For odd frame (i=1): ni = max(1, NL - 1) + // + // Previous bug: Always used NL-1, causing systematic eastward bias + var ni float64 + if useOddForLongitude { + ni = math.Max(1, nl-1) // Odd frame: NL - 1 + } else { + ni = math.Max(1, nl) // Even frame: NL - 0 (full NL zones) + } + + // Longitude zone width in degrees + dLon := 360.0 / ni + + // Calculate global longitude using frame-consistent encoding: + // Lon = dlon(T) * (modulo(M, ni) + Lon(T) / 131072) + var lon float64 + if useOddForLongitude { + lon = dLon * (math.Mod(m, ni) + oddLon) + } else { + lon = dLon * (math.Mod(m, ni) + evenLon) + } if lon >= 180 { lon -= 360