# More on Modeling the Last Flight of MH370 with Monte Carlo - Part 3/3¶

## 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 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

## Next Step - Adding the Other Pings¶

#### 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.¶

In [25]:
#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])


#### 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]

## 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 [26]:
#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 [68]:
#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')

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

#Show below
plt.show()


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

#### 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 [69]:
"""
a function which given a list of discrete probabilities for each destination point,
it will choose one of those points.

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%.

replace "dist_from_sat" with "ping_distance" since that's changing. run 6 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

for i in xrange(len(plane_lat)):
#new_circle gives up possible coords for diff headings

raw_weights = np.zeros(len(new_circle))
for j in xrange(len(new_circle)):

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

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

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 [70]:
"""
a function which given a list of discrete probabilities for each destination point,
it will choose one of those points.

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.

"""

#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

for i in xrange(len(plane_lat)):
#new_circle gives up possible coords for diff headings

raw_weights = np.zeros(len(new_circle))
for j in xrange(len(new_circle)):

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

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

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 [71]:
"""
a function which given a list of discrete probabilities for each destination point,
it will choose one of those points.

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.

"""

#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

for i in xrange(len(plane_lat)):
#new_circle gives up possible coords for diff headings

raw_weights = np.zeros(len(new_circle))
for j in xrange(len(new_circle)):

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

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

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 [72]:
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


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

In [73]:
percenterror1,percenterror2 = 0.05, 0.025

std_dev = 30

In [74]:
plane_hops_5per = []
plane_hops_2per = []

for i in xrange(N):

In [75]:
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 [76]:
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 [77]:
#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)

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)

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)

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')