Merge pull request 'Fix CPR zone ambiguity causing systematic eastward position bias' (#20) from fix/cpr-zone-ambiguity into main

Reviewed-on: #20

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 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)
+	
+	// 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 calculation
-	nl := d.nlFunction(aircraft.Latitude)
-	ni := math.Max(nl-1, 1)
+	// Longitude zone width in degrees
 	dLon := 360.0 / ni
-	m := math.Floor(evenLon*(nl-1) - oddLon*nl + 0.5)
-	lon := dLon * (math.Mod(m, ni) + oddLon)
+	
+	// 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