# Nav
# Program to calculate the great-circle distance between any two points on earth
#
# Jos van der Woude
# V 1.0 3 Jan 1997
# V 2.0 1 Aug 1998
# V 2.1 5 Sep 2012

save set 'settings.ini'

# Define starting and ending locations

# Example from Amsterdam ...
set label "ams +" at 4.8, 52.3 right
frlat = 52.303818
frlon =  4.763002

# ... to Sydney
set label "+ syd" at  151.1663208,-33.9360542		# Sydney, Australia
tolat = -33.9360542
tolon = 151.1663208

# Define number of way points
nwp = 100

# Settings 
set parametric
unset key
set trange [0:nwp]
set samples nwp

# **********************************
# * Solves law of cosines for side *
# **********************************
locside(side1,side2,angle) = acos(cos(side1)*cos(side2) + sin(side1)*sin(side2)*cos(angle))


# ***********************************
# * Solves law of cosines for angle *
# ***********************************
locangle(side1, side2, side3) = acos((cos(side1) - cos(side2) * cos(side3)) / (sin(side2) * sin(side3)))



# Define variables (in gnuplot longitudes East are positive) 
rd    = pi/180.0
a     = 0.5 * pi - frlat * rd
beta1 = frlon * rd
c     = 0.5 * pi - tolat * rd
beta2 = tolon * rd

# Use law of cosines to solve for b
b = locside(a,c,(beta1 - beta2))
	
# Use laws of cosines to solve for alpha
alpha = locangle(a,c,b)
	
# Calculate new latitude
newlat(t) = (0.5 * pi - locside(b * t / nwp,c,alpha)) / rd

# Calculate new longitude
newlon(t) = (beta2 + (sin(beta1 - beta2) < 0 ? -1 : 1) * \
	     locangle(b * t / nwp,locside(b * t / nwp,c,alpha),c)) / rd

#
# Plot world map and route
set xrange [-180:180]
set yrange [-90:90]
plot 'world-cil.dat' with lines ls 2, newlon(t),newlat(t) with lines ls 4 

#
# Save numbers to an output file
# You'll have to get rid of the "i" column yourself ...
set table "results.dat"
set format "%.5f"
plot newlon(t),newlat(t)

# Clean up
unset table
load 'settings.ini'