Deriving the Haversine Formula
Date: 04/20/99 at 14:27:24
From: Dena Neff
Subject: Longitude and Latitude Calculations
I'm an SAS programmer at The Coca-Cola Company. It's been a long, long
time since I have had to use trigonometry. I need to write a program
module to calculate distances given longitude and latitude data. I am
trying to find an object within a mile's radius of its location. Would
you have the equation for this scenario?
Thank you. I hope I am not too old for a response. :)
Date: 04/20/99 at 17:16:15
From: Doctor Rick
Subject: Re: Longitude and Latitude Calculations
Hello, Dena, welcome to Ask Dr. Math. We're allowed to answer "old"
people now and then, if the question is interesting enough! :-)
You can find several expositions of the "Law of Cosines" for spherical
trigonometry in our Dr. Math Archives. For example:
http://forum.swarthmore.edu/dr.math/problems/longandlat.html
http://forum.swarthmore.edu/dr.math/problems/reed12.31.97.html
But a correspondent recently directed our attention to a Web site,
The Geographic Information Systems FAQ
http://www.census.gov/cgi-bin/geo/gisfaq?Q5.1
This site points out that the cosine formula does not give accurate
results for small distances, and presents a better formula, the
Haversine Formula. The form of the Haversine Formula that uses the
atan2() function is probably best for programming, assuming your
programming language has an atan2() function. I will copy some
excerpts here for your convenience, but I recommend looking at the
site for the other information you will find there.
Note that in the description that follows, atan2(y, x) = arctan(y/x).
However, there does seem to be some confusion about the order of the
arguments. C/C++ puts the arguments in the order that I assume, as
indicated by the Unix man page as well as the Visual C++ Programmer's
Guide. On the other hand, Excel describes the function (not too
clearly) as "ATAN2(x_num, y_num) returns the arctangent of the
specified x- and y-coordinates, in radians." I tested it, and indeed it
takes the denominator x first.
=============================
Presuming a spherical Earth with radius R (see below), and that the
locations of the two points in spherical coordinates (longitude and
latitude) are lon1,lat1 and lon2,lat2, then the Haversine Formula
(from R. W. Sinnott, "Virtues of the Haversine," Sky and Telescope,
vol. 68, no. 2, 1984, p. 159):
dlon = lon2 - lon1
dlat = lat2 - lat1
a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2
c = 2 * atan2(sqrt(a), sqrt(1-a))
d = R * c
will give mathematically and computationally exact results. The
intermediate result c is the great circle distance in radians. The
great circle distance d will be in the same units as R.
Most computers require the arguments of trigonometric functions to be
expressed in radians. To convert lon1,lat1 and lon2,lat2 from degrees,
minutes, and seconds to radians, first convert them to decimal
degrees. To convert decimal degrees to radians, multiply the number of
degrees by pi/180 = 0.017453293 radians/degree.
Inverse trigonometric functions return results expressed in radians.
To express c in decimal degrees, multiply the number of radians by
180/pi = 57.295780 degrees/radian. (But be sure to multiply the number
of RADIANS by R to get d.)
The historical definition of a "nautical mile" is "one minute of arc
of a great circle of the earth." Since the earth is not a perfect
sphere, that definition is ambiguous. However, the internationally
accepted (SI) value for the length of a nautical mile is (exactly, by
definition) 1.852 km or exactly 1.852/1.609344 international miles
(that is, approximately 1.15078 miles - either "international" or
"U.S. statute"). Thus, the implied "official" circumference is 360
degrees times 60 minutes/degree times 1.852 km/minute = 40003.2 km.
The implied radius is the circumference divided by 2 pi:
R = 6367 km = 3956 mi
=============================
I'll just add that the intermediate result "a" is the square of half
of the straight-line distance (chord length) between the two points.
a = (AB/2)^2
= (sin^2(dlat/2)) + cos(lat1) * cos(lat2) * sin^2(dlon/2)
If you are interested in a derivation of the formula, I could provide
it, but the figure would be hard to draw with "character graphics."
- Doctors Rick and Peterson, The Math Forum
http://forum.swarthmore.edu/dr.math/
Date: 04/21/99 at 09:14:55
From: Dena Neff
Subject: Longitude and Latitude Calculations
Doctor Rick,
Thank you very much for your quick response. I'm looking forward to
visiting the Web sites you provided to review the expositions
of the "Law of Cosines" and the Haversine Formula.
I think I will always have a nature of curiosity about me. Therefore,
I would LOVE to see a derivation of the formula.
Dena
Date: 04/21/99 at 22:42:45
From: Doctor Rick
Subject: Re: Longitude and Latitude Calculations
Hello again, Dena.
You called my bluff! Good - now I have a reason to write up my
derivation of the haversine formula. I had derived the cosine formula
for great-circle distance, but not the haversine formula. I was
reluctant to send the formula to you without having checked it out, so
I derived it before I answered you. It turned out to be relatively
straightforward.
First, a word about the meaning of "haversine." This is the name of
one of several lesser-known trigonometric functions that were once
familiar to navigators and the like. The VERSINE (or versed sine) of
angle A is 1-cos(A). The HAVERSINE is half the versine, or
(1-cos(A))/2. With a little math, you can show that
hav(A) = (1-cos(A))/2 = sin^2(A/2)
where ^2 means squared. You see the form on the right twice in the
haversine formula.
In proving the formula, we will work with a sphere of radius 1. The
distance we get should be multiplied by R, the radius of the earth.
Let's start with a formula for the length of a chord subtending an
angle theta in the unit circle (circle of radius 1).
O
/|\
/ | \
1 / | \ 1
/ | \
/ | \
/_____|_____\
A C B
O is the center of the circle, and A and B are on the circle, so AB is
the chord. If angle AOB = theta, then angle AOC = theta/2 where OC is
perpendicular to the chord AB. C bisects AB, since triangle AOB is
isosceles. Segment AC = sin(theta/2), so the length of chord AB is
2*sin(theta/2).
Now for the 3-dimensional figure. Here is a sphere of radius 1, with a
wedge removed to show the center, labeled O. The points we are
interested in are A (lat1,lon1) and B (lat2,lon2). The lines of
longitude lon1 and lon2 are shown as curved single lines, meeting at
the north pole, N; the equator is shown as a double line (equal
signs).
***********
**** N ****
**** / | \ ****
*** / | \ ***
** / | \ **
** / | \ **
* / | \ *
* / | \ *
* A----------------------D *
* /|\ | \ *
* | | \ | | *
* | | \ | | *
* | | \ | | *
* | | \ | | *
* C----------\-----------------B *
* | | - _\ | | *
* | | _O_ | *
* = | | _ _ | = *
* === | G_ _ | === *
* ====E______________________________F ===== *
* ========= ========= *
* ============== *
* *
* *
* *
** **
** **
*** ***
**** ****
**** ****
*************
A and B are opposite vertices of an isosceles trapezoid, ACBD, with
additional vertices C (lat2,lon1) and D (lat1,lon2). I haven't tried
to draw the straight lines (chords) connecting A to C and B to D, but
I show the chords connecting A to D and B to C. I won't bother to
prove that the trapezoid ACBD is planar.
If you join A to O and C to O, the angle AOC is the difference between
lat1 and lat2, or dlat. The length of chord AC is therefore
2*sin(dlat/2), using the chord-length formula derived above. Chord BD
is the same length.
Points E and F are the points where the latitude lines lat1 and lat2
meet the equator. The length of EF is 2*sin(dlon/2) by the chord-
length formula.
Points A and D are on a circle of constant latitude lat1. The radius
of this circle is cos(lat1). You can see this by dropping a
perpendicular from A to OE, meeting OE at G. The angle EOA is lat1,
so OG = cos(lat1). But this is equal to the radius of the latitude
circle at A.
Therefore the length of chord AD is 2*sin(dlon/2)*cos(lat1).
Similarly, the length of chord CB is 2*sin(dlon/2)*cos(lat2).
Now let's go back to 2 dimensions, and find the length of a diagonal
of an isosceles trapezoid.
A_________________D
/| * \
/ | * \
/ | * \
/ | * \
/____|_________________*__\
C H B
The length of CH, where AH is perpendicular to CB, is (CB-AD)/2. By
the Pythagorean theorem,
AH^2 = AC^2 - CH^2
= AC^2 - (CB-AD)^2/4
The length of HB is (CB+AD)/2. Using Pythagoras again, we find that
AB^2 = AH^2 + HB^2
= AC^2 - (CB-AD)^2/4 + (CB+AD)^2/4
= AC^2 + CB*AD
We're ready to plug in the lengths of chords AC, AD, and CB from our
work with the sphere:
AB^2 = 4*(sin^2(dlat/2) + cos(lat1)*cos(lat2)*sin^2(dlon/2))
The intermediate result a is the square of half the chord length AB:
a = (AB/2)^2
= (sin^2(dlat/2)) + cos(lat1) * cos(lat2) * sin^2(dlon/2)
The final step is to find the central angle AOB that corresponds to
this chord length. The arctan formula can be found from this figure:
O
/|\
/ | \
/ | \
/ | \
1 / | \ 1
/ | \
/ | \
/ | \
/ | \
/ | \
/ sqrt(1-a)| \
/ | \
/____________|____________\
A sqrt(a) C B
If AC = sqrt(a), then Pythagoras says that
OC = sqrt(OA^2 - AC^2)
= sqrt(1-a)
Therefore tan(<AOC) = AC/OC = sqrt(a)/sqrt(1-a), or
c = 2 * arctan(sqrt(a)/sqrt(1-a))
where c is the angle AOB.
|