#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

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

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

#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

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

#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

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

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

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

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]

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

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

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

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

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

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

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


# 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

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

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

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

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


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

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

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.

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

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

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

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

percenterror = 0.025

std_dev = 30

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

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

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

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