*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

In [1]:

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

In [2]:

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

In [3]:

```
# 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¶ms=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]
```

In [4]:

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

In [5]:

```
#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()
```

In [6]:

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

In [7]:

```
"""
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
```

In [8]:

```
#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()
```

In [9]:

```
#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()
```

In [10]:

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

In [11]:

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

In [12]:

```
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])
```

In [13]:

```
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])
```

In [14]:

```
#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()
```