Using Longitude and Latitude to Determine Distance

Date: 17 Apr 1995 17:47:06 -0400
From: Clarence Warren
Subject: (none)

I've been looking for the equation for finding the distance between two
cities, given the latitude and longitude of both cities. For example: what 
is the formula one would use to find the distance between San Francisco 
(N37 37' 00" latitude, W122 22' 00" longitude) and Paris (N48 44' 00" 
latitude, E02 23' 00" longitude.)? 

I'm trying to write a program in Visual Basic that will include this distance
calculation in it.

Any help that members of "The Swat Team" can offer would be greatly 
appreciated.

Date: 17 Apr 1995 21:38:43 -0400
From: Dr. Ken
Subject: Re: distance on a sphere

Hello there!

This problem can be most easily solved by using spherical coordinates on the
earth.  Have you dealt with those before?  Here's the transformation from
spherical coordinates to normal rectangular coordinates, where a=latitude
and b=longitude, and r is the radius of the earth:

x = r Cos[a] Cos[b]
y = r Cos[a] Sin[b]
z = r Sin[a]

Then we'll use the following property of the dot product (notated [p,q]):

[p,q] = Length[p] * Length[q] * Cos[angle between p & q]

Now, any vector that points to a point on the surface of the earth will have
length r.  So the right side we have r^2 * Cos[angle between p & q].  On the
left side, we can pull the r's out of the dot product, and cancel them with
the r's on the right side.  Let t represent the angle between p and q.  Then
if the latitude and longitude of our two cities, p and q, are (a1,b1) and
(a2,b2), we have 

Cos[a1] Cos[b1] Cos[a2] Cos[b2] + Cos[a1] Sin[b1] Cos[a2] Sin[b2] 
    + Sin[a1] Sin[a2]  =  Cos[t]

So you can compute the angle t as a function of a1, b1, a2, b2, which are
the latitudes and longitudes of our cities p and q.  Then visualize what
you've got:  draw a great circle through the points p and q.  This is just a
plain old Joe-Schmoe circle of radius r, and the angle t is the angle of the
arc that subtends p and q.  The problem from here on out is just figuring
out what the arc length between p and q is.  The relevant formula is 
Arc length = t/360 * 2Pi* r.  So that's your formula.  By substitution, we
have 

Arccos[Cos[a1] Cos[b1] Cos[a2] Cos[b2] + Cos[a1] Sin[b1] Cos[a2] Sin[b2]
    + Sin[a1] Sin[a2]]/360 * 2Pi * r

Oh, by the way, West longitude means negative values of b, and South
latitude means negative values of a.

Enjoy, and let us know if I've made mistakes, or something's not clear!

-Ken "Dr." Math

From: Dr. Ken
Date: October 30, 1997

I thought I'd add a couple of remarks that may help some people use the above
formula.  First, of all, you need to make sure your calculator or computer
is using degrees, not radians, to figure out the Sine, Cosine, and ArcCosine
functions.  If you're using a calculator or computer that uses radians, then
use a different version of the formula:

Arccos[Cos[a1] Cos[b1] Cos[a2] Cos[b2] + Cos[a1] Sin[b1] Cos[a2] Sin[b2]
    + Sin[a1] Sin[a2]] * r

Also, keep in mind that these formulae don't take into account the squashed
nature of the earth.  As you may know, the earth is kind of fat around the
equator, as a result of the centrifugal force it gets from spinning on its
axis.  So that will throw a little error into these calculations.  I've
never tried to come up with a formula that takes the squashing into account,
and I suspect it might be hard.

-Dr. Ken
 The Math Forum

Date: 11/21/2001 at 09:28:12
From: Nat Keller
Subject: Angular earth distance

I saw Chris Michels' Latitude/Longitude Distance Calculation page at

   http://jan.ucc.nau.edu/~cvm/latlongdist.html   

which he says is based on the above answer. I recently needed to program a 
latitude, longitude, distance calculation and found a formula in the SAS tech 
page (as well as a few other places).  I think that your formula can be 
reduced to the following (and they give the same results):

   A = LAT1, B = LONG1
   C = LAT2, D = LONG2 (all converted to radians: degree/57.29577951)

   IF A = C AND B = D THEN DISTANCE = 0; 
   ELSE

   IF [SIN(A)SIN(C)+COS(A)COS(C)COS(B-D)] > 1 THEN DISTANCE = 
   3963.1*ARCOS[1]; /* solved a prob I ran into.  I haven't fully 
   analyzed it yet */                   
   ELSE

   DISTANCE=3963.1*ARCOS[SIN(A)SIN(C)+COS(A)COS(C)COS(B-D)];

With more time I would solidify and show that it equals your formula.

NK

Date: 11/21/2001 at 10:38:54
From: Doctor Rick
Subject: Re: Angular earth distance

Hi, Nat.

If you go to our Dr. Math Search page at

  http://mathforum.org/mathgrepform.html   

and look for items with the words latitude, longitude, and distance, you 
will find several formulas, including yours. Chris Michels must have found 
the formula this way himself; it was not written in response to him. He did 
not find the simplest formula in our archives. He did not find the most 
reliable formula, either; another, the Haversine formula, is supposed to be 
less prone to rounding errors under certain conditions.

You're correct that the formula he found can be simplified to yours. It's 
easy to start with yours and obtain his, by applying the angle-difference 
identity

  cos(b-d) = cos(b)*cos(d) + sin(b)*sin(d)

You can also follow through the derivation of the more complicated formula 
and modify it by rotating the cartesian coordinate system so that point 1 
is in the x-z plane. Then the y coordinate of point 1 is 0, which causes 
terms to vanish. The x coordinate of point 2 is r*cos(a1)*cos(b2-b1), and 
the y coordinate of point 2 is r*cos(a1)*sin(b2-b1).

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