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:

   Using Longitude and Latitude to Determine Distance
   http://mathforum.org/library/drmath/view/51711.html 

   Distance using Latitude and Longitude
   http://mathforum.org/library/drmath/view/54680.html   

But a correspondent recently directed our attention to a Web site,

  The Geographic Information Systems FAQ (See Q5.1)
  http://www.faqs.org/faqs/geography/infosystems-faq/   

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://mathforum.org/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. 

- Doctor Rick, The Math Forum
  http://mathforum.org/dr.math/