cos(D) = sin(latB)·sin(latA) + cos(latB)·cos(latA)·cos(lonB - lonA)
sin(latB) - sin(latA)·cos(D)
cos(C) = ---------------------------
cos(latA)·sin(D)
A - origin
B - destination
D = angular distance along path
C = true bearing from the origin to the destination measured from north.
If the value for sin(lonb - lona) is positive;
otherwise, the true bearing is 360° - C.
In applying the above formula, south latitudes and west longitudes are treated as negative angles. Note: When the origin is exactly at a pole, the bearing or course to the destination cannot be determined, why? The calculated bearing angle or course is not a constant. As one moves along the great circle path, the bearing angle to the destination will change continuously. There are two exceptions to this, what are they?
The formulas above assume that the Earth is spherical. The Earth's shape can more accurately be described as an oblate spheroid, a flattened sphere. If a spherical shape is assumed, differences of several tens of kilometers in the calculated versus the true distance result. To obtain distance calculations that are more precise, the Zeno program implements the method given in chapter 10 of reference 1. Additionally, the program uses WGS 1984 parameters for the Earth's equatorial radius and ellipsoid flattening given in reference 2.
Zeno 1.2 is an interpreter for the Zeno programming language. It is an easy to learn and is suitable for educational purposes.
const PI : real := 2*arcsin(1)
program
var sgn, deg, min, sec : int
var lat1, lat2, lat3, lat4 : real
var lon1, lon2, lon3, lon4 : real
var dist, crse : real
% Sydney, AS 33°53' S 151°13' E
lat1 := latitude( -1, 33, 53, 0 )
lon1 := longitude( +1, 151, 13, 0 )
% London, UK, 51°30' N 0°07' W
lat2 := latitude( +1, 51, 30, 0 )
lon2 := longitude( -1, 0, 7, 0 )
% Boston, MA US 42°21' N 71°04' W
lat3 := latitude( +1, 42, 21, 0 )
lon3 := longitude( -1, 71, 4, 0 )
% Tokyo, JA, 35°41' N 139°45' E
lat4 := latitude( +1, 35, 41, 0 )
lon4 := longitude( +1, 139, 45, 0 )
put "Sydney, AS to London, UK"
dist := distance(lat1, lon1, lat2, lon2)
crse := course(lat1, lon1, lat2, lon2)
put "distance = ", dist:8:2, " kilometers"
put "heading = ", crse:8:2, "° true (initially)\n"
put "London, UK to Boston, MA US"
dist := distance(lat2, lon2, lat3, lon3)
crse := course(lat2, lon2, lat3, lon3)
put "distance = ", dist:8:2, " kilometers"
put "heading = ", crse:8:2, "° true (initially)\n"
put "Boston, MA US to Tokyo, JA"
dist := distance(lat3, lon3, lat4, lon4)
crse := course(lat3, lon3, lat4, lon4)
put "distance = ", dist:8:2, " kilometers"
put "heading = ", crse:8:2, "° true (initially)\n"
put "Tokyo, JA to Sydney, AS"
dist := distance(lat4, lon4, lat1, lon1)
crse := course(lat4, lon4, lat1, lon1)
put "distance = ", dist:8:2, " kilometers"
put "heading = ", crse:8:2, "° true (initially)"
end program
% distance using Meeus approximation
function distance( lat1, lon1, lat2, lon2 : real ) : real
var F : real := (lat1 + lat2)/2.0
var G : real := (lat1 - lat2)/2.0
var L : real := (lon1 - lon2)/2.0
var sinG2 : real := sin(G)^2
var cosG2 : real := cos(G)^2
var sinF2 : real := sin(F)^2
var cosF2 : real := cos(F)^2
var sinL2 : real := sin(L)^2
var cosL2 : real := cos(L)^2
var S : real := sinG2*cosL2 + cosF2*sinL2
var C : real := cosG2*cosL2 + sinF2*sinL2
var w : real := arctan(sqrt(S/C))
var R : real := sqrt(S*C)/w
var a : real := 6378.137 % WGS-84 equatorial radius
var f : real := 1.0/298.257223563 % WGS-84 ellipsoid flattening factor
var D : real := 2*w*a
var H1 : real := (3*R - 1)/(2*C)
var H2 : real := (3*R + 1)/(2*S)
var dist : real := D*(1 + f*H1*sinF2*cosG2 - f*H2*cosF2*sinG2)
return round(100*dist)/100.0
end function
% course using spherical trigonometry
% does not work if one latitude is polar!!!
function course( lat1, lon1, lat2, lon2 : real ) : real
var L, D, cosD, C, cosC : real
L := lon2 - lon1
cosD := sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(L)
D := arccos(cosD)
cosC := (sin(lat2) - cosD*sin(lat1))/(sin(D)*cos(lat1))
% numerical error can result in |cosC| slightly > 1.0
if cosC > 1.0 then
cosC := 1.0
end if
if cosC < -1.0 then
cosC := -1.0
end if
C := 180.0*arccos(cosC)/PI
if sin(L) < 0.0 then
C := 360.0 - C
end if
return round(100*C)/100.0
end function
% convert to radians
function latitude( sgn, deg, min, sec : int ) : real
var lat : real := 0
if (0.0 <= sec) and (sec < 60.0) then
lat := lat + sec/3600.0
else
put "check latitude seconds!"
end if
if (0.0 <= min) and (min < 60.0) then
lat := lat + min/60.0
else
put "check latitude minutes!"
end if
if (0.0 <= deg) and (deg <= 90.0) then
lat := lat + deg
else
put "check latitude degrees!"
end if
if (0.0 <= lat) and (lat <= 90.0) then
lat := PI*lat/180.0
else
put "latitude range error!"
end if
assert (sgn = +1) or (sgn = -1)
return sgn*lat
end function
% convert to radians
function longitude( sgn, deg, min, sec : int ) : real
var lon : real := 0
if (0.0 <= sec) and (sec < 60.0) then
lon := lon + sec / 3600.0
else
put "check longitude seconds!"
end if
if (0.0 <= min) and (min < 60.0) then
lon := lon + min / 60.0
else
put "check longitude minutes!"
end if
if (0.0 <= deg) and (deg <= 180.0) then
lon := lon + deg
else
put "check longitude degrees!"
end if
if (0.0 <= lon) and (lon <= 180.0) then
lon := PI * lon / 180.0
else
put "longitude range error!"
end if
assert (sgn = +1) or (sgn = -1)
return sgn*lon
end function
Sydney, AS to London, UK distance = 16977.12 kilometers heading = 319.14° true (initially) London, UK to Boston, MA US distance = 5280.81 kilometers heading = 288.27° true (initially) Boston, MA US to Tokyo, JA distance = 10814.72 kilometers heading = 335.21° true (initially) Tokyo, JA to Sydney, AS distance = 7780.17 kilometers heading = 169.91° true (initially)
AbCd Classics - free on-line books