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.