FGFS-4.1/gr-air-modes
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Big update to UHD 3.14, Gnuradio 3.8, Python 3.6. Not fully tested.

Browse Source
This commit is contained in:
Nick Foster
2019-09-17 14:13:51 -07:00
parent 0b6c383506
commit 9b17824c49
36 changed files with 544 additions and 1684 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
python/CMakeLists.txt
View File

@@ -40,6 +40,7 @@ GR_PYTHON_INSTALL(
flightgear.py
gui_model.py
kml.py
modes_types.py
parse.py
msprint.py
radio.py
@@ -47,7 +48,6 @@ GR_PYTHON_INSTALL(
rx_path.py
sbs1.py
sql.py
types.py
zmq_socket.py
Quaternion.py
DESTINATION ${GR_PYTHON_DIR}/air_modes

60
python/__init__.py
View File

@@ -28,26 +28,10 @@ KML, and PlanePlotter-compliant SBS-1 emulation output. mlat.py provides
an experimental implementation of a multilateration solver.
'''
# ----------------------------------------------------------------
# Temporary workaround for ticket:181 (swig+python problem)
import sys
_RTLD_GLOBAL = 0
try:
from dl import RTLD_GLOBAL as _RTLD_GLOBAL
except ImportError:
try:
from DLFCN import RTLD_GLOBAL as _RTLD_GLOBAL
except ImportError:
pass
if _RTLD_GLOBAL != 0:
_dlopenflags = sys.getdlopenflags()
sys.setdlopenflags(_dlopenflags|_RTLD_GLOBAL)
# ----------------------------------------------------------------
from __future__ import unicode_literals
# import swig generated symbols into the gr-air-modes namespace
from air_modes_swig import *
from .air_modes_swig import *
# import any pure python here
#
@@ -57,30 +41,24 @@ try:
except ImportError:
raise RuntimeError("PyZMQ not found! Please install libzmq and PyZMQ to run gr-air-modes")
from rx_path import rx_path
from zmq_socket import zmq_pubsub_iface
from parse import *
from msprint import output_print
from sql import output_sql
from sbs1 import output_sbs1
from kml import output_kml, output_jsonp
from raw_server import raw_server
from radio import modes_radio
from exceptions import *
from types import *
from altitude import *
from cpr import cpr_decoder
from html_template import html_template
from .rx_path import rx_path
from .zmq_socket import zmq_pubsub_iface
from .parse import *
from .msprint import output_print
from .sql import output_sql
from .sbs1 import output_sbs1
from .kml import output_kml, output_jsonp
from .raw_server import raw_server
from .radio import modes_radio
from .exceptions import *
from .modes_types import *
from .altitude import *
from .cpr import cpr_decoder
from .html_template import html_template
#this is try/excepted in case the user doesn't have numpy installed
try:
from flightgear import output_flightgear
from Quaternion import *
from .flightgear import output_flightgear
from .Quaternion import *
except ImportError:
print "gr-air-modes warning: numpy+scipy not installed, FlightGear interface not supported"
print("gr-air-modes warning: numpy+scipy not installed, FlightGear interface not supported")
pass
# ----------------------------------------------------------------
# Tail of workaround
if _RTLD_GLOBAL != 0:
sys.setdlopenflags(_dlopenflags) # Restore original flags
# ----------------------------------------------------------------

189
python/altitude.py
View File

@@ -26,120 +26,119 @@
from air_modes.exceptions import *
def decode_alt(alt, bit13):
mbit = alt & 0x0040
qbit = alt & 0x0010
if mbit and bit13:
#nobody uses metric altitude: AFAIK, it's an orphaned part of
#the spec. haven't seen it in three years. as a result, replies
#with mbit set can be considered spurious, and so we discard them here.
#bits 20-25, 27-31 encode alt in meters
#remember that bits are LSB (bit 20 is MSB)
#meters_alt = 0
#for (shift, bit) in enumerate(range(31,26,-1)+range(25,19,-1)):
# meters_alt += ((alt & (1<<bit)) != 0) << shift
#decoded_alt = meters_alt / 0.3048
raise MetricAltError
mbit = alt & 0x0040
qbit = alt & 0x0010
if mbit and bit13:
#nobody uses metric altitude: AFAIK, it's an orphaned part of
#the spec. haven't seen it in three years. as a result, replies
#with mbit set can be considered spurious, and so we discard them here.
#bits 20-25, 27-31 encode alt in meters
#remember that bits are LSB (bit 20 is MSB)
#meters_alt = 0
#for (shift, bit) in enumerate(range(31,26,-1)+range(25,19,-1)):
# meters_alt += ((alt & (1<<bit)) != 0) << shift
#decoded_alt = meters_alt / 0.3048
raise MetricAltError
if qbit: #a mode S-style reply
#bit13 is false for BDS0,5 ADS-B squitters, and is true otherwise
if bit13:
#in this representation, the altitude bits are as follows:
# 12 11 10 9 8 7 (6) 5 (4) 3 2 1 0
# so bits 6 and 4 are the M and Q bits, respectively.
tmp1 = (alt & 0x3F80) >> 2
tmp2 = (alt & 0x0020) >> 1
else:
tmp1 = (alt & 0x1FE0) >> 1
tmp2 = 0
if qbit: #a mode S-style reply
#bit13 is false for BDS0,5 ADS-B squitters, and is true otherwise
if bit13:
#in this representation, the altitude bits are as follows:
# 12 11 10 9 8 7 (6) 5 (4) 3 2 1 0
# so bits 6 and 4 are the M and Q bits, respectively.
tmp1 = (alt & 0x3F80) >> 2
tmp2 = (alt & 0x0020) >> 1
else:
tmp1 = (alt & 0x1FE0) >> 1
tmp2 = 0
decoded_alt = ((alt & 0x0F) | tmp1 | tmp2) * 25 - 1000
decoded_alt = ((alt & 0x0F) | tmp1 | tmp2) * 25 - 1000
else: #a mode C-style reply
#okay, the order they come in is:
#C1 A1 C2 A2 C4 A4 X B1 D1 B2 D2 B4 D4
#the order we want them in is:
#D2 D4 A1 A2 A4 B1 B2 B4
#so we'll reassemble into a Gray-coded representation
else: #a mode C-style reply
#okay, the order they come in is:
#C1 A1 C2 A2 C4 A4 X B1 D1 B2 D2 B4 D4
#the order we want them in is:
#D2 D4 A1 A2 A4 B1 B2 B4
#so we'll reassemble into a Gray-coded representation
if bit13 is False:
alt = (alt & 0x003F) | (alt & 0x0FC0 << 1)
if bit13 is False:
alt = (alt & 0x003F) | (alt & 0x0FC0 << 1)
C1 = 0x1000
A1 = 0x0800
C2 = 0x0400
A2 = 0x0200 #this represents the order in which the bits come
C4 = 0x0100
A4 = 0x0080
B1 = 0x0020
D1 = 0x0010
B2 = 0x0008
D2 = 0x0004
B4 = 0x0002
D4 = 0x0001
C1 = 0x1000
A1 = 0x0800
C2 = 0x0400
A2 = 0x0200 #this represents the order in which the bits come
C4 = 0x0100
A4 = 0x0080
B1 = 0x0020
D1 = 0x0010
B2 = 0x0008
D2 = 0x0004
B4 = 0x0002
D4 = 0x0001
bigpart = ((alt & B4) >> 1) \
+ ((alt & B2) >> 2) \
+ ((alt & B1) >> 3) \
+ ((alt & A4) >> 4) \
+ ((alt & A2) >> 5) \
+ ((alt & A1) >> 6) \
+ ((alt & D4) << 6) \
+ ((alt & D2) << 5)
bigpart = ((alt & B4) >> 1) \
+ ((alt & B2) >> 2) \
+ ((alt & B1) >> 3) \
+ ((alt & A4) >> 4) \
+ ((alt & A2) >> 5) \
+ ((alt & A1) >> 6) \
+ ((alt & D4) << 6) \
+ ((alt & D2) << 5)
#bigpart is now the 500-foot-resolution Gray-coded binary part
decoded_alt = gray2bin(bigpart)
#real_alt is now the 500-foot-per-tick altitude
#bigpart is now the 500-foot-resolution Gray-coded binary part
decoded_alt = gray2bin(bigpart)
#real_alt is now the 500-foot-per-tick altitude
cbits = ((alt & C4) >> 8) + ((alt & C2) >> 9) + ((alt & C1) >> 10)
cval = gray2bin(cbits) #turn them into a real number
cbits = ((alt & C4) >> 8) + ((alt & C2) >> 9) + ((alt & C1) >> 10)
cval = gray2bin(cbits) #turn them into a real number
if cval == 7:
cval = 5 #not a real gray code after all
if cval == 7:
cval = 5 #not a real gray code after all
if decoded_alt % 2:
cval = 6 - cval #since the code is symmetric this unwraps it to see whether to subtract the C bits or add them
if decoded_alt % 2:
cval = 6 - cval #since the code is symmetric this unwraps it to see whether to subtract the C bits or add them
decoded_alt *= 500 #take care of the A,B,D data
decoded_alt += cval * 100 #factor in the C data
decoded_alt -= 1300 #subtract the offset
decoded_alt *= 500 #take care of the A,B,D data
decoded_alt += cval * 100 #factor in the C data
decoded_alt -= 1300 #subtract the offset
return decoded_alt
return decoded_alt
def gray2bin(gray):
i = gray >> 1
i = gray >> 1
while i != 0:
gray ^= i
i >>= 1
while i != 0:
gray ^= i
i >>= 1
return gray
return gray
def encode_alt_modes(alt, bit13):
mbit = False
qbit = True
encalt = (int(alt) + 1000) / 25
mbit = False
qbit = True
encalt = (int(alt) + 1000) / 25
if bit13 is True:
tmp1 = (encalt & 0xfe0) << 2
tmp2 = (encalt & 0x010) << 1
else:
tmp1 = (encalt & 0xff8) << 1
tmp2 = 0
if bit13 is True:
tmp1 = (encalt & 0xfe0) << 2
tmp2 = (encalt & 0x010) << 1
else:
tmp1 = (encalt & 0xff8) << 1
tmp2 = 0
return (encalt & 0x0F) | tmp1 | tmp2 | (mbit << 6) | (qbit << 4)
return (encalt & 0x0F) | tmp1 | tmp2 | (mbit << 6) | (qbit << 4)
if __name__ == "__main__":
try:
for alt in range(-1000, 101400, 25):
dec = decode_alt(encode_alt_modes(alt, False), False)
if dec != alt:
print "Failure at %i with bit13 clear (got %s)" % (alt, dec)
for alt in range(-1000, 101400, 25):
dec = decode_alt(encode_alt_modes(alt, True), True)
if dec != alt:
print "Failure at %i with bit13 set (got %s)" % (alt, dec)
except MetricAltError:
print "Failure at %i due to metric alt bit" % alt
try:
for alt in range(-1000, 101400, 25):
dec = decode_alt(encode_alt_modes(alt, False), False)
if dec != alt:
print("Failure at %i with bit13 clear (got %s)" % (alt, dec))
for alt in range(-1000, 101400, 25):
dec = decode_alt(encode_alt_modes(alt, True), True)
if dec != alt:
print("Failure at %i with bit13 set (got %s)" % (alt, dec))
except MetricAltError:
print("Failure at %i due to metric alt bit" % alt)

1
python/az_map.py
View File

@@ -25,7 +25,6 @@
from PyQt4 import QtCore, QtGui
import threading
import math
import air_modes
from air_modes.exceptions import *
import numpy as np

442
python/cpr.py
View File

@@ -31,302 +31,302 @@ from air_modes.exceptions import *
latz = 15
def nz(ctype):
return 4 * latz - ctype
return 4 * latz - ctype
def dlat(ctype, surface):
if surface == 1:
tmp = 90.0
else:
tmp = 360.0
if surface == 1:
tmp = 90.0
else:
tmp = 360.0
nzcalc = nz(ctype)
if nzcalc == 0:
return tmp
else:
return tmp / nzcalc
nzcalc = nz(ctype)
if nzcalc == 0:
return tmp
else:
return tmp / nzcalc
def nl(declat_in):
if abs(declat_in) >= 87.0:
return 1.0
return math.floor( (2.0*math.pi) * math.acos(1.0- (1.0-math.cos(math.pi/(2.0*latz))) / math.cos( (math.pi/180.0)*abs(declat_in) )**2 )**-1)
if abs(declat_in) >= 87.0:
return 1.0
return math.floor( (2.0*math.pi) * math.acos(1.0- (1.0-math.cos(math.pi/(2.0*latz))) / math.cos( (math.pi/180.0)*abs(declat_in) )**2 )**-1)
def dlon(declat_in, ctype, surface):
if surface:
tmp = 90.0
else:
tmp = 360.0
nlcalc = max(nl(declat_in)-ctype, 1)
return tmp / nlcalc
if surface:
tmp = 90.0
else:
tmp = 360.0
nlcalc = max(nl(declat_in)-ctype, 1)
return tmp / nlcalc
def decode_lat(enclat, ctype, my_lat, surface):
tmp1 = dlat(ctype, surface)
tmp2 = float(enclat) / (2**17)
j = math.floor(my_lat/tmp1) + math.floor(0.5 + ((my_lat % tmp1) / tmp1) - tmp2)
tmp1 = dlat(ctype, surface)
tmp2 = float(enclat) / (2**17)
j = math.floor(my_lat/tmp1) + math.floor(0.5 + ((my_lat % tmp1) / tmp1) - tmp2)
return tmp1 * (j + tmp2)
return tmp1 * (j + tmp2)
def decode_lon(declat, enclon, ctype, my_lon, surface):
tmp1 = dlon(declat, ctype, surface)
tmp2 = float(enclon) / (2**17)
m = math.floor(my_lon / tmp1) + math.floor(0.5 + ((my_lon % tmp1) / tmp1) - tmp2)
tmp1 = dlon(declat, ctype, surface)
tmp2 = float(enclon) / (2**17)
m = math.floor(my_lon / tmp1) + math.floor(0.5 + ((my_lon % tmp1) / tmp1) - tmp2)
return tmp1 * (m + tmp2)
return tmp1 * (m + tmp2)
def cpr_resolve_local(my_location, encoded_location, ctype, surface):
[my_lat, my_lon] = my_location
[enclat, enclon] = encoded_location
[my_lat, my_lon] = my_location
[enclat, enclon] = encoded_location
decoded_lat = decode_lat(enclat, ctype, my_lat, surface)
decoded_lon = decode_lon(decoded_lat, enclon, ctype, my_lon, surface)
decoded_lat = decode_lat(enclat, ctype, my_lat, surface)
decoded_lon = decode_lon(decoded_lat, enclon, ctype, my_lon, surface)
return [decoded_lat, decoded_lon]
return [decoded_lat, decoded_lon]
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)
#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)
evenpos = [float(evenpos[0]), float(evenpos[1])]
oddpos = [float(oddpos[0]), float(oddpos[1])]
j = math.floor(((nz(1)*evenpos[0] - nz(0)*oddpos[0])/2**17) + 0.5) #latitude index
evenpos = [float(evenpos[0]), float(evenpos[1])]
oddpos = [float(oddpos[0]), float(oddpos[1])]
j = math.floor(((nz(1)*evenpos[0] - nz(0)*oddpos[0])/2**17) + 0.5) #latitude index
rlateven = dlateven * ((j % nz(0))+evenpos[0]/2**17)
rlatodd = dlatodd * ((j % nz(1))+ oddpos[0]/2**17)
rlateven = dlateven * ((j % nz(0))+evenpos[0]/2**17)
rlatodd = dlatodd * ((j % nz(1))+ oddpos[0]/2**17)
#limit to -90, 90
if rlateven > 270.0:
rlateven -= 360.0
if rlatodd > 270.0:
rlatodd -= 360.0
#limit to -90, 90
if rlateven > 270.0:
rlateven -= 360.0
if rlatodd > 270.0:
rlatodd -= 360.0
#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):
raise CPRBoundaryStraddleError
#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):
raise CPRBoundaryStraddleError
if mostrecent == 0:
rlat = rlateven
else:
rlat = rlatodd
if mostrecent == 0:
rlat = rlateven
else:
rlat = rlatodd
#disambiguate latitude
if surface:
if mypos[0] < 0:
rlat -= 90
#disambiguate latitude
if surface:
if mypos[0] < 0:
rlat -= 90
dl = dlon(rlat, mostrecent, surface)
nl_rlat = nl(rlat)
dl = dlon(rlat, mostrecent, surface)
nl_rlat = nl(rlat)
m = math.floor(((evenpos[1]*(nl_rlat-1)-oddpos[1]*nl_rlat)/2**17)+0.5) #longitude index
#when surface positions straddle a disambiguation boundary (90 degrees),
#surface decoding will fail. this might never be a problem in real life, but it'll fail in the
#test case. the documentation doesn't mention it.
m = math.floor(((evenpos[1]*(nl_rlat-1)-oddpos[1]*nl_rlat)/2**17)+0.5) #longitude index
#when surface positions straddle a disambiguation boundary (90 degrees),
#surface decoding will fail. this might never be a problem in real life, but it'll fail in the
#test case. the documentation doesn't mention it.
if mostrecent == 0:
enclon = evenpos[1]
else:
enclon = oddpos[1]
if mostrecent == 0:
enclon = evenpos[1]
else:
enclon = oddpos[1]
rlon = dl * ((m % max(nl_rlat-mostrecent,1)) + enclon/2.**17)
rlon = dl * ((m % max(nl_rlat-mostrecent,1)) + enclon/2.**17)
#print "DL: %f nl: %f m: %f rlon: %f" % (dl, nl_rlat, m, rlon)
#print "evenpos: %x, oddpos: %x, mostrecent: %i" % (evenpos[1], oddpos[1], mostrecent)
#print "DL: %f nl: %f m: %f rlon: %f" % (dl, nl_rlat, m, rlon)
#print "evenpos: %x, oddpos: %x, mostrecent: %i" % (evenpos[1], oddpos[1], mostrecent)
if surface:
#longitudes need to be resolved to the nearest 90 degree segment to the receiver.
wat = mypos[1]
if wat < 0:
wat += 360
zone = lambda lon: 90 * (int(lon) / 90)
rlon += (zone(wat) - zone(rlon))
if surface:
#longitudes need to be resolved to the nearest 90 degree segment to the receiver.
wat = mypos[1]
if wat < 0:
wat += 360
zone = lambda lon: 90 * (int(lon) / 90)
rlon += (zone(wat) - zone(rlon))
#limit to (-180, 180)
if rlon > 180:
rlon -= 360.0
#limit to (-180, 180)
if rlon > 180:
rlon -= 360.0
return [rlat, rlon]
return [rlat, rlon]
#calculate range and bearing between two lat/lon points
#should probably throw this in the mlat py somewhere or make another lib
def range_bearing(loc_a, loc_b):
[a_lat, a_lon] = loc_a
[b_lat, b_lon] = loc_b
[a_lat, a_lon] = loc_a
[b_lat, b_lon] = loc_b
esquared = (1/298.257223563)*(2-(1/298.257223563))
earth_radius_mi = 3963.19059 * (math.pi / 180)
esquared = (1/298.257223563)*(2-(1/298.257223563))
earth_radius_mi = 3963.19059 * (math.pi / 180)
delta_lat = b_lat - a_lat
delta_lon = b_lon - a_lon
delta_lat = b_lat - a_lat
delta_lon = b_lon - a_lon
avg_lat = ((a_lat + b_lat) / 2.0) * math.pi / 180
avg_lat = ((a_lat + b_lat) / 2.0) * math.pi / 180
R1 = earth_radius_mi*(1.0-esquared)/pow((1.0-esquared*pow(math.sin(avg_lat),2)),1.5)
R1 = earth_radius_mi*(1.0-esquared)/pow((1.0-esquared*pow(math.sin(avg_lat),2)),1.5)
R2 = earth_radius_mi/math.sqrt(1.0-esquared*pow(math.sin(avg_lat),2))
R2 = earth_radius_mi/math.sqrt(1.0-esquared*pow(math.sin(avg_lat),2))
distance_North = R1*delta_lat
distance_East = R2*math.cos(avg_lat)*delta_lon
distance_North = R1*delta_lat
distance_East = R2*math.cos(avg_lat)*delta_lon
bearing = math.atan2(distance_East,distance_North) * (180.0 / math.pi)
if bearing < 0.0:
bearing += 360.0
bearing = math.atan2(distance_East,distance_North) * (180.0 / math.pi)
if bearing < 0.0:
bearing += 360.0
rnge = math.hypot(distance_East,distance_North)
return [rnge, bearing]
rnge = math.hypot(distance_East,distance_North)
return [rnge, bearing]
class cpr_decoder:
def __init__(self, my_location):
self.my_location = my_location
self.evenlist = {}
self.oddlist = {}
self.evenlist_sfc = {}
self.oddlist_sfc = {}
def __init__(self, my_location):
self.my_location = my_location
self.evenlist = {}
self.oddlist = {}
self.evenlist_sfc = {}
self.oddlist_sfc = {}
def set_location(self, new_location):
self.my_location = new_location
def set_location(self, new_location):
self.my_location = new_location
def weed_poslists(self):
for poslist in [self.evenlist, self.oddlist]:
for key, item in poslist.items():
if time.time() - item[2] > 10:
del poslist[key]
for poslist in [self.evenlist_sfc, self.oddlist_sfc]:
for key, item in poslist.items():
if time.time() - item[2] > 25:
del poslist[key]
def weed_poslists(self):
for poslist in [self.evenlist, self.oddlist]:
for key, item in poslist.items():
if time.time() - item[2] > 10:
del poslist[key]
for poslist in [self.evenlist_sfc, self.oddlist_sfc]:
for key, item in poslist.items():
if time.time() - item[2] > 25:
del poslist[key]
def decode(self, icao24, encoded_lat, encoded_lon, cpr_format, surface):
if surface:
oddlist = self.oddlist_sfc
evenlist = self.evenlist_sfc
else:
oddlist = self.oddlist
evenlist = self.evenlist
def decode(self, icao24, encoded_lat, encoded_lon, cpr_format, surface):
if surface:
oddlist = self.oddlist_sfc
evenlist = self.evenlist_sfc
else:
oddlist = self.oddlist
evenlist = self.evenlist
#add the info to the position reports list for global decoding
if cpr_format==1:
oddlist[icao24] = [encoded_lat, encoded_lon, time.time()]
else:
evenlist[icao24] = [encoded_lat, encoded_lon, time.time()]
#add the info to the position reports list for global decoding
if cpr_format==1:
oddlist[icao24] = [encoded_lat, encoded_lon, time.time()]
else:
evenlist[icao24] = [encoded_lat, encoded_lon, time.time()]
[decoded_lat, decoded_lon] = [None, None]
[decoded_lat, decoded_lon] = [None, None]
#okay, let's traverse the lists and weed out those entries that are older than 10 seconds
self.weed_poslists()
#okay, let's traverse the lists and weed out those entries that are older than 10 seconds
self.weed_poslists()
if (icao24 in evenlist) \
and (icao24 in oddlist):
newer = (oddlist[icao24][2] - evenlist[icao24][2]) > 0 #figure out which report is newer
[decoded_lat, decoded_lon] = cpr_resolve_global(evenlist[icao24][0:2], oddlist[icao24][0:2], self.my_location, newer, surface) #do a global decode
else:
raise CPRNoPositionError
if (icao24 in evenlist) \
and (icao24 in oddlist):
newer = (oddlist[icao24][2] - evenlist[icao24][2]) > 0 #figure out which report is newer
[decoded_lat, decoded_lon] = cpr_resolve_global(evenlist[icao24][0:2], oddlist[icao24][0:2], self.my_location, newer, surface) #do a global decode
else:
raise CPRNoPositionError
if self.my_location is not None:
[rnge, bearing] = range_bearing(self.my_location, [decoded_lat, decoded_lon])
else:
rnge = None
bearing = None
if self.my_location is not None:
[rnge, bearing] = range_bearing(self.my_location, [decoded_lat, decoded_lon])
else:
rnge = None
bearing = None
return [decoded_lat, decoded_lon, rnge, bearing]
return [decoded_lat, decoded_lon, rnge, bearing]
#encode CPR position
def cpr_encode(lat, lon, ctype, surface):
if surface is True:
scalar = 2.**19
else:
scalar = 2.**17
if surface is True:
scalar = 2.**19
else:
scalar = 2.**17
#encode using 360 constant for segment size.
dlati = dlat(ctype, False)
yz = math.floor(scalar * ((lat % dlati)/dlati) + 0.5)
rlat = dlati * ((yz / scalar) + math.floor(lat / dlati))
#encode using 360 constant for segment size.
dlati = dlat(ctype, False)
yz = math.floor(scalar * ((lat % dlati)/dlati) + 0.5)
rlat = dlati * ((yz / scalar) + math.floor(lat / dlati))
#encode using 360 constant for segment size.
dloni = dlon(lat, ctype, False)
xz = math.floor(scalar * ((lon % dloni)/dloni) + 0.5)
#encode using 360 constant for segment size.
dloni = dlon(lat, ctype, False)
xz = math.floor(scalar * ((lon % dloni)/dloni) + 0.5)
yz = int(yz) & (2**17-1)
xz = int(xz) & (2**17-1)
yz = int(yz) & (2**17-1)
xz = int(xz) & (2**17-1)
return (yz, xz) #lat, lon
return (yz, xz) #lat, lon
if __name__ == '__main__':
import sys, random
rounds = 10001
threshold = 1e-3 #0.001 deg lat/lon
#this accuracy is highly dependent on latitude, since at high
#latitudes the corresponding error in longitude is greater
import sys, random
rounds = 10001
threshold = 1e-3 #0.001 deg lat/lon
#this accuracy is highly dependent on latitude, since at high
#latitudes the corresponding error in longitude is greater
bs = 0
surface = False
bs = 0
surface = False
lats = [i/(rounds/170.)-85 for i in range(0,rounds)]
lons = [i/(rounds/360.)-180 for i in range(0,rounds)]
lats = [i/(rounds/170.)-85 for i in range(0,rounds)]
lons = [i/(rounds/360.)-180 for i in range(0,rounds)]
for i in range(0, rounds):
even_lat = lats[i]
#even_lat = random.uniform(-85, 85)
even_lon = lons[i]
#even_lon = random.uniform(-180, 180)
odd_lat = even_lat + 1e-3
odd_lon = min(even_lon + 1e-3, 180)
decoder = cpr_decoder([odd_lat, odd_lon])
for i in range(0, rounds):
even_lat = lats[i]
#even_lat = random.uniform(-85, 85)
even_lon = lons[i]
#even_lon = random.uniform(-180, 180)
odd_lat = even_lat + 1e-3
odd_lon = min(even_lon + 1e-3, 180)
decoder = cpr_decoder([odd_lat, odd_lon])
#encode that position
(evenenclat, evenenclon) = cpr_encode(even_lat, even_lon, False, surface)
(oddenclat, oddenclon) = cpr_encode(odd_lat, odd_lon, True, surface)
#encode that position
(evenenclat, evenenclon) = cpr_encode(even_lat, even_lon, False, surface)
(oddenclat, oddenclon) = cpr_encode(odd_lat, odd_lon, True, surface)
#try to perform a global decode -- this should fail since the decoder
#only has heard one position. need two for global decoding.
icao = random.randint(0, 0xffffff)
try:
evenpos = decoder.decode(icao, evenenclat, evenenclon, False, surface)
raise Exception("CPR test failure: global decode with only one report")
except CPRNoPositionError:
pass
#try to perform a global decode -- this should fail since the decoder
#only has heard one position. need two for global decoding.
icao = random.randint(0, 0xffffff)
try:
evenpos = decoder.decode(icao, evenenclat, evenenclon, False, surface)
raise Exception("CPR test failure: global decode with only one report")
except CPRNoPositionError:
pass
#now try to do a real decode with the last packet's odd complement
#watch for a boundary straddle -- this isn't fatal, it just indicates
#that the even and odd reports lie on either side of a longitudinal boundary
#and so you can't get a position
try:
(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")
#now try to do a real decode with the last packet's odd complement
#watch for a boundary straddle -- this isn't fatal, it just indicates
#that the even and odd reports lie on either side of a longitudinal boundary
#and so you can't get a position
try:
(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")
if abs(odddeclat - odd_lat) > threshold or abs(odddeclon - odd_lon) > threshold:
print "F odddeclat: %f odd_lat: %f" % (odddeclat, odd_lat)
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)
if abs(odddeclat - odd_lat) > threshold or abs(odddeclon - odd_lon) > threshold:
print("F odddeclat: %f odd_lat: %f" % (odddeclat, odd_lat))
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-3
nexteven_lon = min(odd_lon + 1e-3, 180)
nexteven_lat = odd_lat + 1e-3
nexteven_lon = min(odd_lon + 1e-3, 180)
(nexteven_enclat, nexteven_enclon) = cpr_encode(nexteven_lat, nexteven_lon, False, surface)
(nexteven_enclat, nexteven_enclon) = cpr_encode(nexteven_lat, nexteven_lon, False, surface)
#try a locally-referenced decode
try:
(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")
#try a locally-referenced decode
try:
(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")
#check to see if the positions were valid
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")
#check to see if the positions were valid
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)
print("CPR test successful. There were %i boundary straddles over %i rounds." % (bs, rounds))

8
python/flightgear.py
View File

@@ -24,7 +24,7 @@ class output_flightgear:
self._cpr = cprdec
self.sock = socket.socket(socket.AF_INET, socket.SOCK_DGRAM)
self.sock.connect((self.hostname, self.port))
# self.sock.connect((self.hostname, self.port))
pub.subscribe("type17_dl", self.output)
def output(self, msg):
@@ -69,7 +69,7 @@ class output_flightgear:
and (icao24 in self.velocities)\
and (icao24 in self.callsigns)
if complete:
print "FG update: %s" % (self.callsigns[icao24][0])
print("FG update: %s" % (self.callsigns[icao24][0]))
msg = fg_posmsg(self.callsigns[icao24][0],
self.callsigns[icao24][1],
self.positions[icao24][0],
@@ -80,7 +80,7 @@ class output_flightgear:
self.velocities[icao24][2],
self.velocities[icao24][3]).pack()
self.sock.send(msg)
self.sock.sendto(msg, (self.hostname, self.port))
class fg_header:
def __init__(self):
@@ -119,7 +119,7 @@ modelmap = { None: 'Aircraft/777-200/Models/777-200ER.xml'
"SMALL": 'Aircraft/CitationX/Models/Citation-X.xml',
"LARGE": 'Aircraft/CRJ700-family/Models/CRJ700.xml',
"LARGE HIGH VORTEX": 'Aircraft/757-200/Models/757-200.xml',
"HEAVY": 'Aircraft/747-200/Models/boeing747-200.xml',
"HEAVY": 'Aircraft/IDG-A32X/Models/A320neo-CFM.xml',
"HIGH PERFORMANCE": 'Aircraft/SR71-BlackBird/Models/Blackbird-SR71B.xml', #yeah i know
"ROTORCRAFT": 'Aircraft/ec130/Models/ec130b4.xml',
"GLIDER": 'Aircraft/ASK21-MI/Models/ask21mi.xml',

19
python/mlat.py
View File

@@ -55,7 +55,8 @@ wgs84_b2 = wgs84_b**2
#convert ECEF to lat/lon/alt without geoid correction
#returns alt in meters
def ecef2llh((x,y,z)):
def ecef2llh(ecef):
x, y, z = ecef
ep = math.sqrt((wgs84_a2 - wgs84_b2) / wgs84_b2)
p = math.sqrt(x**2+y**2)
th = math.atan2(wgs84_a*z, wgs84_b*p)
@@ -71,7 +72,8 @@ def ecef2llh((x,y,z)):
#convert lat/lon/alt coords to ECEF without geoid correction, WGS84 model
#remember that alt is in meters
def llh2ecef((lat, lon, alt)):
def llh2ecef(lla):
lat, lon, alt = lla
lat *= (math.pi / 180.0)
lon *= (math.pi / 180.0)
@@ -84,7 +86,8 @@ def llh2ecef((lat, lon, alt)):
return [x,y,z]
#do both of the above to get a geoid-corrected x,y,z position
def llh2geoid((lat, lon, alt)):
def llh2geoid(lla):
lat, lon, alt = lla
(x,y,z) = llh2ecef((lat, lon, alt + wgs84_height(lat, lon)))
return [x,y,z]
@@ -185,7 +188,7 @@ if __name__ == '__main__':
10 + numpy.linalg.norm(testplane-numpy.array(llh2geoid(teststations[3]))) / c,
]
print teststamps
print(teststamps)
replies = []
for i in range(0, len(teststations)):
@@ -193,7 +196,7 @@ if __name__ == '__main__':
ans = mlat(replies, testalt)
error = numpy.linalg.norm(numpy.array(llh2ecef(ans))-numpy.array(testplane))
range = numpy.linalg.norm(llh2geoid(ans)-numpy.array(testme))
print testplane-testme
print ans
print "Error: %.2fm" % (error)
print "Range: %.2fkm (from first station in list)" % (range/1000)
print(testplane-testme)
print(ans)
print("Error: %.2fm" % (error))
print("Range: %.2fkm (from first station in list)" % (range/1000))

1
python/mlat_types.py
View File

@@ -1 +0,0 @@

0
python/types.py → python/modes_types.py
View File

3
python/msprint.py
View File

@@ -20,7 +20,6 @@
#
import time, os, sys
from string import split, join
import air_modes
from air_modes.exceptions import *
import math
@@ -44,7 +43,7 @@ class output_print:
def _print(self, msg):
if self._callback is None:
print msg
print(msg)
else:
self._callback(msg)

25
python/parse.py
View File

@@ -20,8 +20,7 @@
#
import time, os, sys
from string import split, join
from altitude import decode_alt
from air_modes.altitude import decode_alt
import math
import air_modes
from air_modes.exceptions import *
@@ -32,14 +31,14 @@ class data_field:
self.data = data
self.fields = self.parse()
types = { }
dtypes = { }
offset = 1 #field offset applied to all fields. used for offsetting
#subtypes to reconcile with the spec. Really just for readability.
#get a particular field from the data
def __getitem__(self, fieldname):
mytype = self.get_type()
if mytype in self.types:
if mytype in self.dtypes:
if fieldname in self.fields: #verify it exists in this packet type
return self.fields[fieldname]
else:
@@ -52,9 +51,9 @@ class data_field:
def parse(self):
fields = {}
mytype = self.get_type()
if mytype in self.types:
for field in self.types[mytype]:
bits = self.types[self.get_type()][field]
if mytype in self.dtypes:
for field in self.dtypes[mytype]:
bits = self.dtypes[self.get_type()][field]
if len(bits) == 3:
obj = bits[2](self.get_bits(bits[0], bits[1]))
fields.update(obj.parse())
@@ -93,7 +92,7 @@ class data_field:
class bds09_reply(data_field):
offset = 6
types = { #BDS0,9 subtype 0
dtypes = { #BDS0,9 subtype 0
0: {"sub": (6,3), "dew": (10,1), "vew": (11,11), "dns": (22,1),
"vns": (23,11), "str": (34,1), "tr": (35,6), "dvr": (41,1),
"vr": (42,9)},
@@ -123,7 +122,7 @@ class bds09_reply(data_field):
class me_reply(data_field):
#types in this format are listed by BDS register
#TODO: add comments explaining these fields
types = { 0x05: {"ftc": (1,5), "ss": (6,2), "saf": (8,1), "alt": (9,12), "time": (21,1), "cpr": (22,1), "lat": (23,17), "lon": (40,17)}, #airborne position
dtypes = { 0x05: {"ftc": (1,5), "ss": (6,2), "saf": (8,1), "alt": (9,12), "time": (21,1), "cpr": (22,1), "lat": (23,17), "lon": (40,17)}, #airborne position
0x06: {"ftc": (1,5), "mvt": (6,7), "gts": (13,1), "gtk": (14,7), "time": (21,1), "cpr": (22,1), "lat": (23,17), "lon": (40,17)}, #surface position
0x07: {"ftc": (1,5),}, #TODO extended squitter status
0x08: {"ftc": (1,5), "cat": (6,3), "ident": (9,48)}, #extended squitter identification and type
@@ -157,7 +156,7 @@ class me_reply(data_field):
#resolves the TCAS reply types from TTI info
class tcas_reply(data_field):
offset = 61
types = { 0: {"tti": (61,2)}, #UNKNOWN
dtypes = { 0: {"tti": (61,2)}, #UNKNOWN
1: {"tti": (61,2), "tid": (63,26)},
2: {"tti": (61,2), "tida": (63,13), "tidr": (76,7), "tidb": (83,6)}
}
@@ -171,7 +170,7 @@ class tcas_reply(data_field):
class mb_reply(data_field):
offset = 33 #fields offset by 33 to match documentation
#types are based on bds1 subfield
types = { 0: {"bds1": (33,4), "bds2": (37,4)}, #TODO
dtypes = { 0: {"bds1": (33,4), "bds2": (37,4)}, #TODO
1: {"bds1": (33,4), "bds2": (37,4), "cfs": (41,4), "acs": (45,20), "bcs": (65,16), "ecs": (81,8)},
2: {"bds1": (33,4), "bds2": (37,4), "ais": (41,48)},
3: {"bds1": (33,4), "bds2": (37,4), "ara": (41,14), "rac": (55,4), "rat": (59,1),
@@ -195,7 +194,7 @@ class mb_reply(data_field):
class mv_reply(data_field):
offset = 33
types = { "ara": (41,14), "mte": (60,1), "rac": (55,4), "rat": (59,1),
dtypes = { "ara": (41,14), "mte": (60,1), "rac": (55,4), "rat": (59,1),
"vds": (33,8), "vds1": (33,4), "vds2": (37,4)
}
@@ -211,7 +210,7 @@ class mv_reply(data_field):
#the whole Mode S packet type
class modes_reply(data_field):
types = { 0: {"df": (1,5), "vs": (6,1), "cc": (7,1), "sl": (9,3), "ri": (14,4), "ac": (20,13), "ap": (33,24)},
dtypes = { 0: {"df": (1,5), "vs": (6,1), "cc": (7,1), "sl": (9,3), "ri": (14,4), "ac": (20,13), "ap": (33,24)},
4: {"df": (1,5), "fs": (6,3), "dr": (9,5), "um": (14,6), "ac": (20,13), "ap": (33,24)},
5: {"df": (1,5), "fs": (6,3), "dr": (9,5), "um": (14,6), "id": (20,13), "ap": (33,24)},
11: {"df": (1,5), "ca": (6,3), "aa": (9,24), "pi": (33,24)},

59
python/qa_gr-air-modes.py
View File

@@ -1,59 +0,0 @@
#!/usr/bin/env python
#
# Copyright 2004,2007 Free Software Foundation, Inc.
#
# This file is part of GNU Radio
#
# GNU Radio is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 3, or (at your option)
# any later version.
#
# GNU Radio is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with GNU Radio; see the file COPYING. If not, write to
# the Free Software Foundation, Inc., 51 Franklin Street,
# Boston, MA 02110-1301, USA.
#
from gnuradio import gr, gr_unittest
import gr-air-modes_swig
class qa_gr-air-modes (gr_unittest.TestCase):
def setUp (self):
self.tb = gr.top_block ()
def tearDown (self):
self.tb = None
def test_001_square_ff (self):
src_data = (-3, 4, -5.5, 2, 3)
expected_result = (9, 16, 30.25, 4, 9)
src = gr.vector_source_f (src_data)
sqr = gr-air-modes_swig.square_ff ()
dst = gr.vector_sink_f ()
self.tb.connect (src, sqr)
self.tb.connect (sqr, dst)
self.tb.run ()
result_data = dst.data ()
self.assertFloatTuplesAlmostEqual (expected_result, result_data, 6)
def test_002_square2_ff (self):
src_data = (-3, 4, -5.5, 2, 3)
expected_result = (9, 16, 30.25, 4, 9)
src = gr.vector_source_f (src_data)
sqr = gr-air-modes_swig.square2_ff ()
dst = gr.vector_sink_f ()
self.tb.connect (src, sqr)
self.tb.connect (sqr, dst)
self.tb.run ()
result_data = dst.data ()
self.assertFloatTuplesAlmostEqual (expected_result, result_data, 6)
if __name__ == '__main__':
gr_unittest.main ()

26
python/radio.py
View File

@@ -129,7 +129,7 @@ class modes_radio (gr.top_block, pubsub):
def set_gain(self, gain):
if self.live_source():
self._u.set_gain(gain)
print "Gain is %f" % self.get_gain()
print("Gain is %f" % self.get_gain())
return self.get_gain()
def set_rate(self, rate):
@@ -164,13 +164,19 @@ class modes_radio (gr.top_block, pubsub):
if options.source == "uhd":
#UHD source by default
from gnuradio import uhd
self._u = uhd.single_usrp_source(options.args, uhd.io_type_t.COMPLEX_FLOAT32, 1)
self._u = uhd.usrp_source(
options.args,
uhd.stream_args(
cpu_format="fc32",
channels=range(1),
),
)
if(options.subdev):
self._u.set_subdev_spec(options.subdev, 0)
if not self._u.set_center_freq(options.freq):
print "Failed to set initial frequency"
print("Failed to set initial frequency")
#check for GPSDO
#if you have a GPSDO, UHD will automatically set the timestamp to UTC time
@@ -188,9 +194,9 @@ class modes_radio (gr.top_block, pubsub):
g = self._u.get_gain_range()
options.gain = (g.start()+g.stop()) / 2.0
print "Setting gain to %i" % options.gain
print("Setting gain to %i" % options.gain)
self._u.set_gain(options.gain)
print "Gain is %i" % self._u.get_gain()
print("Gain is %i" % self._u.get_gain())
#TODO: detect if you're using an RTLSDR or Jawbreaker
#and set up accordingly.
@@ -200,13 +206,13 @@ class modes_radio (gr.top_block, pubsub):
# self._u.set_sample_rate(3.2e6) #fixed for RTL dongles
self._u.set_sample_rate(options.rate)
if not self._u.set_center_freq(options.freq):
print "Failed to set initial frequency"
print("Failed to set initial frequency")
# self._u.set_gain_mode(0) #manual gain mode
if options.gain is None:
options.gain = 34
self._u.set_gain(options.gain)
print "Gain is %i" % self._u.get_gain()
print("Gain is %i" % self._u.get_gain())
#Note: this should only come into play if using an RTLSDR.
# lpfiltcoeffs = gr.firdes.low_pass(1, 5*3.2e6, 1.6e6, 300e3)
@@ -220,12 +226,12 @@ class modes_radio (gr.top_block, pubsub):
except:
raise Exception("Please input UDP source e.g. 192.168.10.1:12345")
self._u = blocks.udp_source(gr.sizeof_gr_complex, ip, int(port))
print "Using UDP source %s:%s" % (ip, port)
print("Using UDP source %s:%s" % (ip, port))
else:
self._u = blocks.file_source(gr.sizeof_gr_complex, options.source)
print "Using file source %s" % options.source
print("Using file source %s" % options.source)
print "Rate is %i" % (options.rate,)
print("Rate is %i" % (options.rate,))
def close(self):
self.stop()

5
python/raw_server.py
View File

@@ -21,7 +21,6 @@
import time, os, sys, socket
from string import split, join
from datetime import *
class raw_server:
@@ -41,12 +40,12 @@ class raw_server:
conn.send(msg)
except socket.error:
self._conns.remove(conn)
print "Connections: ", len(self._conns)
print("Connections: ", len(self._conns))
def add_pending_conns(self):
try:
conn, addr = self._s.accept()
self._conns.append(conn)
print "Connections: ", len(self._conns)
print("Connections: ", len(self._conns))
except socket.error:
pass

1
python/rx_path.py
View File

@@ -20,7 +20,6 @@
#
from gnuradio import gr, blocks, filter
import air_modes_swig
class rx_path(gr.hier_block2):

5
python/sbs1.py
View File

@@ -21,7 +21,6 @@
import time, os, sys, socket
from string import split, join
import air_modes
import datetime
from air_modes.exceptions import *
@@ -98,7 +97,7 @@ class output_sbs1:
conn.send(sbs1_msg)
except socket.error:
self._conns.remove(conn)
print "Connections: ", len(self._conns)
print("Connections: ", len(self._conns))
except ADSBError:
pass
@@ -106,7 +105,7 @@ class output_sbs1:
try:
conn, addr = self._s.accept()
self._conns.append(conn)
print "Connections: ", len(self._conns)
print("Connections: ", len(self._conns))
except socket.error:
pass

1
python/sql.py
View File

@@ -20,7 +20,6 @@
#
import time, os, sys, threading
from string import split, join
import air_modes
import sqlite3
from air_modes.exceptions import *

6
python/zmq_socket.py
View File

@@ -27,13 +27,13 @@ import time
import threading
import zmq
from gnuradio.gr.pubsub import pubsub
import Queue
import queue
class zmq_pubsub_iface(threading.Thread):
def __init__(self, context, subaddr=None, pubaddr=None):
threading.Thread.__init__(self)
#private data
self._queue = Queue.Queue()
self._queue = queue.queue()
self._subsocket = context.socket(zmq.SUB)
self._pubsocket = context.socket(zmq.PUB)
self._subaddr = subaddr
@@ -114,7 +114,7 @@ class zmq_pubsub_iface(threading.Thread):
self.finished.wait(0.2)
def pr(x):
print x
print(x)
if __name__ == "__main__":
#create socket pair