Notebook

More on Modeling the Last Flight of MH370 with Monte Carlo¶

Conor L. Myhrvold¶

Harvard University, SEAS, IACS , Computational Science & Engineering (CSE)

AM207 : Advanced Scientific Computing: Stochastic Optimization Methods. Monte Carlo Methods for Inference and Data Analysis

Final Project, Spring 2014. Current as of March 31.

For the first version of the model, with runways, see: IPython Notebook Link ¶

**For articles about this first work, see Fast Company Labs' :
How I Narrowed Down The Location Of Malaysia Air Using "Monte Carlo" Data Models

*More About Our Methodology: Tracking MH370 With Monte Carlo Data Models***
¶

Contact Info:

conor.myhrvold@gmail.com

Twitter: @conormyhrvold

LinkedIn: www.linkedin.com/pub/conor-myhrvold/37/583/a7/

Fast Company Labs Contributor Page

Update Explanation: I'm updating from version one, incorporating new info released since I worked on the model. Since the goal is to make this version standalone as well, I've inserted the new material seamlessly without calling too much attention to what came first. I've also eliminated parts which are no longer relevant. For example, I took out plotting all of the runways in the region that a 777 could land on; while cool, the most likely scenario for MH370 now, according to my model and Inmarsat's analysis, is that it crashed and is at the bottom of the Indian Ocean. And, I fixed a few bugs.¶

Goal: The overall goal of this project is to use publically available information to pinpoint possible locations of flight MH370, a Malaysian Airlines Boeing 777-200ER which went missing in flight on March 8, 2014.¶

Let's get started. Here are the libraries we need. Note that after installing the Anaconda distribution of Python 2.7, which has IPython, you need to install seaborn and Basemap, which is part of matplotlib but is not included in the default installation. You need Basemap to make the maps -- there's really no other way around it.¶

In [5]:

#using Anaconda 64 bit distro of Python 2.7, Windows 7 environment (other specs will work though)
#import libraries

#numpy
import numpy as np

#scipy
import scipy as sp
from scipy import stats
from scipy.stats import norm
from scipy.stats import expon

#plotting
%matplotlib inline 
# ^^ make plots appear in IPython notebook instead of separate window
import matplotlib.pyplot as plt
import seaborn as sns
from IPython import display
from mpl_toolkits.basemap import Basemap # use for plotting lat lon & world maps

Overview - MH370¶

Malaysia Airlines 370 disappeared. A variety of information sources have been obtained, which can clue us into its possible locations. I have read a ton of articles and also reached out to some of the journalists (which I have a background in) and organizations involved, in an effort to get the best data possible to put into our Monte Carlo (MC) simulation. In particular, the New York Times Hong Kong Bureau has been of great assistance to me. Why Monte Carlo? MC is an ideal approach to model where the airplane could be since unlike other air disasters, we don't know a lot about what happened yet. Several weeks have passed and no definitive wreckage has been found, let alone the airplane crash site or how it got there. This is despite a determined search effort.¶

Getting the Data - Flight MH370 Known Flight Coordinates + Inmarsat Coords¶

Here are some useful constants we need for converting purposes:¶

In [1]:

eq_deg_km = 111.32 # number of km/degree at eq Source: http://en.wikipedia.org/wiki/Decimal_degrees
earth_radius = 6371 #in km, http://en.wikipedia.org/wiki/Earth

#Inmarsat satellite information
sat_height = 42170 #Inmarsat satellite height off ground, in meters
elevation_angle = np.radians(40) #elevation angle of satellite; convert degrees to radians. Source: NYT Hong Kong Bureau

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-1-509c27f0aa21> in <module>()
      4 #Inmarsat satellite information
      5 sat_height = 42170 #Inmarsat satellite height off ground, in meters
----> 6 elevation_angle = np.radians(40) #elevation angle of satellite; convert degrees to radians. Source: NYT Hong Kong Bureau

NameError: name 'np' is not defined

Now we'll input the coordinates relevant to us:¶

In [2]:

# The Inmarsat satellite is at 0,64.5 -- it's geostationary.
inmarsat = [0, 64.5]

#Now we plot the plane's known positions

#Kuala Lumpur International Airport Coordinates: http://www.distancesfrom.com/my/Kuala-Lumpur-Airport-(KUL)-Malaysia-latitude-longitude-Kuala-Lumpur-Airport-(KUL)-Malaysia-latitude-/LatLongHistory/3308940.aspx
kualalumpur = [2.7544829, 101.7011363]

#Pulau Perak coordinates: http://tools.wmflabs.org/geohack/geohack.php?pagename=Pulau_Perak&params=5_40_50_N_98_56_27_E_type:isle_region:MY
# http://en.wikipedia.org/wiki/Perak_Island -> Indonesia military radar detected near island
pulauperak = [5.680556, 98.940833]

# Igari Waypoint. Source: # http://www.fallingrain.com/waypoint/SN/IGARI.html Given in hours,minutes,seconds.
igariwaypoint = [6. + 56./60. + 12./3600., 103. + 35./60. + 6./3600.] 

print "inmarsat lat/lon:", inmarsat[0], inmarsat[1]
print "kuala lumpur int'l airport coord lat/lon:", kualalumpur[0],kualalumpur[1]
print "pulau perak lat/lon:", pulauperak[0],pulauperak[1]
print "igari waypoint lat/lon:", igariwaypoint[0],igariwaypoint[1]

inmarsat lat/lon: 0 64.5
kuala lumpur int'l airport coord lat/lon: 2.7544829 101.7011363
pulau perak lat/lon: 5.680556 98.940833
igari waypoint lat/lon: 6.93666666667 103.585

Note that we don't use the Kuala Lumpur airport coordinates in our actual simulations. We do use the Igari Waypoint coordinate, however, in conjunction with the Pulau Perak coordinates to generate our initial heading. (Which was done in a Mathematica notebook, separate to this file, but could be easily replicable given that I provide the coordinates for both.) By changing the distribution see how sensitive this initial heading will be to our final results. However, I'll point out that irrespective of the sensitivity, this is our best guess of where the plane was last headed.¶

Plotting the Data¶

Before we plot the data, let's formally define our latitude and longitude bounds for our map, up front. That way, we can more readily change them across the document without manually editing one by one. These will double duty to make our probabilities as well:¶

In [3]:

#create lat/long grid of distance

lat_min = -50 #50 S
lat_max = 50  #50 N
lon_min = 50  #50 E
lon_max = 140 #130 E
lat_space = 5 #spacing for plotting latitudes and longitudes
lon_space = 5

Now we'll make a map with those points on it, using Basemap, an extension of Matplotlib:¶

In [4]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=13000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

# plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='lower left',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('Initial Map & Data', fontsize=30)

#Show below
plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-4-6a350768797f> in <module>()
      1 #Set figure size
----> 2 fig = plt.figure(figsize=[30,20])
      3 
      4 #Setup Basemap
      5 fig = Basemap(width=10000000,height=13000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

NameError: name 'plt' is not defined

This map is the basis for more complicated maps showing the plane's possible position at a given point in time, along with its flight path. The middle of the South Indian Ocean in this map is basically the middle of nowhere, but from our results this region of the South Indian Ocean, especially further down, is a very interesting area to search.¶

Satellite Ping Data -- the 5th Ping¶

In [5]:

#To get the satellite calculate we have to solve several equations which were computed on Mathematica
"""
Computes the ping arc distance from the satellite to the plane.
Returns the angle in degrees.
"""
def satellite_calc(radius,orbit,angle):
    interim = (np.sqrt(-radius**2+orbit**2*(1./np.cos(angle)**2))-orbit*np.tan(angle))/np.float(orbit+radius)
    return np.degrees(2*np.arctan(interim))

ping_arc_dist = satellite_calc(earth_radius,sat_height,elevation_angle)
print "ping arc distance in degrees:", ping_arc_dist

dist_from_sat = earth_radius*np.radians(satellite_calc(earth_radius,sat_height,elevation_angle))
print "distance from satellite", dist_from_sat

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-5-7bd537eaffc2> in <module>()
      8     return np.degrees(2*np.arctan(interim))
      9 
---> 10 ping_arc_dist = satellite_calc(earth_radius,sat_height,elevation_angle)
     11 print "ping arc distance in degrees:", ping_arc_dist
     12 

NameError: name 'elevation_angle' is not defined

There are two types of approximations we can make, for a given radius of the geostationary satellite or the airplane. The first is a simple circle, which is pretty good. However, Earth is not a sphere, and so at large radii this approximation is not ideal. So I write a more complicated function to reflect the ellipse-like shape the arc will really be like. I don't include the spherical trigonometry / mathematical derivations (done in a Mathematica notebook), since they are not in the scope of finding MH370.¶

In [6]:

"""
write circle function. generates x&y pairs in a circle, 360 degrees
angle in degrees, number of points to put in circle, satellite location lat/lon
"""
def make_circle1(radius,pts,lon_loc,lat_loc): 
    coords_array = np.zeros((pts, 2))
    for i in xrange(pts):
        coords_array[i][0] =  radius*np.cos(np.radians(i)) + lat_loc 
        coords_array[i][1] =  radius*np.sin(np.radians(i)) + lon_loc 
    return coords_array

"""
write ellipse function
since across the earth it's not actually a circle
derived from spherical trigonometry
"""
def make_circle(radius,pts,lon_loc,lat_loc):
    coords_array = np.zeros((pts, 2))
    for i in xrange(pts):
        coords_array[i][0] =  radius*np.cos(np.radians(i)) + lat_loc 
        coords_array[i][1] =  radius*np.sin(np.radians(i))/(np.cos(np.cos(np.radians(radius))*(lat_loc+radius*np.cos(np.radians(i)))*(np.pi/180.))) + lon_loc 
    return coords_array

Now we test the two approximations to better understand how they compare to each other:¶

In [7]:

#test near the equator, -5 (5 S)

test = make_circle(8.12971613367,360,90,-5)

test_lat = []
test_lon = []

for i in xrange(len(test)):
    test_lat.append(test[i][0])
    test_lon.append(test[i][1])

test1 = make_circle1(8.12971613367,360,90,-5)

test_lat1 = []
test_lon1 = []

for i in xrange(len(test1)):
    test_lat1.append(test1[i][0])
    test_lon1.append(test1[i][1])

#create figure
plt.figure()
plt.plot(test_lon1,test_lat1,color='red',label='simple circle')
plt.plot(test_lon,test_lat,color='blue',label='ellipsoid')
legend = plt.legend(loc='lower left',fontsize=8,frameon=True,markerscale=1,prop={'size':5})
plt.title('comparing circle approximations to each other')
plt.legend()
plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-7-2077abb64802> in <module>()
      1 #test near the equator, -5 (5 S)
      2 
----> 3 test = make_circle(8.12971613367,360,90,-5)
      4 
      5 test_lat = []

<ipython-input-6-17c869eba246> in make_circle(radius, pts, lon_loc, lat_loc)
     16 """
     17 def make_circle(radius,pts,lon_loc,lat_loc):
---> 18     coords_array = np.zeros((pts, 2))
     19     for i in xrange(pts):
     20         coords_array[i][0] =  radius*np.cos(np.radians(i)) + lat_loc

NameError: global name 'np' is not defined

In [8]:

#test at a low latitude, -40 (40 S)

test = make_circle(8.12971613367,360,90,-40)

test_lat = []
test_lon = []

for i in xrange(len(test)):
    test_lat.append(test[i][0])
    test_lon.append(test[i][1])

test1 = make_circle1(8.12971613367,360,90,-40)

test_lat1 = []
test_lon1 = []

for i in xrange(len(test1)):
    test_lat1.append(test1[i][0])
    test_lon1.append(test1[i][1])

#create figure
plt.figure()
plt.plot(test_lon1,test_lat1,color='red',label='simple circle')
plt.plot(test_lon,test_lat,color='blue',label='ellipsoid')
legend = plt.legend(loc='lower left',fontsize=8,frameon=True,markerscale=1,prop={'size':5})
plt.title('comparing circle approximations to each other')
plt.legend()
plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-8-21b10bbc9ff0> in <module>()
      1 #test at a low latitude, -40 (40 S)
      2 
----> 3 test = make_circle(8.12971613367,360,90,-40)
      4 
      5 test_lat = []

<ipython-input-6-17c869eba246> in make_circle(radius, pts, lon_loc, lat_loc)
     16 """
     17 def make_circle(radius,pts,lon_loc,lat_loc):
---> 18     coords_array = np.zeros((pts, 2))
     19     for i in xrange(pts):
     20         coords_array[i][0] =  radius*np.cos(np.radians(i)) + lat_loc

NameError: global name 'np' is not defined

Above we see the difference between the two circles. Remember that at the equator, 1 degree of latitude and longitude stay the same, but at the poles, the degrees of latitude contract tremendously. The ellipse we have created is slightly better at the equator, but much better at lower latitudes. And that's where it matters most to find our plane, if our first version is any guide (either one -- even the one with bugs, which basically had the lower latitude plot, near the equator, which caused that lower latitude plot to be out of line for a small fraction of the make_vector dots.)¶

In [9]:

circle_pts = make_circle(ping_arc_dist,360,64.5,0)
#print circle_pts

circle_lat = []
for i in xrange(len(circle_pts)):
    circle_lat.append(circle_pts[i][0])

circle_lon = []
for i in xrange(len(circle_pts)):
    circle_lon.append(circle_pts[i][1])

print circle_lat[0:10]
print "-------------------"
print circle_lon[0:10]

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-9-936e41cc3e84> in <module>()
----> 1 circle_pts = make_circle(ping_arc_dist,360,64.5,0)
      2 #print circle_pts
      3 
      4 circle_lat = []
      5 for i in xrange(len(circle_pts)):

NameError: name 'ping_arc_dist' is not defined

Incorporating New Info on Ping Distance Error¶

Right after my model was done, on March 24, excellent analysis from the TMF Associates article Understanding the "satellite ping" conclusion... was published, pointed out that the ping distance error would be 1 or 2 degrees. Essentially, the biggest uncertainty for the ping arc is the lag time, which is apparently on the order of 300 microseconds, which round trip is 600 microseconds according to the linked calculation (recall that the plane "saw" the satellite at a 40 degree angle, which on the graph is 600 microseconds.) It's 1 or 2 degrees because it hasn't been established whether the measurement by the satellite is done one-way or round trip.¶

In [10]:

print "Percent error 1 way is:", (1./40)*100
print "Percent error round trip is:", (2./40)*100

Percent error 1 way is: 2.5
Percent error round trip is: 5.0

Thus, I will use 2.5 and 5% error. (Fortuitously, I'd picked 5% before! So with the advent of this extra information and analysis, that number is spot on, for one ping. Which also means 10 and 20% were reasonable if you wanted to cover an area that included multiple pings, instead of just using the 5th ping.) Now what we'll do is vary the distance from the satellite (radius) by 5% and 2.5% to build in error in where the actually was (the geostationary satellite is assumed to essentially have no error in its position, relative to the error in knowing where the plane is).¶

Inmarsat 5% error¶

In [11]:

err1_5per = 0.95*ping_arc_dist
err2_5per = 1.05*ping_arc_dist

circle_pts_err1_5per = make_circle(err1_5per,360,64.5,0)
circle_pts_err2_5per = make_circle(err2_5per,360,64.5,0)

circle_lon_err1_5per = []
circle_lon_err2_5per = []
circle_lat_err1_5per = []
circle_lat_err2_5per = []

for i in xrange(len(circle_pts_err1_5per)):
    circle_lon_err1_5per.append(circle_pts_err1_5per[i][1])
    
for i in xrange(len(circle_pts_err2_5per)):
    circle_lon_err2_5per.append(circle_pts_err2_5per[i][1])
    
for i in xrange(len(circle_pts_err1_5per)):
    circle_lat_err1_5per.append(circle_pts_err1_5per[i][0])
    
for i in xrange(len(circle_pts_err2_5per)):
    circle_lat_err2_5per.append(circle_pts_err2_5per[i][0])

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-11-7dfca8693360> in <module>()
----> 1 err1_5per = 0.95*ping_arc_dist
      2 err2_5per = 1.05*ping_arc_dist
      3 
      4 circle_pts_err1_5per = make_circle(err1_5per,360,64.5,0)
      5 circle_pts_err2_5per = make_circle(err2_5per,360,64.5,0)

NameError: name 'ping_arc_dist' is not defined

Inmarsat 2.5% error¶

In [12]:

err1_2per = 0.975*ping_arc_dist
err2_2per = 1.025*ping_arc_dist

circle_pts_err1_2per = make_circle(err1_2per,360,64.5,0)
circle_pts_err2_2per = make_circle(err2_2per,360,64.5,0)

circle_lon_err1_2per = []
circle_lon_err2_2per = []
circle_lat_err1_2per = []
circle_lat_err2_2per = []

for i in xrange(len(circle_pts_err1_2per)):
    circle_lon_err1_2per.append(circle_pts_err1_2per[i][1])
    
for i in xrange(len(circle_pts_err2_2per)):
    circle_lon_err2_2per.append(circle_pts_err2_2per[i][1])
    
for i in xrange(len(circle_pts_err1_2per)):
    circle_lat_err1_2per.append(circle_pts_err1_2per[i][0])
    
for i in xrange(len(circle_pts_err2_2per)):
    circle_lat_err2_2per.append(circle_pts_err2_2per[i][0])

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-12-9b1f4fce56fa> in <module>()
----> 1 err1_2per = 0.975*ping_arc_dist
      2 err2_2per = 1.025*ping_arc_dist
      3 
      4 circle_pts_err1_2per = make_circle(err1_2per,360,64.5,0)
      5 circle_pts_err2_2per = make_circle(err2_2per,360,64.5,0)

NameError: name 'ping_arc_dist' is not defined

Now we make a map plot showing these circles. The dashed lines away from the solid red line is the error, 2.5%.¶

In [13]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=13000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- 2.5% error
x5,y5 = fig(circle_lon,circle_lat)
x6,y6 = fig(circle_lon_err1_2per,circle_lat_err1_2per)
x7,y7 = fig(circle_lon_err2_2per,circle_lat_err2_2per)
x8,y8 = fig(circle_lon_err1_5per,circle_lat_err1_5per)
x9,y9 = fig(circle_lon_err2_5per,circle_lat_err2_5per)

#Draw circle showing extent of Inmarsat sat radar detection
fig.plot(x5,y5,'r-',markersize=5,label='Inferred MH370 Location from Satellite, 5th Ping')
fig.plot(x6,y6,'r--',markersize=5,label='with 2.5 & 5% error')
fig.plot(x7,y7,'r--',markersize=5)
fig.plot(x8,y8,'r--',markersize=5)
fig.plot(x9,y9,'r--',markersize=5)

# plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='lower left',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('Inmarsat Ping Estimation', fontsize=30)

#Show below
plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-13-eaaee86c608e> in <module>()
      1 #Set figure size
----> 2 fig = plt.figure(figsize=[30,20])
      3 
      4 #Setup Basemap
      5 fig = Basemap(width=10000000,height=13000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

NameError: name 'plt' is not defined

Now we have the starting plot that we saw in so many newspaper websites all around the world in March. We have gathered an impressive data set which forms the foundation of the second version of our Monte Carlo model.

Making the Monte Carlo Model¶

Let's keep our Haversine formula from our airplane runways section of the first version of this model; we need this for computing distances below. Apparently, there's a Python module haversine which has this, but let's write this function ourselves:¶

In [14]:

"""
Haversine equation.
Computes the great circle distance between two pairs of longitude/latitude.
Returns the distance in m or km depending on input (I use meters.) 
"""
def haversine(r,lat1,lon1,lat2,lon2):
    d = 2.0*r*np.arcsin(np.sqrt(np.sin(np.radians(lat2-lat1)/2.0)**2+np.cos(np.radians(lat1))*np.cos(np.radians(lat2))*np.sin(np.radians(lon2-lon1)/2.0)**2))
    return d   

Convert the Inmarsat line location and error bounds into probabilities¶

In [15]:

"""
ping_prob_normal is specifically the 5th ping probability, which we have.
we center a normal probability distribution upon the location of the radius line.

d = dist_from_sat in my case and is more abstractly 'd', a distance
r = earth_radius in my case and is is more abstractly 'r', a radius
lat1 = sat latitude in my case
lon1 = sat longitude in my case
lat2, lon2 we iterate through for our function
err is the 5,10,or 20% error Inmarsat error
"""
def ping_prob_normal(lat1,lon1,lat2,lon2,err,d,r): 
    return np.exp(-0.5*((haversine(r,lat1,lon1,lat2,lon2)-d)/(err*d))**2)/(d*np.sqrt(2*np.pi)) 

Next we make a probability grid of where the plane could be located, across our latitudes and longitudes (which we defined at the top of the document). We'll initialize this probability from the 5th ping perspective, since we've already calculated appropriate latitudes and longitudes from that. We'll keep this at a normal distribution, since the standard deviation as a percentange of the radius is small (std dev = 2.5 and 5% of the radius.)¶

In [16]:

lat_range = np.abs(lat_max)+np.abs(lat_min)
lon_range = np.abs(lon_max)-np.abs(lon_min)
print "latitude range in degrees:", lat_range
print "longitude range in degrees:", lon_range

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-16-5ca86b89f636> in <module>()
----> 1 lat_range = np.abs(lat_max)+np.abs(lat_min)
      2 lon_range = np.abs(lon_max)-np.abs(lon_min)
      3 print "latitude range in degrees:", lat_range
      4 print "longitude range in degrees:", lon_range

NameError: name 'np' is not defined

In [17]:

# x coord = longitude
# y coord = latitude
mult_fac = 10 #Use to make a finer grid, so 10 per 1 degree lat/lon , and fill with probabilities.

prob_grid = np.zeros((lat_range*mult_fac,lon_range*mult_fac)) #initialize grid as numpy array

for i in xrange(lat_range*mult_fac): #across all lat grid-lets
    for j in xrange(lon_range*mult_fac): #across all lon grid-lets (so now we're iterating through rows + columns)
        prob_grid[i][j] = ping_prob_normal(inmarsat[0],inmarsat[1],(lat_min+i/np.float(mult_fac)),(lon_min+j/np.float(mult_fac)),0.05,dist_from_sat,earth_radius) 
        #assuming 5% error ... 5% == 0.05

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-17-f3e9ed508a74> in <module>()
      3 mult_fac = 10 #Use to make a finer grid, so 10 per 1 degree lat/lon , and fill with probabilities.
      4 
----> 5 prob_grid = np.zeros((lat_range*mult_fac,lon_range*mult_fac)) #initialize grid as numpy array
      6 
      7 for i in xrange(lat_range*mult_fac): #across all lat grid-lets

NameError: name 'np' is not defined

In [18]:

plt.figure()
plt.plot(prob_grid)
plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-18-ddd942a63958> in <module>()
----> 1 plt.figure()
      2 plt.plot(prob_grid)
      3 plt.show()

NameError: name 'plt' is not defined

Note that I'm now making a grid with a multiplication factor of 10 instead of 1. This fancier option was not exercised in the first version of the model, as I just kept it to 1 for simplicity's sake (and to show that the "coarseness" of the probability grid did not have a large effect on results -- results which initially didn't have independent corroboration, but that right around the time I finished, Inmarsat released their analysis showing the plane to be near where I said it would be. So they looked reasonable to me...and then soon afterward an independent analysis -- the official one, I might add -- confirmed what I came up with, during more or less the same time interval.)¶

In [19]:

#test out with just longitude at 0 latitude
prob_grid_lon =  np.zeros((lon_range*mult_fac))

for j in xrange(lon_range*mult_fac): #across all lon grid-lets (so now we're iterating through rows + columns)
    prob_grid_lon[j] = ping_prob_normal(inmarsat[0],inmarsat[1],0,(lon_min+j/np.float(mult_fac)),0.05,dist_from_sat,earth_radius)

plt.figure()
plt.plot(prob_grid_lon)
plt.title('longitude ping probability across the equator (0 deg lat)')
plt.xlabel('array index -- to get longitude need to add 500, which is minimum longitude')
plt.ylabel('probability')
plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-19-03cd57230985> in <module>()
      1 #test out with just longitude at 0 latitude
----> 2 prob_grid_lon =  np.zeros((lon_range*mult_fac))
      3 
      4 for j in xrange(lon_range*mult_fac): #across all lon grid-lets (so now we're iterating through rows + columns)
      5     prob_grid_lon[j] = ping_prob_normal(inmarsat[0],inmarsat[1],0,(lon_min+j/np.float(mult_fac)),0.05,dist_from_sat,earth_radius)

NameError: name 'np' is not defined

Determining MH370's Final Location -- (up to) The Last Hour¶

After this, I write a function that will make same-heading vectors from the 5th ping, to whereever the plane may finally be located. Note that the largely-offsetting bugs from the first version have been fixed!¶

In [20]:

"""
write function that plots final destination from plane, from the point of the last ping.
this will be some period of time in between 0 minutes and an hour -- or else it would have pinged again.
make a point at a distance on a heading to see where the plane would be if it continued on a straight line,
from the 5th ping.
"""
def make_vector(ang_radius,heading,lon_loc,lat_loc):
    vec_lat = ang_radius*np.cos(np.radians(heading)) + lat_loc
    vec_lon = ang_radius*np.sin(np.radians(heading))/(np.cos(np.cos(np.radians(ang_radius))*(lat_loc+ang_radius*np.cos(np.radians(heading)))*(np.pi/180.))) + lon_loc
    return vec_lat,vec_lon

Initial Location Uncertainty and Randomness¶

In [21]:

"""
takes in a heading, a starting location, and a std dev for picking a new heading
a std dev of 30, for example, with a normal distribution, means that -30 to 30 degrees will be
where most of the new heading draws will come out of (a std dev applies to both sides of the distribution.) 
and the highest probability will be straight ahead. this is of course not including the ping probability.
"""
def normal_prob_step(old_heading, std_dev, start_lon, start_lat, ang_dist):
    
    angle_list = range(0,360) # ; print angle_list
    
    #we create a radius of 360 points that would be the heading for 360 possible degrees
    circle = make_circle(ang_dist,len(angle_list),start_lon,start_lat)
    
    weights = np.zeros(len(angle_list)) #make 360 array
    for i in xrange(len(weights)):
        weights[i] = np.exp(-(((np.float(angle_list[i])-180.)/std_dev)**2) /2.) / (std_dev*np.sqrt(2*np.pi)) 
        #makes array of normally distributed weights as if heading was 180 degrees. Sort of a hack to make it periodic.  
    
    #Now we have to translate it back to whatever the heading should be, instead of 180 degrees
    #Fortunately we can use numpy's roll. Implementing our own would be a pain.
    s_weights = np.roll(weights,-180+np.int(old_heading))
    
    #initialize new possible coordinates within an hr's distance, new weights for the odds, and new angles
    #(depending on whether the plane would go off the grid or not)
    new_circle = []
    new_weights = []
    new_angles = []
    
    #make sure lat & lon are in bounds
    for i in xrange(len(circle)):
        if circle[i][0] >= lat_min and circle[i][0] <= lat_max and circle[i][1] >= lon_min and circle[i][1] <= lon_max:
            new_circle.append(circle[i])       
            new_weights.append(s_weights[i])
            new_angles.append(angle_list[i])

    return new_circle,new_weights,new_angles

We incorporate the above function our 5 step model -- 1 for each hour, since that's when the pings were received. Afterward, the plane went down sometime in an hr, so I plot where the final contact spot could be at 1/4,1/2,3/4, and 59/60 hrs.¶

In [22]:

"""
a function which given a list of discrete probabilities for each destination point, 
it will choose one of those points.

heading_init -- initial direction was headed at last known point
lon_init,lat_init -- last known point of plane in longitude and latitude
km_hop -- how far the plane went in the time interval, 1 hr. So in simplest case, the 777's cruising speed/hr.
std_dev -- the standard deviation of the heading, based on a normal distribution from the current heading (0 deg).
ping_percent_err -- what % error you assume in the Inmarsat 5th ping. either 2.5 or 5%.

uses normal distribution for heading

"""
def five_hop_model_normal(heading_init,lon_init,lat_init,km_hop,std_dev,ping_percent_err):   
    
    #initialize
    plane_lat = np.zeros(5) #initialize plane location after each hop (assumed to be 1 hr for now)
    plane_lon = np.zeros(5)  
    lat = lat_init
    lon = lon_init
    heading = heading_init
        
    for i in xrange(len(plane_lat)):
        new_circle,new_weights,new_angles = normal_prob_step(heading,std_dev,lon,lat,km_hop/eq_deg_km)
        #new_circle gives up possible coords for diff headings
        
        raw_weights = np.zeros(len(new_circle))
        for j in xrange(len(new_circle)):
            raw_weights[j] = new_weights[j]*ping_prob_normal(inmarsat[0],inmarsat[1],new_circle[j][0],new_circle[j][1],ping_percent_err,dist_from_sat,earth_radius) 
        
        probs = raw_weights / np.sum(raw_weights) #normalize
        
        index = range(len(new_circle))
        chosen = np.random.choice(index,size=None,p=probs)
        #print "chosen",chosen
        
        heading = new_angles[chosen] #update heading
        
        plane_lat[i],plane_lon[i] = new_circle[chosen] #update position
        lat = plane_lat[i]
        lon = plane_lon[i]
    
    #at end of simulation, run the last location & heading for plane for 4 different times
    route1 = make_vector(0.25*km_hop/eq_deg_km,heading,lon,lat)
    route2 = make_vector(0.5*km_hop/eq_deg_km,heading,lon,lat)
    route3 = make_vector(0.75*km_hop/eq_deg_km,heading,lon,lat)
    route4 = make_vector((59./60.)*km_hop/eq_deg_km,heading,lon,lat)

    new_plane_lat = np.zeros(9)
    new_plane_lon = np.zeros(9)
    
    for i in xrange(len(plane_lat)):
        new_plane_lat[i] = plane_lat[i]
        new_plane_lon[i] = plane_lon[i]
    
    new_plane_lat[5] = route1[0]
    new_plane_lat[6] = route2[0]
    new_plane_lat[7] = route3[0]
    new_plane_lat[8] = route4[0]
    new_plane_lon[5] = route1[1]
    new_plane_lon[6] = route2[1]
    new_plane_lon[7] = route3[1]
    new_plane_lon[8] = route4[1]
    
    return new_plane_lat,new_plane_lon

Modeling Constants (can vary)¶

Now we test out using the last known heading, a 777's typical cruising speed, and the last known coordinates, below:¶

In [23]:

last_known_heading = 255.136 #calculated in Mathematica from MH370's two last publically known locations:
                             #when it deviated from its flight path, and when it was last detected by Malyasian military radar
                             #0 degrees is due north, so this is basically to the west (270 degrees), but slightly south

km_hop = 905 #assuming 1 hr intervals, at 905 km/hr which is 777 cruising speed -- use for test case
             # max speed of a Boeing 777 is 950 km/hr FYI

N = 1000 #define number of simulations to run

5% Inmarsat error¶

In [24]:

percenterror = 0.05

std_dev = 30 #the standard deviation is the only arbitrary choice for us, along with the N=1000 simulations.
             #but it fits the notion that the plane is likely to continue on the same course, allowing for
             #some turns / heading change over the course of an hour. I stand by this number, although it's easy
             #to change to 15 or 45 and see what the difference is. The smaller the number the narrower the
             #ability to change heading is. The larger the number the more like a true "random walk" the plane's
             #direction and by consequence, location will be, at each time step.

In [25]:

plane_hops = []

for i in xrange(N):
    plane_hops.append(five_hop_model_normal(last_known_heading,pulauperak[1],pulauperak[0],km_hop,std_dev,percenterror))

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-25-7605bc96ecd9> in <module>()
      2 
      3 for i in xrange(N):
----> 4     plane_hops.append(five_hop_model_normal(last_known_heading,pulauperak[1],pulauperak[0],km_hop,std_dev,percenterror))

<ipython-input-22-5bd89d04f2e3> in five_hop_model_normal(heading_init, lon_init, lat_init, km_hop, std_dev, ping_percent_err)
     15 
     16     #initialize
---> 17     plane_lat = np.zeros(5) #initialize plane location after each hop (assumed to be 1 hr for now)
     18     plane_lon = np.zeros(5)
     19     lat = lat_init

NameError: global name 'np' is not defined

In [26]:

first_lat = []
two_lat = []
three_lat = []
four_lat = []
final_lat = []

first_lon = []
two_lon = []
three_lon = []
four_lon = []
final_lon = []

route1_lat = []
route2_lat = []
route3_lat = []
route4_lat = []

route1_lon = []
route2_lon = []
route3_lon = []
route4_lon = []

for i in xrange(len(plane_hops)):
    first_lat.append(plane_hops[i][0][0])
    first_lon.append(plane_hops[i][1][0])
    two_lat.append(plane_hops[i][0][1])
    two_lon.append(plane_hops[i][1][1])
    three_lat.append(plane_hops[i][0][2])
    three_lon.append(plane_hops[i][1][2])
    four_lat.append(plane_hops[i][0][3])
    four_lon.append(plane_hops[i][1][3])
    final_lat.append(plane_hops[i][0][4])
    final_lon.append(plane_hops[i][1][4])
    
    route1_lat.append(plane_hops[i][0][5])
    route1_lon.append(plane_hops[i][1][5])
    route2_lat.append(plane_hops[i][0][6])
    route2_lon.append(plane_hops[i][1][6])
    route3_lat.append(plane_hops[i][0][7])
    route3_lon.append(plane_hops[i][1][7])
    route4_lat.append(plane_hops[i][0][8])
    route4_lon.append(plane_hops[i][1][8])

In [27]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=13000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- 5% error
x5,y5 = fig(circle_lon,circle_lat)
x6,y6 = fig(circle_lon_err1_5per,circle_lat_err1_5per)
x7,y7 = fig(circle_lon_err2_5per,circle_lat_err2_5per)

#Add points after each hr
x9,y9 = fig(first_lon,first_lat)
x10,y10 = fig(two_lon,two_lat)
x11,y11 = fig(three_lon,three_lat)
x12,y12 = fig(four_lon,four_lat)
x13,y13 = fig(final_lon,final_lat)

#Draw circle showing extent of Inmarsat sat radar detection
fig.plot(x5,y5,'r-',markersize=5,label='5th Ping')
fig.plot(x6,y6,'r--',markersize=5,label='with 5% error')
fig.plot(x7,y7,'r--',markersize=5)

#Plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Add monte carlo points
fig.plot(x9,y9,'yo',markersize=5,label='after 1 hr')
fig.plot(x10,y10,'co',markersize=5,label='after 2 hrs')
fig.plot(x11,y11,'mo',markersize=5,label= 'after 3 hrs')
fig.plot(x12,y12,'wo',markersize=5,label='after 4 hrs')
fig.plot(x13,y13,'ro',markersize=7,label='after 5 hrs')

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='lower left',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('MH370 Position Progression Over Time', fontsize=30)

#Show below
plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-27-515feb35531b> in <module>()
      1 #Set figure size
----> 2 fig = plt.figure(figsize=[30,20])
      3 
      4 #Setup Basemap
      5 fig = Basemap(width=10000000,height=13000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

NameError: name 'plt' is not defined

In [28]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- 5% error
x5,y5 = fig(circle_lon,circle_lat)
x6,y6 = fig(circle_lon_err1_5per,circle_lat_err1_5per)
x7,y7 = fig(circle_lon_err2_5per,circle_lat_err2_5per)

#Add points after each hr
x9,y9 = fig(first_lon,first_lat)
x10,y10 = fig(two_lon,two_lat)
x11,y11 = fig(three_lon,three_lat)
x12,y12 = fig(four_lon,four_lat)
x13,y13 = fig(final_lon,final_lat)

#Add ultimate locations of MH370
x14,y14 = fig(route1_lon,route1_lat)
x15,y15 = fig(route2_lon,route2_lat)
x16,y16 = fig(route3_lon,route3_lat)
x17,y17 = fig(route4_lon,route4_lat)

#Draw circle showing extent of Inmarsat sat radar detection
fig.plot(x5,y5,'r-',markersize=5,label='5th Ping')
fig.plot(x6,y6,'r--',markersize=5,label='with 5% error')
fig.plot(x7,y7,'r--',markersize=5)

#Plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Add monte carlo points
fig.plot(x9,y9,'yo',markersize=5,label='after 1 hr')
fig.plot(x10,y10,'co',markersize=5,label='after 2 hrs')
fig.plot(x11,y11,'mo',markersize=5,label= 'after 3 hrs')
fig.plot(x12,y12,'wo',markersize=5,label='after 4 hrs')
fig.plot(x13,y13,'ro',markersize=7,label='after 5 hrs')

#Plot ultimate locations of MH370
fig.plot(x14,y14,'bo',markersize=5,label='in final hr')
fig.plot(x15,y15,'bo',markersize=5)
fig.plot(x16,y16,'bo',markersize=5)
fig.plot(x17,y17,'bo',markersize=5)

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='upper right',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('MH370 Position Progression Over Time', fontsize=30)

#Show below
plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-28-831074dafd8f> in <module>()
      1 #Set figure size
----> 2 fig = plt.figure(figsize=[30,20])
      3 
      4 #Setup Basemap
      5 fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

NameError: name 'plt' is not defined

2.5% Inmarsat error¶

In [ ]:

percenterror = 0.025

std_dev = 30

In [ ]:

plane_hops = []

for i in xrange(N):
    plane_hops.append(five_hop_model_normal(last_known_heading,pulauperak[1],pulauperak[0],km_hop,std_dev,percenterror))

In [ ]:

first_lat = []
two_lat = []
three_lat = []
four_lat = []
final_lat = []

first_lon = []
two_lon = []
three_lon = []
four_lon = []
final_lon = []

route1_lat = []
route2_lat = []
route3_lat = []
route4_lat = []

route1_lon = []
route2_lon = []
route3_lon = []
route4_lon = []

for i in xrange(len(plane_hops)):
    first_lat.append(plane_hops[i][0][0])
    first_lon.append(plane_hops[i][1][0])
    two_lat.append(plane_hops[i][0][1])
    two_lon.append(plane_hops[i][1][1])
    three_lat.append(plane_hops[i][0][2])
    three_lon.append(plane_hops[i][1][2])
    four_lat.append(plane_hops[i][0][3])
    four_lon.append(plane_hops[i][1][3])
    final_lat.append(plane_hops[i][0][4])
    final_lon.append(plane_hops[i][1][4])
    
    route1_lat.append(plane_hops[i][0][5])
    route1_lon.append(plane_hops[i][1][5])
    route2_lat.append(plane_hops[i][0][6])
    route2_lon.append(plane_hops[i][1][6])
    route3_lat.append(plane_hops[i][0][7])
    route3_lon.append(plane_hops[i][1][7])
    route4_lat.append(plane_hops[i][0][8])
    route4_lon.append(plane_hops[i][1][8])

In [ ]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=13000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- 2.5% error
x5,y5 = fig(circle_lon,circle_lat)
x6,y6 = fig(circle_lon_err1_2per,circle_lat_err1_2per)
x7,y7 = fig(circle_lon_err2_2per,circle_lat_err2_2per)

#Add points after each hr
x9,y9 = fig(first_lon,first_lat)
x10,y10 = fig(two_lon,two_lat)
x11,y11 = fig(three_lon,three_lat)
x12,y12 = fig(four_lon,four_lat)
x13,y13 = fig(final_lon,final_lat)

#Draw circle showing extent of Inmarsat sat radar detection
fig.plot(x5,y5,'r-',markersize=5,label='5th Ping')
fig.plot(x6,y6,'r--',markersize=5,label='with 2.5% error')
fig.plot(x7,y7,'r--',markersize=5)

#Plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Add monte carlo points
fig.plot(x9,y9,'yo',markersize=5,label='after 1 hr')
fig.plot(x10,y10,'co',markersize=5,label='after 2 hrs')
fig.plot(x11,y11,'mo',markersize=5,label= 'after 3 hrs')
fig.plot(x12,y12,'wo',markersize=5,label='after 4 hrs')
fig.plot(x13,y13,'ro',markersize=7,label='after 5 hrs')

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='lower left',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('MH370 Position Progression Over Time', fontsize=30)

#Show below
plt.show()

In [ ]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- 2.5% error
x5,y5 = fig(circle_lon,circle_lat)
x6,y6 = fig(circle_lon_err1_2per,circle_lat_err1_2per)
x7,y7 = fig(circle_lon_err2_2per,circle_lat_err2_2per)

#Add points after each hr
x9,y9 = fig(first_lon,first_lat)
x10,y10 = fig(two_lon,two_lat)
x11,y11 = fig(three_lon,three_lat)
x12,y12 = fig(four_lon,four_lat)
x13,y13 = fig(final_lon,final_lat)

#Add ultimate locations of MH370
x14,y14 = fig(route1_lon,route1_lat)
x15,y15 = fig(route2_lon,route2_lat)
x16,y16 = fig(route3_lon,route3_lat)
x17,y17 = fig(route4_lon,route4_lat)

#Draw circle showing extent of Inmarsat sat radar detection
fig.plot(x5,y5,'r-',markersize=5,label='5th Ping')
fig.plot(x6,y6,'r--',markersize=5,label='with 2.5% error')
fig.plot(x7,y7,'r--',markersize=5)

#Plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Add monte carlo points
fig.plot(x9,y9,'yo',markersize=5,label='after 1 hr')
fig.plot(x10,y10,'co',markersize=5,label='after 2 hrs')
fig.plot(x11,y11,'mo',markersize=5,label= 'after 3 hrs')
fig.plot(x12,y12,'wo',markersize=5,label='after 4 hrs')
fig.plot(x13,y13,'ro',markersize=7,label='after 5 hrs')

#Plot ultimate locations of MH370
fig.plot(x14,y14,'bo',markersize=5,label='in final hr')
fig.plot(x15,y15,'bo',markersize=5)
fig.plot(x16,y16,'bo',markersize=5)
fig.plot(x17,y17,'bo',markersize=5)

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='upper right',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('MH370 Position Progression Over Time', fontsize=30)

#Show below
plt.show()

Different Distributions Besides Normal for the Heading -- Von Mises and Wrapped Cauchy¶

Our normal heading distribution was great for a lot of reasons (including its ability to act like a uniform distribution if you pump the standard deviation up.) But, now it's time to consider a few other distributions as well, ones more exotic than our 'plain vanilla' one (well, the normal one that people use!) I'll pick ones that are popular to do along a unit circle; since our heading choice will be along a circle.¶

First, let's try a Von Mises distribution for our heading. Here's a general function that can do that:¶

In [ ]:

"""
takes in angle and heading in degrees (note conversion to radians in function)
k goes from 0 to infinity and is a point at infinity, and a uniform distribution at 0
"""
def von_mises(angle,heading,k):
    return np.exp(k*np.cos(np.radians(angle)-np.radians(heading)))/(2*np.pi*np.i0(k))

Let's test this function to see if it works. We'll test in with a heading of 150 degrees (which means that the curves should be centered at 150, i.e. the mean should be 150):¶

In [ ]:

circle_test = []

for i in xrange(360):
    circle_test.append(von_mises(i,150,0))
    
plt.figure()
plt.plot(circle_test)
plt.title('k=0 -- should be uniform')
plt.show()

circle_test2 = []

for i in xrange(360):
    circle_test2.append(von_mises(i,150,1))

plt.figure()
plt.plot(circle_test2)
plt.title('k=1 -- should be bell curve')
plt.show()

circle_test3 = []

for i in xrange(360):
    circle_test3.append(von_mises(i,150,10))
    
plt.figure()
plt.plot(circle_test3)
plt.title('k=10 -- should be contracting/converging toward a point')
plt.show()

It works! The other distribution we can try is a Wrapped Cauchy distribution for our heading. It's also a popular one to use on a unit circle (so you can say I used 'circular logic' when deciding to pick these two other distributions!)¶

In [ ]:

"""
takes in angle and heading in degrees (note conversion to radians in function)
k goes from 0-1 and is a point at 1, and a uniform distribution at 0
"""
def wrapped_cauchy(angle,heading,k):
    return np.float((1-k**2))/(1+k**2-2*k*np.cos(np.radians(angle)-np.radians(heading)))

We'll test this function in a similar manner:¶

In [ ]:

circle_test = []

for i in xrange(360):
    circle_test.append(wrapped_cauchy(i,150,0))
    
plt.figure()
plt.plot(circle_test)
plt.title('k=0 -- should be uniform')
plt.show()

circle_test2 = []

for i in xrange(360):
    circle_test2.append(wrapped_cauchy(i,150,0.5))

plt.figure()
plt.plot(circle_test2)
plt.title('k=0.5 -- should be bell curve')
plt.show()

circle_test3 = []

for i in xrange(360):
    circle_test3.append(wrapped_cauchy(i,150,0.99))
    
plt.figure()
plt.plot(circle_test3)
plt.title('k=0.99 -- should be contracting/converging toward a point')
plt.show()

It works as well. :0) :0) :0) Right, so now we have to update the prob_step functions to have Von Mises and Wrapped Cauchy for the heading, instead of a normal distribution. Note that k is effectively standard deviation, in that it varies the distribution to be point-like to uniform at the extreme ends:¶

In [ ]:

"""
takes in a heading, a starting location, and a k for picking a new heading based on the Von Mises distribution.
k = 0 -> uniform distribution. jump all over the place, true random walk.
k = 1,10 -> similar to before. tends to go straight ahead or near-straight ahead.
            1 is like a normal distribution with a large standard deviation from the mean.
            10 is like a normal distribution with a small standard deviation from the mean.
"""
def von_mises_prob_step(heading, k, start_lon, start_lat, ang_dist):
    
    angle_list = range(0,360)
    
    #we create a radius of 360 points that would be the heading for 360 possible degrees
    circle = make_circle(ang_dist,len(angle_list),start_lon,start_lat)
    
    weights = np.zeros(len(angle_list)) #make 360 array
    for i in xrange(len(weights)):
        weights[i] = von_mises(angle_list[i],heading,k)
        #makes array of von mises distributed weights.  
    
    #initialize new possible coordinates within an hr's distance, new weights for the odds, and new angles
    new_circle = []
    new_weights = []
    new_angles = []
    
    #make sure lat & lon are in bounds
    for i in xrange(len(circle)):
        if circle[i][0] >= lat_min and circle[i][0] <= lat_max and circle[i][1] >= lon_min and circle[i][1] <= lon_max:
            new_circle.append(circle[i])       
            new_weights.append(weights[i])
            new_angles.append(angle_list[i])
    
    return new_circle,new_weights,new_angles

In [ ]:

"""
takes in a heading, a starting location, and a k for picking a new heading based on the Wrapped Cauchy distribution.
k = 0 -> uniform distribution. jump all over the place, true random walk.
k = 0.5 -> similar to before. tends to go straight ahead or near-straight ahead.
k = 0.95 -> plane makes a bee-line for straight ahead. there's just no "convincing" it otherwise. 
            almost no deviation from 0 degrees relative to previous heading.
"""
def wrapped_cauchy_prob_step(heading, k, start_lon, start_lat, ang_dist):
    
    angle_list = range(0,360)
    
    #we create a radius of 360 points that would be the heading for 360 possible degrees
    circle = make_circle(ang_dist,len(angle_list),start_lon,start_lat)
    
    weights = np.zeros(len(angle_list)) #make 360 array
    for i in xrange(len(weights)):
        weights[i] = wrapped_cauchy(angle_list[i],heading,k)

    #initialize new possible coordinates within an hr's distance, new weights for the odds, and new angles
    new_circle = []
    new_weights = []
    new_angles = []
    
    #make sure lat & lon are in bounds
    for i in xrange(len(circle)):
        if circle[i][0] >= lat_min and circle[i][0] <= lat_max and circle[i][1] >= lon_min and circle[i][1] <= lon_max:
            new_circle.append(circle[i])       
            new_weights.append(weights[i])
            new_angles.append(angle_list[i])
    
    return new_circle,new_weights,new_angles

Finally we update our five hop function with the other distributions. There really isn't any math involved here -- just replacing with the appropriate functions and modifying the names of the inputs to better reflect what's going on in the Wrapped Cauchy and Von Mises distributions for headings:¶

In [ ]:

"""
a function which given a list of discrete probabilities for each destination point, 
it will choose one of those points.

heading_init -- initial direction was headed at last known point
lon_init,lat_init -- last known point of plane in longitude and latitude
km_hop -- how far the plane went in the time interval, 1 hr. So in simplest case, the 777's cruising speed/hr.
k -- affects the heading distribution, based on a Von Mises distribution from the current heading (0 deg).
ping_percent_err -- what % error you assume in the Inmarsat 5th ping. either 2.5 or 5%.

uses Von Mises distribution for heading

"""
def five_hop_model_von_mises(heading_init,lon_init,lat_init,km_hop,k,ping_percent_err):   
    
    #initialize
    plane_lat = np.zeros(5) #initialize plane location after each hop (assumed to be 1 hr for now)
    plane_lon = np.zeros(5)  
    lat = lat_init
    lon = lon_init
    heading = heading_init
        
    for i in xrange(len(plane_lat)):
        new_circle,new_weights,new_angles = von_mises_prob_step(heading,k,lon,lat,km_hop/eq_deg_km)
        #new_circle gives up possible coords for diff headings
        
        raw_weights = np.zeros(len(new_circle))
        for j in xrange(len(new_circle)):
            raw_weights[j] = new_weights[j]*ping_prob_normal(inmarsat[0],inmarsat[1],new_circle[j][0],new_circle[j][1],ping_percent_err,dist_from_sat,earth_radius) 
        
        probs = raw_weights / np.sum(raw_weights) #normalize
        
        index = range(len(new_circle))
        chosen = np.random.choice(index,size=None,p=probs)
        #print "chosen",chosen
        
        heading = new_angles[chosen] #update heading
        
        plane_lat[i],plane_lon[i] = new_circle[chosen] #update position
        lat = plane_lat[i]
        lon = plane_lon[i]
    
    #at end of simulation, run the last location & heading for plane for 4 different times
    route1 = make_vector(0.25*km_hop/eq_deg_km,heading,lon,lat)
    route2 = make_vector(0.5*km_hop/eq_deg_km,heading,lon,lat)
    route3 = make_vector(0.75*km_hop/eq_deg_km,heading,lon,lat)
    route4 = make_vector((59./60.)*km_hop/eq_deg_km,heading,lon,lat)

    new_plane_lat = np.zeros(9)
    new_plane_lon = np.zeros(9)
    
    for i in xrange(len(plane_lat)):
        new_plane_lat[i] = plane_lat[i]
        new_plane_lon[i] = plane_lon[i]
    
    new_plane_lat[5] = route1[0]
    new_plane_lat[6] = route2[0]
    new_plane_lat[7] = route3[0]
    new_plane_lat[8] = route4[0]
    new_plane_lon[5] = route1[1]
    new_plane_lon[6] = route2[1]
    new_plane_lon[7] = route3[1]
    new_plane_lon[8] = route4[1]
    
    return new_plane_lat,new_plane_lon

In [ ]:

"""
a function which given a list of discrete probabilities for each destination point, 
it will choose one of those points.

heading_init -- initial direction was headed at last known point
lon_init,lat_init -- last known point of plane in longitude and latitude
km_hop -- how far the plane went in the time interval, 1 hr. So in simplest case, the 777's cruising speed/hr.
k -- affects the heading distribution, based on a Wrapped Cauchy distribution from the current heading (0 deg).
ping_percent_err -- what % error you assume in the Inmarsat 5th ping. either 2.5 or 5%.

uses Wrapped Cauchy distribution for heading

"""
def five_hop_model_wrapped_cauchy(heading_init,lon_init,lat_init,km_hop,k,ping_percent_err):   
    
    #initialize
    plane_lat = np.zeros(5) #initialize plane location after each hop (assumed to be 1 hr for now)
    plane_lon = np.zeros(5)  
    lat = lat_init
    lon = lon_init
    heading = heading_init
        
    for i in xrange(len(plane_lat)):
        new_circle,new_weights,new_angles = wrapped_cauchy_prob_step(heading,k,lon,lat,km_hop/eq_deg_km)
        #new_circle gives up possible coords for diff headings
        
        raw_weights = np.zeros(len(new_circle))
        for j in xrange(len(new_circle)):
            raw_weights[j] = new_weights[j]*ping_prob_normal(inmarsat[0],inmarsat[1],new_circle[j][0],new_circle[j][1],ping_percent_err,dist_from_sat,earth_radius) 
        
        probs = raw_weights / np.sum(raw_weights) #normalize
        
        index = range(len(new_circle))
        chosen = np.random.choice(index,size=None,p=probs)
        #print "chosen",chosen
        
        heading = new_angles[chosen] #update heading
        
        plane_lat[i],plane_lon[i] = new_circle[chosen] #update position
        lat = plane_lat[i]
        lon = plane_lon[i]
    
    #at end of simulation, run the last location & heading for plane for 4 different times
    route1 = make_vector(0.25*km_hop/eq_deg_km,heading,lon,lat)
    route2 = make_vector(0.5*km_hop/eq_deg_km,heading,lon,lat)
    route3 = make_vector(0.75*km_hop/eq_deg_km,heading,lon,lat)
    route4 = make_vector((59./60.)*km_hop/eq_deg_km,heading,lon,lat)

    new_plane_lat = np.zeros(9)
    new_plane_lon = np.zeros(9)
    
    for i in xrange(len(plane_lat)):
        new_plane_lat[i] = plane_lat[i]
        new_plane_lon[i] = plane_lon[i]
    
    new_plane_lat[5] = route1[0]
    new_plane_lat[6] = route2[0]
    new_plane_lat[7] = route3[0]
    new_plane_lat[8] = route4[0]
    new_plane_lon[5] = route1[1]
    new_plane_lon[6] = route2[1]
    new_plane_lon[7] = route3[1]
    new_plane_lon[8] = route4[1]
    
    return new_plane_lat,new_plane_lon

Running Code with New Heading Distributions¶

Now we'll run the code with 5% and 2.5% to see if there are any major changes. We've made our probability grid, so all we need now is to rereun the simulations substituting in our different heading distribution. I'll vary these by more than what I got for the normal distribution with a std dev of 30; otherwise, it follows from the above functions that we could pick parameters for all 3 distributions (normal, von mises, wrapped cauchy) that would all result in the similar looking curves.¶

In [ ]:

last_known_heading = 255.136

km_hop = 905 #assuming 1 hr intervals, at 905 km/hr which is 777 cruising speed -- use for test case
             # max speed of a Boeing 777 is 950 km/hr FYI

N = 1000 #define number of simulations to run -- set to 1000, but set low number to make sure it's working first.

Von Mises¶

To show what a uniform distribution (random walk) to a point (straight ahead scenario), along with the realistic in between (tends to go straight, but can turn), I'll use the same k values that I used to test the Von Mises function. And I'll do these for both the 2.5% std dev error and the 5% std dev error for the ping probability distribution, which remains normal.¶

In [ ]:

percenterror1, percenterror2 = 0.05, 0.025

k1, k2, k3 = 0, 1, 10 #Von Mises constants -- uniform, bell, point-like

In [ ]:

#5 percent error w/ different ks
plane_hops_5per_k1 = []
plane_hops_5per_k2 = []
plane_hops_5per_k3 = []

#2.5 percent error w/ different ks
plane_hops_2per_k1 = []
plane_hops_2per_k2 = []
plane_hops_2per_k3 = []

for i in xrange(N):
    #5% error runs
    plane_hops_5per_k1.append(five_hop_model_von_mises(last_known_heading,pulauperak[1],pulauperak[0],km_hop,k1,percenterror1))
    plane_hops_5per_k2.append(five_hop_model_von_mises(last_known_heading,pulauperak[1],pulauperak[0],km_hop,k2,percenterror1))
    plane_hops_5per_k3.append(five_hop_model_von_mises(last_known_heading,pulauperak[1],pulauperak[0],km_hop,k3,percenterror1))
    #2.5% error runs
    plane_hops_2per_k1.append(five_hop_model_von_mises(last_known_heading,pulauperak[1],pulauperak[0],km_hop,k1,percenterror2))
    plane_hops_2per_k2.append(five_hop_model_von_mises(last_known_heading,pulauperak[1],pulauperak[0],km_hop,k2,percenterror2))
    plane_hops_2per_k3.append(five_hop_model_von_mises(last_known_heading,pulauperak[1],pulauperak[0],km_hop,k3,percenterror2))

The above loop should take a little while to run. Unfortunately, given the way I've written the code by appending to lists, the below will also be a little long (in terms of space); but just run the below cells and it should quickly (in terms of run time) gather all of the data we've generated into lists for plotting purposes:¶

In [ ]:

# 5per_k1 run

first_lat_5per_k1 = []
two_lat_5per_k1 = []
three_lat_5per_k1 = []
four_lat_5per_k1 = []
final_lat_5per_k1 = []

first_lon_5per_k1 = []
two_lon_5per_k1 = []
three_lon_5per_k1 = []
four_lon_5per_k1 = []
final_lon_5per_k1 = []

route1_lat_5per_k1 = []
route2_lat_5per_k1 = []
route3_lat_5per_k1 = []
route4_lat_5per_k1 = []

route1_lon_5per_k1 = []
route2_lon_5per_k1 = []
route3_lon_5per_k1 = []
route4_lon_5per_k1 = []

for i in xrange(len(plane_hops_5per_k1)):
    first_lat_5per_k1.append(plane_hops_5per_k1[i][0][0])
    first_lon_5per_k1.append(plane_hops_5per_k1[i][1][0])
    two_lat_5per_k1.append(plane_hops_5per_k1[i][0][1])
    two_lon_5per_k1.append(plane_hops_5per_k1[i][1][1])
    three_lat_5per_k1.append(plane_hops_5per_k1[i][0][2])
    three_lon_5per_k1.append(plane_hops_5per_k1[i][1][2])
    four_lat_5per_k1.append(plane_hops_5per_k1[i][0][3])
    four_lon_5per_k1.append(plane_hops_5per_k1[i][1][3])
    final_lat_5per_k1.append(plane_hops_5per_k1[i][0][4])
    final_lon_5per_k1.append(plane_hops_5per_k1[i][1][4])
    
    route1_lat_5per_k1.append(plane_hops_5per_k1[i][0][5])
    route1_lon_5per_k1.append(plane_hops_5per_k1[i][1][5])
    route2_lat_5per_k1.append(plane_hops_5per_k1[i][0][6])
    route2_lon_5per_k1.append(plane_hops_5per_k1[i][1][6])
    route3_lat_5per_k1.append(plane_hops_5per_k1[i][0][7])
    route3_lon_5per_k1.append(plane_hops_5per_k1[i][1][7])
    route4_lat_5per_k1.append(plane_hops_5per_k1[i][0][8])
    route4_lon_5per_k1.append(plane_hops_5per_k1[i][1][8])

In [ ]:

# 5per_k2 run

first_lat_5per_k2 = []
two_lat_5per_k2 = []
three_lat_5per_k2 = []
four_lat_5per_k2 = []
final_lat_5per_k2 = []

first_lon_5per_k2 = []
two_lon_5per_k2 = []
three_lon_5per_k2 = []
four_lon_5per_k2 = []
final_lon_5per_k2 = []

route1_lat_5per_k2 = []
route2_lat_5per_k2 = []
route3_lat_5per_k2 = []
route4_lat_5per_k2 = []

route1_lon_5per_k2 = []
route2_lon_5per_k2 = []
route3_lon_5per_k2 = []
route4_lon_5per_k2 = []

for i in xrange(len(plane_hops_5per_k2)):
    first_lat_5per_k2.append(plane_hops_5per_k2[i][0][0])
    first_lon_5per_k2.append(plane_hops_5per_k2[i][1][0])
    two_lat_5per_k2.append(plane_hops_5per_k2[i][0][1])
    two_lon_5per_k2.append(plane_hops_5per_k2[i][1][1])
    three_lat_5per_k2.append(plane_hops_5per_k2[i][0][2])
    three_lon_5per_k2.append(plane_hops_5per_k2[i][1][2])
    four_lat_5per_k2.append(plane_hops_5per_k2[i][0][3])
    four_lon_5per_k2.append(plane_hops_5per_k2[i][1][3])
    final_lat_5per_k2.append(plane_hops_5per_k2[i][0][4])
    final_lon_5per_k2.append(plane_hops_5per_k2[i][1][4])
    
    route1_lat_5per_k2.append(plane_hops_5per_k2[i][0][5])
    route1_lon_5per_k2.append(plane_hops_5per_k2[i][1][5])
    route2_lat_5per_k2.append(plane_hops_5per_k2[i][0][6])
    route2_lon_5per_k2.append(plane_hops_5per_k2[i][1][6])
    route3_lat_5per_k2.append(plane_hops_5per_k2[i][0][7])
    route3_lon_5per_k2.append(plane_hops_5per_k2[i][1][7])
    route4_lat_5per_k2.append(plane_hops_5per_k2[i][0][8])
    route4_lon_5per_k2.append(plane_hops_5per_k2[i][1][8])

In [ ]:

# 5per_k3 run

first_lat_5per_k3 = []
two_lat_5per_k3 = []
three_lat_5per_k3 = []
four_lat_5per_k3 = []
final_lat_5per_k3 = []

first_lon_5per_k3 = []
two_lon_5per_k3 = []
three_lon_5per_k3 = []
four_lon_5per_k3 = []
final_lon_5per_k3 = []

route1_lat_5per_k3 = []
route2_lat_5per_k3 = []
route3_lat_5per_k3 = []
route4_lat_5per_k3 = []

route1_lon_5per_k3 = []
route2_lon_5per_k3 = []
route3_lon_5per_k3 = []
route4_lon_5per_k3 = []

for i in xrange(len(plane_hops_5per_k3)):
    first_lat_5per_k3.append(plane_hops_5per_k3[i][0][0])
    first_lon_5per_k3.append(plane_hops_5per_k3[i][1][0])
    two_lat_5per_k3.append(plane_hops_5per_k3[i][0][1])
    two_lon_5per_k3.append(plane_hops_5per_k3[i][1][1])
    three_lat_5per_k3.append(plane_hops_5per_k3[i][0][2])
    three_lon_5per_k3.append(plane_hops_5per_k3[i][1][2])
    four_lat_5per_k3.append(plane_hops_5per_k3[i][0][3])
    four_lon_5per_k3.append(plane_hops_5per_k3[i][1][3])
    final_lat_5per_k3.append(plane_hops_5per_k3[i][0][4])
    final_lon_5per_k3.append(plane_hops_5per_k3[i][1][4])
    
    route1_lat_5per_k3.append(plane_hops_5per_k3[i][0][5])
    route1_lon_5per_k3.append(plane_hops_5per_k3[i][1][5])
    route2_lat_5per_k3.append(plane_hops_5per_k3[i][0][6])
    route2_lon_5per_k3.append(plane_hops_5per_k3[i][1][6])
    route3_lat_5per_k3.append(plane_hops_5per_k3[i][0][7])
    route3_lon_5per_k3.append(plane_hops_5per_k3[i][1][7])
    route4_lat_5per_k3.append(plane_hops_5per_k3[i][0][8])
    route4_lon_5per_k3.append(plane_hops_5per_k3[i][1][8])

In [29]:

# 2per_k1 run

first_lat_2per_k1 = []
two_lat_2per_k1 = []
three_lat_2per_k1 = []
four_lat_2per_k1 = []
final_lat_2per_k1 = []

first_lon_2per_k1 = []
two_lon_2per_k1 = []
three_lon_2per_k1 = []
four_lon_2per_k1 = []
final_lon_2per_k1 = []

route1_lat_2per_k1 = []
route2_lat_2per_k1 = []
route3_lat_2per_k1 = []
route4_lat_2per_k1 = []

route1_lon_2per_k1 = []
route2_lon_2per_k1 = []
route3_lon_2per_k1 = []
route4_lon_2per_k1 = []

for i in xrange(len(plane_hops_2per_k1)):
    first_lat_2per_k1.append(plane_hops_2per_k1[i][0][0])
    first_lon_2per_k1.append(plane_hops_2per_k1[i][1][0])
    two_lat_2per_k1.append(plane_hops_2per_k1[i][0][1])
    two_lon_2per_k1.append(plane_hops_2per_k1[i][1][1])
    three_lat_2per_k1.append(plane_hops_2per_k1[i][0][2])
    three_lon_2per_k1.append(plane_hops_2per_k1[i][1][2])
    four_lat_2per_k1.append(plane_hops_2per_k1[i][0][3])
    four_lon_2per_k1.append(plane_hops_2per_k1[i][1][3])
    final_lat_2per_k1.append(plane_hops_2per_k1[i][0][4])
    final_lon_2per_k1.append(plane_hops_2per_k1[i][1][4])
    
    route1_lat_2per_k1.append(plane_hops_2per_k1[i][0][5])
    route1_lon_2per_k1.append(plane_hops_2per_k1[i][1][5])
    route2_lat_2per_k1.append(plane_hops_2per_k1[i][0][6])
    route2_lon_2per_k1.append(plane_hops_2per_k1[i][1][6])
    route3_lat_2per_k1.append(plane_hops_2per_k1[i][0][7])
    route3_lon_2per_k1.append(plane_hops_2per_k1[i][1][7])
    route4_lat_2per_k1.append(plane_hops_2per_k1[i][0][8])
    route4_lon_2per_k1.append(plane_hops_2per_k1[i][1][8])

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-29-b87991f6d41f> in <module>()
     23 route4_lon_2per_k1 = []
     24 
---> 25 for i in xrange(len(plane_hops_2per_k1)):
     26     first_lat_2per_k1.append(plane_hops_2per_k1[i][0][0])
     27     first_lon_2per_k1.append(plane_hops_2per_k1[i][1][0])

NameError: name 'plane_hops_2per_k1' is not defined

In [ ]:

# 2per_k2 run

first_lat_2per_k2 = []
two_lat_2per_k2 = []
three_lat_2per_k2 = []
four_lat_2per_k2 = []
final_lat_2per_k2 = []

first_lon_2per_k2 = []
two_lon_2per_k2 = []
three_lon_2per_k2 = []
four_lon_2per_k2 = []
final_lon_2per_k2 = []

route1_lat_2per_k2 = []
route2_lat_2per_k2 = []
route3_lat_2per_k2 = []
route4_lat_2per_k2 = []

route1_lon_2per_k2 = []
route2_lon_2per_k2 = []
route3_lon_2per_k2 = []
route4_lon_2per_k2 = []

for i in xrange(len(plane_hops_2per_k2)):
    first_lat_2per_k2.append(plane_hops_2per_k2[i][0][0])
    first_lon_2per_k2.append(plane_hops_2per_k2[i][1][0])
    two_lat_2per_k2.append(plane_hops_2per_k2[i][0][1])
    two_lon_2per_k2.append(plane_hops_2per_k2[i][1][1])
    three_lat_2per_k2.append(plane_hops_2per_k2[i][0][2])
    three_lon_2per_k2.append(plane_hops_2per_k2[i][1][2])
    four_lat_2per_k2.append(plane_hops_2per_k2[i][0][3])
    four_lon_2per_k2.append(plane_hops_2per_k2[i][1][3])
    final_lat_2per_k2.append(plane_hops_2per_k2[i][0][4])
    final_lon_2per_k2.append(plane_hops_2per_k2[i][1][4])
    
    route1_lat_2per_k2.append(plane_hops_2per_k2[i][0][5])
    route1_lon_2per_k2.append(plane_hops_2per_k2[i][1][5])
    route2_lat_2per_k2.append(plane_hops_2per_k2[i][0][6])
    route2_lon_2per_k2.append(plane_hops_2per_k2[i][1][6])
    route3_lat_2per_k2.append(plane_hops_2per_k2[i][0][7])
    route3_lon_2per_k2.append(plane_hops_2per_k2[i][1][7])
    route4_lat_2per_k2.append(plane_hops_2per_k2[i][0][8])
    route4_lon_2per_k2.append(plane_hops_2per_k2[i][1][8])

In [ ]:

# 2per_k3 run

first_lat_2per_k3 = []
two_lat_2per_k3 = []
three_lat_2per_k3 = []
four_lat_2per_k3 = []
final_lat_2per_k3 = []

first_lon_2per_k3 = []
two_lon_2per_k3 = []
three_lon_2per_k3 = []
four_lon_2per_k3 = []
final_lon_2per_k3 = []

route1_lat_2per_k3 = []
route2_lat_2per_k3 = []
route3_lat_2per_k3 = []
route4_lat_2per_k3 = []

route1_lon_2per_k3 = []
route2_lon_2per_k3 = []
route3_lon_2per_k3 = []
route4_lon_2per_k3 = []

for i in xrange(len(plane_hops_2per_k3)):
    first_lat_2per_k3.append(plane_hops_2per_k3[i][0][0])
    first_lon_2per_k3.append(plane_hops_2per_k3[i][1][0])
    two_lat_2per_k3.append(plane_hops_2per_k3[i][0][1])
    two_lon_2per_k3.append(plane_hops_2per_k3[i][1][1])
    three_lat_2per_k3.append(plane_hops_2per_k3[i][0][2])
    three_lon_2per_k3.append(plane_hops_2per_k3[i][1][2])
    four_lat_2per_k3.append(plane_hops_2per_k3[i][0][3])
    four_lon_2per_k3.append(plane_hops_2per_k3[i][1][3])
    final_lat_2per_k3.append(plane_hops_2per_k3[i][0][4])
    final_lon_2per_k3.append(plane_hops_2per_k3[i][1][4])
    
    route1_lat_2per_k3.append(plane_hops_2per_k3[i][0][5])
    route1_lon_2per_k3.append(plane_hops_2per_k3[i][1][5])
    route2_lat_2per_k3.append(plane_hops_2per_k3[i][0][6])
    route2_lon_2per_k3.append(plane_hops_2per_k3[i][1][6])
    route3_lat_2per_k3.append(plane_hops_2per_k3[i][0][7])
    route3_lon_2per_k3.append(plane_hops_2per_k3[i][1][7])
    route4_lat_2per_k3.append(plane_hops_2per_k3[i][0][8])
    route4_lon_2per_k3.append(plane_hops_2per_k3[i][1][8])

Now we plot. To conserve space I am plotting all of the points at once. I will not make a separate plot just showing up to the 5th ping:¶

5per_k1 plot - Von Mises Heading, k=0, 5% ping error¶

In [ ]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- 5% error
x5,y5 = fig(circle_lon,circle_lat)
x6,y6 = fig(circle_lon_err1_5per,circle_lat_err1_5per)
x7,y7 = fig(circle_lon_err2_5per,circle_lat_err2_5per)

#Add points after each hr
x9,y9 = fig(first_lon_5per_k1,first_lat_5per_k1)
x10,y10 = fig(two_lon_5per_k1,two_lat_5per_k1)
x11,y11 = fig(three_lon_5per_k1,three_lat_5per_k1)
x12,y12 = fig(four_lon_5per_k1,four_lat_5per_k1)
x13,y13 = fig(final_lon_5per_k1,final_lat_5per_k1)

#Add ultimate locations of MH370
x14,y14 = fig(route1_lon_5per_k1,route1_lat_5per_k1)
x15,y15 = fig(route2_lon_5per_k1,route2_lat_5per_k1)
x16,y16 = fig(route3_lon_5per_k1,route3_lat_5per_k1)
x17,y17 = fig(route4_lon_5per_k1,route4_lat_5per_k1)

#Draw circle showing extent of Inmarsat sat radar detection
fig.plot(x5,y5,'r-',markersize=5,label='5th Ping')
fig.plot(x6,y6,'r--',markersize=5,label='with 5% error')
fig.plot(x7,y7,'r--',markersize=5)

#Plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Add monte carlo points
fig.plot(x9,y9,'yo',markersize=5,label='after 1 hr')
fig.plot(x10,y10,'co',markersize=5,label='after 2 hrs')
fig.plot(x11,y11,'mo',markersize=5,label= 'after 3 hrs')
fig.plot(x12,y12,'wo',markersize=5,label='after 4 hrs')
fig.plot(x13,y13,'ro',markersize=7,label='after 5 hrs')

#Plot ultimate locations of MH370
fig.plot(x14,y14,'bo',markersize=5,label='in final hr')
fig.plot(x15,y15,'bo',markersize=5)
fig.plot(x16,y16,'bo',markersize=5)
fig.plot(x17,y17,'bo',markersize=5)

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='upper right',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('MH370 Position Progression Over Time', fontsize=30)

#Show below
plt.show()

5per_k2 plot - Von Mises Heading, k=1, 5% ping error¶

In [ ]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- 5% error
x5,y5 = fig(circle_lon,circle_lat)
x6,y6 = fig(circle_lon_err1_5per,circle_lat_err1_5per)
x7,y7 = fig(circle_lon_err2_5per,circle_lat_err2_5per)

#Add points after each hr
x9,y9 = fig(first_lon_5per_k2,first_lat_5per_k2)
x10,y10 = fig(two_lon_5per_k2,two_lat_5per_k2)
x11,y11 = fig(three_lon_5per_k2,three_lat_5per_k2)
x12,y12 = fig(four_lon_5per_k2,four_lat_5per_k2)
x13,y13 = fig(final_lon_5per_k2,final_lat_5per_k2)

#Add ultimate locations of MH370
x14,y14 = fig(route1_lon_5per_k2,route1_lat_5per_k2)
x15,y15 = fig(route2_lon_5per_k2,route2_lat_5per_k2)
x16,y16 = fig(route3_lon_5per_k2,route3_lat_5per_k2)
x17,y17 = fig(route4_lon_5per_k2,route4_lat_5per_k2)

#Draw circle showing extent of Inmarsat sat radar detection
fig.plot(x5,y5,'r-',markersize=5,label='5th Ping')
fig.plot(x6,y6,'r--',markersize=5,label='with 5% error')
fig.plot(x7,y7,'r--',markersize=5)

#Plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Add monte carlo points
fig.plot(x9,y9,'yo',markersize=5,label='after 1 hr')
fig.plot(x10,y10,'co',markersize=5,label='after 2 hrs')
fig.plot(x11,y11,'mo',markersize=5,label= 'after 3 hrs')
fig.plot(x12,y12,'wo',markersize=5,label='after 4 hrs')
fig.plot(x13,y13,'ro',markersize=7,label='after 5 hrs')

#Plot ultimate locations of MH370
fig.plot(x14,y14,'bo',markersize=5,label='in final hr')
fig.plot(x15,y15,'bo',markersize=5)
fig.plot(x16,y16,'bo',markersize=5)
fig.plot(x17,y17,'bo',markersize=5)

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='upper right',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('MH370 Position Progression Over Time', fontsize=30)

#Show below
plt.show()

5per_k3 plot - Von Mises Heading, k=10, 5% ping error¶

In [ ]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- 5% error
x5,y5 = fig(circle_lon,circle_lat)
x6,y6 = fig(circle_lon_err1_5per,circle_lat_err1_5per)
x7,y7 = fig(circle_lon_err2_5per,circle_lat_err2_5per)

#Add points after each hr
x9,y9 = fig(first_lon_5per_k3,first_lat_5per_k3)
x10,y10 = fig(two_lon_5per_k3,two_lat_5per_k3)
x11,y11 = fig(three_lon_5per_k3,three_lat_5per_k3)
x12,y12 = fig(four_lon_5per_k3,four_lat_5per_k3)
x13,y13 = fig(final_lon_5per_k3,final_lat_5per_k3)

#Add ultimate locations of MH370
x14,y14 = fig(route1_lon_5per_k3,route1_lat_5per_k3)
x15,y15 = fig(route2_lon_5per_k3,route2_lat_5per_k3)
x16,y16 = fig(route3_lon_5per_k3,route3_lat_5per_k3)
x17,y17 = fig(route4_lon_5per_k3,route4_lat_5per_k3)

#Draw circle showing extent of Inmarsat sat radar detection
fig.plot(x5,y5,'r-',markersize=5,label='5th Ping')
fig.plot(x6,y6,'r--',markersize=5,label='with 5% error')
fig.plot(x7,y7,'r--',markersize=5)

#Plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Add monte carlo points
fig.plot(x9,y9,'yo',markersize=5,label='after 1 hr')
fig.plot(x10,y10,'co',markersize=5,label='after 2 hrs')
fig.plot(x11,y11,'mo',markersize=5,label= 'after 3 hrs')
fig.plot(x12,y12,'wo',markersize=5,label='after 4 hrs')
fig.plot(x13,y13,'ro',markersize=7,label='after 5 hrs')

#Plot ultimate locations of MH370
fig.plot(x14,y14,'bo',markersize=5,label='in final hr')
fig.plot(x15,y15,'bo',markersize=5)
fig.plot(x16,y16,'bo',markersize=5)
fig.plot(x17,y17,'bo',markersize=5)

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='upper right',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('MH370 Position Progression Over Time', fontsize=30)

#Show below
plt.show()

2per_k1 plot - Von Mises Heading, k=0, 2.5% ping error¶

In [ ]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- 2.5% error
x5,y5 = fig(circle_lon,circle_lat)
x6,y6 = fig(circle_lon_err1_2per,circle_lat_err1_2per)
x7,y7 = fig(circle_lon_err2_2per,circle_lat_err2_2per)

#Add points after each hr
x9,y9 = fig(first_lon_5per_k1,first_lat_5per_k1)
x10,y10 = fig(two_lon_5per_k1,two_lat_5per_k1)
x11,y11 = fig(three_lon_5per_k1,three_lat_5per_k1)
x12,y12 = fig(four_lon_5per_k1,four_lat_5per_k1)
x13,y13 = fig(final_lon_5per_k1,final_lat_5per_k1)

#Add ultimate locations of MH370
x14,y14 = fig(route1_lon_5per_k1,route1_lat_5per_k1)
x15,y15 = fig(route2_lon_5per_k1,route2_lat_5per_k1)
x16,y16 = fig(route3_lon_5per_k1,route3_lat_5per_k1)
x17,y17 = fig(route4_lon_5per_k1,route4_lat_5per_k1)

#Draw circle showing extent of Inmarsat sat radar detection
fig.plot(x5,y5,'r-',markersize=5,label='5th Ping')
fig.plot(x6,y6,'r--',markersize=5,label='with 5% error')
fig.plot(x7,y7,'r--',markersize=5)

#Plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Add monte carlo points
fig.plot(x9,y9,'yo',markersize=5,label='after 1 hr')
fig.plot(x10,y10,'co',markersize=5,label='after 2 hrs')
fig.plot(x11,y11,'mo',markersize=5,label= 'after 3 hrs')
fig.plot(x12,y12,'wo',markersize=5,label='after 4 hrs')
fig.plot(x13,y13,'ro',markersize=7,label='after 5 hrs')

#Plot ultimate locations of MH370
fig.plot(x14,y14,'bo',markersize=5,label='in final hr')
fig.plot(x15,y15,'bo',markersize=5)
fig.plot(x16,y16,'bo',markersize=5)
fig.plot(x17,y17,'bo',markersize=5)

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='upper right',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('MH370 Position Progression Over Time', fontsize=30)

#Show below
plt.show()

2per_k2 plot - Von Mises Heading, k=1, 2.5% ping error¶

In [ ]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- 2.5% error
x5,y5 = fig(circle_lon,circle_lat)
x6,y6 = fig(circle_lon_err1_2per,circle_lat_err1_2per)
x7,y7 = fig(circle_lon_err2_2per,circle_lat_err2_2per)

#Add points after each hr
x9,y9 = fig(first_lon_5per_k2,first_lat_5per_k2)
x10,y10 = fig(two_lon_5per_k2,two_lat_5per_k2)
x11,y11 = fig(three_lon_5per_k2,three_lat_5per_k2)
x12,y12 = fig(four_lon_5per_k2,four_lat_5per_k2)
x13,y13 = fig(final_lon_5per_k2,final_lat_5per_k2)

#Add ultimate locations of MH370
x14,y14 = fig(route1_lon_5per_k2,route1_lat_5per_k2)
x15,y15 = fig(route2_lon_5per_k2,route2_lat_5per_k2)
x16,y16 = fig(route3_lon_5per_k2,route3_lat_5per_k2)
x17,y17 = fig(route4_lon_5per_k2,route4_lat_5per_k2)

#Draw circle showing extent of Inmarsat sat radar detection
fig.plot(x5,y5,'r-',markersize=5,label='5th Ping')
fig.plot(x6,y6,'r--',markersize=5,label='with 5% error')
fig.plot(x7,y7,'r--',markersize=5)

#Plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Add monte carlo points
fig.plot(x9,y9,'yo',markersize=5,label='after 1 hr')
fig.plot(x10,y10,'co',markersize=5,label='after 2 hrs')
fig.plot(x11,y11,'mo',markersize=5,label= 'after 3 hrs')
fig.plot(x12,y12,'wo',markersize=5,label='after 4 hrs')
fig.plot(x13,y13,'ro',markersize=7,label='after 5 hrs')

#Plot ultimate locations of MH370
fig.plot(x14,y14,'bo',markersize=5,label='in final hr')
fig.plot(x15,y15,'bo',markersize=5)
fig.plot(x16,y16,'bo',markersize=5)
fig.plot(x17,y17,'bo',markersize=5)

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='upper right',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('MH370 Position Progression Over Time', fontsize=30)

#Show below
plt.show()

2per_k3 plot - Von Mises Heading, k=10, 2.5% ping error¶

In [ ]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- 2.5% error
x5,y5 = fig(circle_lon,circle_lat)
x6,y6 = fig(circle_lon_err1_2per,circle_lat_err1_2per)
x7,y7 = fig(circle_lon_err2_2per,circle_lat_err2_2per)

#Add points after each hr
x9,y9 = fig(first_lon_5per_k3,first_lat_5per_k3)
x10,y10 = fig(two_lon_5per_k3,two_lat_5per_k3)
x11,y11 = fig(three_lon_5per_k3,three_lat_5per_k3)
x12,y12 = fig(four_lon_5per_k3,four_lat_5per_k3)
x13,y13 = fig(final_lon_5per_k3,final_lat_5per_k3)

#Add ultimate locations of MH370
x14,y14 = fig(route1_lon_5per_k3,route1_lat_5per_k3)
x15,y15 = fig(route2_lon_5per_k3,route2_lat_5per_k3)
x16,y16 = fig(route3_lon_5per_k3,route3_lat_5per_k3)
x17,y17 = fig(route4_lon_5per_k3,route4_lat_5per_k3)

#Draw circle showing extent of Inmarsat sat radar detection
fig.plot(x5,y5,'r-',markersize=5,label='5th Ping')
fig.plot(x6,y6,'r--',markersize=5,label='with 5% error')
fig.plot(x7,y7,'r--',markersize=5)

#Plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Add monte carlo points
fig.plot(x9,y9,'yo',markersize=5,label='after 1 hr')
fig.plot(x10,y10,'co',markersize=5,label='after 2 hrs')
fig.plot(x11,y11,'mo',markersize=5,label= 'after 3 hrs')
fig.plot(x12,y12,'wo',markersize=5,label='after 4 hrs')
fig.plot(x13,y13,'ro',markersize=7,label='after 5 hrs')

#Plot ultimate locations of MH370
fig.plot(x14,y14,'bo',markersize=5,label='in final hr')
fig.plot(x15,y15,'bo',markersize=5)
fig.plot(x16,y16,'bo',markersize=5)
fig.plot(x17,y17,'bo',markersize=5)

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='upper right',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('MH370 Position Progression Over Time', fontsize=30)

#Show below
plt.show()

Okay, let's take a minute to see what these plots mean!¶

......... ........... ........ ......... ........¶

If you look at the k=0 , it follows that a uniform distribution will be a true random walk; so the plane will tend to move any which way, and influenced by the ping. Unfortunately, this does not tell us much about where along the arc the plane is most likely to be located. On reddit (MH370 subbreddit), me not using a uniform distribution was a criticism which got lots of upvotes; well, I've implemented it now, and the results are not what we're looking for! They're not helpful because they do not constrain the position of MH370; and the assumption inherent in a uniform distribution is not an accurate representation of how an airplane moves when it's trying to get somewhere.¶

As you increase k, the distribution deviating from the heading becomes narrower. And as that becomes apparent, we do in fact see that the plane tends to the south instead of the north. And from looking at a uniform distribution, this is very likely due to not only the geometry of the arc, but also the plane's last known heading -- 255 degrees, slight south and to the west (due west is 270), as opposed to slightly to the north.¶

Finally, in theory 2.5% should have stronger results than 5%, in terms of restricting the search area location at the end...but in practice both are relatively small errors; so we're looking at the same patch of ocean west of Australia.¶

Wrapped Cauchy¶

Similar to the Von Mises runs, I'll use the test k values for Wrapped Cauchy at 2.5 and 5% error for the ping probability distribution.¶

In [ ]:

percenterror1, percenterror2 = 0.05, 0.025

k1, k2, k3 = 0, 0.5, 0.99 #Wrapped Cauchy constants -- uniform, bell, point-like

In [ ]:

#5 percent error w/ different ks
plane_hops_5per_k1 = []
plane_hops_5per_k2 = []
plane_hops_5per_k3 = []

#2.5 percent error w/ different ks
plane_hops_2per_k1 = []
plane_hops_2per_k2 = []
plane_hops_2per_k3 = []

for i in xrange(N):
    #5% error runs
    plane_hops_5per_k1.append(five_hop_model_wrapped_cauchy(last_known_heading,pulauperak[1],pulauperak[0],km_hop,k1,percenterror1))
    plane_hops_5per_k2.append(five_hop_model_wrapped_cauchy(last_known_heading,pulauperak[1],pulauperak[0],km_hop,k2,percenterror1))
    plane_hops_5per_k3.append(five_hop_model_wrapped_cauchy(last_known_heading,pulauperak[1],pulauperak[0],km_hop,k3,percenterror1))
    #2.5% error runs
    plane_hops_2per_k1.append(five_hop_model_wrapped_cauchy(last_known_heading,pulauperak[1],pulauperak[0],km_hop,k1,percenterror2))
    plane_hops_2per_k2.append(five_hop_model_wrapped_cauchy(last_known_heading,pulauperak[1],pulauperak[0],km_hop,k2,percenterror2))
    plane_hops_2per_k3.append(five_hop_model_wrapped_cauchy(last_known_heading,pulauperak[1],pulauperak[0],km_hop,k3,percenterror2))

As they would say in Ye Olde English, "the loops runneth over." Once that's done, proceed with the below:¶

In [ ]:

# 5per_k1 run

first_lat_5per_k1 = []
two_lat_5per_k1 = []
three_lat_5per_k1 = []
four_lat_5per_k1 = []
final_lat_5per_k1 = []

first_lon_5per_k1 = []
two_lon_5per_k1 = []
three_lon_5per_k1 = []
four_lon_5per_k1 = []
final_lon_5per_k1 = []

route1_lat_5per_k1 = []
route2_lat_5per_k1 = []
route3_lat_5per_k1 = []
route4_lat_5per_k1 = []

route1_lon_5per_k1 = []
route2_lon_5per_k1 = []
route3_lon_5per_k1 = []
route4_lon_5per_k1 = []

for i in xrange(len(plane_hops_5per_k1)):
    first_lat_5per_k1.append(plane_hops_5per_k1[i][0][0])
    first_lon_5per_k1.append(plane_hops_5per_k1[i][1][0])
    two_lat_5per_k1.append(plane_hops_5per_k1[i][0][1])
    two_lon_5per_k1.append(plane_hops_5per_k1[i][1][1])
    three_lat_5per_k1.append(plane_hops_5per_k1[i][0][2])
    three_lon_5per_k1.append(plane_hops_5per_k1[i][1][2])
    four_lat_5per_k1.append(plane_hops_5per_k1[i][0][3])
    four_lon_5per_k1.append(plane_hops_5per_k1[i][1][3])
    final_lat_5per_k1.append(plane_hops_5per_k1[i][0][4])
    final_lon_5per_k1.append(plane_hops_5per_k1[i][1][4])
    
    route1_lat_5per_k1.append(plane_hops_5per_k1[i][0][5])
    route1_lon_5per_k1.append(plane_hops_5per_k1[i][1][5])
    route2_lat_5per_k1.append(plane_hops_5per_k1[i][0][6])
    route2_lon_5per_k1.append(plane_hops_5per_k1[i][1][6])
    route3_lat_5per_k1.append(plane_hops_5per_k1[i][0][7])
    route3_lon_5per_k1.append(plane_hops_5per_k1[i][1][7])
    route4_lat_5per_k1.append(plane_hops_5per_k1[i][0][8])
    route4_lon_5per_k1.append(plane_hops_5per_k1[i][1][8])

In [ ]:

# 5per_k2 run

first_lat_5per_k2 = []
two_lat_5per_k2 = []
three_lat_5per_k2 = []
four_lat_5per_k2 = []
final_lat_5per_k2 = []

first_lon_5per_k2 = []
two_lon_5per_k2 = []
three_lon_5per_k2 = []
four_lon_5per_k2 = []
final_lon_5per_k2 = []

route1_lat_5per_k2 = []
route2_lat_5per_k2 = []
route3_lat_5per_k2 = []
route4_lat_5per_k2 = []

route1_lon_5per_k2 = []
route2_lon_5per_k2 = []
route3_lon_5per_k2 = []
route4_lon_5per_k2 = []

for i in xrange(len(plane_hops_5per_k2)):
    first_lat_5per_k2.append(plane_hops_5per_k2[i][0][0])
    first_lon_5per_k2.append(plane_hops_5per_k2[i][1][0])
    two_lat_5per_k2.append(plane_hops_5per_k2[i][0][1])
    two_lon_5per_k2.append(plane_hops_5per_k2[i][1][1])
    three_lat_5per_k2.append(plane_hops_5per_k2[i][0][2])
    three_lon_5per_k2.append(plane_hops_5per_k2[i][1][2])
    four_lat_5per_k2.append(plane_hops_5per_k2[i][0][3])
    four_lon_5per_k2.append(plane_hops_5per_k2[i][1][3])
    final_lat_5per_k2.append(plane_hops_5per_k2[i][0][4])
    final_lon_5per_k2.append(plane_hops_5per_k2[i][1][4])
    
    route1_lat_5per_k2.append(plane_hops_5per_k2[i][0][5])
    route1_lon_5per_k2.append(plane_hops_5per_k2[i][1][5])
    route2_lat_5per_k2.append(plane_hops_5per_k2[i][0][6])
    route2_lon_5per_k2.append(plane_hops_5per_k2[i][1][6])
    route3_lat_5per_k2.append(plane_hops_5per_k2[i][0][7])
    route3_lon_5per_k2.append(plane_hops_5per_k2[i][1][7])
    route4_lat_5per_k2.append(plane_hops_5per_k2[i][0][8])
    route4_lon_5per_k2.append(plane_hops_5per_k2[i][1][8])

In [30]:

# 5per_k3 run

first_lat_5per_k3 = []
two_lat_5per_k3 = []
three_lat_5per_k3 = []
four_lat_5per_k3 = []
final_lat_5per_k3 = []

first_lon_5per_k3 = []
two_lon_5per_k3 = []
three_lon_5per_k3 = []
four_lon_5per_k3 = []
final_lon_5per_k3 = []

route1_lat_5per_k3 = []
route2_lat_5per_k3 = []
route3_lat_5per_k3 = []
route4_lat_5per_k3 = []

route1_lon_5per_k3 = []
route2_lon_5per_k3 = []
route3_lon_5per_k3 = []
route4_lon_5per_k3 = []

for i in xrange(len(plane_hops_5per_k3)):
    first_lat_5per_k3.append(plane_hops_5per_k3[i][0][0])
    first_lon_5per_k3.append(plane_hops_5per_k3[i][1][0])
    two_lat_5per_k3.append(plane_hops_5per_k3[i][0][1])
    two_lon_5per_k3.append(plane_hops_5per_k3[i][1][1])
    three_lat_5per_k3.append(plane_hops_5per_k3[i][0][2])
    three_lon_5per_k3.append(plane_hops_5per_k3[i][1][2])
    four_lat_5per_k3.append(plane_hops_5per_k3[i][0][3])
    four_lon_5per_k3.append(plane_hops_5per_k3[i][1][3])
    final_lat_5per_k3.append(plane_hops_5per_k3[i][0][4])
    final_lon_5per_k3.append(plane_hops_5per_k3[i][1][4])
    
    route1_lat_5per_k3.append(plane_hops_5per_k3[i][0][5])
    route1_lon_5per_k3.append(plane_hops_5per_k3[i][1][5])
    route2_lat_5per_k3.append(plane_hops_5per_k3[i][0][6])
    route2_lon_5per_k3.append(plane_hops_5per_k3[i][1][6])
    route3_lat_5per_k3.append(plane_hops_5per_k3[i][0][7])
    route3_lon_5per_k3.append(plane_hops_5per_k3[i][1][7])
    route4_lat_5per_k3.append(plane_hops_5per_k3[i][0][8])
    route4_lon_5per_k3.append(plane_hops_5per_k3[i][1][8])

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-30-9b461cf8f44c> in <module>()
     23 route4_lon_5per_k3 = []
     24 
---> 25 for i in xrange(len(plane_hops_5per_k3)):
     26     first_lat_5per_k3.append(plane_hops_5per_k3[i][0][0])
     27     first_lon_5per_k3.append(plane_hops_5per_k3[i][1][0])

NameError: name 'plane_hops_5per_k3' is not defined

In [ ]:

# 2per_k1 run

first_lat_2per_k1 = []
two_lat_2per_k1 = []
three_lat_2per_k1 = []
four_lat_2per_k1 = []
final_lat_2per_k1 = []

first_lon_2per_k1 = []
two_lon_2per_k1 = []
three_lon_2per_k1 = []
four_lon_2per_k1 = []
final_lon_2per_k1 = []

route1_lat_2per_k1 = []
route2_lat_2per_k1 = []
route3_lat_2per_k1 = []
route4_lat_2per_k1 = []

route1_lon_2per_k1 = []
route2_lon_2per_k1 = []
route3_lon_2per_k1 = []
route4_lon_2per_k1 = []

for i in xrange(len(plane_hops_2per_k1)):
    first_lat_2per_k1.append(plane_hops_2per_k1[i][0][0])
    first_lon_2per_k1.append(plane_hops_2per_k1[i][1][0])
    two_lat_2per_k1.append(plane_hops_2per_k1[i][0][1])
    two_lon_2per_k1.append(plane_hops_2per_k1[i][1][1])
    three_lat_2per_k1.append(plane_hops_2per_k1[i][0][2])
    three_lon_2per_k1.append(plane_hops_2per_k1[i][1][2])
    four_lat_2per_k1.append(plane_hops_2per_k1[i][0][3])
    four_lon_2per_k1.append(plane_hops_2per_k1[i][1][3])
    final_lat_2per_k1.append(plane_hops_2per_k1[i][0][4])
    final_lon_2per_k1.append(plane_hops_2per_k1[i][1][4])
    
    route1_lat_2per_k1.append(plane_hops_2per_k1[i][0][5])
    route1_lon_2per_k1.append(plane_hops_2per_k1[i][1][5])
    route2_lat_2per_k1.append(plane_hops_2per_k1[i][0][6])
    route2_lon_2per_k1.append(plane_hops_2per_k1[i][1][6])
    route3_lat_2per_k1.append(plane_hops_2per_k1[i][0][7])
    route3_lon_2per_k1.append(plane_hops_2per_k1[i][1][7])
    route4_lat_2per_k1.append(plane_hops_2per_k1[i][0][8])
    route4_lon_2per_k1.append(plane_hops_2per_k1[i][1][8])

In [ ]:

# 2per_k2 run

first_lat_2per_k2 = []
two_lat_2per_k2 = []
three_lat_2per_k2 = []
four_lat_2per_k2 = []
final_lat_2per_k2 = []

first_lon_2per_k2 = []
two_lon_2per_k2 = []
three_lon_2per_k2 = []
four_lon_2per_k2 = []
final_lon_2per_k2 = []

route1_lat_2per_k2 = []
route2_lat_2per_k2 = []
route3_lat_2per_k2 = []
route4_lat_2per_k2 = []

route1_lon_2per_k2 = []
route2_lon_2per_k2 = []
route3_lon_2per_k2 = []
route4_lon_2per_k2 = []

for i in xrange(len(plane_hops_2per_k2)):
    first_lat_2per_k2.append(plane_hops_2per_k2[i][0][0])
    first_lon_2per_k2.append(plane_hops_2per_k2[i][1][0])
    two_lat_2per_k2.append(plane_hops_2per_k2[i][0][1])
    two_lon_2per_k2.append(plane_hops_2per_k2[i][1][1])
    three_lat_2per_k2.append(plane_hops_2per_k2[i][0][2])
    three_lon_2per_k2.append(plane_hops_2per_k2[i][1][2])
    four_lat_2per_k2.append(plane_hops_2per_k2[i][0][3])
    four_lon_2per_k2.append(plane_hops_2per_k2[i][1][3])
    final_lat_2per_k2.append(plane_hops_2per_k2[i][0][4])
    final_lon_2per_k2.append(plane_hops_2per_k2[i][1][4])
    
    route1_lat_2per_k2.append(plane_hops_2per_k2[i][0][5])
    route1_lon_2per_k2.append(plane_hops_2per_k2[i][1][5])
    route2_lat_2per_k2.append(plane_hops_2per_k2[i][0][6])
    route2_lon_2per_k2.append(plane_hops_2per_k2[i][1][6])
    route3_lat_2per_k2.append(plane_hops_2per_k2[i][0][7])
    route3_lon_2per_k2.append(plane_hops_2per_k2[i][1][7])
    route4_lat_2per_k2.append(plane_hops_2per_k2[i][0][8])
    route4_lon_2per_k2.append(plane_hops_2per_k2[i][1][8])

In [ ]:

# 2per_k3 run

first_lat_2per_k3 = []
two_lat_2per_k3 = []
three_lat_2per_k3 = []
four_lat_2per_k3 = []
final_lat_2per_k3 = []

first_lon_2per_k3 = []
two_lon_2per_k3 = []
three_lon_2per_k3 = []
four_lon_2per_k3 = []
final_lon_2per_k3 = []

route1_lat_2per_k3 = []
route2_lat_2per_k3 = []
route3_lat_2per_k3 = []
route4_lat_2per_k3 = []

route1_lon_2per_k3 = []
route2_lon_2per_k3 = []
route3_lon_2per_k3 = []
route4_lon_2per_k3 = []

for i in xrange(len(plane_hops_2per_k3)):
    first_lat_2per_k3.append(plane_hops_2per_k3[i][0][0])
    first_lon_2per_k3.append(plane_hops_2per_k3[i][1][0])
    two_lat_2per_k3.append(plane_hops_2per_k3[i][0][1])
    two_lon_2per_k3.append(plane_hops_2per_k3[i][1][1])
    three_lat_2per_k3.append(plane_hops_2per_k3[i][0][2])
    three_lon_2per_k3.append(plane_hops_2per_k3[i][1][2])
    four_lat_2per_k3.append(plane_hops_2per_k3[i][0][3])
    four_lon_2per_k3.append(plane_hops_2per_k3[i][1][3])
    final_lat_2per_k3.append(plane_hops_2per_k3[i][0][4])
    final_lon_2per_k3.append(plane_hops_2per_k3[i][1][4])
    
    route1_lat_2per_k3.append(plane_hops_2per_k3[i][0][5])
    route1_lon_2per_k3.append(plane_hops_2per_k3[i][1][5])
    route2_lat_2per_k3.append(plane_hops_2per_k3[i][0][6])
    route2_lon_2per_k3.append(plane_hops_2per_k3[i][1][6])
    route3_lat_2per_k3.append(plane_hops_2per_k3[i][0][7])
    route3_lon_2per_k3.append(plane_hops_2per_k3[i][1][7])
    route4_lat_2per_k3.append(plane_hops_2per_k3[i][0][8])
    route4_lon_2per_k3.append(plane_hops_2per_k3[i][1][8])

Now we plot once again:¶

5per_k1 plot - Wrapped Cauchy Heading, k=0, 5% ping error¶

In [ ]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- 5% error
x5,y5 = fig(circle_lon,circle_lat)
x6,y6 = fig(circle_lon_err1_5per,circle_lat_err1_5per)
x7,y7 = fig(circle_lon_err2_5per,circle_lat_err2_5per)

#Add points after each hr
x9,y9 = fig(first_lon_5per_k1,first_lat_5per_k1)
x10,y10 = fig(two_lon_5per_k1,two_lat_5per_k1)
x11,y11 = fig(three_lon_5per_k1,three_lat_5per_k1)
x12,y12 = fig(four_lon_5per_k1,four_lat_5per_k1)
x13,y13 = fig(final_lon_5per_k1,final_lat_5per_k1)

#Add ultimate locations of MH370
x14,y14 = fig(route1_lon_5per_k1,route1_lat_5per_k1)
x15,y15 = fig(route2_lon_5per_k1,route2_lat_5per_k1)
x16,y16 = fig(route3_lon_5per_k1,route3_lat_5per_k1)
x17,y17 = fig(route4_lon_5per_k1,route4_lat_5per_k1)

#Draw circle showing extent of Inmarsat sat radar detection
fig.plot(x5,y5,'r-',markersize=5,label='5th Ping')
fig.plot(x6,y6,'r--',markersize=5,label='with 5% error')
fig.plot(x7,y7,'r--',markersize=5)

#Plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Add monte carlo points
fig.plot(x9,y9,'yo',markersize=5,label='after 1 hr')
fig.plot(x10,y10,'co',markersize=5,label='after 2 hrs')
fig.plot(x11,y11,'mo',markersize=5,label= 'after 3 hrs')
fig.plot(x12,y12,'wo',markersize=5,label='after 4 hrs')
fig.plot(x13,y13,'ro',markersize=7,label='after 5 hrs')

#Plot ultimate locations of MH370
fig.plot(x14,y14,'bo',markersize=5,label='in final hr')
fig.plot(x15,y15,'bo',markersize=5)
fig.plot(x16,y16,'bo',markersize=5)
fig.plot(x17,y17,'bo',markersize=5)

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='upper right',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('MH370 Position Progression Over Time', fontsize=30)

#Show below
plt.show()

5per_k2 plot - Wrapped Cauchy Heading, k=0.5, 5% ping error¶

In [ ]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- 5% error
x5,y5 = fig(circle_lon,circle_lat)
x6,y6 = fig(circle_lon_err1_5per,circle_lat_err1_5per)
x7,y7 = fig(circle_lon_err2_5per,circle_lat_err2_5per)

#Add points after each hr
x9,y9 = fig(first_lon_5per_k2,first_lat_5per_k2)
x10,y10 = fig(two_lon_5per_k2,two_lat_5per_k2)
x11,y11 = fig(three_lon_5per_k2,three_lat_5per_k2)
x12,y12 = fig(four_lon_5per_k2,four_lat_5per_k2)
x13,y13 = fig(final_lon_5per_k2,final_lat_5per_k2)

#Add ultimate locations of MH370
x14,y14 = fig(route1_lon_5per_k2,route1_lat_5per_k2)
x15,y15 = fig(route2_lon_5per_k2,route2_lat_5per_k2)
x16,y16 = fig(route3_lon_5per_k2,route3_lat_5per_k2)
x17,y17 = fig(route4_lon_5per_k2,route4_lat_5per_k2)

#Draw circle showing extent of Inmarsat sat radar detection
fig.plot(x5,y5,'r-',markersize=5,label='5th Ping')
fig.plot(x6,y6,'r--',markersize=5,label='with 5% error')
fig.plot(x7,y7,'r--',markersize=5)

#Plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Add monte carlo points
fig.plot(x9,y9,'yo',markersize=5,label='after 1 hr')
fig.plot(x10,y10,'co',markersize=5,label='after 2 hrs')
fig.plot(x11,y11,'mo',markersize=5,label= 'after 3 hrs')
fig.plot(x12,y12,'wo',markersize=5,label='after 4 hrs')
fig.plot(x13,y13,'ro',markersize=7,label='after 5 hrs')

#Plot ultimate locations of MH370
fig.plot(x14,y14,'bo',markersize=5,label='in final hr')
fig.plot(x15,y15,'bo',markersize=5)
fig.plot(x16,y16,'bo',markersize=5)
fig.plot(x17,y17,'bo',markersize=5)

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='upper right',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('MH370 Position Progression Over Time', fontsize=30)

#Show below
plt.show()

5per_k3 plot - Wrapped Cauchy Heading, k=0.99, 5% ping error¶

In [ ]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- 5% error
x5,y5 = fig(circle_lon,circle_lat)
x6,y6 = fig(circle_lon_err1_5per,circle_lat_err1_5per)
x7,y7 = fig(circle_lon_err2_5per,circle_lat_err2_5per)

#Add points after each hr
x9,y9 = fig(first_lon_5per_k3,first_lat_5per_k3)
x10,y10 = fig(two_lon_5per_k3,two_lat_5per_k3)
x11,y11 = fig(three_lon_5per_k3,three_lat_5per_k3)
x12,y12 = fig(four_lon_5per_k3,four_lat_5per_k3)
x13,y13 = fig(final_lon_5per_k3,final_lat_5per_k3)

#Add ultimate locations of MH370
x14,y14 = fig(route1_lon_5per_k3,route1_lat_5per_k3)
x15,y15 = fig(route2_lon_5per_k3,route2_lat_5per_k3)
x16,y16 = fig(route3_lon_5per_k3,route3_lat_5per_k3)
x17,y17 = fig(route4_lon_5per_k3,route4_lat_5per_k3)

#Draw circle showing extent of Inmarsat sat radar detection
fig.plot(x5,y5,'r-',markersize=5,label='5th Ping')
fig.plot(x6,y6,'r--',markersize=5,label='with 5% error')
fig.plot(x7,y7,'r--',markersize=5)

#Plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Add monte carlo points
fig.plot(x9,y9,'yo',markersize=5,label='after 1 hr')
fig.plot(x10,y10,'co',markersize=5,label='after 2 hrs')
fig.plot(x11,y11,'mo',markersize=5,label= 'after 3 hrs')
fig.plot(x12,y12,'wo',markersize=5,label='after 4 hrs')
fig.plot(x13,y13,'ro',markersize=7,label='after 5 hrs')

#Plot ultimate locations of MH370
fig.plot(x14,y14,'bo',markersize=5,label='in final hr')
fig.plot(x15,y15,'bo',markersize=5)
fig.plot(x16,y16,'bo',markersize=5)
fig.plot(x17,y17,'bo',markersize=5)

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='upper right',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('MH370 Position Progression Over Time', fontsize=30)

#Show below
plt.show()

2per_k1 plot - Wrapped Cauchy Heading, k=0, 2.5% ping error¶

In [ ]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])


#Add circle coords -- 2.5% error
x5,y5 = fig(circle_lon,circle_lat)
x6,y6 = fig(circle_lon_err1_2per,circle_lat_err1_2per)
x7,y7 = fig(circle_lon_err2_2per,circle_lat_err2_2per)

#Add points after each hr
x9,y9 = fig(first_lon_5per_k1,first_lat_5per_k1)
x10,y10 = fig(two_lon_5per_k1,two_lat_5per_k1)
x11,y11 = fig(three_lon_5per_k1,three_lat_5per_k1)
x12,y12 = fig(four_lon_5per_k1,four_lat_5per_k1)
x13,y13 = fig(final_lon_5per_k1,final_lat_5per_k1)

#Add ultimate locations of MH370
x14,y14 = fig(route1_lon_5per_k1,route1_lat_5per_k1)
x15,y15 = fig(route2_lon_5per_k1,route2_lat_5per_k1)
x16,y16 = fig(route3_lon_5per_k1,route3_lat_5per_k1)
x17,y17 = fig(route4_lon_5per_k1,route4_lat_5per_k1)

#Draw circle showing extent of Inmarsat sat radar detection
fig.plot(x5,y5,'r-',markersize=5,label='5th Ping')
fig.plot(x6,y6,'r--',markersize=5,label='with 5% error')
fig.plot(x7,y7,'r--',markersize=5)

#Plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Add monte carlo points
fig.plot(x9,y9,'yo',markersize=5,label='after 1 hr')
fig.plot(x10,y10,'co',markersize=5,label='after 2 hrs')
fig.plot(x11,y11,'mo',markersize=5,label= 'after 3 hrs')
fig.plot(x12,y12,'wo',markersize=5,label='after 4 hrs')
fig.plot(x13,y13,'ro',markersize=7,label='after 5 hrs')

#Plot ultimate locations of MH370
fig.plot(x14,y14,'bo',markersize=5,label='in final hr')
fig.plot(x15,y15,'bo',markersize=5)
fig.plot(x16,y16,'bo',markersize=5)
fig.plot(x17,y17,'bo',markersize=5)

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='upper right',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('MH370 Position Progression Over Time', fontsize=30)

#Show below
plt.show()

2per_k2 plot - Wrapped Cauchy Heading, k=0.5, 2.5% ping error¶

In [ ]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- 2.5% error
x5,y5 = fig(circle_lon,circle_lat)
x6,y6 = fig(circle_lon_err1_2per,circle_lat_err1_2per)
x7,y7 = fig(circle_lon_err2_2per,circle_lat_err2_2per)

#Add points after each hr
x9,y9 = fig(first_lon_5per_k2,first_lat_5per_k2)
x10,y10 = fig(two_lon_5per_k2,two_lat_5per_k2)
x11,y11 = fig(three_lon_5per_k2,three_lat_5per_k2)
x12,y12 = fig(four_lon_5per_k2,four_lat_5per_k2)
x13,y13 = fig(final_lon_5per_k2,final_lat_5per_k2)

#Add ultimate locations of MH370
x14,y14 = fig(route1_lon_5per_k2,route1_lat_5per_k2)
x15,y15 = fig(route2_lon_5per_k2,route2_lat_5per_k2)
x16,y16 = fig(route3_lon_5per_k2,route3_lat_5per_k2)
x17,y17 = fig(route4_lon_5per_k2,route4_lat_5per_k2)

#Draw circle showing extent of Inmarsat sat radar detection
fig.plot(x5,y5,'r-',markersize=5,label='5th Ping')
fig.plot(x6,y6,'r--',markersize=5,label='with 5% error')
fig.plot(x7,y7,'r--',markersize=5)

#Plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Add monte carlo points
fig.plot(x9,y9,'yo',markersize=5,label='after 1 hr')
fig.plot(x10,y10,'co',markersize=5,label='after 2 hrs')
fig.plot(x11,y11,'mo',markersize=5,label= 'after 3 hrs')
fig.plot(x12,y12,'wo',markersize=5,label='after 4 hrs')
fig.plot(x13,y13,'ro',markersize=7,label='after 5 hrs')

#Plot ultimate locations of MH370
fig.plot(x14,y14,'bo',markersize=5,label='in final hr')
fig.plot(x15,y15,'bo',markersize=5)
fig.plot(x16,y16,'bo',markersize=5)
fig.plot(x17,y17,'bo',markersize=5)

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='upper right',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('MH370 Position Progression Over Time', fontsize=30)

#Show below
plt.show()

2per_k3 plot - Wrapped Cauchy Heading, k=0.99, 2.5% ping error¶

In [ ]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- 2.5% error
x5,y5 = fig(circle_lon,circle_lat)
x6,y6 = fig(circle_lon_err1_2per,circle_lat_err1_2per)
x7,y7 = fig(circle_lon_err2_2per,circle_lat_err2_2per)

#Add points after each hr
x9,y9 = fig(first_lon_5per_k3,first_lat_5per_k3)
x10,y10 = fig(two_lon_5per_k3,two_lat_5per_k3)
x11,y11 = fig(three_lon_5per_k3,three_lat_5per_k3)
x12,y12 = fig(four_lon_5per_k3,four_lat_5per_k3)
x13,y13 = fig(final_lon_5per_k3,final_lat_5per_k3)

#Add ultimate locations of MH370
x14,y14 = fig(route1_lon_5per_k3,route1_lat_5per_k3)
x15,y15 = fig(route2_lon_5per_k3,route2_lat_5per_k3)
x16,y16 = fig(route3_lon_5per_k3,route3_lat_5per_k3)
x17,y17 = fig(route4_lon_5per_k3,route4_lat_5per_k3)

#Draw circle showing extent of Inmarsat sat radar detection
fig.plot(x5,y5,'r-',markersize=5,label='5th Ping')
fig.plot(x6,y6,'r--',markersize=5,label='with 5% error')
fig.plot(x7,y7,'r--',markersize=5)

#Plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Add monte carlo points
fig.plot(x9,y9,'yo',markersize=5,label='after 1 hr')
fig.plot(x10,y10,'co',markersize=5,label='after 2 hrs')
fig.plot(x11,y11,'mo',markersize=5,label= 'after 3 hrs')
fig.plot(x12,y12,'wo',markersize=5,label='after 4 hrs')
fig.plot(x13,y13,'ro',markersize=7,label='after 5 hrs')

#Plot ultimate locations of MH370
fig.plot(x14,y14,'bo',markersize=5,label='in final hr')
fig.plot(x15,y15,'bo',markersize=5)
fig.plot(x16,y16,'bo',markersize=5)
fig.plot(x17,y17,'bo',markersize=5)

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='upper right',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('MH370 Position Progression Over Time', fontsize=30)

#Show below
plt.show()

Intriguingly, even with an extremely narrow distribution the Wrapped Cauchy is of little use in determining where plane is. We'll try it again with all of the pings. With this heading distribution we can now say, "it's a wrap!"¶

Conclusions and Comparisons to First Version of the Monte Carlo Model:¶

I cleaned up the code, fixed a few (minor to the results) bugs, used more realistic Inmarsat ping error values based on new information (reducing the error), and used two different distributions for picking the heading besides normal. Remember, the heading is my main state, so changing how the heading is chosen can affect the model outcome, as we saw from the plots.¶

The most realistic intermediary k value cases for Von Mises show the same results as the normal distribution option -- MH370 ends up off the coast of Australia in most cases. The Wrapped Cauchy, only just using the 5th ping arc, does not make meaningful inferences about where the plane is.¶

So, what we see here is accurate for the final location of where MH370 is most likely to be located. But, one of the advantages of having a 10% standard deviation error (or even 20%) is that with our normal probability distribution established from a grid with that amount of ping arc error, we encompassed many pings. So, the plane didn't just travel along the arc. The smaller the error, the more the plane will travel along the arc. Now, we know it was on the arc at the 5th ping, so this is not the end of the world. But even better would be to try to recreate the overall trajectory of the aircraft at all points in time. In that way, if we looked at varying the speed, we could just look further to the north along the southern corridor instead of varying the speed (in other words, we'd do the calculation by hand, and say, at a high speed, the plane would just go down where it was at the 5th ping -- and would run out of fuel almost immediately thereafter.)¶

Next Step - Adding the Other Pings¶

We've now incorporated the reasonable amount of error Inmarsat would expect from their ping estimate -- just 2.5 or 5% -- but we haven't incorporated some idea of where the other pings have been located. So, let's do that.¶

Unfortunately we don't have the actual numbers. But, from scrutinizing infographics, some of which may have had that information, in addition to graphics produced by people "in the know" with that info (such as the NTSB, National Transporatation Safety Board), we can infer a pretty good estimate of where the other pings are. These estimates ultimately can't vary by much, as "what you see is what you get" on an accurate map, which most major media organizations have with multiple pings. Let's do that and see how it affects the trajectory of MH370.¶

These were actually determined in a separate Mathematica 9 notebook (like the spherical trig derivations, as well as the double-checking of the Python implementation of this model. Same results and output, by the way.) However, I'll grab the radius values of the other pings and plot them, and them include links to infographic maps showing the ping locations, so you can see for yourself that these are close to what the actual ping values must be:¶

In [ ]:

#add other ping distances, and then replot

ping_distances = np.array([4036.99, 4194.65, 4352.32, 4509.99, 4667.65, 4825.32])
ping_times = np.array([0.9333, 1, 1, 1, 1, 1]) #make 1st hop slightly smaller owing to time difference
ping_arcs = np.array([34.8485, 36.2649, 37.6812, 39.0976, 40.5139, 41.9303, 43.3466])

WHAT?! You say. What about the fact that you've repeatedly said there were 5 pings afterward. Now there are 6?! Well, this is one of many contradictory pieces of information which has emerged both to the media at large, and to myself as I continue to search for the best sources of information. If you look at the times, apparently, 6 out of the 7 pings happened after the plane disappeared, not 5. A Wall Street Journal article which was picked up in numerous other newspapers and media outlets, such as Yahoo News and several in Malaysia (who you'd think be close to the matter with easier access to officials speaking off the record, for instance), said: "Malaysia Airlines' missing jet transmitted its location repeatedly to satellites over the course of five hours after it disappeared from radar, people briefed on the matter said...". This has never been corrected or further addressed, even though other corrections have been run for facts that are in dispute (such as whether the Rolls Royce engine communication system was responsible for the signal routing through Inmarsat, versus the "the idling satellite communications system" which now appears to be the case, for example.) (Of course, here's a 4 hour claim to further complicate the matter.)¶

Here's an article, published right after the first version of my model finished, with a recap of what we know MH370's location over time -- essentially that the first ping was at 2:11 am -- which was before MH370's disappearance from radar -- and that the last was at 8:11 am.¶

Nevertheless, the last ping that was plotted originally was the last ping. It's just that the timing must have been off, somehow, given the current timeline:¶

[2:11 am], 3:11 am, 4:11 am, 5:11 am, 6:11 am, 7:11 am, 8:11 am, [plus a "partial handshake" which isn't understood]

That's 6 after disappearance, and that's the final count that I'll deal with in my model; the previous iterations of the model used the best information known at the time. So I'll leave those alone as is -- because effectively what I am doing over the course of 5 hrs, instead of 6 hours, is using 5/6 of the cruising speed which is possible. The last ping's location, and the ping arc, remains the same. None of that changes. But in being in the air longer to achieve that arc, the plane can fly further south.¶

An inconsistency this leaves us with, by the way, is that the search is now taken place further north, because the "airplane was going faster and ran out of fuel earlier" -- but that does not reconcile with amount of time it would be in the air for (6+ hrs instead of 5+ hrs.) So I still do not understand the reasoning behind that argument.¶

Replotting the Ping on the Map¶

We'll use the ping distances from the array and Mathematica above, to confirm that they look sensible compared to the infographics out there. We could use a 'for' loop but there are only 6, so I'll do them one-by-one:¶

In [ ]:

#make points for 6 circles -- opt not to use for loop
ping_circle_211am = make_circle(ping_arcs[0],360,64.5,0)
ping_circle_311am = make_circle(ping_arcs[1],360,64.5,0)
ping_circle_411am = make_circle(ping_arcs[2],360,64.5,0)
ping_circle_511am = make_circle(ping_arcs[3],360,64.5,0)
ping_circle_611am = make_circle(ping_arcs[4],360,64.5,0)
ping_circle_711am = make_circle(ping_arcs[5],360,64.5,0)
ping_circle_811am = make_circle(ping_arcs[6],360,64.5,0)

#initialize lat & lon lists
circle_lon_211am = []
circle_lat_211am = []
circle_lat_311am = []
circle_lon_311am = []
circle_lat_411am = []
circle_lon_411am = []
circle_lat_511am = []
circle_lon_511am = []
circle_lat_611am = []
circle_lon_611am = []
circle_lat_711am = []
circle_lon_711am = []
circle_lat_811am = []
circle_lon_811am = []

for i in xrange(len(ping_circle_211am)): #they're all the same length so just do it once
    circle_lat_211am.append(ping_circle_211am[i][0])
    circle_lon_211am.append(ping_circle_211am[i][1])
    
    circle_lat_311am.append(ping_circle_311am[i][0])
    circle_lon_311am.append(ping_circle_311am[i][1])

    circle_lat_411am.append(ping_circle_411am[i][0])
    circle_lon_411am.append(ping_circle_411am[i][1])

    circle_lat_511am.append(ping_circle_511am[i][0])
    circle_lon_511am.append(ping_circle_511am[i][1])

    circle_lat_611am.append(ping_circle_611am[i][0])
    circle_lon_611am.append(ping_circle_611am[i][1])

    circle_lat_711am.append(ping_circle_711am[i][0])
    circle_lon_711am.append(ping_circle_711am[i][1])

    circle_lat_811am.append(ping_circle_811am[i][0])
    circle_lon_811am.append(ping_circle_811am[i][1])

In [31]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Draw great circle to show path autopilot would have taken
fig.drawgreatcircle(pulauperak[1],pulauperak[0],85,-40,linewidth=3,color='black',label='Great Circle Path')

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- these are for each ping. will not plot the 2.5 and 5% error
x5,y5 = fig(circle_lon_211am,circle_lat_211am)
x6,y6 = fig(circle_lon_311am,circle_lat_311am)
x7,y7 = fig(circle_lon_411am,circle_lat_411am)
x8,y8 = fig(circle_lon_511am,circle_lat_511am)
x9,y9 = fig(circle_lon_611am,circle_lat_611am)
x10,y10 = fig(circle_lon_711am,circle_lat_711am)
x11,y11 = fig(circle_lon_811am,circle_lat_811am)

#Draw circle showing extent of Inmarsat sat radar detection for each of the pings
fig.plot(x5,y5,'r--',markersize=5,label='1st Ping Arc')
fig.plot(x6,y6,'r-',markersize=5, label='Ping Arcs After Disappearance')
fig.plot(x7,y7,'r-',markersize=5)
fig.plot(x8,y8,'r-',markersize=5)
fig.plot(x9,y9,'r-',markersize=5)
fig.plot(x10,y10,'r-',markersize=5)
fig.plot(x11,y11,'r-',markersize=5)

# plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='upper right',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('Inmarsat Ping Estimation -- Individual Pings', fontsize=30)

#Show below
plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-31-47b9268a98e6> in <module>()
      1 #Set figure size
----> 2 fig = plt.figure(figsize=[30,20])
      3 
      4 #Setup Basemap
      5 fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

NameError: name 'plt' is not defined

Notice how the very first one goes through the last point. Technically it's ~59 miles before -- but we can't see at that resolution. The first ping is dashed since it's before the plane disappeared, and it's right where we saw the last ping. For the remaining plots I've plotted a great circle line from the plane's last known location, to 40 S 85 E. The fact that it does not quite meet the last ping is not relevant to the conlusions from the below analysis (I reran the notebook to be able to type comments such as this.)¶

Running All of the Above with 5 -- (No Wait is it 6?!) -- Individual Satellite Pings¶

The Ping Probability Grid is kept at a 10x multiplication factor. I re-run all of the heading options with the most likely standard deviation or k value that represents a plane that is trying to get somewhere, but can change direction to some degree. (Again, if you take issue with this, this code is downloadable and re-runnable with different parameters. It's easy to plug and play. You only might have to wait a while for the simulations to finish running.)¶

I will run a normal heading distribution with std dev 30 (wide), von mises distribution with k=10 (medium, but narrower allowance for turns), and finally wrapped cauchy distribution with k=0.99 (almost a point; exceedingly small probability of turning from its previous heading) to capture a range of flying behaviors for our missing airplane, MH370. This is the final analysis undertaken for this model, with the most complete information as of March 30, 2014.¶

Note that one assumption I make which is a slight oversimplification, is that I assume 2.5-5% error for each of the pings even though technically there should be less standard deviation error for the earlier pings because they are closer to the satellite. I do not believe this is significant to the results (or else I would rectify it), but is worth stating. [Update: after running the results, it's actually a stronger argument for them. For my assumption introduces a sligthly wider range of paths MH370 could travel over; and there is still no way it can take a Great Circle route or travel by magnetic bearing -- which the 777 autopilot does by default in between set waypoints -- to locations off the coast of Australia which my model and Inmarsat believe it to be near.)¶

In [ ]:

"""
a function which given a list of discrete probabilities for each destination point, 
it will choose one of those points.

heading_init -- initial direction was headed at last known point
lon_init,lat_init -- last known point of plane in longitude and latitude
km_hop -- how far the plane went in the time interval, 1 hr. So in simplest case, the 777's cruising speed/hr.
std_dev -- the standard deviation of the heading, based on a normal distribution from the current heading (0 deg).
ping_percent_err -- what % error you assume in the Inmarsat 5th ping. either 2.5 or 5%.

uses normal distribution for heading

replace "dist_from_sat" with "ping_distance" since that's changing. run 6 times.

"""
def six_hop_model_normal(heading_init,lon_init,lat_init,km_hop,std_dev,ping_percent_err,ping_distances,ping_times):   
    
    #initialize
    plane_lat = np.zeros(6) #initialize plane location after each hop (assumed to be 1 hr for now)
    plane_lon = np.zeros(6)  
    lat = lat_init
    lon = lon_init
    heading = heading_init
    
    for i in xrange(len(plane_lat)):
        new_circle,new_weights,new_angles = normal_prob_step(heading,std_dev,lon,lat,(km_hop/eq_deg_km)*ping_times[i])
        #new_circle gives up possible coords for diff headings
        
        raw_weights = np.zeros(len(new_circle))
        for j in xrange(len(new_circle)):
            raw_weights[j] = new_weights[j]*ping_prob_normal(inmarsat[0],inmarsat[1],new_circle[j][0],new_circle[j][1],ping_percent_err,ping_distances[i],earth_radius) 
        
        probs = raw_weights / np.sum(raw_weights) #normalize    
        
        index = range(len(new_circle))
        chosen = np.random.choice(index,size=None,p=probs)
        #print "chosen",chosen
        
        heading = new_angles[chosen] #update heading
        
        plane_lat[i],plane_lon[i] = new_circle[chosen] #update position
        lat = plane_lat[i]
        lon = plane_lon[i]
    
    #at end of simulation, run the last location & heading for plane for 4 different times
    route1 = make_vector(0.25*km_hop/eq_deg_km,heading,lon,lat)
    route2 = make_vector(0.5*km_hop/eq_deg_km,heading,lon,lat)
    route3 = make_vector(0.75*km_hop/eq_deg_km,heading,lon,lat)
    route4 = make_vector((59./60.)*km_hop/eq_deg_km,heading,lon,lat)

    new_plane_lat = np.zeros(10)
    new_plane_lon = np.zeros(10)
    
    for i in xrange(len(plane_lat)):
        new_plane_lat[i] = plane_lat[i]
        new_plane_lon[i] = plane_lon[i]
    
    new_plane_lat[6] = route1[0] # add 1 for 6 hops instead of 5
    new_plane_lat[7] = route2[0] # add 1 for 6 hops instead of 5
    new_plane_lat[8] = route3[0] # add 1 for 6 hops instead of 5
    new_plane_lat[9] = route4[0] # add 1 for 6 hops instead of 5
    new_plane_lon[6] = route1[1] # add 1 for 6 hops instead of 5
    new_plane_lon[7] = route2[1] # add 1 for 6 hops instead of 5
    new_plane_lon[8] = route3[1] # add 1 for 6 hops instead of 5
    new_plane_lon[9] = route4[1] # add 1 for 6 hops instead of 5
    
    return new_plane_lat,new_plane_lon

In [ ]:

"""
a function which given a list of discrete probabilities for each destination point, 
it will choose one of those points.

heading_init -- initial direction was headed at last known point
lon_init,lat_init -- last known point of plane in longitude and latitude
km_hop -- how far the plane went in the time interval, 1 hr. So in simplest case, the 777's cruising speed/hr.
k -- affects the heading distribution, based on a Von Mises distribution from the current heading (0 deg).
ping_percent_err -- what % error you assume in the Inmarsat 5th ping. either 2.5 or 5%.

uses Von Mises distribution for heading

replace "dist_from_sat" with "ping_distance" since that's changing. run 6 times.

"""
def six_hop_model_von_mises(heading_init,lon_init,lat_init,km_hop,k,ping_percent_err,ping_distances,ping_times):   
    
    #initialize
    plane_lat = np.zeros(6) #initialize plane location after each hop (assumed to be 1 hr for now)
    plane_lon = np.zeros(6)  
    lat = lat_init
    lon = lon_init
    heading = heading_init
    
    for i in xrange(len(plane_lat)):
        new_circle,new_weights,new_angles = von_mises_prob_step(heading,k,lon,lat,(km_hop/eq_deg_km)*ping_times[i])
        #new_circle gives up possible coords for diff headings
        
        raw_weights = np.zeros(len(new_circle))
        for j in xrange(len(new_circle)):
            raw_weights[j] = new_weights[j]*ping_prob_normal(inmarsat[0],inmarsat[1],new_circle[j][0],new_circle[j][1],ping_percent_err,ping_distances[i],earth_radius) 
        
        probs = raw_weights / np.sum(raw_weights) #normalize    
        
        index = range(len(new_circle))
        chosen = np.random.choice(index,size=None,p=probs)
        #print "chosen",chosen
        
        heading = new_angles[chosen] #update heading
        
        plane_lat[i],plane_lon[i] = new_circle[chosen] #update position
        lat = plane_lat[i]
        lon = plane_lon[i]
    
    #at end of simulation, run the last location & heading for plane for 4 different times
    route1 = make_vector(0.25*km_hop/eq_deg_km,heading,lon,lat)
    route2 = make_vector(0.5*km_hop/eq_deg_km,heading,lon,lat)
    route3 = make_vector(0.75*km_hop/eq_deg_km,heading,lon,lat)
    route4 = make_vector((59./60.)*km_hop/eq_deg_km,heading,lon,lat)

    new_plane_lat = np.zeros(10)
    new_plane_lon = np.zeros(10)
    
    for i in xrange(len(plane_lat)):
        new_plane_lat[i] = plane_lat[i]
        new_plane_lon[i] = plane_lon[i]
    
    new_plane_lat[6] = route1[0] # add 1 for 6 hops instead of 5
    new_plane_lat[7] = route2[0] # add 1 for 6 hops instead of 5
    new_plane_lat[8] = route3[0] # add 1 for 6 hops instead of 5
    new_plane_lat[9] = route4[0] # add 1 for 6 hops instead of 5
    new_plane_lon[6] = route1[1] # add 1 for 6 hops instead of 5
    new_plane_lon[7] = route2[1] # add 1 for 6 hops instead of 5
    new_plane_lon[8] = route3[1] # add 1 for 6 hops instead of 5
    new_plane_lon[9] = route4[1] # add 1 for 6 hops instead of 5
    
    return new_plane_lat,new_plane_lon

In [ ]:

"""
a function which given a list of discrete probabilities for each destination point, 
it will choose one of those points.

heading_init -- initial direction was headed at last known point
lon_init,lat_init -- last known point of plane in longitude and latitude
km_hop -- how far the plane went in the time interval, 1 hr. So in simplest case, the 777's cruising speed/hr.
k -- affects the heading distribution, based on a Wrapped Cauchy distribution from the current heading (0 deg).
ping_percent_err -- what % error you assume in the Inmarsat 5th ping. either 2.5 or 5%.

uses Wrapped Cauchy distribution for heading

replace "dist_from_sat" with "ping_distance" since that's changing. run 6 times.

"""
def six_hop_model_wrapped_cauchy(heading_init,lon_init,lat_init,km_hop,k,ping_percent_err,ping_distances,ping_times):   
    
    #initialize
    plane_lat = np.zeros(6) #initialize plane location after each hop (assumed to be 1 hr for now)
    plane_lon = np.zeros(6)  
    lat = lat_init
    lon = lon_init
    heading = heading_init
    
    for i in xrange(len(plane_lat)):
        new_circle,new_weights,new_angles = wrapped_cauchy_prob_step(heading,k,lon,lat,(km_hop/eq_deg_km)*ping_times[i])
        #new_circle gives up possible coords for diff headings
        
        raw_weights = np.zeros(len(new_circle))
        for j in xrange(len(new_circle)):
            raw_weights[j] = new_weights[j]*ping_prob_normal(inmarsat[0],inmarsat[1],new_circle[j][0],new_circle[j][1],ping_percent_err,ping_distances[i],earth_radius) 
        
        probs = raw_weights / np.sum(raw_weights) #normalize    
        
        index = range(len(new_circle))
        chosen = np.random.choice(index,size=None,p=probs)
        #print "chosen",chosen
        
        heading = new_angles[chosen] #update heading
        
        plane_lat[i],plane_lon[i] = new_circle[chosen] #update position
        lat = plane_lat[i]
        lon = plane_lon[i]
    
    #at end of simulation, run the last location & heading for plane for 4 different times
    route1 = make_vector(0.25*km_hop/eq_deg_km,heading,lon,lat)
    route2 = make_vector(0.5*km_hop/eq_deg_km,heading,lon,lat)
    route3 = make_vector(0.75*km_hop/eq_deg_km,heading,lon,lat)
    route4 = make_vector((59./60.)*km_hop/eq_deg_km,heading,lon,lat)

    new_plane_lat = np.zeros(10)
    new_plane_lon = np.zeros(10)
    
    for i in xrange(len(plane_lat)):
        new_plane_lat[i] = plane_lat[i]
        new_plane_lon[i] = plane_lon[i]
    
    new_plane_lat[6] = route1[0] # add 1 for 6 hops instead of 5
    new_plane_lat[7] = route2[0] # add 1 for 6 hops instead of 5
    new_plane_lat[8] = route3[0] # add 1 for 6 hops instead of 5
    new_plane_lat[9] = route4[0] # add 1 for 6 hops instead of 5
    new_plane_lon[6] = route1[1] # add 1 for 6 hops instead of 5
    new_plane_lon[7] = route2[1] # add 1 for 6 hops instead of 5
    new_plane_lon[8] = route3[1] # add 1 for 6 hops instead of 5
    new_plane_lon[9] = route4[1] # add 1 for 6 hops instead of 5
    
    return new_plane_lat,new_plane_lon

In [ ]:

last_known_heading = 255.136 #calculated in Mathematica from MH370's two last publically known locations:
                             #when it deviated from its flight path, and when it was last detected by Malyasian military radar
                             #0 degrees is due north, so this is basically to the west (270 degrees), but slightly south

km_hop = 905 #assuming 1 hr intervals, at 905 km/hr which is 777 cruising speed -- use for test case
             # max speed of a Boeing 777 is 950 km/hr FYI

N = 1000 #define number of simulations to run

Normal Distribution Heading¶

Standard Deviation of 30 degrees (allows for some turning each time)¶

In [ ]:

percenterror1,percenterror2 = 0.05, 0.025

std_dev = 30

In [ ]:

plane_hops_5per = []
plane_hops_2per = []

for i in xrange(N):
    plane_hops_5per.append(six_hop_model_normal(last_known_heading,pulauperak[1],pulauperak[0],km_hop,std_dev,percenterror1,ping_distances,ping_times))
    plane_hops_2per.append(six_hop_model_normal(last_known_heading,pulauperak[1],pulauperak[0],km_hop,std_dev,percenterror2,ping_distances,ping_times))

In [ ]:

first_lat_5per = []
two_lat_5per = []
three_lat_5per = []
four_lat_5per = []
five_lat_5per = []
final_lat_5per = []

first_lon_5per = []
two_lon_5per = []
three_lon_5per = []
four_lon_5per = []
five_lon_5per = []
final_lon_5per = []

route1_lat_5per = []
route2_lat_5per = []
route3_lat_5per = []
route4_lat_5per = []

route1_lon_5per = []
route2_lon_5per = []
route3_lon_5per = []
route4_lon_5per = []

for i in xrange(len(plane_hops_5per)):
    first_lat_5per.append(plane_hops_5per[i][0][0])
    first_lon_5per.append(plane_hops_5per[i][1][0])
    two_lat_5per.append(plane_hops_5per[i][0][1])
    two_lon_5per.append(plane_hops_5per[i][1][1])
    three_lat_5per.append(plane_hops_5per[i][0][2])
    three_lon_5per.append(plane_hops_5per[i][1][2])
    four_lat_5per.append(plane_hops_5per[i][0][3])
    four_lon_5per.append(plane_hops_5per[i][1][3])
    five_lat_5per.append(plane_hops_5per[i][0][4])
    five_lon_5per.append(plane_hops_5per[i][1][4])
    final_lat_5per.append(plane_hops_5per[i][0][5])
    final_lon_5per.append(plane_hops_5per[i][1][5])
    
    route1_lat_5per.append(plane_hops_5per[i][0][6])
    route1_lon_5per.append(plane_hops_5per[i][1][6])
    route2_lat_5per.append(plane_hops_5per[i][0][7])
    route2_lon_5per.append(plane_hops_5per[i][1][7])
    route3_lat_5per.append(plane_hops_5per[i][0][8])
    route3_lon_5per.append(plane_hops_5per[i][1][8])
    route4_lat_5per.append(plane_hops_5per[i][0][9])
    route4_lon_5per.append(plane_hops_5per[i][1][9])

In [ ]:

first_lat_2per = []
two_lat_2per = []
three_lat_2per = []
four_lat_2per = []
five_lat_2per = []
final_lat_2per = []

first_lon_2per = []
two_lon_2per = []
three_lon_2per = []
four_lon_2per = []
five_lon_2per = []
final_lon_2per = []

route1_lat_2per = []
route2_lat_2per = []
route3_lat_2per = []
route4_lat_2per = []

route1_lon_2per = []
route2_lon_2per = []
route3_lon_2per = []
route4_lon_2per = []

for i in xrange(len(plane_hops_2per)):
    first_lat_2per.append(plane_hops_2per[i][0][0])
    first_lon_2per.append(plane_hops_2per[i][1][0])
    two_lat_2per.append(plane_hops_2per[i][0][1])
    two_lon_2per.append(plane_hops_2per[i][1][1])
    three_lat_2per.append(plane_hops_2per[i][0][2])
    three_lon_2per.append(plane_hops_2per[i][1][2])
    four_lat_2per.append(plane_hops_2per[i][0][3])
    four_lon_2per.append(plane_hops_2per[i][1][3])
    five_lat_2per.append(plane_hops_2per[i][0][4])
    five_lon_2per.append(plane_hops_2per[i][1][4])
    final_lat_2per.append(plane_hops_2per[i][0][5])
    final_lon_2per.append(plane_hops_2per[i][1][5])
    
    route1_lat_2per.append(plane_hops_2per[i][0][6])
    route1_lon_2per.append(plane_hops_2per[i][1][6])
    route2_lat_2per.append(plane_hops_2per[i][0][7])
    route2_lon_2per.append(plane_hops_2per[i][1][7])
    route3_lat_2per.append(plane_hops_2per[i][0][8])
    route3_lon_2per.append(plane_hops_2per[i][1][8])
    route4_lat_2per.append(plane_hops_2per[i][0][9])
    route4_lon_2per.append(plane_hops_2per[i][1][9])

5% results:¶

In [ ]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Draw great circle to show path autopilot would have taken
fig.drawgreatcircle(pulauperak[1],pulauperak[0],85,-40,linewidth=3,color='black',label='Great Circle Path')

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- these are for each ping. will not plot the 2.5 and 5% error
x5,y5 = fig(circle_lon_211am,circle_lat_211am)
x6,y6 = fig(circle_lon_311am,circle_lat_311am)
x7,y7 = fig(circle_lon_411am,circle_lat_411am)
x8,y8 = fig(circle_lon_511am,circle_lat_511am)
x9,y9 = fig(circle_lon_611am,circle_lat_611am)
x10,y10 = fig(circle_lon_711am,circle_lat_711am)
x11,y11 = fig(circle_lon_811am,circle_lat_811am)

#Add points after each hr
x12,y12 = fig(first_lon_5per,first_lat_5per)
x13,y13 = fig(two_lon_5per,two_lat_5per)
x14,y14 = fig(three_lon_5per,three_lat_5per)
x15,y15 = fig(four_lon_5per,four_lat_5per)
x16,y16 = fig(five_lon_5per,five_lat_5per)
x17,y17 = fig(final_lon_5per,final_lat_5per)

#Add ultimate locations of MH370
x18,y18 = fig(route1_lon_5per,route1_lat_5per)
x19,y19 = fig(route2_lon_5per,route2_lat_5per)
x20,y20 = fig(route3_lon_5per,route3_lat_5per)
x21,y21 = fig(route4_lon_5per,route4_lat_5per)

# plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Draw circle showing extent of Inmarsat sat radar detection for each of the pings
fig.plot(x5,y5,'r--',markersize=5,label='1st Ping Arc')
fig.plot(x6,y6,'r-',markersize=5, label='Ping Arcs After Disappearance')
fig.plot(x7,y7,'r-',markersize=5)
fig.plot(x8,y8,'r-',markersize=5)
fig.plot(x9,y9,'r-',markersize=5)
fig.plot(x10,y10,'r-',markersize=5)
fig.plot(x11,y11,'r-',markersize=5)

#Add monte carlo points
fig.plot(x12,y12,'yo',markersize=5,label='after 1 hrs')
fig.plot(x13,y13,'mo',markersize=5,label= 'after 2 hrs')
fig.plot(x14,y14,'wo',markersize=5,label='after 3 hrs')
fig.plot(x15,y15,'bo',markersize=5,label='after 4 hrs')
fig.plot(x16,y16,'go',markersize=5,label='after 5 hrs')
fig.plot(x17,y17,'ro',markersize=7,label='after 6 hrs')

#Plot ultimate locations of MH370
fig.plot(x18,y18,'bo',markersize=5,label='in final hr')
fig.plot(x19,y19,'bo',markersize=5)
fig.plot(x20,y20,'bo',markersize=5)
fig.plot(x21,y21,'bo',markersize=5)

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='upper right',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('Inmarsat Ping Estimation -- Individual Pings', fontsize=30)

#Show below
plt.show()

2.5% results:¶

In [ ]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Draw great circle to show path autopilot would have taken
fig.drawgreatcircle(pulauperak[1],pulauperak[0],85,-40,linewidth=3,color='black',label='Great Circle Path')

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- these are for each ping. will not plot the 2.5 and 5% error
x5,y5 = fig(circle_lon_211am,circle_lat_211am)
x6,y6 = fig(circle_lon_311am,circle_lat_311am)
x7,y7 = fig(circle_lon_411am,circle_lat_411am)
x8,y8 = fig(circle_lon_511am,circle_lat_511am)
x9,y9 = fig(circle_lon_611am,circle_lat_611am)
x10,y10 = fig(circle_lon_711am,circle_lat_711am)
x11,y11 = fig(circle_lon_811am,circle_lat_811am)

#Add points after each hr
x12,y12 = fig(first_lon_2per,first_lat_2per)
x13,y13 = fig(two_lon_2per,two_lat_2per)
x14,y14 = fig(three_lon_2per,three_lat_2per)
x15,y15 = fig(four_lon_2per,four_lat_2per)
x16,y16 = fig(five_lon_2per,five_lat_2per)
x17,y17 = fig(final_lon_2per,final_lat_2per)

#Add ultimate locations of MH370
x18,y18 = fig(route1_lon_2per,route1_lat_2per)
x19,y19 = fig(route2_lon_2per,route2_lat_2per)
x20,y20 = fig(route3_lon_2per,route3_lat_2per)
x21,y21 = fig(route4_lon_2per,route4_lat_2per)

# plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Draw circle showing extent of Inmarsat sat radar detection for each of the pings
fig.plot(x5,y5,'r--',markersize=5,label='1st Ping Arc')
fig.plot(x6,y6,'r-',markersize=5, label='Ping Arcs After Disappearance')
fig.plot(x7,y7,'r-',markersize=5)
fig.plot(x8,y8,'r-',markersize=5)
fig.plot(x9,y9,'r-',markersize=5)
fig.plot(x10,y10,'r-',markersize=5)
fig.plot(x11,y11,'r-',markersize=5)

#Add monte carlo points
fig.plot(x12,y12,'yo',markersize=5,label='after 1 hrs')
fig.plot(x13,y13,'mo',markersize=5,label= 'after 2 hrs')
fig.plot(x14,y14,'wo',markersize=5,label='after 3 hrs')
fig.plot(x15,y15,'bo',markersize=5,label='after 4 hrs')
fig.plot(x16,y16,'go',markersize=5,label='after 5 hrs')
fig.plot(x17,y17,'ro',markersize=7,label='after 6 hrs')

#Plot ultimate locations of MH370
fig.plot(x18,y18,'bo',markersize=5,label='in final hr')
fig.plot(x19,y19,'bo',markersize=5)
fig.plot(x20,y20,'bo',markersize=5)
fig.plot(x21,y21,'bo',markersize=5)

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='upper right',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('Inmarsat Ping Estimation -- Individual Pings', fontsize=30)

#Show below
plt.show()

The normal distribution, with individual pings and 2.5-5% std dev error for each ping, demonstrates a narrow range of possible flight paths to an eventual endpoint in the Southern Ocean, roughly 40-45 degrees south and between 75-90 E. These results are consistent, but further refined, than Inmarsat's best guess (as it incorporates much more than two possibilities -- the benefit of using a Monte Carlo model), fitting to the constraints of the pings and their errors. Note that a speed variance would merely move the plane to the red dots, which would be the plane's location at the last ping, or of the plane was super slow, perhaps toward the green dots along the same arc. So they would not change this plot per se; just how you would interpret it. Also note the preferential movement of the plane to the south.¶

An important condition for later, is that the plane needs to progressively get farther away at each ping (which we have been told is the case). This means that the plane cannot continue in a straight line down south, which would rule out autopilot (which on a 777 defaults to a Great Circle path). None of the existing infographics show all of the pings, which is why they haven't noticed this.¶

And note that with the error margins already built in, you have no room on the map to move the pings between the first and the last one around with each other, and have the plane go in a straight line. Note that there are a few possibilities for a Great Circle arc flight path in the 5% error case...but this would mean that ping locations are substantially greater than what has been reported. So this is extremely unlikely.¶

Von Mises Distribution Heading¶

k = 10 (make pointed)¶

In [ ]:

percenterror1,percenterror2 = 0.05, 0.025

k = 10

In [ ]:

plane_hops_5per = []
plane_hops_2per = []

for i in xrange(N):
    plane_hops_5per.append(six_hop_model_von_mises(last_known_heading,pulauperak[1],pulauperak[0],km_hop,k,percenterror1,ping_distances,ping_times))
    plane_hops_2per.append(six_hop_model_von_mises(last_known_heading,pulauperak[1],pulauperak[0],km_hop,k,percenterror2,ping_distances,ping_times))

In [ ]:

first_lat_5per = []
two_lat_5per = []
three_lat_5per = []
four_lat_5per = []
five_lat_5per = []
final_lat_5per = []

first_lon_5per = []
two_lon_5per = []
three_lon_5per = []
four_lon_5per = []
five_lon_5per = []
final_lon_5per = []

route1_lat_5per = []
route2_lat_5per = []
route3_lat_5per = []
route4_lat_5per = []

route1_lon_5per = []
route2_lon_5per = []
route3_lon_5per = []
route4_lon_5per = []

for i in xrange(len(plane_hops_5per)):
    first_lat_5per.append(plane_hops_5per[i][0][0])
    first_lon_5per.append(plane_hops_5per[i][1][0])
    two_lat_5per.append(plane_hops_5per[i][0][1])
    two_lon_5per.append(plane_hops_5per[i][1][1])
    three_lat_5per.append(plane_hops_5per[i][0][2])
    three_lon_5per.append(plane_hops_5per[i][1][2])
    four_lat_5per.append(plane_hops_5per[i][0][3])
    four_lon_5per.append(plane_hops_5per[i][1][3])
    five_lat_5per.append(plane_hops_5per[i][0][4])
    five_lon_5per.append(plane_hops_5per[i][1][4])
    final_lat_5per.append(plane_hops_5per[i][0][5])
    final_lon_5per.append(plane_hops_5per[i][1][5])
    
    route1_lat_5per.append(plane_hops_5per[i][0][6])
    route1_lon_5per.append(plane_hops_5per[i][1][6])
    route2_lat_5per.append(plane_hops_5per[i][0][7])
    route2_lon_5per.append(plane_hops_5per[i][1][7])
    route3_lat_5per.append(plane_hops_5per[i][0][8])
    route3_lon_5per.append(plane_hops_5per[i][1][8])
    route4_lat_5per.append(plane_hops_5per[i][0][9])
    route4_lon_5per.append(plane_hops_5per[i][1][9])

In [ ]:

first_lat_2per = []
two_lat_2per = []
three_lat_2per = []
four_lat_2per = []
five_lat_2per = []
final_lat_2per = []

first_lon_2per = []
two_lon_2per = []
three_lon_2per = []
four_lon_2per = []
five_lon_2per = []
final_lon_2per = []

route1_lat_2per = []
route2_lat_2per = []
route3_lat_2per = []
route4_lat_2per = []

route1_lon_2per = []
route2_lon_2per = []
route3_lon_2per = []
route4_lon_2per = []

for i in xrange(len(plane_hops_2per)):
    first_lat_2per.append(plane_hops_2per[i][0][0])
    first_lon_2per.append(plane_hops_2per[i][1][0])
    two_lat_2per.append(plane_hops_2per[i][0][1])
    two_lon_2per.append(plane_hops_2per[i][1][1])
    three_lat_2per.append(plane_hops_2per[i][0][2])
    three_lon_2per.append(plane_hops_2per[i][1][2])
    four_lat_2per.append(plane_hops_2per[i][0][3])
    four_lon_2per.append(plane_hops_2per[i][1][3])
    five_lat_2per.append(plane_hops_2per[i][0][4])
    five_lon_2per.append(plane_hops_2per[i][1][4])
    final_lat_2per.append(plane_hops_2per[i][0][5])
    final_lon_2per.append(plane_hops_2per[i][1][5])
    
    route1_lat_2per.append(plane_hops_2per[i][0][6])
    route1_lon_2per.append(plane_hops_2per[i][1][6])
    route2_lat_2per.append(plane_hops_2per[i][0][7])
    route2_lon_2per.append(plane_hops_2per[i][1][7])
    route3_lat_2per.append(plane_hops_2per[i][0][8])
    route3_lon_2per.append(plane_hops_2per[i][1][8])
    route4_lat_2per.append(plane_hops_2per[i][0][9])
    route4_lon_2per.append(plane_hops_2per[i][1][9])

5% results:¶

In [32]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Draw great circle to show path autopilot would have taken
fig.drawgreatcircle(pulauperak[1],pulauperak[0],85,-40,linewidth=3,color='black',label='Great Circle Path')

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- these are for each ping. will not plot the 2.5 and 5% error
x5,y5 = fig(circle_lon_211am,circle_lat_211am)
x6,y6 = fig(circle_lon_311am,circle_lat_311am)
x7,y7 = fig(circle_lon_411am,circle_lat_411am)
x8,y8 = fig(circle_lon_511am,circle_lat_511am)
x9,y9 = fig(circle_lon_611am,circle_lat_611am)
x10,y10 = fig(circle_lon_711am,circle_lat_711am)
x11,y11 = fig(circle_lon_811am,circle_lat_811am)

#Add points after each hr
x12,y12 = fig(first_lon_5per,first_lat_5per)
x13,y13 = fig(two_lon_5per,two_lat_5per)
x14,y14 = fig(three_lon_5per,three_lat_5per)
x15,y15 = fig(four_lon_5per,four_lat_5per)
x16,y16 = fig(five_lon_5per,five_lat_5per)
x17,y17 = fig(final_lon_5per,final_lat_5per)

#Add ultimate locations of MH370
x18,y18 = fig(route1_lon_5per,route1_lat_5per)
x19,y19 = fig(route2_lon_5per,route2_lat_5per)
x20,y20 = fig(route3_lon_5per,route3_lat_5per)
x21,y21 = fig(route4_lon_5per,route4_lat_5per)

# plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Draw circle showing extent of Inmarsat sat radar detection for each of the pings
fig.plot(x5,y5,'r--',markersize=5,label='1st Ping Arc')
fig.plot(x6,y6,'r-',markersize=5, label='Ping Arcs After Disappearance')
fig.plot(x7,y7,'r-',markersize=5)
fig.plot(x8,y8,'r-',markersize=5)
fig.plot(x9,y9,'r-',markersize=5)
fig.plot(x10,y10,'r-',markersize=5)
fig.plot(x11,y11,'r-',markersize=5)

#Add monte carlo points
fig.plot(x12,y12,'yo',markersize=5,label='after 1 hrs')
fig.plot(x13,y13,'mo',markersize=5,label= 'after 2 hrs')
fig.plot(x14,y14,'wo',markersize=5,label='after 3 hrs')
fig.plot(x15,y15,'bo',markersize=5,label='after 4 hrs')
fig.plot(x16,y16,'go',markersize=5,label='after 5 hrs')
fig.plot(x17,y17,'ro',markersize=7,label='after 6 hrs')

#Plot ultimate locations of MH370
fig.plot(x18,y18,'bo',markersize=5,label='in final hr')
fig.plot(x19,y19,'bo',markersize=5)
fig.plot(x20,y20,'bo',markersize=5)
fig.plot(x21,y21,'bo',markersize=5)

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='upper right',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('Inmarsat Ping Estimation -- Individual Pings', fontsize=30)

#Show below
plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-32-2f2ff7fd3d11> in <module>()
      1 #Set figure size
----> 2 fig = plt.figure(figsize=[30,20])
      3 
      4 #Setup Basemap
      5 fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

NameError: name 'plt' is not defined

2.5% results:¶

In [ ]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Draw great circle to show path autopilot would have taken
fig.drawgreatcircle(pulauperak[1],pulauperak[0],85,-40,linewidth=3,color='black',label='Great Circle Path')

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- these are for each ping. will not plot the 2.5 and 5% error
x5,y5 = fig(circle_lon_211am,circle_lat_211am)
x6,y6 = fig(circle_lon_311am,circle_lat_311am)
x7,y7 = fig(circle_lon_411am,circle_lat_411am)
x8,y8 = fig(circle_lon_511am,circle_lat_511am)
x9,y9 = fig(circle_lon_611am,circle_lat_611am)
x10,y10 = fig(circle_lon_711am,circle_lat_711am)
x11,y11 = fig(circle_lon_811am,circle_lat_811am)

#Add points after each hr
x12,y12 = fig(first_lon_2per,first_lat_2per)
x13,y13 = fig(two_lon_2per,two_lat_2per)
x14,y14 = fig(three_lon_2per,three_lat_2per)
x15,y15 = fig(four_lon_2per,four_lat_2per)
x16,y16 = fig(five_lon_2per,five_lat_2per)
x17,y17 = fig(final_lon_2per,final_lat_2per)

#Add ultimate locations of MH370
x18,y18 = fig(route1_lon_2per,route1_lat_2per)
x19,y19 = fig(route2_lon_2per,route2_lat_2per)
x20,y20 = fig(route3_lon_2per,route3_lat_2per)
x21,y21 = fig(route4_lon_2per,route4_lat_2per)

# plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Draw circle showing extent of Inmarsat sat radar detection for each of the pings
fig.plot(x5,y5,'r--',markersize=5,label='1st Ping Arc')
fig.plot(x6,y6,'r-',markersize=5, label='Ping Arcs After Disappearance')
fig.plot(x7,y7,'r-',markersize=5)
fig.plot(x8,y8,'r-',markersize=5)
fig.plot(x9,y9,'r-',markersize=5)
fig.plot(x10,y10,'r-',markersize=5)
fig.plot(x11,y11,'r-',markersize=5)

#Add monte carlo points
fig.plot(x12,y12,'yo',markersize=5,label='after 1 hrs')
fig.plot(x13,y13,'mo',markersize=5,label= 'after 2 hrs')
fig.plot(x14,y14,'wo',markersize=5,label='after 3 hrs')
fig.plot(x15,y15,'bo',markersize=5,label='after 4 hrs')
fig.plot(x16,y16,'go',markersize=5,label='after 5 hrs')
fig.plot(x17,y17,'ro',markersize=7,label='after 6 hrs')

#Plot ultimate locations of MH370
fig.plot(x18,y18,'bo',markersize=5,label='in final hr')
fig.plot(x19,y19,'bo',markersize=5)
fig.plot(x20,y20,'bo',markersize=5)
fig.plot(x21,y21,'bo',markersize=5)

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='upper right',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('Inmarsat Ping Estimation -- Individual Pings', fontsize=30)

#Show below
plt.show()

Wrapped Cauchy Distribution Heading¶

k = 0.99 (make pointed if it was going straight ahead)¶

In [ ]:

percenterror1,percenterror2 = 0.05, 0.025

k = 0.99

In [ ]:

plane_hops_5per = []
plane_hops_2per = []

for i in xrange(N):
    plane_hops_5per.append(six_hop_model_wrapped_cauchy(last_known_heading,pulauperak[1],pulauperak[0],km_hop,k,percenterror1,ping_distances,ping_times))
    plane_hops_2per.append(six_hop_model_wrapped_cauchy(last_known_heading,pulauperak[1],pulauperak[0],km_hop,k,percenterror2,ping_distances,ping_times))

In [ ]:

first_lat_5per = []
two_lat_5per = []
three_lat_5per = []
four_lat_5per = []
five_lat_5per = []
final_lat_5per = []

first_lon_5per = []
two_lon_5per = []
three_lon_5per = []
four_lon_5per = []
five_lon_5per = []
final_lon_5per = []

route1_lat_5per = []
route2_lat_5per = []
route3_lat_5per = []
route4_lat_5per = []

route1_lon_5per = []
route2_lon_5per = []
route3_lon_5per = []
route4_lon_5per = []

for i in xrange(len(plane_hops_5per)):
    first_lat_5per.append(plane_hops_5per[i][0][0])
    first_lon_5per.append(plane_hops_5per[i][1][0])
    two_lat_5per.append(plane_hops_5per[i][0][1])
    two_lon_5per.append(plane_hops_5per[i][1][1])
    three_lat_5per.append(plane_hops_5per[i][0][2])
    three_lon_5per.append(plane_hops_5per[i][1][2])
    four_lat_5per.append(plane_hops_5per[i][0][3])
    four_lon_5per.append(plane_hops_5per[i][1][3])
    five_lat_5per.append(plane_hops_5per[i][0][4])
    five_lon_5per.append(plane_hops_5per[i][1][4])
    final_lat_5per.append(plane_hops_5per[i][0][5])
    final_lon_5per.append(plane_hops_5per[i][1][5])
    
    route1_lat_5per.append(plane_hops_5per[i][0][6])
    route1_lon_5per.append(plane_hops_5per[i][1][6])
    route2_lat_5per.append(plane_hops_5per[i][0][7])
    route2_lon_5per.append(plane_hops_5per[i][1][7])
    route3_lat_5per.append(plane_hops_5per[i][0][8])
    route3_lon_5per.append(plane_hops_5per[i][1][8])
    route4_lat_5per.append(plane_hops_5per[i][0][9])
    route4_lon_5per.append(plane_hops_5per[i][1][9])

In [ ]:

first_lat_2per = []
two_lat_2per = []
three_lat_2per = []
four_lat_2per = []
five_lat_2per = []
final_lat_2per = []

first_lon_2per = []
two_lon_2per = []
three_lon_2per = []
four_lon_2per = []
five_lon_2per = []
final_lon_2per = []

route1_lat_2per = []
route2_lat_2per = []
route3_lat_2per = []
route4_lat_2per = []

route1_lon_2per = []
route2_lon_2per = []
route3_lon_2per = []
route4_lon_2per = []

for i in xrange(len(plane_hops_2per)):
    first_lat_2per.append(plane_hops_2per[i][0][0])
    first_lon_2per.append(plane_hops_2per[i][1][0])
    two_lat_2per.append(plane_hops_2per[i][0][1])
    two_lon_2per.append(plane_hops_2per[i][1][1])
    three_lat_2per.append(plane_hops_2per[i][0][2])
    three_lon_2per.append(plane_hops_2per[i][1][2])
    four_lat_2per.append(plane_hops_2per[i][0][3])
    four_lon_2per.append(plane_hops_2per[i][1][3])
    five_lat_2per.append(plane_hops_2per[i][0][4])
    five_lon_2per.append(plane_hops_2per[i][1][4])
    final_lat_2per.append(plane_hops_2per[i][0][5])
    final_lon_2per.append(plane_hops_2per[i][1][5])
    
    route1_lat_2per.append(plane_hops_2per[i][0][6])
    route1_lon_2per.append(plane_hops_2per[i][1][6])
    route2_lat_2per.append(plane_hops_2per[i][0][7])
    route2_lon_2per.append(plane_hops_2per[i][1][7])
    route3_lat_2per.append(plane_hops_2per[i][0][8])
    route3_lon_2per.append(plane_hops_2per[i][1][8])
    route4_lat_2per.append(plane_hops_2per[i][0][9])
    route4_lon_2per.append(plane_hops_2per[i][1][9])

5% results:¶

In [ ]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Draw great circle to show path autopilot would have taken
fig.drawgreatcircle(pulauperak[1],pulauperak[0],85,-40,linewidth=3,color='black',label='Great Circle Path')

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- these are for each ping. will not plot the 2.5 and 5% error
x5,y5 = fig(circle_lon_211am,circle_lat_211am)
x6,y6 = fig(circle_lon_311am,circle_lat_311am)
x7,y7 = fig(circle_lon_411am,circle_lat_411am)
x8,y8 = fig(circle_lon_511am,circle_lat_511am)
x9,y9 = fig(circle_lon_611am,circle_lat_611am)
x10,y10 = fig(circle_lon_711am,circle_lat_711am)
x11,y11 = fig(circle_lon_811am,circle_lat_811am)

#Add points after each hr
x12,y12 = fig(first_lon_5per,first_lat_5per)
x13,y13 = fig(two_lon_5per,two_lat_5per)
x14,y14 = fig(three_lon_5per,three_lat_5per)
x15,y15 = fig(four_lon_5per,four_lat_5per)
x16,y16 = fig(five_lon_5per,five_lat_5per)
x17,y17 = fig(final_lon_5per,final_lat_5per)

#Add ultimate locations of MH370
x18,y18 = fig(route1_lon_5per,route1_lat_5per)
x19,y19 = fig(route2_lon_5per,route2_lat_5per)
x20,y20 = fig(route3_lon_5per,route3_lat_5per)
x21,y21 = fig(route4_lon_5per,route4_lat_5per)

# plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Draw circle showing extent of Inmarsat sat radar detection for each of the pings
fig.plot(x5,y5,'r--',markersize=5,label='1st Ping Arc')
fig.plot(x6,y6,'r-',markersize=5, label='Ping Arcs After Disappearance')
fig.plot(x7,y7,'r-',markersize=5)
fig.plot(x8,y8,'r-',markersize=5)
fig.plot(x9,y9,'r-',markersize=5)
fig.plot(x10,y10,'r-',markersize=5)
fig.plot(x11,y11,'r-',markersize=5)

#Add monte carlo points
fig.plot(x12,y12,'yo',markersize=5,label='after 1 hrs')
fig.plot(x13,y13,'mo',markersize=5,label= 'after 2 hrs')
fig.plot(x14,y14,'wo',markersize=5,label='after 3 hrs')
fig.plot(x15,y15,'bo',markersize=5,label='after 4 hrs')
fig.plot(x16,y16,'go',markersize=5,label='after 5 hrs')
fig.plot(x17,y17,'ro',markersize=7,label='after 6 hrs')

#Plot ultimate locations of MH370
fig.plot(x18,y18,'bo',markersize=5,label='in final hr')
fig.plot(x19,y19,'bo',markersize=5)
fig.plot(x20,y20,'bo',markersize=5)
fig.plot(x21,y21,'bo',markersize=5)

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='upper right',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('Inmarsat Ping Estimation -- Individual Pings', fontsize=30)

#Show below
plt.show()

2.5% results:¶

In [ ]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Draw great circle to show path autopilot would have taken
fig.drawgreatcircle(pulauperak[1],pulauperak[0],85,-40,linewidth=3,color='black',label='Great Circle Path')

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- these are for each ping. will not plot the 2.5 and 5% error
x5,y5 = fig(circle_lon_211am,circle_lat_211am)
x6,y6 = fig(circle_lon_311am,circle_lat_311am)
x7,y7 = fig(circle_lon_411am,circle_lat_411am)
x8,y8 = fig(circle_lon_511am,circle_lat_511am)
x9,y9 = fig(circle_lon_611am,circle_lat_611am)
x10,y10 = fig(circle_lon_711am,circle_lat_711am)
x11,y11 = fig(circle_lon_811am,circle_lat_811am)

#Add points after each hr
x12,y12 = fig(first_lon_2per,first_lat_2per)
x13,y13 = fig(two_lon_2per,two_lat_2per)
x14,y14 = fig(three_lon_2per,three_lat_2per)
x15,y15 = fig(four_lon_2per,four_lat_2per)
x16,y16 = fig(five_lon_2per,five_lat_2per)
x17,y17 = fig(final_lon_2per,final_lat_2per)

#Add ultimate locations of MH370
x18,y18 = fig(route1_lon_2per,route1_lat_2per)
x19,y19 = fig(route2_lon_2per,route2_lat_2per)
x20,y20 = fig(route3_lon_2per,route3_lat_2per)
x21,y21 = fig(route4_lon_2per,route4_lat_2per)

# plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Draw circle showing extent of Inmarsat sat radar detection for each of the pings
fig.plot(x5,y5,'r--',markersize=5,label='1st Ping Arc')
fig.plot(x6,y6,'r-',markersize=5, label='Ping Arcs After Disappearance')
fig.plot(x7,y7,'r-',markersize=5)
fig.plot(x8,y8,'r-',markersize=5)
fig.plot(x9,y9,'r-',markersize=5)
fig.plot(x10,y10,'r-',markersize=5)
fig.plot(x11,y11,'r-',markersize=5)

#Add monte carlo points
fig.plot(x12,y12,'yo',markersize=5,label='after 1 hrs')
fig.plot(x13,y13,'mo',markersize=5,label= 'after 2 hrs')
fig.plot(x14,y14,'wo',markersize=5,label='after 3 hrs')
fig.plot(x15,y15,'bo',markersize=5,label='after 4 hrs')
fig.plot(x16,y16,'go',markersize=5,label='after 5 hrs')
fig.plot(x17,y17,'ro',markersize=7,label='after 6 hrs')

#Plot ultimate locations of MH370
fig.plot(x18,y18,'bo',markersize=5,label='in final hr')
fig.plot(x19,y19,'bo',markersize=5)
fig.plot(x20,y20,'bo',markersize=5)
fig.plot(x21,y21,'bo',markersize=5)

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='upper right',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('Inmarsat Ping Estimation -- Individual Pings', fontsize=30)

#Show below
plt.show()

Okay, so the wrapped cauchy results are worth an addition comment here: with all of the pings we do see some propensity for the MH370 to switch directions midway, and go back to the north when it goes to the south; so the plane does wind up in a variety of places on the last ping arc (what used to be called the 5th ping, and is now the 6th ping.) However, note that in the slight majority of cases, MH370 winds up to the south instead of the north (you can judge this by the blue dots, and the red dots to some extent.) So even under whacky unit circle heading distribution choices the plane still winds up to the south more than the nort.h¶

Overall Conclusions from Monte Carlo Modeling MH370's Last Flight¶

What I've also plotted in these final versions of the model is a Great Circle line (using Basemap's built-in function, with a default of 100 pts created to make the line) of the Palau Perak coordinates -- MH370's final location -- to 40 S, 85 E, which in all of my 6 hop simulations was the arc location where the plane wound up near.¶

Arguably the most stringent scenario for the plane moving west, is the wrapped cauchy at k=0.99 distribution, with a 2.5% standard deviation error for the ping probability distributions. We see that the plane only really starts to move to the west when the last heading dictates up to 1 hr's flight after the last ping; which after the wrapped cauchy condition ends. (Recall the last step takes the last heading and projects it out up to 1 hr, in increments of 1/4, 1/2, 3/4, and 59/60 hrs.)¶

In [ ]:

#Set figure size
fig = plt.figure(figsize=[30,20])

#Setup Basemap
fig = Basemap(width=10000000,height=18000000,projection='lcc',resolution='c',lat_0=10,lon_0=90,suppress_ticks=True)

#Draw coasts
fig.drawcoastlines()

#Draw boundary
fig.drawmapboundary(fill_color='lightblue')

#Fill background
fig.fillcontinents(color='#FFD39B',lake_color='lightblue')

#Draw parallels
parallels = np.arange(lat_min,lat_max,lat_space)
fig.drawparallels(np.arange(lat_min,lat_max,lat_space),labels=[1,1,0,1], fontsize=15)

#Draw meridians
meridians = np.arange(lon_min,lon_max,lon_space)
fig.drawmeridians(np.arange(lon_min,lon_max,lon_space),labels=[1,1,0,1], fontsize=15)

#Draw great circle to show path autopilot would have taken
fig.drawgreatcircle(pulauperak[1],pulauperak[0],85,-40,linewidth=3,color='black',label='Great Circle Path')

#Translate coords into map coord system to plot

#Known 777 Locs
x,y = fig(kualalumpur[1],kualalumpur[0]) #plotted as lon,lat NOT lat,lon -- watch out!!
x2,y2 = fig(igariwaypoint[1],igariwaypoint[0])
x3,y3 = fig(pulauperak[1],pulauperak[0])

#Inmarsat Satellite Loc
x4,y4 = fig(inmarsat[1],inmarsat[0])

#Add circle coords -- these are for each ping. will not plot the 2.5 and 5% error
x5,y5 = fig(circle_lon_211am,circle_lat_211am)
x6,y6 = fig(circle_lon_311am,circle_lat_311am)
x7,y7 = fig(circle_lon_411am,circle_lat_411am)
x8,y8 = fig(circle_lon_511am,circle_lat_511am)
x9,y9 = fig(circle_lon_611am,circle_lat_611am)
x10,y10 = fig(circle_lon_711am,circle_lat_711am)
x11,y11 = fig(circle_lon_811am,circle_lat_811am)

#Add points after each hr
x12,y12 = fig(first_lon_2per,first_lat_2per)
x13,y13 = fig(two_lon_2per,two_lat_2per)
x14,y14 = fig(three_lon_2per,three_lat_2per)
x15,y15 = fig(four_lon_2per,four_lat_2per)
x16,y16 = fig(five_lon_2per,five_lat_2per)
x17,y17 = fig(final_lon_2per,final_lat_2per)

#Add ultimate locations of MH370
x18,y18 = fig(route1_lon_2per,route1_lat_2per)
x19,y19 = fig(route2_lon_2per,route2_lat_2per)
x20,y20 = fig(route3_lon_2per,route3_lat_2per)
x21,y21 = fig(route4_lon_2per,route4_lat_2per)

# plot coords w/ filled circles
fig.plot(x,y,'bo',markersize=10,label='MH370 Flight Path')
fig.plot(x2,y2,'bo',markersize=10)
fig.plot(x3,y3,'go',markersize=10,label='MH370 Last Known Coords')
fig.plot(x4,y4,'ro',markersize=10,label='Inmarsat 3-F1')

#Draw circle showing extent of Inmarsat sat radar detection for each of the pings
fig.plot(x5,y5,'r--',markersize=5,label='1st Ping Arc')
fig.plot(x6,y6,'r-',markersize=5, label='Ping Arcs After Disappearance')
fig.plot(x7,y7,'r-',markersize=5)
fig.plot(x8,y8,'r-',markersize=5)
fig.plot(x9,y9,'r-',markersize=5)
fig.plot(x10,y10,'r-',markersize=5)
fig.plot(x11,y11,'r-',markersize=5)

#Add monte carlo points
fig.plot(x12,y12,'yo',markersize=5,label='after 1 hrs')
fig.plot(x13,y13,'mo',markersize=5,label= 'after 2 hrs')
fig.plot(x14,y14,'wo',markersize=5,label='after 3 hrs')
fig.plot(x15,y15,'bo',markersize=5,label='after 4 hrs')
fig.plot(x16,y16,'go',markersize=5,label='after 5 hrs')
fig.plot(x17,y17,'ro',markersize=7,label='after 6 hrs')

#Plot ultimate locations of MH370
fig.plot(x18,y18,'bo',markersize=5,label='in final hr')
fig.plot(x19,y19,'bo',markersize=5)
fig.plot(x20,y20,'bo',markersize=5)
fig.plot(x21,y21,'bo',markersize=5)

#Draw arrows showing flight path
arrow1 = plt.arrow(x,y,x2-x,y2-y,linewidth=3,color='blue',linestyle='dashed',label='flight path')
arrow2 = plt.arrow(x2,y2,x3-x2,y3-y2,linewidth=3,color='blue',linestyle='dashed',label='flight path')

#Make legend
legend = plt.legend(loc='upper right',fontsize=10,frameon=True,title='Legend',markerscale=1,prop={'size':15})
legend.get_title().set_fontsize('20')

#Add title
plt.title('Great Circle Path on Most Stringent Scenario to Constrain Moving West', fontsize=15)

#Show below
plt.show()

Monte Carlo with Constraints and "Priors" is Better Than Cherry-Picking Flight Paths¶

Here's an important point which hasn't been brought up by the media much -- at the present time of this writing, NTSB has provided 2 possible paths. Here are those same 2, and then 2 by Malayia. Why 2? Surely, putting your eggs in 2 baskets is better than 1, to resurrect this cliche. But what a Monte Carlo model can do is put your eggs in thousands of baskets (I have run 1000 for now, but could reduce the dot size on the plot and increase to 10000 trials, which would take an overnight to run given the number of permutations and the fact that I haven't vectorized my code and the main functions have double 'for' loops which is O(n^2) complexity). Then, from there, you can establish from this data what you think the most likely scenarios are -- either by taking an average if you think error is normally distributed (like Galton's Ox), or by making assumptions ("priors", but from an informal sense as opposed to a strict Bayesian calculation included in the simulation, which I have not had time to implement) which would rule out a good fraction of the plotted locations. So you get a much better idea of where the plane is most likely to be, instead of trying to cherry pick what you think is possible in advance. It's a more scientific approach for how you go about searching for the plane.¶

MH370 Location Is Sensitive to Heading¶

I show this through picking uniform like distributions, and Wrapped Cauchy at k=0.99 which is almost a point; yet we see scattered results. However, the ~255 degree heading is, to the best of my knowledge, the last known heading we have for MH370.¶

MH370 is Likely Off Australia, Still, 5 or 6 Pings Regardless¶

MH370 is likely off the coast of Australia, in the general vicinity of where searches have been taken place over the past week (and starting right around the time when I finished the first version of my model, without that information). I initially thought it was further north -- when I assumed 5 pings (or, equivalently, 6 pings with a 5/6 cruising speed, or ~750 km/hr instead of 905 km/hr cruising speed). Those simulations are all still valid, but would involve a lower speed for the plane than a 777's cruising speed (which is possible, but I haven't seen evidence suggesting this is the case.)¶

The blue dot locations in the scenarios that involve a limited standard devation or high k value is where additional search effort (particularly trolling with ships with a black box ping detector) should be concentrated for the next two weeks. And even after this, scanning the seafloor for MH370. 777s are big, so although the ocean is big they should see it with the proper equipment if they're in the vicinity at the sea surface.¶

Disagree? Choose-Your-Own-Adventure Features Enabled¶

This code is written such that if better speed parameters, or more accurate ping information in released or uncovered, the prediction of where MH370 is most likely to be, can be updated accordingly with minimal effort. I've documented the code and made it available for precisely this purpose. And I will point out bugs and errors I find, etc., whether or not they affect the analysis. As I have done in the past.¶

MH370 Could Not Have Flown a Great Circle Route or via Magnetic Bearing -- Autopilot Ruled Out?¶

Finally, based on the most charitable scenario of the many simulations I ran, with the ping locations inherently constrained by the 1st and last ping, I show almost no possible route to the south that MH370 takes, satisfies a Great Circle requirement for auto pilot (it could appear to, up to 4 hours, which might be accomplished if MH370 is going really slow, but there are 2+ unexplained hours left). If you look closely at the above paths, if I ran 10k or 100k trials there's a chance some might work but this is an exceedingly low probability. So, from a Monte Carlo standpoint, is not a likely conclusion.¶

If the plane continued to fly southwest, and then deviated -- then the Great Circle path seems even more unlikely to reconcile with the pings, with the most likely locations of MH370 according to my model results. Furthermore, autopilot and a "ghost flight" scenario like Helios Airways Flight 522 from Cyprus to Greece in 2005, would be unlikely to vary the speed in between pings if on autopilot. And so far it seems like regardless of whether you move the pings in irregular time intervals, there isn't enough variance from that to account for straight Great Circle path scenarios.¶

MH370's drifting with the pings over time seems to be an inevitable consequence of the geometry of the ping arc constraints. I am not sure how to reconcile this, except that it would mean that one of the points along that route was programmed into the flight management system, or the plane was flown manually. Therefore, I conclude based on my Monte Carlo modeling and analysis that whereever MH370 may be...someone wanted it there.¶

More on Modeling the Last Flight of MH370 with Monte Carlo¶

Conor L. Myhrvold¶

For the first version of the model, with runways, see: IPython Notebook Link¶

For articles about this first work, see Fast Company Labs' : How I Narrowed Down The Location Of Malaysia Air Using "Monte Carlo" Data Models More About Our Methodology: Tracking MH370 With Monte Carlo Data Models ¶

Goal: The overall goal of this project is to use publically available information to pinpoint possible locations of flight MH370, a Malaysian Airlines Boeing 777-200ER which went missing in flight on March 8, 2014.¶

Overview - MH370¶

Getting the Data - Flight MH370 Known Flight Coordinates + Inmarsat Coords¶

Here are some useful constants we need for converting purposes:¶

Now we'll input the coordinates relevant to us:¶

Plotting the Data¶

Before we plot the data, let's formally define our latitude and longitude bounds for our map, up front. That way, we can more readily change them across the document without manually editing one by one. These will double duty to make our probabilities as well:¶

Now we'll make a map with those points on it, using Basemap, an extension of Matplotlib:¶

Satellite Ping Data -- the 5th Ping¶

Now we test the two approximations to better understand how they compare to each other:¶

Incorporating New Info on Ping Distance Error¶

Inmarsat 5% error¶

Inmarsat 2.5% error¶

Now we make a map plot showing these circles. The dashed lines away from the solid red line is the error, 2.5%.¶

Now we have the starting plot that we saw in so many newspaper websites all around the world in March. We have gathered an impressive data set which forms the foundation of the second version of our Monte Carlo model.

Making the Monte Carlo Model¶

Let's keep our Haversine formula from our airplane runways section of the first version of this model; we need this for computing distances below. Apparently, there's a Python module haversine which has this, but let's write this function ourselves:¶

Convert the Inmarsat line location and error bounds into probabilities¶

Determining MH370's Final Location -- (up to) The Last Hour¶

After this, I write a function that will make same-heading vectors from the 5th ping, to whereever the plane may finally be located. Note that the largely-offsetting bugs from the first version have been fixed!¶

Initial Location Uncertainty and Randomness¶

We incorporate the above function our 5 step model -- 1 for each hour, since that's when the pings were received. Afterward, the plane went down sometime in an hr, so I plot where the final contact spot could be at 1/4,1/2,3/4, and 59/60 hrs.¶

Modeling Constants (can vary)¶

Now we test out using the last known heading, a 777's typical cruising speed, and the last known coordinates, below:¶

5% Inmarsat error¶

2.5% Inmarsat error¶

Different Distributions Besides Normal for the Heading -- Von Mises and Wrapped Cauchy¶

First, let's try a Von Mises distribution for our heading. Here's a general function that can do that:¶

Let's test this function to see if it works. We'll test in with a heading of 150 degrees (which means that the curves should be centered at 150, i.e. the mean should be 150):¶

It works! The other distribution we can try is a Wrapped Cauchy distribution for our heading. It's also a popular one to use on a unit circle (so you can say I used 'circular logic' when deciding to pick these two other distributions!)¶

We'll test this function in a similar manner:¶

Finally we update our five hop function with the other distributions. There really isn't any math involved here -- just replacing with the appropriate functions and modifying the names of the inputs to better reflect what's going on in the Wrapped Cauchy and Von Mises distributions for headings:¶

Running Code with New Heading Distributions¶

Von Mises¶

Now we plot. To conserve space I am plotting all of the points at once. I will not make a separate plot just showing up to the 5th ping:¶

5per_k1 plot - Von Mises Heading, k=0, 5% ping error¶

5per_k2 plot - Von Mises Heading, k=1, 5% ping error¶

5per_k3 plot - Von Mises Heading, k=10, 5% ping error¶

2per_k1 plot - Von Mises Heading, k=0, 2.5% ping error¶

2per_k2 plot - Von Mises Heading, k=1, 2.5% ping error¶

2per_k3 plot - Von Mises Heading, k=10, 2.5% ping error¶

Okay, let's take a minute to see what these plots mean!¶

......... ........... ........ ......... ........¶

Finally, in theory 2.5% should have stronger results than 5%, in terms of restricting the search area location at the end...but in practice both are relatively small errors; so we're looking at the same patch of ocean west of Australia.¶

Wrapped Cauchy¶

Similar to the Von Mises runs, I'll use the test k values for Wrapped Cauchy at 2.5 and 5% error for the ping probability distribution.¶

As they would say in Ye Olde English, "the loops runneth over." Once that's done, proceed with the below:¶

Now we plot once again:¶

5per_k1 plot - Wrapped Cauchy Heading, k=0, 5% ping error¶

5per_k2 plot - Wrapped Cauchy Heading, k=0.5, 5% ping error¶

5per_k3 plot - Wrapped Cauchy Heading, k=0.99, 5% ping error¶

2per_k1 plot - Wrapped Cauchy Heading, k=0, 2.5% ping error¶

2per_k2 plot - Wrapped Cauchy Heading, k=0.5, 2.5% ping error¶

2per_k3 plot - Wrapped Cauchy Heading, k=0.99, 2.5% ping error¶

Intriguingly, even with an extremely narrow distribution the Wrapped Cauchy is of little use in determining where plane is. We'll try it again with all of the pings. With this heading distribution we can now say, "it's a wrap!"¶

Conclusions and Comparisons to First Version of the Monte Carlo Model:¶

The most realistic intermediary k value cases for Von Mises show the same results as the normal distribution option -- MH370 ends up off the coast of Australia in most cases. The Wrapped Cauchy, only just using the 5th ping arc, does not make meaningful inferences about where the plane is.¶

Next Step - Adding the Other Pings¶

We've now incorporated the reasonable amount of error Inmarsat would expect from their ping estimate -- just 2.5 or 5% -- but we haven't incorporated some idea of where the other pings have been located. So, let's do that.¶

Here's an article, published right after the first version of my model finished, with a recap of what we know MH370's location over time -- essentially that the first ping was at 2:11 am -- which was before MH370's disappearance from radar -- and that the last was at 8:11 am.¶

Nevertheless, the last ping that was plotted originally was the last ping. It's just that the timing must have been off, somehow, given the current timeline:¶

Replotting the Ping on the Map¶

We'll use the ping distances from the array and Mathematica above, to confirm that they look sensible compared to the infographics out there. We could use a 'for' loop but there are only 6, so I'll do them one-by-one:¶

Running All of the Above with 5 -- (No Wait is it 6?!) -- Individual Satellite Pings¶

Normal Distribution Heading¶

Standard Deviation of 30 degrees (allows for some turning each time)¶

5% results:¶

2.5% results:¶

Von Mises Distribution Heading¶

k = 10 (make pointed)¶

5% results:¶

2.5% results:¶

Wrapped Cauchy Distribution Heading¶

k = 0.99 (make pointed if it was going straight ahead)¶

5% results:¶

2.5% results:¶

For the first version of the model, with runways, see: IPython Notebook Link ¶

**For articles about this first work, see Fast Company Labs' :
How I Narrowed Down The Location Of Malaysia Air Using "Monte Carlo" Data Models

*More About Our Methodology: Tracking MH370 With Monte Carlo Data Models***
¶