Here is an illustration of input that is quite useless to the group
but a subject that all naval officers are required to know as part of
the qualifying syllabus.
If you are interested in the subject matter please take it to:
sci.navigation.romp
Haversine Formula and Half Log Haversine Formula:
To calculate distance and bearing between two Latitude/Longitude
points without the use of tables.
This script calculates great-circle distances between the two points –
that is, the shortest distance over the earth’s surface – using the
‘Haversine’ formula.
It assumes a spherical earth, ignoring ellipsoidal effects – which is
accurate enough* for most purposes… – giving an ‘as-the-crow-flies’
distance between the two points (ignoring any hills!).
Enter the co-ordinates into the text boxes to try it out. It accepts a
variety of formats:
deg-min-sec suffixed with N/S/E/W (e.g. 40°44'55?N, 73 59 11W), or
signed decimal degrees without compass direction, where negative
indicates west/south (e.g. 40.7486, -73.9864):
Lat 1: Long 1:
Lat 2: Long 2:
And you can see it on a map (thanks to the nice guys at Google Maps)
Haversine formula:
R = earth’s radius (mean radius = 6,371km)
?lat = lat2- lat1
?long = long2- long1
a = sin²(?lat/2) + cos(lat1).cos(lat2).sin²(?long/2)
c = 2.atan2(va, v(1-a))
d = R.c
(Note that angles need to be in radians to pass to trig functions).
The Haversine formula remains particularly well-conditioned for
numerical computation even at small distances – unlike calculations
based on the spherical law of cosines. (It was published by R W
Sinnott in Sky and Telescope, 1984; the ‘half-versed-sine’ is
(1-cos?)/2, or sin²(?/2) – don’t ask, I’m not a mathematician).
The JavaScript implementation is shown below.
Note (July 2005): on further investigation, I’ve found that the
numeric precision of JavaScript is so good (15 significant figures
using IEEE 754 floating-point numbers) that the simple spherical law
of cosines formula gives well-conditioned results down to distances as
small as around 1 metre. In view of this it is probably worth, in most
situations, using either the simpler law of cosines or the more
accurate ellipsoidal Vincenty formula in preference to Haversine! (See
notes below on the limitations in accuracy of the spherical model).
Spherical law of cosines: d =
acos(sin(lat1).sin(lat2)+cos(lat1).cos(lat2).cos(l ong2-long1)).R
Formula: ? = atan2( sin(?long).cos(lat2),
cos(lat1).sin(lat2) - sin(lat1).cos(lat2).cos(?long) )
Since atan2 returns values in the range -p ... +p, to normalise the
result to a compass bearing, multiply ? by 180/p then use (?+360) %
360, where % is modulo.
This is the initial bearing which if followed in a straight line along
a great-circle arc (orthodrome) will take you from the start point to
the end point; in general, the bearing you are following will have
varied by the time you get to the end point (if you were to go from
say 35°N,45°E (Baghdad) to 35°N,135°E (Osaka), you would start on a
bearing of 60° and end up on a bearing of 120°!).
For final bearing, take the initial bearing from the end point to the
start point and reverse it (using ? = (?+180) % 360).
Formula: Bx = cos(lat2).cos(?long)
By = cos(lat2).sin(?long)
latm = atan2(sin(lat1) + sin(lat2), v((cos(lat1)+Bx)² + By²))
lonm = lon1 + atan2(By, cos(lat1)+Bx)
Just as the initial bearing may vary from the final bearing, the
midpoint may not be located half-way between latitudes/longitudes; the
midpoint between 35°N,45°E and 35°N,135°E is around 45°N,90°E.
Destination point given distance and bearing from start point
This page is steadily growing! Given a start point, initial bearing,
and distance, this will calculate the destination point and final
bearing travelling along a (shortest distance) great circle arc.
Start Lat: Start Long:
Bearing (deg): Distance (km):
Formula: lat2 = asin(sin(lat1)*cos(d/R) +
cos(lat1)*sin(d/R)*cos(brng))
lon2 = lon1 + atan2(sin(brng)*sin(d/R)*cos(lat1),
cos(d/R)-sin(lat1)*sin(lat2))
d/R is the angular distance (in radians), where d is the distance
travelled and R is the earth’s radius
For final bearing, take the initial bearing from the end point to the
start point and reverse it (using ? = (?+p) % 2p).
Rhumb lines
A ‘rhumb line’ (or loxodrome) is a path of constant bearing, which
crosses all meridians at the same angle.
Sailors used to navigate along rhumb lines since it is easier to
follow a constant compass bearing than to constantly adjust the
bearing as is needed to follow a great circle. Rhumb lines are
straight lines on a Mercator Projection map.
Lat 1: Long 1:
Lat 2: Long 2:
Formula: ?f = ln(tan(lat2/2+p/4)/tan(lat1/2+p/4))
if E:W line q = cos(lat1) (length of a circle of latitude)
otherwise q = ?lat/?f
d = v[?lat² + q².?lon²].R
? = atan2(?lon, ?f)
where ln is natural log, ?lon is taking shortest route (<180º), and
R is the earth’s radius
Given a start point and a distance d along constant bearing ?, this
will calculate the destination point. If you maintain a constant
bearing along a rhumb line, you will gradually spiral in towards one
of the poles.
Start Lat: Start Long:
Bearing (deg): Distance (km):
Formula: a = d/R (angular distance)
lat2 = lat1 + a.cos(?)
?f = ln(tan(lat2/2+p/4)/tan(lat1/2+p/4))
if E:W line q = cos(lat1) (length of a circle of latitude)
otherwise q = ?lat/?f
?lon = a.sin(?)/q
lon2 = (lon1+?lon+p) % 2.p - p
where ln is natural log and % is modulo, ?lon is taking shortest
route (<180º), and R is the earth’s radius
If you use Ordnance Survey Grid References, I have implemented a
script for converting between Lat/Long & OS Grid References.
--------------------------------------------------------------------------------
Convert between degrees-minutes-seconds & decimal degrees
Latitude Longitude 1° ˜ 111 km (110.57 eq’l — 111.70 polar)
1' ˜ 1.85 km (= 1 nm) 0.001° ˜ 111 m
1? ˜ 30.9 m 0.00001° ˜ 1 m
--------------------------------------------------------------------------------
*Since the earth is not quite a sphere, there are small errors in
using spherical geometry; the earth is actually roughly ellipsoidal
(or more precisely, oblate spheroidal) with a radius varying between
about 6,378km (equatorial) and 6,357km (polar), and local radius of
curvature varying from 6,336km (equatorial meridian) to 6,399km
(polar). This means that errors from assuming spherical geometry might
be up to 0.55% crossing the equator, though generally below 0.3%,
depending on latitude and direction of travel. An accuracy of better
than 3m in 1km is good enough for me, but if you want greater
accuracy, you could refine the result by using the local radius of
curvature, as explained in the US Census Bureau GIS FAQ.
Take equatorial radius a = 6378 km
polar radius b = 6357 km
eccentricity e = v(1-b²/a²)
if midpoint latitude f = as given above
local meridional (N/S) radius of curvature is ? = a.(1 - e²) / (1 -
e².sin²f)3/2
local transverse radius of curvature is ? = a / v(1 - e².sin²f)
if bearing a = as given above, taken from midpoint to endpoint
then local radius of curvature R' = ?.? / (?.cos²a+?.sin²a)
Errors would then be below 0.01% over smaller distances, up to about
0.1% for trans-continental distances. Thanks to Jason Eisner at Johns
Hopkins University for helping with this.
Alternatively, you could use the Vincenty formula for calculating
geodesic distances on ellipsoids, which gives results accurate to
within 1mm. Out of sheer perversity (I’ve never needed such accuracy),
I looked up this formula and discovered the JavaScript implementation
was simpler than I expected.
--------------------------------------------------------------------------------
Notes: trig functions take arguments in radians, so latitude,
longitude, and bearings in degrees (either decimal or
degrees/minutes/seconds) need to be converted to radians, rad =
p.deg/180. When converting radians back to degrees (deg = 180.rad/p),
West is negative if using signed decimal degrees. For bearings, values
in the range -p to +p (-180° to +180°) need to be converted to 0 to
+2p (0°–360°); this can be done by (brng+2.p)%2.p where % is the
modulo operator. View page source to see JavaScript functions to
handle these conversions.
The atan2() function widely used here takes two arguments, atan2(y,
x), and computes the arc tangent of the ratio y/x. It is more flexible
than atan(y/x), since it handles x=0, and it also returns values in
all 4 quadrants -p to +p (the atan function returns values in the
range -p/2 to +p/2).
If you implement any formula involving atan2 in Microsoft Excel, you
will need to reverse the arguments, as Excel has them the opposite way
around from JavaScript – conventional order is atan2(y, x), but Excel
uses atan2(x, y)
For miles, divide km by 1.609344
For nautical miles, divide km by 1.852
Thanks to Ed Williams’ Aviation Formulary for many of the formulae
Principal JavaScript functions for distance calculation are shown
below; use ‘View Source’ to see the full JavaScript implementation
including all other formulae shown above. These should be simple to
translate into other languages if required. If you can’t or don’t want
to use JavaScript constructors, functions could accept separate values
for the latitude and longitude of a point as arguments, and return a
point as an array (e.g. return [lat, lon]
. You are welcome to re-use
these scripts [without any warranty express or implied] provided you
retain my copyright notice and when possible a link to my website. If
you have any queries or find any problems, please contact me.
© 2002-2006 Chris Veness
Posting in order to be considered academically brilliant! 
Linear Mode
