CPR changes. Fixed a lot of the surface position errors. At least one bug still remains.

python/cpr.py

@@ -28,12 +28,6 @@ from air_modes.exceptions import *
 
 latz = 15
 
-def nbits(surface):
-	if surface == 1:
-		return 19
-	else:
-		return 17
-
 def nz(ctype):
 	return 4 * latz - ctype
 
@@ -49,9 +43,6 @@ def dlat(ctype, surface):
 	else:
 		return tmp / nzcalc
 
-def nl_eo(declat_in, ctype):
-	return nl(declat_in) - ctype
-
 def nl(declat_in):
 	if abs(declat_in) >= 87.0:
 		return 1.0
@@ -62,25 +53,23 @@ def dlon(declat_in, ctype, surface):
 		tmp = 90.0
 	else:
 		tmp = 360.0
-	nlcalc = nl_eo(declat_in, ctype)
+	nlcalc = nl(declat_in)-ctype
 	if nlcalc == 0:
 		return tmp
 	else:
-		return tmp / max(nlcalc, 1.0)
+		return tmp / nlcalc
 
 def decode_lat(enclat, ctype, my_lat, surface):
 	tmp1 = dlat(ctype, surface)
-	tmp2 = float(enclat) / (2**nbits(surface))
+	tmp2 = float(enclat) / (2**17)
 	j = math.floor(my_lat/tmp1) + math.floor(0.5 + ((my_lat % tmp1) / tmp1) - tmp2)
-#	print "dlat gives " + "%.6f " % tmp1 + "with j = " + "%.6f " % j + " and tmp2 = " + "%.6f" % tmp2 + " given enclat " + "%x" % enclat
 
 	return tmp1 * (j + tmp2)
 
 def decode_lon(declat, enclon, ctype, my_lon, surface):
 	tmp1 = dlon(declat, ctype, surface)
-	tmp2 = float(enclon) / (2.0**nbits(surface))
+	tmp2 = float(enclon) / (2**17)
 	m = math.floor(my_lon / tmp1) + math.floor(0.5 + ((my_lon % tmp1) / tmp1) - tmp2)
-#	print "dlon gives " + "%.6f " % tmp1 + "with m = " + "%.6f " % m + " and tmp2 = " + "%.6f" % tmp2 + " given enclon " + "%x" % enclon
 
 	return tmp1 * (m + tmp2)
 
@@ -93,7 +82,11 @@ def cpr_resolve_local(my_location, encoded_location, ctype, surface):
 
 	return [decoded_lat, decoded_lon]
 
-def cpr_resolve_global(evenpos, oddpos, mostrecent, surface):
+def cpr_resolve_global(evenpos, oddpos, mypos, mostrecent, surface):
+	#cannot resolve surface positions unambiguously without knowing receiver position
+	if surface and mypos is None:
+		raise CPRNoPositionError
+	
 	dlateven = dlat(0, surface)
 	dlatodd  = dlat(1, surface)
 
@@ -114,7 +107,6 @@ def cpr_resolve_global(evenpos, oddpos, mostrecent, surface):
 	#This checks to see if the latitudes of the reports straddle a transition boundary
 	#If so, you can't get a globally-resolvable location.
 	if nl(rlateven) != nl(rlatodd):
-		#print "Boundary straddle!"
 		raise CPRBoundaryStraddleError
 
 	if mostrecent == 0:
@@ -122,21 +114,41 @@ def cpr_resolve_global(evenpos, oddpos, mostrecent, surface):
 	else:
 		rlat = rlatodd
 
+	#disambiguate latitude
+	if surface:
+		if mypos[0] < 0:
+			rlat -= 90
+
 	dl = dlon(rlat, mostrecent, surface)
 	nlthing = nl(rlat)
-	ni = max(nlthing - mostrecent, 1)
+	ni = nlthing - mostrecent
 
-	m =  math.floor(((evenpos[1]*(nlthing-1)-oddpos[1]*(nlthing))/2**17)+0.5) #longitude index
+	#THIS COLLAPSES WHEN M IS NOT A FACTOR OF FOUR
+	m = math.floor(((evenpos[1]*(nlthing-1)-oddpos[1]*nlthing)/2**17)+0.5) #longitude index
+	if m%4:
+		print "IM GONNA DIE"
+
+	print "DL: %f nlthing: %f ni: %f m: %f" % (dl, nlthing, ni, m)
 
 	if mostrecent == 0:
 		enclon = evenpos[1]
 	else:
 		enclon = oddpos[1]
 
-	rlon = dl * (((ni+m) % ni)+enclon/2**17)
+	rlon = dl * ((m % ni)+enclon/2**17)
 
+	if surface:
+		#longitudes need to be resolved to the nearest 90 degree segment to the receiver.
+		wat = mypos[1]
+		if wat < 0:
+			wat += 360
+		print "mylon: %f mylon unwrapped: %f rlon raw: %f" % (mypos[1], wat, rlon)
+		zone = lambda lon: 90 * (int(lon) / 90)
+		rlon += (zone(wat) - zone(rlon))
+
+	#limit to (-180, 180)
 	if rlon > 180:
-		rlon = rlon - 360.0
+		rlon -= 360.0
 
 	return [rlat, rlon]
 
@@ -199,8 +211,7 @@ class cpr_decoder:
 		if (icao24 in self.evenlist) \
 		  and (icao24 in self.oddlist):
 			newer = (self.oddlist[icao24][2] - self.evenlist[icao24][2]) > 0 #figure out which report is newer
-   			[decoded_lat, decoded_lon] = cpr_resolve_global(self.evenlist[icao24][0:2], self.oddlist[icao24][0:2], newer, surface) #do a global decode
+   			[decoded_lat, decoded_lon] = cpr_resolve_global(self.evenlist[icao24][0:2], self.oddlist[icao24][0:2], self.my_location, newer, surface) #do a global decode
-
 		else:
 			raise CPRNoPositionError
 
@@ -223,11 +234,11 @@ def cpr_encode(lat, lon, ctype, surface):
 	yz = math.floor(scalar * ((lat % dlati)/dlati) + 0.5)
 	rlat = dlati * ((yz / scalar) + math.floor(lat / dlati))
 
-	nleo = nl_eo(rlat, ctype)
+	nleo = nl(rlat)-ctype
 	if nleo == 0:
 		dloni = 360.0
 	else:
-		dloni = 360.0 / nl_eo(rlat, ctype)
+		dloni = 360.0 / nleo
 
 	xz = math.floor(scalar * ((lon % dloni)/dloni) + 0.5)
 
@@ -245,53 +256,57 @@ if __name__ == '__main__':
 	#latitudes the corresponding error in longitude is greater
 
 	bs = 0
+	surface = True
 
 	for i in range(0, rounds):
-		decoder = cpr_decoder(None)
 		even_lat = random.uniform(-85, 85)
 		even_lon = random.uniform(-180, 180)
 		odd_lat = even_lat + 1e-2
 		odd_lon = min(even_lon + 1e-2, 180)
+		decoder = cpr_decoder([odd_lat, odd_lon])
 
 		#encode that position
-		(evenenclat, evenenclon) = cpr_encode(even_lat, even_lon, False, False)
+		(evenenclat, evenenclon) = cpr_encode(even_lat, even_lon, False, surface)
-		(oddenclat, oddenclon)   = cpr_encode(odd_lat, odd_lon, True, False)
+		(oddenclat, oddenclon)   = cpr_encode(odd_lat, odd_lon, True, surface)
 
 		#perform a global decode
 		icao = random.randint(0, 0xffffff)
 		try:
-			evenpos = decoder.decode(icao, evenenclat, evenenclon, False, False)
+			evenpos = decoder.decode(icao, evenenclat, evenenclon, False, surface)
 			#print "CPR global decode with only one report: %f %f" % (evenpos[0], evenpos[1])
 			raise Exception("CPR test failure: global decode with only one report")
 		except CPRNoPositionError:
 			pass
 
 		try:
-			(odddeclat, odddeclon, rng, brg) = decoder.decode(icao, oddenclat, oddenclon, True, False)
+			(odddeclat, odddeclon, rng, brg) = decoder.decode(icao, oddenclat, oddenclon, True, surface)
 		except CPRBoundaryStraddleError:
 			bs += 1
 			continue
 		except CPRNoPositionError:
 			raise Exception("CPR test failure: no decode after even/odd inputs")
 
-		#print "Lat: %f Lon: %f" % (ac_lat, ac_lon)
-
 		if abs(odddeclat - odd_lat) > threshold or abs(odddeclon - odd_lon) > threshold:
-			print "odddeclat: %f odd_lat: %f" % (odddeclat, odd_lat)
+			print "F odddeclat: %f odd_lat: %f" % (odddeclat, odd_lat)
-			print "odddeclon: %f odd_lon: %f" % (odddeclon, odd_lon)
+			print "F odddeclon: %f odd_lon: %f" % (odddeclon, odd_lon)
 			raise Exception("CPR test failure: global decode error greater than threshold")
+		else:
+			print "S odddeclat: %f odd_lat: %f" % (odddeclat, odd_lat)
+			print "S odddeclon: %f odd_lon: %f" % (odddeclon, odd_lon)
 
 		nexteven_lat = odd_lat + 1e-2
 		nexteven_lon = min(odd_lon + 1e-2, 180)
 
-		(nexteven_enclat, nexteven_enclon) = cpr_encode(nexteven_lat, nexteven_lon, False, False)
+		(nexteven_enclat, nexteven_enclon) = cpr_encode(nexteven_lat, nexteven_lon, False, surface)
 
 		try:
-			(evendeclat, evendeclon) = cpr_resolve_local([even_lat, even_lon], [nexteven_enclat, nexteven_enclon], False, False)
+			(evendeclat, evendeclon) = cpr_resolve_local([even_lat, even_lon], [nexteven_enclat, nexteven_enclon], False, surface)
 		except CPRNoPositionError:
 			raise Exception("CPR test failure: local decode failure to resolve")
 		
 		if abs(evendeclat - nexteven_lat) > threshold or abs(evendeclon - nexteven_lon) > threshold:
+			print "F evendeclat: %f nexteven_lat: %f evenlat: %f" % (evendeclat, nexteven_lat, even_lat)
+			print "F evendeclon: %f nexteven_lon: %f evenlon: %f" % (evendeclon, nexteven_lon, even_lon)
 			raise Exception("CPR test failure: local decode error greater than threshold")
 
 	print "CPR test successful. There were %i boundary straddles over %i rounds." % (bs, rounds)