Initial commit

2010-09-14 22:01:56 -07:00
commit e39021ff8a
189 changed files with 21343 additions and 0 deletions
--- a/src/python/Makefile.am
+++ b/src/python/Makefile.am
@@ -0,0 +1,32 @@
+#
+# Copyright 2004 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.
+# 
+
+include $(top_srcdir)/Makefile.common
+
+EXTRA_DIST = run_tests.in
+
+
+TESTS = 				\
+	run_tests
+
+
+noinst_PYTHON = 			\
+	qa_howto.py			
--- a/src/python/altitude.py
+++ b/src/python/altitude.py
@@ -0,0 +1,71 @@
+#!/usr/bin/env python
+#from string import split, join
+
+#betcha this would be faster if you used a table for mode C
+#you could strip out D1 since it's never used, that leaves 11 bits (table is 2048 entries)
+def decode_alt(alt, bit13):
+	if alt & 0x40 and bit13 is True:
+		return "METRIC ERROR"
+
+	if alt & 0x10: #a mode S-style reply
+		if bit13 is True:
+			tmp1 = (alt & 0x1F80) >> 2 #first 6 bits get shifted 2 down
+			tmp2 = (alt & 0x20) >> 1 #that bit gets shifted 1 down
+		else:
+			tmp1 = (alt & 0x0FE0) >> 1 #first 7 bits get shifted 1 down but there are only 12 bits in the representation
+			tmp2 = 0
+
+		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
+
+		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
+
+		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
+
+		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 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
+
+	return decoded_alt
+
+def gray2bin(gray):
+	i = gray >> 1
+
+	while i != 0:
+		gray ^= i
+		i >>= 1
+
+	return gray
--- a/src/python/altitude.pyc
+++ b/src/python/altitude.pyc
--- a/src/python/cpr-test.py
+++ b/src/python/cpr-test.py
@@ -0,0 +1,38 @@
+#!/usr/bin/env python
+
+from modes_parse import *
+from cpr import *
+import sys
+
+my_location = [37.76225, -122.44254]
+
+shortdata = long(sys.argv[1], 16)
+longdata = long(sys.argv[2], 16)
+parity = long(sys.argv[3], 16)
+ecc = long(sys.argv[4], 16)
+
+[altitude, decoded_lat, decoded_lon, rnge, bearing] = parseBDS05(shortdata, longdata, parity, ecc)
+
+if decoded_lat is not None:
+	print "Altitude: %i\nLatitude: %.6f\nLongitude: %.6f\nRange: %.2f\nBearing: %i\n" % (altitude, decoded_lat, decoded_lon, rnge, bearing,)
+
+print "Decomposing...\n"
+
+subtype = (longdata >> 51) & 0x1F;
+
+encoded_lon = longdata & 0x1FFFF
+encoded_lat = (longdata >> 17) & 0x1FFFF
+cpr_format = (longdata >> 34) & 1
+
+enc_alt = (longdata >> 36) & 0x0FFF
+
+print "Subtype: %i\nEncoded longitude: %x\nEncoded latitude: %x\nCPR format: %i\nEncoded altitude: %x\n" % (subtype, encoded_lon, encoded_lat, cpr_format, enc_alt,)
+
+#print "First argument is order %i, second %i" % ((evendata >> 34) & 1, (odddata >> 34) & 1,)
+
+#evenencpos = [(evendata >> 17) & 0x1FFFF, evendata & 0x1FFFF]
+#oddencpos = [(odddata >> 17) & 0x1FFFF, odddata & 0x1FFFF]
+
+#[declat, declon] = cpr_decode_global(evenencpos, oddencpos, newer)
+
+#print "Global latitude: %.6f\nGlobal longitude: %.6f" % (declat, declon,)
--- a/src/python/cpr.py
+++ b/src/python/cpr.py
@@ -0,0 +1,246 @@
+#!/usr/bin/env python
+#from string import split, join
+#from math import pi, floor, cos, acos
+import math, time
+#this implements CPR position decoding. local only for now.
+
+latz = 15
+nbits = 17
+my_lat = 37.76225 #update these later!
+my_lon = -122.44254
+
+
+def nz(ctype):
+	return 4 * latz - ctype
+
+def dlat(ctype, surface):
+	if surface == 1:
+		tmp = 90.0
+	else:
+		tmp = 360.0
+
+	nzcalc = nz(ctype)
+	if nzcalc == 0:
+		return tmp
+	else:
+		return tmp / nzcalc
+
+def nl_eo(declat_in, ctype):
+	return nl(declat_in) - ctype
+
+def nl(declat_in):
+	return math.floor( (2.0*math.pi) * pow(math.acos(1.0- (1.0-math.cos(math.pi/(2.0*latz))) / pow( math.cos( (math.pi/180.0)*abs(declat_in) ) ,2.0) ),-1.0))
+
+def dlon(declat_in, ctype, surface):
+	if surface == 1:
+		tmp = 90.0
+	else:
+		tmp = 360.0
+	nlcalc = nl_eo(declat_in, ctype)
+	if nlcalc == 0:
+		return tmp
+	else:
+		return tmp / nlcalc
+
+def decode_lat(enclat, ctype, my_lat, surface):
+	tmp1 = dlat(ctype, surface)
+	tmp2 = float(enclat) / (2**nbits)
+	j = math.floor(my_lat/tmp1) + math.floor(0.5 + (mod(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)
+	m = math.floor(my_lon / tmp1) + math.floor(0.5 + (mod(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)
+
+
+def mod(a, b):
+	if a < 0:
+		a += 360.0
+
+	return a - b * math.floor(a / b)
+
+def cpr_resolve_local(my_location, encoded_location, ctype, surface):
+	[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)
+
+	return [decoded_lat, decoded_lon]
+
+def cpr_resolve_global(evenpos, oddpos, mostrecent, surface): #ok this is considered working, tentatively
+	dlateven = dlat(0, surface);
+	dlatodd  = dlat(1, surface);
+
+#	print dlateven;
+#	print dlatodd;
+
+	evenpos = [float(evenpos[0]), float(evenpos[1])]
+	oddpos = [float(oddpos[0]), float(oddpos[1])]
+
+	#print "Even position: %x, %x\nOdd position: %x, %x" % (evenpos[0], evenpos[1], oddpos[0], oddpos[1],)
+
+	j = math.floor(((59*evenpos[0] - 60*oddpos[0])/2**nbits) + 0.5) #latitude index
+
+	#print "Latitude index: %i" % j #should be 6, getting 5?
+
+	rlateven = dlateven * (mod(j, 60)+evenpos[0]/2**nbits)
+	rlatodd  = dlatodd  * (mod(j, 59)+ oddpos[0]/2**nbits)
+
+	#print "Rlateven: %f\nRlatodd: %f" % (rlateven, rlatodd,)
+
+	if nl(rlateven) != nl(rlatodd):
+		#print "Boundary straddle!"
+		return (None, None,)
+
+	if mostrecent == 0:
+		rlat = rlateven
+	else:
+		rlat = rlatodd
+
+	if rlat > 90:
+		rlat = rlat - 180.0
+
+	dl = dlon(rlat, mostrecent, surface)
+	nlthing = nl(rlat)
+	ni = nlthing - mostrecent
+
+	#print "ni is %i" % ni
+
+	m =  math.floor(((evenpos[1]*(nlthing-1)-oddpos[1]*(nlthing))/2**nbits)+0.5) #longitude index, THIS LINE IS CORRECT
+	#print "m is %f" % m #should be -16
+
+
+	if mostrecent == 0:
+		enclon = evenpos[1]
+	else:
+		enclon = oddpos[1]
+
+	rlon = dl * (mod(ni+m, ni)+enclon/2**nbits)
+
+	if rlon > 180:
+		rlon = rlon - 360.0
+
+	return [rlat, rlon]
+
+
+def cpr_decode(icao24, encoded_lat, encoded_lon, cpr_format, evenlist, oddlist, lkplist, surface, longdata):
+	#this is a stopgap measure to catch those packets which aren't really position packets. what gives?
+#	if encoded_lat == 0 or encoded_lon == 0:
+		#print "debug: lat or lon zero for longdata %x" % (longdata,)
+#		return [None, None, None, None]
+
+	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]
+
+	#okay, let's traverse the lists and weed out those entries that are older than 15 minutes, as they're unlikely to be useful.
+	for key, item in lkplist.items():
+		if time.time() - item[2] > 900:
+			del lkplist[key]
+
+	for key, item in evenlist.items():
+		if time.time() - item[2] > 900:
+			del evenlist[key]
+
+	for key, item in oddlist.items():
+		if time.time() - item[2] > 900:
+			del oddlist[key]
+
+	#here we perform global/emitter-centered CPR decoding as follows:
+	#first, check for the ICAO number in the planelist. if there is a decoded position in there, use that for emitter-centered decoding and be done with it.
+
+	if surface==1:
+		validrange = 45
+	else:
+		validrange = 180
+
+	if icao24 in lkplist:
+		#print "debug: icao found in LKP table. EC decoding with local position list %s" % str(lkplist[icao24][0:2])
+		[decoded_lat, decoded_lon] = cpr_resolve_local(lkplist[icao24][0:2], [encoded_lat, encoded_lon], cpr_format, surface) #do emitter-centered local decoding
+		lkplist[icao24] = [decoded_lat, decoded_lon, time.time()] #update the local position for next time
+
+############debug info for plotting strange position reports###############
+#		[lkprange, lkpbearing] = range_bearing(lkplist[icao24][0:2], [decoded_lat, decoded_lon])
+#		lkpdeltat = time.time() - lkplist[icao24][2]
+#		#the units are now mi/sec
+#		#an SR-71 can move at 0.6 miles per second, so let's say if it's over 1.0mi/s it's probably a bug
+#		if lkprange / lkpdeltat > 1.0:
+#			print "debug: buggy position packet detected from icao %x, encoded lat %x, encoded lon %x, CPR format %i, longdata %x." % (icao24, encoded_lat, encoded_lon, cpr_format, longdata)
+#		
+############debug info for plotting strange position reports###############
+
+	else: #no LKP available
+		#print "debug: icao %x not found. attempting local decode." % icao24
+		[local_lat, local_lon] = cpr_resolve_local([my_lat, my_lon], [encoded_lat, encoded_lon], cpr_format, surface) #try local decoding
+#		print "debug: local resolve gives %.6f, %.6f" % (local_lat, local_lon)
+		[rnge, bearing] = range_bearing([my_lat, my_lon], [local_lat, local_lon])
+		if rnge < validrange: #if the local decoding can be guaranteed valid
+			#print "debug: range < 180nm, position valid."
+			lkplist[icao24] = [local_lat, local_lon, time.time()] #update the local position for next time
+			[decoded_lat, decoded_lon] = [local_lat, local_lon]
+		else: #if the local decoding can't be guaranteed valid AND you couldn't find an LKP
+#			print "debug: range > %inm, attempting global decode." % validrange
+			#attempt global decode
+			if (icao24 in evenlist) and (icao24 in oddlist):
+#				print "debug: ICAOs found in both lists."
+				if abs(evenlist[icao24][2] - oddlist[icao24][2]) < 10: #if there's less than 10 seconds of time difference between the reports
+#					print "debug: valid even/odd positions, performing global decode."
+					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], newer, surface) #do a global decode
+					if decoded_lat is not None:
+						lkplist[icao24] = [decoded_lat, decoded_lon, time.time()]
+#				else:
+#					print "debug: timestamps not close enough to be valid."
+#			else:
+#					print "debug: even/odd information not found."
+
+	#print "settled on position: %.6f, %.6f" % (decoded_lat, decoded_lon,)
+	if decoded_lat is not None:
+		[rnge, bearing] = range_bearing([my_lat, my_lon], [decoded_lat, decoded_lon])
+	else:
+		rnge = None
+		bearing = None
+
+	return [decoded_lat, decoded_lon, rnge, bearing]
+
+
+
+def range_bearing(loc_a, 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)
+
+	delta_lat = b_lat - a_lat
+	delta_lon = b_lon - a_lon
+
+	avg_lat = (a_lat + b_lat) / 2.0
+
+	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))
+
+	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
+
+	rnge = math.hypot(distance_East,distance_North)
+
+
+	return [rnge, bearing]
--- a/src/python/cpr.pyc
+++ b/src/python/cpr.pyc
--- a/src/python/modes_parse.py
+++ b/src/python/modes_parse.py
@@ -0,0 +1,247 @@
+#!/usr/bin/env python
+
+import time, os, sys
+from string import split, join
+from altitude import decode_alt
+from cpr import cpr_decode
+import math
+
+def parse0(shortdata, parity, ecc):
+#	shortdata = long(shortdata, 16)
+	#parity = long(parity)
+
+	vs = bool(shortdata >> 26 & 0x1)    #ground sensor -- airborne when 0
+	cc = bool(shortdata >> 25 & 0x1)    #crosslink capability, binary
+	sl = shortdata >> 21 & 0x07  #operating sensitivity of onboard TCAS system. 0 means no TCAS sensitivity reported, 1-7 give TCAS sensitivity
+	ri = shortdata >> 15 & 0x0F #speed coding: 0 = no onboard TCAS, 1 = NA, 2 = TCAS w/inhib res, 3 = TCAS w/vert only, 4 = TCAS w/vert+horiz, 5-7 = NA, 8 = no max A/S avail,
+								  #9 = A/S <= 75kt, 10 = A/S (75-150]kt, 11 = (150-300]kt, 12 = (300-600]kt, 13 = (600-1200]kt, 14 = >1200kt, 15 = NA
+
+	altitude = decode_alt(shortdata & 0x1FFF, True) #bit 13 is set for type 0
+
+	return [vs, cc, sl, ri, altitude]
+
+def parse4(shortdata, parity, ecc):
+#	shortdata = long(shortdata, 16)
+	fs = shortdata >> 24 & 0x07 #flight status: 0 is airborne normal, 1 is ground normal, 2 is airborne alert, 3 is ground alert, 4 is alert SPI, 5 is normal SPI
+	dr = shortdata >> 19 & 0x1F #downlink request: 0 means no req, bit 0 is Comm-B msg rdy bit, bit 1 is TCAS info msg rdy, bit 2 is Comm-B bcast #1 msg rdy, bit2+bit0 is Comm-B bcast #2 msg rdy,
+								#bit2+bit1 is TCAS info and Comm-B bcast #1 msg rdy, bit2+bit1+bit0 is TCAS info and Comm-B bcast #2 msg rdy, 8-15 N/A, 16-31 req to send N-15 segments
+	um = shortdata >> 13 & 0x3F #transponder status readouts, no decoding information available
+
+	altitude = decode_alt(shortdata & 0x1FFF, True)
+
+	return [fs, dr, um, altitude]
+
+
+
+def parse5(shortdata, parity, ecc):
+#	shortdata = long(shortdata, 16)
+	fs = shortdata >> 24 & 0x07 #flight status: 0 is airborne normal, 1 is ground normal, 2 is airborne alert, 3 is ground alert, 4 is alert SPI, 5 is normal SPI
+	dr = shortdata >> 19 & 0x1F #downlink request: 0 means no req, bit 0 is Comm-B msg rdy bit, bit 1 is TCAS info msg rdy, bit 2 is Comm-B bcast #1 msg rdy, bit2+bit0 is Comm-B bcast #2 msg rdy,
+								#bit2+bit1 is TCAS info and Comm-B bcast #1 msg rdy, bit2+bit1+bit0 is TCAS info and Comm-B bcast #2 msg rdy, 8-15 N/A, 16-31 req to send N-15 segments
+	um = shortdata >> 13 & 0x3F #transponder status readouts, no decoding information available
+
+	return [fs, dr, um]
+
+def parse11(shortdata, parity, ecc):
+#	shortdata = long(shortdata, 16)
+	interrogator = ecc & 0x0F
+	
+	ca = shortdata >> 13 & 0x3F #capability
+	icao24 = shortdata & 0xFFFFFF
+
+	return [icao24, interrogator, ca]
+
+#def parse17(shortdata, longdata, parity, ecc):
+#	shortdata = long(shortdata, 16)
+#	longdata = long(longdata, 16)
+#	parity = long(parity, 16)
+#	ecc = long(ecc, 16)
+
+#	subtype = (longdata >> 51) & 0x1F;
+
+	#the subtypes are:
+	#0: No position information
+	#1: Identification (Category set D)
+	#2: Identification (Category set C)
+	#3: "" (B)
+	#4: "" (A)
+	#5: Surface position accurate to 7.5m
+	#6: "" to 25m
+	#7: "" to 185.2m (0.1nm)
+	#8: "" above 185.2m
+	#9: Airborne position to 7.5m
+	#10-18: Same with less accuracy
+	#19: Airborne velocity
+	#20: Airborne position w/GNSS height above earth
+	#21: same to 25m
+	#22: same above 25m
+	#23: Reserved
+	#24: Reserved for surface system status
+	#25-27: Reserved
+	#28: Extended squitter aircraft status
+	#29: Current/next trajectory change point
+	#30: Aircraft operational coordination
+	#31: Aircraft operational status
+
+
+#	if subtype == 4:
+#		retstr = parseBDS08(shortdata, longdata, parity, ecc)
+#	elif subtype >= 9 and subtype <= 18:
+#		retstr = parseBDS05(shortdata, longdata, parity, ecc)
+#	elif subtype == 19:
+#		subsubtype = (longdata >> 48) & 0x07
+#		if subsubtype == 0:
+#			retstr = parseBDS09_0(shortdata, longdata, parity, ecc)
+#		elif subsubtype == 1:
+#			retstr = parseBDS09_1(shortdata, longdata, parity, ecc)
+#		else:
+#			retstr = "BDS09 subtype " + str(subsubtype) + " not implemented"
+#	else:
+#		retstr = "Type 17, subtype " + str(subtype) + " not implemented"
+
+#	return retstr
+
+def parseBDS08(shortdata, longdata, parity, ecc):
+	icao24 = shortdata & 0xFFFFFF
+
+	msg = ""
+	for i in range(0, 8):
+		msg += charmap( longdata >> (42-6*i) & 0x3F)
+
+	#retstr = "Type 17 subtype 04 (ident) from " + "%x" % icao24 + " with data " + msg
+
+	return msg
+
+def charmap(d):
+	if d > 0 and d < 27:
+		retval = chr(ord("A")+d-1)
+	elif d == 32:
+		retval = " "
+	elif d > 47 and d < 58:
+		retval = chr(ord("0")+d-48)
+	else:
+		retval = " "
+
+	return retval
+
+
+#lkplist is the last known position, for emitter-centered decoding. evenlist and oddlist are the last 
+#received encoded position data for each reporting type. all dictionaries indexed by ICAO number.
+lkplist = {}
+evenlist = {}
+oddlist = {}
+evenlist_ground = {}
+oddlist_ground = {}
+
+#the above dictionaries are all in the format [lat, lon, time].
+
+def parseBDS05(shortdata, longdata, parity, ecc):
+	icao24 = shortdata & 0xFFFFFF
+
+	encoded_lon = longdata & 0x1FFFF
+	encoded_lat = (longdata >> 17) & 0x1FFFF
+	cpr_format = (longdata >> 34) & 1
+
+	enc_alt = (longdata >> 36) & 0x0FFF
+
+	altitude = decode_alt(enc_alt, False)
+
+	[decoded_lat, decoded_lon, rnge, bearing] = cpr_decode(icao24, encoded_lat, encoded_lon, cpr_format, evenlist, oddlist, lkplist, 0, longdata)
+
+	return [altitude, decoded_lat, decoded_lon, rnge, bearing]
+
+
+#welp turns out it looks like there's only 17 bits in the BDS0,6 ground packet after all. fuck.
+def parseBDS06(shortdata, longdata, parity, ecc):
+	icao24 = shortdata & 0xFFFFFF
+
+	encoded_lon = longdata & 0x1FFFF
+	encoded_lat = (longdata >> 17) & 0x1FFFF
+	cpr_format = (longdata >> 34) & 1
+
+#	enc_alt = (longdata >> 36) & 0x0FFF
+
+	altitude = 0
+
+	[decoded_lat, decoded_lon, rnge, bearing] = cpr_decode(icao24, encoded_lat, encoded_lon, cpr_format, evenlist_ground, oddlist_ground, lkplist, 1, longdata)
+
+	return [altitude, decoded_lat, decoded_lon, rnge, bearing]
+
+
+def parseBDS09_0(shortdata, longdata, parity, ecc):
+	icao24 = shortdata & 0xFFFFFF
+	vert_spd = ((longdata >> 6) & 0x1FF) * 32
+	ud = bool((longdata >> 15) & 1)
+	if ud:
+		vert_spd = 0 - vert_spd
+	turn_rate = (longdata >> 16) & 0x3F
+	turn_rate = turn_rate * 15/62
+	rl = bool((longdata >> 22) & 1)
+	if rl:
+		turn_rate = 0 - turn_rate
+	ns_vel = (longdata >> 23) & 0x7FF - 1
+	ns = bool((longdata >> 34) & 1)
+	ew_vel = (longdata >> 35) & 0x7FF - 1
+	ew = bool((longdata >> 46) & 1)
+	subtype = (longdata >> 48) & 0x07
+
+	velocity = math.hypot(ns_vel, ew_vel)
+	if ew:
+		ew_vel = 0 - ew_vel
+	if ns:
+		ns_vel = 0 - ns_vel
+	heading = math.atan2(ew_vel, ns_vel) * (180.0 / math.pi)
+	if heading < 0:
+		heading += 360
+
+	#retstr = "Type 17 subtype 09-0 (track report) from " + "%x" % icao24 + " with velocity " + "%.0f" % velocity + "kt heading " + "%.0f" % heading + " VS " + "%.0f" % vert_spd
+
+	return [velocity, heading, vert_spd]
+
+def parseBDS09_1(shortdata, longdata, parity, ecc):
+	icao24 = shortdata & 0xFFFFFF
+	alt_geo_diff = longdata & 0x7F - 1
+	above_below = bool((longdata >> 7) & 1)
+	if above_below:
+		alt_geo_diff = 0 - alt_geo_diff;
+	vert_spd = float((longdata >> 10) & 0x1FF - 1)
+	ud = bool((longdata >> 19) & 1)
+	if ud:
+		vert_spd = 0 - vert_spd
+	vert_src = bool((longdata >> 20) & 1)
+	ns_vel = float((longdata >> 21) & 0x3FF - 1)
+	ns = bool((longdata >> 31) & 1)
+	ew_vel = float((longdata >> 32) & 0x3FF - 1)
+	ew = bool((longdata >> 42) & 1)
+	subtype = (longdata >> 48) & 0x07
+
+
+	if subtype == 0x02:
+		ns_vel *= 4
+		ew_vel *= 4
+
+
+	vert_spd *= 64
+	alt_geo_diff *= 25
+	
+	velocity = math.hypot(ns_vel, ew_vel)
+	if ew:
+		ew_vel = 0 - ew_vel
+	
+	if ns_vel == 0:
+		heading = 0
+	else:
+		heading = math.atan(float(ew_vel) / float(ns_vel)) * (180.0 / math.pi)
+	if ns:
+		heading = 180 - heading
+	if heading < 0:
+		heading += 360
+
+	#retstr = "Type 17 subtype 09-1 (track report) from " + "%x" % icao24 + " with velocity " + "%.0f" % velocity + "kt heading " + "%.0f" % heading + " VS " + "%.0f" % vert_spd
+
+	return [velocity, heading, vert_spd]
+
+def parse20(shortdata, longdata, parity, ecc):
+	return "Message 20 not yet implemented"
+
+
--- a/src/python/modes_parse.pyc
+++ b/src/python/modes_parse.pyc
--- a/src/python/modes_print.py
+++ b/src/python/modes_print.py
@@ -0,0 +1,142 @@
+#!/usr/bin/env python
+import time, os, sys
+from string import split, join
+from modes_parse import *
+
+def modes_print(message):
+
+	#a mode S parser for all message types. first, split the input into data fields
+	[msgtype, shortdata, longdata, parity, ecc, reference] = message.split()
+
+	shortdata = long(shortdata, 16)
+	longdata = long(longdata, 16)
+	parity = long(parity, 16)
+	ecc = long(ecc, 16)
+	reference = float(reference)
+
+	msgtype = int(msgtype)
+
+	output = str("")
+
+	if msgtype == 0:
+		output = print0(shortdata, parity, ecc)
+	elif msgtype == 4:
+		output = print4(shortdata, parity, ecc)
+	elif msgtype == 5:
+		output = print5(shortdata, parity, ecc)
+	elif msgtype == 11:
+		output = print11(shortdata, parity, ecc)
+	elif msgtype == 17:
+		output = print17(shortdata, longdata, parity, ecc)
+	elif msgtype == 20:
+		output = parse20(shortdata, longdata, parity, ecc)
+	else:
+		output = "No handler for message type " + str(msgtype) + " from " + str(ecc)
+
+	output = "(%.0f) " % float(reference) + output
+
+	return output
+
+def print0(shortdata, parity, ecc):
+	[vs, cc, sl, ri, altitude] = parse0(shortdata, parity, ecc)
+	
+	retstr = "Type 0 (short A-A surveillance) from " + "%x" % ecc + " at " + str(altitude) + "ft"
+	# the ri values below 9 are used for other things. might want to print those someday.
+	if ri == 9:
+		retstr = retstr + " (speed <75kt)"
+	elif ri > 9:
+		retstr = retstr + " (speed " + str(75 * (1 << (ri-10))) + "-" + str(75 * (1 << (ri-9))) + "kt)"
+
+	if vs is True:
+		retstr = retstr + " (aircraft is on the ground)"
+
+	return retstr
+
+def print4(shortdata, parity, ecc):
+
+	[fs, dr, um, altitude] = parse4(shortdata, parity, ecc)
+
+	retstr = "Type 4 (short surveillance altitude reply) from " + "%x" % ecc + " at " + str(altitude) + "ft"
+
+	if fs == 1:
+		retstr = retstr + " (aircraft is on the ground)"
+	elif fs == 2:
+		retstr = retstr + " (AIRBORNE ALERT)"
+	elif fs == 3:
+		retstr = retstr + " (GROUND ALERT)"
+	elif fs == 4:
+		retstr = retstr + " (SPI ALERT)"
+	elif fs == 5:
+		retstr = retstr + " (SPI)"
+
+	return retstr
+
+def print5(shortdata, parity, ecc):
+	[fs, dr, um] = parse5(shortdata, parity, ecc)
+
+	retstr = "Type 5 (short surveillance ident reply) from " + "%x" % ecc + " with ident " + str(shortdata & 0x1FFF)
+
+	if fs == 1:
+		retstr = retstr + " (aircraft is on the ground)"
+	elif fs == 2:
+		retstr = retstr + " (AIRBORNE ALERT)"
+	elif fs == 3:
+		retstr = retstr + " (GROUND ALERT)"
+	elif fs == 4:
+		retstr = retstr + " (SPI ALERT)"
+	elif fs == 5:
+		retstr = retstr + " (SPI)"
+
+	return retstr
+
+def print11(shortdata, parity, ecc):
+	[icao24, interrogator, ca] = parse11(shortdata, parity, ecc)
+
+	retstr = "Type 11 (all call reply) from " + "%x" % icao24 + " in reply to interrogator " + str(interrogator)
+	return retstr
+
+def print17(shortdata, longdata, parity, ecc):
+	icao24 = shortdata & 0xFFFFFF
+	subtype = (longdata >> 51) & 0x1F;
+
+	if subtype == 4:
+		msg = parseBDS08(shortdata, longdata, parity, ecc)
+		retstr = "Type 17 subtype 04 (ident) from " + "%x" % icao24 + " with data " + msg
+
+	elif subtype >= 5 and subtype <= 8:
+		[altitude, decoded_lat, decoded_lon, rnge, bearing] = parseBDS06(shortdata, longdata, parity, ecc)
+		if decoded_lat==0: #no unambiguously valid position available
+			retstr = ""
+		else:
+			#retstr = "INSERT INTO plane_positions (icao, seen, alt, lat, lon) VALUES ('" + "%x" % icao24 + "', now(), " + str(altitude) + ", " + "%.6f" % decoded_lat + ", " + "%.6f" % decoded_lon + ")"
+			retstr = "Type 17 subtype 06 (surface report) from " + "%x" % icao24 + " at (" + "%.6f" % decoded_lat + ", " + "%.6f" % decoded_lon + ") (" + "%.2f" % rnge + " @ " + "%.0f" % bearing + ")"
+
+	elif subtype >= 9 and subtype <= 18:
+		[altitude, decoded_lat, decoded_lon, rnge, bearing] = parseBDS05(shortdata, longdata, parity, ecc)
+
+		retstr = "Type 17 subtype 05 (position report) from " + "%x" % icao24 + " at (" + "%.6f" % decoded_lat + ", " + "%.6f" % decoded_lon + ") (" + "%.2f" % rnge + " @ " + "%.0f" % bearing + ") at " + str(altitude) + "ft"
+
+#	this is a trigger to capture the bizarre BDS0,5 squitters you keep seeing on the map with latitudes all over the place
+#		if icao24 == 0xa1ede9:
+#			print "Buggy squitter with shortdata %s longdata %s parity %s ecc %s" % (str(shortdata), str(longdata), str(parity), str(ecc),)
+
+	elif subtype == 19:
+		subsubtype = (longdata >> 48) & 0x07
+		if subsubtype == 0:
+			[velocity, heading, vert_spd] = parseBDS09_0(shortdata, longdata, parity, ecc)
+			retstr = "Type 17 subtype 09-0 (track report) from " + "%x" % icao24 + " with velocity " + "%.0f" % velocity + "kt heading " + "%.0f" % heading + " VS " + "%.0f" % vert_spd
+
+		elif subsubtype == 1:
+			[velocity, heading, vert_spd] = parseBDS09_1(shortdata, longdata, parity, ecc)
+			retstr = "Type 17 subtype 09-1 (track report) from " + "%x" % icao24 + " with velocity " + "%.0f" % velocity + "kt heading " + "%.0f" % heading + " VS " + "%.0f" % vert_spd
+
+		else:
+			retstr = "BDS09 subtype " + str(subsubtype) + " not implemented"
+	else:
+		retstr = "Type 17, subtype " + str(subtype) + " not implemented"
+
+	return retstr
+
+
+
+
--- a/src/python/modes_sql.py
+++ b/src/python/modes_sql.py
@@ -0,0 +1,85 @@
+#!/usr/bin/env python
+import time, os, sys
+from string import split, join
+from modes_parse import *
+
+def modes_sql(message):
+
+	#assembles a MySQLdb query tailored to Owen's database
+	#this version ignores anything that isn't Type 17 for now, because we just don't care
+
+	[msgtype, shortdata, longdata, parity, ecc, reference] = message.split()
+
+	shortdata = long(shortdata, 16)
+	longdata = long(longdata, 16)
+	parity = long(parity, 16)
+	ecc = long(ecc, 16)
+#	reference = float(reference)
+
+	msgtype = int(msgtype)
+
+	query = None
+
+#	if msgtype == 0:
+#		query = sql0(shortdata, parity, ecc)
+#	elif msgtype == 4:
+#		query = sql4(shortdata, parity, ecc)
+#	elif msgtype == 5:
+#		query = sql5(shortdata, parity, ecc)
+#	elif msgtype == 11:
+#		query = sql11(shortdata, parity, ecc)
+#	elif msgtype == 17:
+	if msgtype == 17:
+		query = sql17(shortdata, longdata, parity, ecc)
+#	elif msgtype == 20:
+#		output = parse20(shortdata, longdata, parity, ecc)
+#	else:
+		#output = "No handler for message type " + str(msgtype) + " from " + str(ecc)
+
+#	output = "(%.0f) " % float(reference) + output
+
+	return query
+
+def sql17(shortdata, longdata, parity, ecc):
+	icao24 = shortdata & 0xFFFFFF	
+	subtype = (longdata >> 51) & 0x1F
+
+	retstr = None
+
+	if subtype == 4:
+		msg = parseBDS08(shortdata, longdata, parity, ecc)
+		retstr = "INSERT INTO plane_metadata (icao, ident) VALUES ('" + "%x" % icao24 + "', '" + msg + "') ON DUPLICATE KEY UPDATE seen=now(), ident=values(ident)"
+
+	elif subtype >= 5 and subtype <= 8:
+		[altitude, decoded_lat, decoded_lon, rnge, bearing] = parseBDS06(shortdata, longdata, parity, ecc)
+		if decoded_lat is None: #no unambiguously valid position available
+			retstr = None
+		else:
+			retstr = "INSERT INTO plane_positions (icao, seen, alt, lat, lon) VALUES ('" + "%x" % icao24 + "', now(), " + str(altitude) + ", " + "%.6f" % decoded_lat + ", " + "%.6f" % decoded_lon + ")"
+
+	elif subtype >= 9 and subtype <= 18 and subtype != 15: #i'm eliminating type 15 records because they don't appear to be valid position reports.
+		[altitude, decoded_lat, decoded_lon, rnge, bearing] = parseBDS05(shortdata, longdata, parity, ecc)
+		if decoded_lat is None: #no unambiguously valid position available
+			retstr = None
+		else:
+			retstr = "INSERT INTO plane_positions (icao, seen, alt, lat, lon) VALUES ('" + "%x" % icao24 + "', now(), " + str(altitude) + ", " + "%.6f" % decoded_lat + ", " + "%.6f" % decoded_lon + ")"
+
+	elif subtype == 19:
+		subsubtype = (longdata >> 48) & 0x07
+		if subsubtype == 0:
+			[velocity, heading, vert_spd] = parseBDS09_0(shortdata, longdata, parity, ecc)
+			retstr = "INSERT INTO plane_metadata (icao, seen, speed, heading, vertical) VALUES ('" + "%x" % icao24 + "', now(), " + "%.0f" % velocity + ", " + "%.0f" % heading + ", " + "%.0f" % vert_spd + ") ON DUPLICATE KEY UPDATE seen=now(), speed=values(speed), heading=values(heading), vertical=values(vertical)"
+
+		elif subsubtype == 1:
+			[velocity, heading, vert_spd] = parseBDS09_1(shortdata, longdata, parity, ecc)
+			retstr = "INSERT INTO plane_metadata (icao, seen, speed, heading, vertical) VALUES ('" + "%x" % icao24 + "', now(), " + "%.0f" % velocity + ", " + "%.0f" % heading + ", " + "%.0f" % vert_spd + ") ON DUPLICATE KEY UPDATE seen=now(), speed=values(speed), heading=values(heading), vertical=values(vertical)"
+
+	else:
+		print "debug (modes_sql): unknown subtype %i with data %x %x %x" % (subtype, shortdata, longdata, parity,)
+
+
+	return retstr
+
+
+
+
--- a/src/python/usrp_modes.py
+++ b/src/python/usrp_modes.py
@@ -0,0 +1,183 @@
+#!/usr/bin/env python
+
+from gnuradio import gr, gru, usrp, optfir, eng_notation, blks2, air
+from gnuradio.eng_option import eng_option
+from optparse import OptionParser
+import time, os, sys
+from string import split, join
+from usrpm import usrp_dbid
+from modes_print import modes_print
+from modes_sql import modes_sql
+import gnuradio.gr.gr_threading as _threading
+import MySQLdb
+
+
+class top_block_runner(_threading.Thread):
+    def __init__(self, tb):
+        _threading.Thread.__init__(self)
+        self.setDaemon(1)
+        self.tb = tb
+        self.done = False
+        self.start()
+
+    def run(self):
+        self.tb.run()
+        self.done = True
+
+
+
+def pick_subdevice(u):
+#this should pick which USRP subdevice if none was specified on the command line
+#since the only thing that will receive ADS-B appears to be the DBSRX, we'll default to that
+
+	return usrp.pick_subdev(u, (usrp_dbid.DBS_RX,
+				    usrp_dbid.DBS_RX_REV_2_1,
+				    usrp_dbid.BASIC_RX))
+"""
+
+The following are optional command line parameters:
+
+-R SUBDEV    Daughter board specification, defaults to first found
+-f FREQ      USRP receive frequency (1090 MHz Default)
+-g GAIN      Daughterboard gain setting. Defaults to mid-range.
+-d DECIM     USRP decimation rate
+-t THRESH    Receiver valid pulse threshold
+-a           Output all frames. Defaults only output frames
+
+Once the program is running, ctrl-break (Ctrl-C) stops operation.
+"""
+
+class adsb_rx_block (gr.top_block):
+
+	def __init__(self, options, args, queue):
+		gr.top_block.__init__(self)
+
+		self.options = options
+		self.args = args
+
+		if options.filename is None:
+			if options.decim < 8:
+				self.fpga_filename="std_4rx_0tx.rbf" #can go down to decim 4
+				self.u = usrp.source_c(fpga_filename=self.fpga_filename)
+			else :
+				self.u = usrp.source_c()
+
+			if options.rx_subdev_spec is None:
+				options.rx_subdev_spec = pick_subdevice(self.u)
+
+			self.u.set_mux(usrp.determine_rx_mux_value(self.u, options.rx_subdev_spec))
+			self.subdev = usrp.selected_subdev(self.u, options.rx_subdev_spec)
+			print "Using RX d'board %s" % (self.subdev.side_and_name(),)
+			self.u.set_decim_rate(options.decim)
+
+			if options.gain is None: #set to halfway
+				g = self.subdev.gain_range()
+				options.gain = (g[0]+g[1]) / 2.0
+
+			if not(self.tune(options.freq)):
+				print "Failed to set initial frequency"
+
+			print "Setting gain to %i" % (options.gain,)
+			self.subdev.set_gain(options.gain)
+			self.subdev.set_bw(self.options.bandwidth) #only for DBSRX
+
+			rate = self.u.adc_rate() / options.decim
+
+		else:
+			rate = int(64e6 / options.decim)
+			self.u = gr.file_source(gr.sizeof_gr_complex, options.filename)
+
+		print "Rate is %i" % (rate,)
+		pass_all = 0
+		if options.output_all :
+			pass_all = 1
+
+
+
+		self.demod = gr.complex_to_mag()
+		self.avg = gr.moving_average_ff(100, 1.0/100, 400);
+		self.preamble = air.modes_preamble(rate, options.threshold)
+		self.framer = air.modes_framer(rate)
+		self.slicer = air.modes_slicer(rate, queue)
+
+		if options.decim < 16:
+				#there's a really nasty spur at 1088 caused by a multiple of the USRP xtal. if you use a decimation of 16, it gets filtered out by the CIC. if not, it really fucks with you unless you filter it out.
+				filter_coeffs = gr.firdes.band_reject(1.0, rate, 1.7e6, 2.3e6, 0.5e6, gr.firdes.WIN_HAMMING)
+				self.filt = gr.fir_filter_ccf(1, filter_coeffs)
+				self.connect(self.u, self.filt)
+		else:
+				self.filt = self.u
+
+		self.connect(self.filt, self.demod)
+		self.connect(self.demod, self.avg)
+		self.connect(self.demod, (self.preamble, 0))
+		self.connect(self.avg, (self.preamble, 1))
+		self.connect(self.demod, (self.framer, 0))
+		self.connect(self.preamble, (self.framer, 1))
+		self.connect(self.demod, (self.slicer, 0))
+		self.connect(self.framer, (self.slicer, 1))
+
+	def tune(self, freq):
+		result = usrp.tune(self.u, 0, self.subdev, freq)
+		return True
+
+
+if __name__ == '__main__':
+	usage = "%prog: [options] output filename"
+	parser = OptionParser(option_class=eng_option, usage=usage)
+	parser.add_option("-R", "--rx-subdev-spec", type="subdev",
+		      help="select USRP Rx side A or B", metavar="SUBDEV")
+	parser.add_option("-f", "--freq", type="eng_float", default=1090e6,
+                      help="set receive frequency in Hz [default=%default]", metavar="FREQ")
+	parser.add_option("-g", "--gain", type="int", default=None,
+                      help="set RF gain", metavar="dB")
+	parser.add_option("-d", "--decim", type="int", default=16,
+                      help="set fgpa decimation rate [default=%default]")
+	parser.add_option("-T", "--threshold", type="eng_float", default=3.0,
+                      help="set pulse detection threshold above noise in dB [default=%default]")
+	parser.add_option("-a","--output-all", action="store_true", default=False,
+                      help="output all frames")
+	parser.add_option("-b","--bandwidth", type="eng_float", default=5e6,
+		      help="set DBSRX front-end bandwidth in Hz [default=5e6]")
+	parser.add_option("-F","--filename", type="string", default=None,
+					  help="read data from file instead of USRP")
+	parser.add_option("-D","--database", action="store_true", default=False,
+											help="send to database instead of printing to screen")
+	(options, args) = parser.parse_args()
+#	if len(args) != 1:
+#		parser.print_help()
+#		sys.exit(1)
+
+#	filename = args[0]
+
+	queue = gr.msg_queue()
+
+	if options.database is True:
+		db = MySQLdb.connect(host="localhost", user="planes", passwd="planes", db="planes")
+
+	fg = adsb_rx_block(options, args, queue)
+	runner = top_block_runner(fg)
+
+	while 1:
+		try:
+			if queue.empty_p() == 0 :
+				while queue.empty_p() == 0 :
+					msg = queue.delete_head() #blocking read
+					if options.database is False:					
+						print modes_print(msg.to_string())
+					else:
+						query = modes_sql(msg.to_string())
+						if query is not None:
+							c = db.cursor()
+							c.execute(query)
+
+			elif runner.done:
+				break
+			else:
+				time.sleep(0.1)
+
+		except KeyboardInterrupt:
+			fg.stop()
+			runner = None
+			break
+