COMPUTING DISTANCES
From: Alex Gleyzer
To: G22_3033_010_fa01
Subject: [G22_3033_010_fa01] another distance formula
I was getting some weird results from the formula given in the
assignment, so I found another one at this location:
http://www.census.gov/cgi-bin/geo/gisfaq?Q5.1
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 * arcsin(min(1,sqrt(a)))
d = R * c
where R is the radius of the Earth.
It appears to generate more believable results, at least for my data.
Also my guess is that latitude/longitude info in Tiger is stored in
decimal degrees. To get radians I multiply them by PI/180.
Alex
Following is a Mirror of
http://www.census.gov/geo/www/faq-index.html
Q5.1: What is the best way to calculate the distance between 2 points?
From: Bob Chamberlain <rgc@jpl.nasa.gov>
PRE-GEOMETRY DISCUSSION OF CALCULATING DISTANCES
ON THE EARTH'S SURFACE
(The professional-level discussion follows this brief summary.)
The Earth is round, but big, so we can consider it flat for short
distances. But, even though the circumference of the Earth is about
25,000 miles (40,000 kilometers), flat-Earth formulas for calculating
the distance between two points start showing noticeable errors when
the distance is more than about 12 miles (20 kilometers). Of course,
how much error is "noticeable" depends on how you are going to use the
result.
If lots of distances are going to be computed, the amount of
calculation has to be considered. I studied this issue because I was
working on a computer program to support U.S. Army training exercises
- our program typically calculated hundreds of thousands of distances
every hour.
Cartesian coordinates express distances in two different directions,
such as north-south for one direction and east-west for the other.
The straight line distance between two points can then be thought of
as the long side of a right triangle with one of the short sides being
the north-south distance between the points and the other being the
east-west distance. (A right triangle is one that has a square
corner.) The usual formula for computing the length of the long side
of a right triangle is the Pythagorean Theorem. Using this formula
from geometry requires knowing about square roots.
Near the North Pole and near the South Pole, the longitude lines,
which go north-south and are called the meridians, approach each other
noticeably - in fact, they meet at the pole. The latitude lines, which
go east-west, are circles around the pole. Treating differences in
locations along these directions as if they were the sides of a right
triangle leads to errors in the computation of distance. Very close
to the pole, the answer could be VERY wrong - but a different
flat-Earth approximation, obtained from plane trigonometry, can be
used for short distances: the Polar Coordinate Flat-Earth Formula.
Using this formula - and the others mentioned below - requires knowing
something about trigonometry.
Latitude and longitude are spherical coordinates, based on recognition
that the Earth is round. Their definition does not require that the
Earth be exactly spherical, but approximating the Earth as a sphere is
satisfactory for most needs. Converting from spherical coordinates to
flat (Cartesian) coordinates is what most map projections are all
about; a lot of calculation is sometimes needed. The easiest is to use
north-south for one direction and east-west for the other - but that
has problems near the poles, as I described above.
When you study spherical trigonometry, you learn a bunch of formulas.
One of them is called the Law of Cosines for Spherical Trigonometry.
It is a perfectly fine formula when it is used for the right purposes.
Solving for short distances on the surface of the Earth is not one of
those purposes. The problem is as follows: Suppose you have a right
triangle with a very small angle. The ratio between the looooooong
short side and the long side is very close to 1.0. The formula
computes that ratio first, then requires the computer to find the angle
that has that ratio. In principle, the computer can do so - after all,
the formula is mathematically correct - but ordinary computers
approximate all numbers to some number of what are called "significant
digits". With 7 or 8 significant digits, the computer cannot
distinguish between the ratios for angles smaller than about a minute
of arc (a minute is 1/60 of a degree). Since the angle being computed
has its apex at the center of the Earth, a minute of arc corresponds to
about a mile on the surface.
Since the formula is mathematically correct, it can be manipulated
into other forms. The Haversine Formula is one result of such
manipulations. It has a similar problem, but it is "poorly
conditioned" when the two points are all the way around the Earth
from each other, rather than when they are close to each other. The
discussion below gives a second version of the haversine formula that
is easier to program on some computers.
It also discusses the fact that the Earth is not quite a sphere, and
gives several references for further reading. There is discussion of
where to look if even the ellipsoid approximation is too coarse. It
also has a pointer to one source of navigation formulas on the
Internet.
******** PROFESSIONAL-LEVEL DISCUSSION ********
CALCULATING DISTANCES ON THE SURFACE OF THE EARTH
What is the best way to calculate the distance between 2 points?
From: Bob Chamberlain <rgc@jpl.nasa.gov>
The distances considered here are along the surface of the Earth,
and deliberately ignore the effect of differences in elevation.
Distances on the surface of the terrain, whether geodesic, on roads,
or cross-country, depend on relief (including elevation differences),
the status of engineering projects, and perhaps even route selection.
Hence, computation is idiosyncratic and not well suited to simple
approximations.
If the distance is less than about 20 km (12 mi) and the locations of
the two points in Cartesian coordinates are X1,Y1 and X2,Y2 then the
Pythagorean Theorem
d = sqrt((X2 - X1)^2 + (Y2 - Y1)^2)
will require the least amount of computation and will be in error by
less than 30 meters (100 ft) for latitudes less than 70 degrees
less than 20 meters ( 66 ft) for latitudes less than 50 degrees
less than 9 meters ( 30 ft) for latitudes less than 30 degrees
(These error statements reflect both the convergence of the meridians
and the curvature of the parallels. The error is non-linear with
distance; shorter distances will have better percentage errors.)
The flat-Earth distance d will be expressed in the same units as the
coordinates.
If the locations are not already in Cartesian coordinates, the
computational cost of converting from spherical coordinates and then
using the flat-Earth model may exceed that of using the more accurate
spherical model.
Otherwise, presuming a spherical Earth with radius R (see below), and
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 * arcsin(min(1,sqrt(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.
When the two points are antipodal (on opposite sides of the Earth),
the Haversine Formula is ill-conditioned (see the discussion below
the Law of Cosines for Spherical Trigonometry), but the error, perhaps
as large as 2 km (1 mi), is in the context of a distance near
20,000 km (12,000 mi). Further, there is a possibility that round-off
errors might cause the value of sqrt(a) to exceed 1.0, which would
cause the inverse sine to crash without the bulletproofing provided by
the min() function.
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 Haversine Formula can be expressed in terms of a two-argument
inverse tangent function, atan2(y,x), instead of an inverse sine
as follows (no bulletproofing is needed for an inverse tangent):
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
In this context, a one-argument inverse tangent function would also
work: replace "atan2( sqrt(a), sqrt(1-a) )" by "arctan( sqrt(a) /
sqrt(1-a) )", but insert a test to ensure you will not divide by
zero if a = 1. (In the case a = 1, c = pi radians = 180 degrees,
and d is halfway around the Earth = 3.14159265... * R.)
The problem of determining the great circle distance on a sphere has
been around for hundreds of years, as have both the Law of Cosines
solution (given below but not recommended) and the Haversine Formula.
Sinnott gets the credit here because he was quoted by Snyder (cited
below). Perhaps someone will provide the truly seminal reference so
the proper attribution can be given?
The Pythagorean flat-Earth approximation assumes that meridians are
parallel, that the parallels of latitude are negligibly different from
great circles, and that great circles are negligibly different from
straight lines. Close to the poles, the parallels of latitude are not
only shorter than great circles, but indispensably curved. Taking this
into account leads to the use of polar coordinates and the planar law
of cosines for computing short distances near the poles: The
Polar Coordinate Flat-Earth Formula
a = pi/2 - lat1
b = pi/2 - lat2
c = sqrt( a^2 + b^2 - 2 * a * b * cos(lon2 - lon1) )
d = R * c
is computationally only a little more expensive than the Pythagorean
Theorem and will give smaller maximum errors for higher latitudes and
greater distances. The maximum errors, which depend upon azimuth in
addition to separation distance, are equal at 80 degrees latitude when
the separation is 33 km (20 mi), 82 degrees at 18 km (11 mi),
84 degrees at 9 km (5.4 mi). But even at 88 degrees the polar error
can be as large as 20 meters (66 ft) when the distance between the
points is 20 km (12 mi).
The latitudes lat1 and lat2 must be expressed in radians (see above);
pi/2 = 1.5707963. Again, the intermediate result c is the distance in
radians and the distance d is in the same units as R.
An UNRELIABLE way to calculate distance on a spherical Earth is the
Law of Cosines for Spherical Trigonometry
** NOT RECOMMENDED **
a = sin(lat1) * sin(lat2)
b = cos(lat1) * cos(lat2) * cos(lon2 - lon1)
c = arccos(a + b)
d = R * c
Although this formula is mathematically exact, it is unreliable for
small distances because the inverse cosine is ill-conditioned. Sinnott
(in the article cited above) offers the following table to illustrate
the point:
cos (5 degrees) = 0.996194698
cos (1 degree) = 0.999847695
cos (1 minute) = 0.9999999577
cos (1 second) = 0.9999999999882
cos (0.05 sec) = 0.999999999999971
A computer carrying seven significant figures cannot distinguish the
cosines of any distances smaller than about one minute of arc.
The function min(1,(a + b)) could replace (a + b) as the argument for
the inverse cosine to guard against possible round-off errors, but
doing so would be to "polish a cannonball".
5.1a: What value should I use for the radius of the Earth, R?
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
5.1b: When is it NOT okay to assume the Earth is a sphere?
A quick (?) test is: Compare the results produced by using the two
extreme values of the radius of curvature for the Earth:
minimum radius of curvature: 6336 km (3937 mi)
maximum radius of curvature: 6399 km (3976 mi)
in your application. If the results are different enough to cause you
to change your action (or your recommendation, or your interpretation
of the implication of the results, etc.), then assuming the Earth is
spherical is NOT okay.
The shape of the Earth is well approximated by an oblate spheroid. The
radius of curvature varies with direction and latitude. According to
formulas given on pages 24 and 25 of the book by Snyder,
"Map Projections - A Working Manual", by John P. Snyder,
U.S. Geological Survey Professional Paper 1395,
United States Government Printing Office, Washington DC, 1987,
the radius of curvature of an ellipsoidal Earth in the plane of the
meridian is given by
R' = a * (1 - e^2) / (1 - e^2 * (sin(lat))^2)^(3/2)
where a is the equatorial radius,
b is the polar radius, and
e is the eccentricity of the ellipsoid = sqrt(1 - b^2/a^2)
and the radius of curvature in a plane perpendicular to the meridian
and perpendicular to a plane tangent to the surface is given by
N = a/sqrt(1-e^2*(sin(lat)^2))
A Swedish book by Ilmar Ussisoo, _Kartprojektioner_ [map projections]
(published by the National Land Survey, Sweden, Professional papers
1977/6) suggests use of the geometric mean of these two radii of
curvature for all azimuths, as it produces errors of order of
magnitude 0.1% for distances within 500 km (300 mi) at 60 degrees
latitude. The formula for that average is no more complicated than
either of its components. That is, r = sqrt(R' * N) becomes
r = a * sqrt(1 - e^2) / (1 - e^2 * (sin(lat))^2)
Using these formulas with
a = 6378 km (3963 mi) Equatorial radius (surface to center distance)
b = 6357 km (3950 mi) Polar radius (surface to center distance)
e = 0.081082 Eccentricity
gives the following table of values for the
Radii of Curvature:
latitude...........r...................R'..................N..........
00 degrees . 6357 km (3950 mi) . 6336 km (3937 mi) . 6378 km (3963 mi)
15 degrees . 6360 km (3952 mi) . 6340 km (3940 mi) . 6379 km (3964 mi)
30 degrees . 6367 km (3957 mi) . 6352 km (3947 mi) . 6383 km (3966 mi)
45 degrees . 6378 km (3963 mi) . 6367 km (3957 mi) . 6389 km (3970 mi)
60 degrees . 6388 km (3970 mi) . 6383 km (3966 mi) . 6394 km (3973 mi)
75 degrees . 6396 km (3974 mi) . 6395 km (3974 mi) . 6398 km (3975 mi)
90 degrees . 6399 km (3976 mi) . 6399 km (3976 mi) . 6399 km (3976 mi)
Note that the radius of curvature for an ellipsoid is not the same as
the distance from the surface of the ellipsoid to the center. In fact,
the radius of curvature increases as the radius decreases. Also, be
aware that a variety of ellipsoids with slightly different parameters
have been fit to the Earth; the preferred ellipsoid may depend on the
region in which you are most interested.
Cliff Tiedemann, <clifft@uic.edu>, at the University of Illinois at
Chicago pointed out that spherical earth computations will provide
underestimates of real world distances measured in the direction of
the equator (and especially for trans-equatorial links) and
overestimates for those measured in the direction of the poles (and
especially for trans-polar ones).
For most purposes, it is quite satisfactory to treat the Earth as a
sphere. If not, an ellipsoid can provide a better approximation.
Some standard textbooks that may be helpful follow (reviews are by
Steve Robertson, <stever@mindlink.bc.ca>, of Tangent Survey Systems in
Canada):
Bomford, Guy 1980 _Geodesy_ Clarendon Press, Oxford
ISBN 0-19-851946-X
Review: For geodetic computations, this is pretty well the standard in
English. It's a cookbook and offers no development, however.
Vanicek, Petr, and Krakiwsky, Edward 1986 _Geodesy, the Concepts_
North-Holland, Amsterdam
ISBN 0-444-87775-4
Review: This offers a great, but quite involved, discussion of the
concepts behind geometrical (and all other) geodesy.
Torge, Wolfgang 1980 _Geodesy_ de Gruyter, Berlin (translated to
English by C. Jekeli)
ISBN 3-11-007232-7
Review: This concentrates mostly on gravimetric geodesy, but has some
geometrical stuff, well explained without too much mathematics.
Software for solving distance and azimuth problems on the ellipsoid
can be obtained (as of 10/10/96) by anonymous ftp from several
sources, two of which are listed below:
The URL of the National Geodetic Survey (of the National Oceanic and
Atmospheric Administration in the US Department of Commerce) is:
ftp://www.ngs.noaa.gov/pub/pcsoft/for_inv.3d/
Review (by Ronald C. McConnell, <rcmcc@cc.bellcore.com>, of Bellcore):
They have Fortran source and PC executable versions of both the normal
"inverse" great circle calculations (two lat/long pairs to distance
and bearing), and the less used "forward" calculation (one lat/long
pair plus bearing and distance to the second lat/long pair). They have
both 2-dimensional and 3-dimensional versions of each. The inverse
program works to within a few seconds or a few minutes, depending on
the Fortran compiler, of the antipodal points. The forward program
seems immune to any and all problem locations and pairs of locations.
You can choose among a couple of dozen ellipsoids.
See the read.me file for explanations. The NGS software directory may
contain other listing of interest. Its URL is:
http://www.ngs.noaa.gov/PC_PROD/pc_prod.shtml/
(Case is relevant in many URLs - e.g.: this one.)
The NGS has a public BBS at 301-713-4181 and provides FTP access at:
ftp://ftp.ngs.noaa.gov/pub/pcsoft
Another anonymous ftp source for ellipsoid software is the US
Geological Survey (of the US Department of the Interior), at:
http://kai.er.usgs.gov/pub/Proj.4/
Again, see the readme.txt file for explanations. The URL for the USGS
home page is:
http://kai.er.usgs.gov/homepage.html
5.1c: When is it NOT okay to assume the Earth is an ellipsoid?
The shape the Earth would assume if it were all measured at mean sea
level is called the geoid. The geoid varies no more than about a
hundred meters above or below a well-fitting ellipsoid, a variation
far less than the ellipsoid varies from the sphere. Terrain relief is
reported relative to the geoid. (Paraphrased from p. 11 of the book by
Snyder cited above.)
Distances on the surface of the geoid are not particularly meaningful.
However, there are applications, such as long-term prediction of
orbits of Earth satellites, that require better approximations than
are provided by an ellipsoid. Astrodynamics texts, such as
Kaula, William M. 1966 _Theory of Satellite Geodesy_ Blaisdell
Publishing Co., Waltham, MA (This book may be out of print.)
Battin, Richard H. 1964 _Astronautical Guidance_ McGraw-Hill,
New York (There may be later editions.)
may be consulted for further information.
5.1d: Where can I find formulas for great circle navigation problems,
like what course to follow and where will I be after following a known
azimuth for a given distance?
See Ed Williams' navigation formulas at
http://www.best.com/~williams/avform.htm
which can also be accessed by clicking on 'FAQ Lists' at
http://vancouver-webpages.com/peter/#faq
5.1e: How can I calculate distance on an ellipsoid.
On May 18, 1999, Steven Michael Robbins <stever@jeff.cs.mcgill.ca>
posted a discussion of this issue. He said:
About two weeks ago, I posted the following question:
> I need to compute the distance between two points on the surface of
> an ellipsoid. I'm looking for the distance *on* the surface. In
> short, I'm after the length of the shortest geodesic curve joining
> the points.
I was naively expecting a nice solution, akin to the "great circles"
on a sphere. I was soon disabused of this notion.
A big thanks to all who replied. Here is a summary of the responses.
I found two approaches to this problem.
One approach involves computing a discretized approximation to a
geodesic. This works on any ellipsoid. The first two software
pointers follow this approach. As far as I can tell, all the other
references follow the second approach.
The second approach is designed for the Earth, exploiting two special
features: (1) it is an ellipsoid of revolution, and (2) it is nearly
spherical. These properties allow one to construct a solution as a
power series in the small "flattening" parameter. One truncates the
series depending upon the accuracy required.
I was looking for a solution that works for any ellipsoid. Most of
what follows is geodesy-centric, so I have only explored a handful of
the suggestions listed below. Also --- as one person noted --- a web
search with key words of "ellipsoid & distance & arc" will generate a
list of hundreds of sites. (All dealing with the Earth, I might add).
What follows is heavily cut'n'pasted from the responses or from web
pages. Don't take these as advice or endorsement from me -- these are
not my words below!
Available Software on the net:
http://www.magic-software.com/gr_surf.htm
The code here computes a polygonal approximation to a geodesic and
then refines it until the geodesic curvature vanishes at each vertex
on the path.
http://www.netlib.org/ode/geodesic/
Similar to the above, except that it actually works on any manifold --
you supply a routine to compute the metric at any point, etc.
ftp://kai.er.usgs.gov/pub/PROJ.4/
Well-proven code, forward and inverse problem.
ftp://ftp.ngs.noaa.gov/pub/pcsoft/for_inv.3d
>From the (U.S.) National Geodetic Survey:
INVERSE/FORWARD3D (Version 1.0)
Comprises four programs - Inverse (Version 2.0) which computes the
geodetic azimuth and distance between two points, given their
geographic positions; Forward (Version 2.0) which computes the
geographic position of a point, given the geodetic azimuth and
distance from a point with known geographic position; and the
three-dimensional versions of these programs . INVERS3D (Version 1.0)
and FORWRD3D (Version 1.0), which include the height component.
ftp://mapping.usgs.gov/pub/software/current_software/w5501
This program will solve either the geodetic position or inverse
problems in any quadrant of the Earth using the parameters of one of
five commonly-used ellipsoids.
http://www.globalserve.net/~nac/products.html
A product called NACNav (60USD/license) which can calculate the
shortest distance between two points anywhere on the earth surface and
the angle between the direction and the magnetic north (determined by
local magnetic declination).
Suggested References
If you need very high accuracy, I suggest reading B.R. Bowring, "Total
Inverse Solution for the Geodesic and Great Elliptic", Survey Review,
33, 261 (July 1996).
If you are satisfied with the relative error of the square of the
Earth's flattening, then see J. Meeus, Astronomical Algorithms, 1991.
Two people recommended: Bomford, G., 1952, Geodesy: London, England,
Oxford University Press.
Clark, D., 1963, Plane and Geodetic Surveying, II: London, England,
Constable & Co., Ltd.
Hosmer, G.L., 1930, Geodesy: New York, New York, John Wiley & Sons, Inc.
Lambert, W.D. and Swick, C.H.,1935, Formulas and Tables for the
Computation of Geodetic Positions on the International Ellipsoid:
U.S. Coast and Geodetic Survey Special Publication No. 200.
U.S.C.& G.S., Formulas and Tables for the Computation of Geodetic
Positions, U.S. Coast and Geodetic Survey Special Publication No. 8.
"Direct and Inverse Solutions of Geodesics on the Ellipsoid with
Application of Nested Equations", T. Vincenty, Survey Review 22, April
1975.
Thomas, Paul D., 1965, Geodesic arc-length on the reference ellipsoid
to second-order terms in the flattening, Journal of Geophysical
Research, vol.70(14), 3331-3340.
========================
This answer was prepared by Robert G. Chamberlain of Caltech (JPL),
<rgc@jpl.nasa.gov>, and reviewed on the comp.infosystems.gis newsgroup
in Oct 1996. It was revised in November 1997 and February 1998.
Formatting was updated and the hyperlinks were checked in January
2000. Minor cosmetic changes were made and hyperlinks were updated
in February 2001.
The inverse tangent formula and the pointer to Ed Williams' formulas
were added in November 1997, and the "min(1,sqrt(a))" bulletproofing
was removed from the haversine formula.
Mikael Rittri, <Mikael.Rittri@carmenta.se>, pointed out that the simple
approximation to the radius of the Earth estimated the distance
between the surface and the center, rather than the radius of
curvature. This approximation (and the associated warning) were
removed in February 1998; the azimuth-averaged radius was offered in
its stead, as he suggested. Mikael also noted that the inverse tangent
formula does not need an adjustment for quadrant, which allowed
simplification of that formula. He also convinced me that it _might_
be possible for the "sqrt(a)" term in the haversine formula to exceed
1, so I reinstated the "min(1,sqrt(a))" bulletproofing, ugly though it
is.
Steven Michael Robbins' discussion of distance computation on an
ellipsoid was added in May 1999.
The pre-geometry summary was added in January 2000 in response to a
request from Karen Kast, who was seeking information that might be
useful to a 6th grade student doing a science project.
Following is a Mirror of
http://home.pipeline.com/~hbaker/FAQ-lat-long.txt.
Draft FAQ: Computing "Great Circle Distances" from latitudes and longitudes.
This problem (and its solutions) are more than 500 years old.
Henry Baker (hbaker1@pipeline.com), May 1995.
Problem: Compute the great circle distance (either angular or linear)
between the two points whose coordinates are given by
(latitude,Longitude) pairs. Specifically, given P1=(l1,L1),
P2=(l2,L2), compute the angular distance d between these two points on
the globe. The linear distance is then r*d, where r is the radius of
the globe (assuming that d is given in radians).
Plan: Transform coordinates into Cartesian coordinates, compute the
"dot" product, and then take the arccosine of the result.
We can either perform laborious trigonometric calculations, or we can
utilize geometric insight to realize that the great circle distance
doesn't depend upon the actual values of the longitudes L1,L2, but
only on their _difference_ dL=L2-L1. We therefore translate the
problem to P1'=(l1,0), P2'=(l2,L2-L1)=(l2,dL).
P1': x1=cos(l1), y1=0, z1=sin(l1).
P2': x2=cos(l2)cos(dL), y2=cos(l2)sin(dL), z3=sin(l2).
P1'.P2' = cos(d) = cos(l1)cos(l2)cos(dL) + sin(l1)sin(l2).
We are now mathematically (but not computationally) done, because we
have the following expression for the angular distance we seek:
d = arccos(cos(l1)cos(l2)cos(dL) + sin(l1)sin(l2)).
Discussion: There is a problem with this formula, however. If dL is
very small, cos(dL) is essentially 1, and arccos(1)=0. In other
words, we lose all accuracy in the (common) case where dL and |l1-l2|
are both small.
The problem is that cosines and arccosines do not have the property
that f(0)=0. What we really want is a formula based on sines/arcsines
or tangents/arctangents, which _do_ have this property.
We therefore use the following centuries-old half-angle trick:
sin(dL/2)^2 = (1-cos(dL))/2, or equivalently, cos(dL) = 1-2sin(dL/2)^2.
cos(d) = cos(l1)cos(l2)(1 - 2 sin(dL/2)^2) + sin(l1)sin(l2)
= cos(l1)cos(l2) - 2 cos(l1)cos(l2)sin(dL/2)^2 + sin(l1)sin(l2)
(recognizing the cosine difference formula)
= cos(l2-l1) - 2 cos(l1)cos(l2)sin(dL/2)^2
= cos(dl) - 2 cos(l1)cos(l2)sin(dL/2)^2 [dl = l1-l2]
(using the half-angle trick again with dl)
= 1 - 2 sin(dl/2)^2 - 2 cos(l1)cos(l2)sin(dL/2)^2
(using the half-angle trick a 3rd time with d)
1 - 2 sin(d/2)^2 = 1 - 2 sin(dl/2)^2 - 2 cos(l1)cos(l2)sin(dL/2)^2
sin(d/2)^2 = sin(dl/2)^2 + cos(l1)cos(l2)sin(dL/2)^2
(This now looks a lot like a Pythagorean formula, except with
half-sines and a correction coefficient of cos(l1)cos(l2).)
sin(d/2) = sqrt(sin(dl/2)^2 + cos(l1)cos(l2)sin(dL/2)^2)
d = 2 asin( sqrt( sin(dl/2)^2 + cos(l1)cos(l2)sin(dL/2)^2 ) )
This formula is now accurate when dl and/or dL are near zero.
[As a check, if l1=l2=l, this formula simplifies to:
d = 2 asin(min(1,cos(l)sin(dL/2))). If l=0, then d = dL.
Furthermore, if dL=0, this formula simplifies to:
d = 2 asin(sin(dl/2)) = 2 (dl/2) = dl.]
Furthermore, sin()^2 >=0, and cos() >=0, since -pi/2 <= l1,l2 <= pi/2,
so the argument to the square root is always >=0. It is also
mathematically <=1, but errors may violate this by small amounts, so
use min to guarantee that the argument to asin is in the range [0,1].
We finally have:
d = 2 asin(min(1,sqrt(sin(dl/2)^2 + cos(l1)cos(l2)sin(dL/2)^2)))
Although in current subroutine libraries this appears to be a more
expensive calculation, classical books tabulate the function
versine(x) = vers(x) = 2 haversine(x) = 2 sin(x/2)^2 and its inverse.
Therefore, if you're going to do this calculation a lot, you should
also have these specialized subroutines for versine and its inverse.
Since vers(x) = 2 sin(x/2)^2 = 1 - cos(x), its Taylor series is
x^2/2! - x^4/4! + x^6/6! - x^8/8! + ... = cos(0) - cos(x),
so versine is no more difficult to compute than cosine. Versine is
more useful, however, because it is much more accurate near zero.
Since accuracy is hard to get and easy to lose, the primitive
trignometric functions found in your subroutine library should be
_versine_ and _arcversine_ instead of _cosine_ and _arccosine_, and
those (usually inappropriately) desiring cosine and arccosine can
easily construct them as the macros:
cos(x) = (1 - vers(x))
acos(x) = (pi/2 - asin(x)).
[Also realize that the Earth is not really a sphere, but an oblate
spheroid, so even these calculations are not accurate for the Earth.]
[Historial note: The ancient mathematicians (astronomers and
astrologers, mostly, who actually had to _compute numerically_ for a
living instead of doing just symbolic algebra) used the function
chord(x) = cd(x) = 2 sin(x/2) instead of sin(x). Using cd(x), the
Pythagoreanness of the formula above is particularly striking:
cd(d)^2 = cd(dl)^2 + cos(l1) cos(l2) cd(dL)^2
= cd(dl)^2 + (1-vers(l1)) (1-vers(l2)) cd(dL)^2
The Arab mathematicians circa 1000 AD used the "versed sine", or
versine, in addition to the sine and chord. The cosine was almost
never used, and apparently became popular much later as a result of
Euler's identity: exp(iy) = cos(y) + i sin(y). In other words,
cosines are important for _algebraic manipulation_, but not for
_numerical calculation_. Apparently, those doing numerical
calculations on computers still have much to learn from the ancient
human computers.
Reference: Smith, David E. History of Mathematics, Vol. II. Dover
Publs, NY, 1925.]
Copyright (c) 1995 by Henry Baker. All rights reserved.