Created
April 15, 2014 09:48
-
-
Save seanmcleod/10718607 to your computer and use it in GitHub Desktop.
Python program to render possible paths for MH-370 backwards from last ping using BFO Doppler data and ping rings.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import math | |
import matplotlib.pyplot as plt | |
from matplotlib.patches import Circle | |
satelliteInfo0011 = { 'Y': 0.0, 'Vz': 159.58, 'LOSSpeed': 250.43, 'PreviousPingTimeOffset': 91.0, 'PingRadius': 2652, 'Color': 'b', 'Elevation': 39.33 } | |
satelliteInfo2240 = { 'Y': 233.91 - 451.20, 'Vz': 123.31, 'LOSSpeed': 175.49, 'PreviousPingTimeOffset': 60.0, 'PingRadius': 2206, 'Color': 'r', 'Elevation': 47.54 } | |
satelliteInfo2140 = { 'Y': 233.91 - 557.62, 'Vz': 88.38, 'LOSSpeed': 118.34, 'PreviousPingTimeOffset': 60.0, 'PingRadius': 1965, 'Color': 'g', 'Elevation': 52.01 } | |
satelliteInfo2040 = { 'Y': 233.91 - 625.88, 'Vz': 47.42, 'LOSSpeed': 57.02, 'PreviousPingTimeOffset': 60.0, 'PingRadius': 1806, 'Color': 'c', 'Elevation': 54.98 } | |
satelliteInfo1940 = { 'Y': 233.91 - 651.30, 'Vz': 3.22, 'LOSSpeed': 4.97, 'PreviousPingTimeOffset': 78.0, 'PingRadius': 1760, 'Color': 'm', 'Elevation': 55.80 } | |
satelliteInfos = [satelliteInfo0011, satelliteInfo2240, satelliteInfo2140, satelliteInfo2040, satelliteInfo1940] | |
LOSSpeeds = [ (250.43, 125.35), (175.49, 100.64), (118.34, 79.85), (57.02, 65.80), (4.97, 39.14) ] | |
aircraftGS = 500 | |
# Setup LOS Speeds | |
LOSSpeedsSelection = 1 | |
for index in range(0, len(satelliteInfos)): | |
first, second = LOSSpeeds[index] | |
if LOSSpeedsSelection == 1: | |
satelliteInfos[index]['LOSSpeed'] = first | |
else: | |
satelliteInfos[index]['LOSSpeed'] = second | |
# Set up plot | |
fig = plt.figure(1, figsize=(8.5, 11)) | |
ax = plt.subplot(111, aspect='equal') | |
#ax.set_xlim(-3000,3000) | |
#ax.set_ylim(3000,-3000) | |
ax.set_xlim(-200,3000) | |
ax.set_ylim(3000,-1500) | |
for satelliteInfo in satelliteInfos: | |
pingCircle = Circle((0,satelliteInfo['Y']), radius=satelliteInfo['PingRadius'], fill=False) | |
ax.add_artist(pingCircle) | |
ax.scatter(0, satelliteInfo['Y'], c=satelliteInfo['Color']) | |
ax.set_title("MH370 Backtrack - %dkt - Last Doppler LOS Speed - %dkt" % (aircraftGS, satelliteInfos[0]['LOSSpeed'])) | |
ax.text(0, 2800, '00:11') | |
ax.text(0, 2150, '22:40') | |
ax.text(0, 1800, '21:40') | |
ax.text(0, 1550, '20:40') | |
ax.text(0, 1300, '19:40') | |
# Last radar position lat - 6.5381 lon - 96.408 | |
ax.scatter(1906.59, -354.75, c='y', s=50, marker='D') | |
aircraftCircle = Circle((1906.59,-354.75), radius=650, fill=False, color='y') | |
ax.add_artist(aircraftCircle) | |
ax.annotate('Northern 19:40 Arc', xy=(1668,-958), xytext=(700,-1100), arrowprops=dict(facecolor='black')) | |
ax.annotate('Southern 19:40 Arc', xy=(1634,235), xytext=(600,100), arrowprops=dict(facecolor='black')) | |
latitudes = [(35, ' '), (333, 'S5'), (632, 'S10'), (931, 'S15'), (1223, 'S20'), (1528, 'S25'), (1828, 'S30'), (2127, 'S35'), (-263, 'N5'), (-562, 'N10'), (-861, 'N15')] | |
for latitude in latitudes: | |
yCart, latText = latitude | |
ax.plot((-3000,3000), (yCart,yCart), color='grey', alpha=0.5) | |
ax.text(2750, yCart, latText) | |
def distance(x1, y1, x2, y2): | |
dx = x1 - x2 | |
dy = y1 - y2 | |
return math.sqrt(dx * dx + dy * dy) | |
def calc(aircraftX, aircraftY, satelliteIndex): | |
if satelliteIndex == len(satelliteInfos) - 1: | |
return | |
# Unpack satellite info | |
satelliteInfo = satelliteInfos[satelliteIndex] | |
satelliteVz = satelliteInfo['Vz'] | |
satelliteElevation = satelliteInfo['Elevation'] | |
satelliteY = satelliteInfo['Y'] | |
LOSSpeed = satelliteInfo['LOSSpeed'] | |
backtrackTime = satelliteInfo['PreviousPingTimeOffset'] | |
distanceFromSatelliteToAircraft = distance(0.0, satelliteY, aircraftX, aircraftY) | |
bearingFromSatelliteToAircaft = math.asin(aircraftX/distanceFromSatelliteToAircraft) | |
aircraftDistanceCoveredBetweenLastTwoPings = (backtrackTime / 60.0) * aircraftGS | |
satelliteCone = satelliteVz * math.cos(math.radians(satelliteElevation)) | |
aircraftCone = aircraftGS * math.cos(math.radians(satelliteElevation)) | |
satelliteLOS = satelliteCone * math.cos(bearingFromSatelliteToAircaft) | |
try: | |
aircraftBearingAngleDifference = math.degrees(math.acos((LOSSpeed + satelliteLOS) / aircraftCone)) | |
except ValueError: | |
return | |
aircraftLOS = aircraftCone * math.cos(math.radians(aircraftBearingAngleDifference)) | |
finalLOS = aircraftLOS - satelliteLOS | |
bearingPositive = 180 - math.degrees(bearingFromSatelliteToAircaft) + aircraftBearingAngleDifference | |
bearingNegative = 180 - math.degrees(bearingFromSatelliteToAircaft) - aircraftBearingAngleDifference | |
matches = [] | |
for bearing in [bearingPositive, bearingNegative]: | |
previousAircraftX = -9999 | |
previousAircraftY = -9999 | |
if bearing > 90.0 and bearing <= 180.0: | |
previousAircraftX = aircraftX - aircraftDistanceCoveredBetweenLastTwoPings * math.sin(math.radians(180.0 - bearing)) | |
previousAircraftY = aircraftY - aircraftDistanceCoveredBetweenLastTwoPings * math.cos(math.radians(180.0 - bearing)) | |
elif bearing > 180.0 and bearing <= 270.0: | |
previousAircraftX = aircraftX + aircraftDistanceCoveredBetweenLastTwoPings * math.sin(math.radians(bearing - 180.0)) | |
previousAircraftY = aircraftY - aircraftDistanceCoveredBetweenLastTwoPings * math.cos(math.radians(bearing - 180.0)) | |
if previousAircraftX != -9999 and previousAircraftY != -9999: | |
distanceFromPreviousSatelitePoint = distance(0.0, satelliteInfos[satelliteIndex + 1]['Y'], previousAircraftX, previousAircraftY) | |
rangeRingDiff = math.fabs(distanceFromPreviousSatelitePoint - satelliteInfos[satelliteIndex + 1]['PingRadius']) | |
if rangeRingDiff < (0.02 + (satelliteIndex * 0.5)/100.0) * satelliteInfos[satelliteIndex + 1]['PingRadius']: | |
if satelliteIndex == 0: | |
ax.scatter(aircraftX, aircraftY, c=satelliteInfo['Color']) | |
ax.scatter(previousAircraftX, previousAircraftY, c=satelliteInfos[satelliteIndex + 1]['Color']) | |
ax.plot([aircraftX, previousAircraftX], [aircraftY, previousAircraftY]) | |
calc(previousAircraftX, previousAircraftY, satelliteIndex + 1) | |
return | |
for bearing in range(1, 89, 1): | |
pingRadius = satelliteInfos[0]['PingRadius'] | |
aircraftX = pingRadius * math.sin(math.radians(bearing)) | |
aircraftY = pingRadius * math.cos(math.radians(bearing)) | |
calc(aircraftX, aircraftY, 0) | |
plt.savefig("MH370 Backtrack - %dkt - Last Doppler LOS Speed - %dkt.png" % (aircraftGS, satelliteInfos[0]['LOSSpeed'])) | |
plt.show() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment