Welcome to the NavList Message Boards.


A Community Devoted to the Preservation and Practice of Celestial Navigation and Other Methods of Traditional Wayfinding

Compose Your Message

Add Images & Files
    Re: Haversine formulae for Great Circles
    From: Brian Whatcott
    Date: 2001 Nov 19, 5:51 AM

    It's rare to find such throw away insightfulness
    (Ah! So that's what haversines are all about...)
    I suspect that the code in a few GPSs could use a little tweak in the
    way George mentioned (taken with his errata...)
    Thanks a lot
    Brian Whatcott
    At 09:04 AM 11/18/01, you wrote:
    >What an interesting contribution from Lu Abel! (copied below). Perhaps I
    >can add something to it.
    >When I first scanned it through, my heart sank slightly. Not haversines
    >again! I said to myself. Surely, those went out with logarithms, as soon as
    >we started to use calculators and computers.
    >Navigators used logs to get round the difficulties they always had doing
    >long-multiplication and long-division. The problem was that logarithms
    >couldn't be used with negative numbers, the log of a negative number being
    >a meaningless concept. As ordinary trig functions ranged over positive and
    >negative values, something had to be done.
    >So a new trig function, the "versine" was invented, which never went
    >negative, but varied between zero and +2. The haversine, as you might
    >guess, is simply half the versine, so it varies, more conveniently, between
    >zero and +1. And then all the trig formulae that navigators used were bent
    >and twisted into forms that used one of these new functions, and nothing
    >that required a log ever took a negative value.
    >Since calculators took over, the need for logs vanished, and haversines
    >vanished with them.
    >Lu has pointed to a difficulty in the use of the standard great-circle
    >distance formula. Using his own notation-
    >"L1, Lo1 are L/Lo of the starting point, L2/Lo2 are the L/Lo of the
    >the usual expression for the distance in miles is-
    >60 * acs ((sin L2 * sin L1) + (cos L2 * cos L1 * cos (Lo2 - Lo1))),
    >where the angles are in degrees and acs means "the angle whose cos is..."
    >or arc.cos.
    >As Lu has said, this presents difficulties in situations where the
    >destination is close to the starting point. His explanation omits an
    >important matter. The real problem arises in that the expression above is
    >trying to derive an angle from its cosine, when that angle is small. In
    >that case, the cosine is very near to +1, and if you plot out a cosine, you
    >will see that it changes very little from +1 as the angle increases. For
    >small angles, then, the cosine is a very "flat" function, so it just isn't
    >possible to derive an angle from its cosine when that angle is very small.
    >It's not just a problem with a computer: you would meet it if you tried to
    >use cosine tables instead.
    >I have tried the above formula with a Casio programmable calculator, which
    >has an internal precision of 12 decimal digits. Even with this high
    >accuracy, distances of 0.1 nautical miles have started to go a bit
    >inaccurate, and distances which should be .01 miles are calculated as zero.
    >With a computer using single-precision, things would be a lot worse. So the
    >problem Lu has raised is a real one, in these special circumstances of
    >close approach.
    >Now for the alternative formula that Lu quotes, to which he gives the name
    >"Haversine" formula. I regret that choice of name, though I can see how it
    >has arisen. It's distinctly off-putting, to a navigator who has never been
    >exposed to the concept of haversines, as it invokes the magic and mystery
    >of the occult craft of navigation. And it's quite unnecessary, as to use
    >the formula that he provides, you need no knowledge about haversines at
    >all. Haversines are not mentioned in the formula. So let's keep things as
    >simple as possible, and leave haversines out of it altogether, shall we?
    >For anyone that would like to use the procedure that Lu quotes as a single
    >line in a program, we can recast his formula slightly.
    >The distance in miles, when angles are in degrees, is-
    >120*asn(sqr(sin((L2-L1)/2)^2 + cosL2*cosL1*(sin((Lo2-Lo1)/2)^2)))
    >where asn means "the angle whose sine is..." or arc.sin.
    >This formula works beautifully for small angles, preserving full accuracy.
    >Does it have any snags? Well just as the cos function becomes flat for
    >angles near zero degrees, the sin function becomes flat near 90 degrees,
    >and this happens when the distance nears 10800 miles, on the other side of
    >the world. It must be rare for an accurate distance to be required to a
    >destination near the observer's antipode, though perhaps an accurate
    >azimuth may occasionally be of interest (see footnote below). So the
    >formula that Lu quotes remains useful over any distance that is of real
    >interest. He has highlighted a significant problem and provided a useful
    >At the end, Lu asks-
    >"Are there equivalent "haversine" formulae for initial direction and the
    >L/Lo of intermediate points?"
    >Consider the problem of calculating the observer's azimuth of a heavenly
    >body, which is almost the same as calculating the initial direction to a
    >Equivalent quantities in these two problems, using Lu's notation, are-
    >L1=lat, L2=dec, (Lo2-Lo1)=hour angle, distance = zenith angle= 90-altitude.
    >The distance here should be expressed in degrees (of 60 miles) rather than
    >in miles.
    >Latitudes and declinations should be given signs + or - if they are N or S.
    >Hour angles are measured westwards.
    >Having calculated the altitude of a body, many navigators seem to obtain
    >the azimuth from-
    >asn ((cos (90 - hour angle)*cos dec) / cos altitude)
    >This is a terrible choice of formula. It gives a completely ambiguous
    >result, for any azimuths that are near 90 degrees (and also near 270). For
    >example, an angle of 80 degrees has exactly the same sine as an angle of
    >100, so the formula is quite unable to distinguish between these two
    >solutions. And I know of no way of distinguishing between these solutions
    >"by inspection". If anyone else knows how to do this, I would be interested
    >to learn. Avoid this method
    >Another option is to obtain azimuth from-
    >acos((sin L2 - (sin L1 cos distance)) / (sin distance cos L1))
    >This has a similar problem that we considered above, that for some azimuths
    >(i.e. directions very near due North and South) this can become rather
    >inaccurate because cos of a small angle varies so little from 1.
    >But why not obtain azimuth from an arc.tan formula, that has no such
    >problems? Obtain azimuth, in degrees, for an observed body, from-
    >atn( sin hour angle / (cos ha sin lat - cos lat tan dec))
    >or, for a great-circle, initial azimuth is-
    >atn( sin (Lo2 - Lo1) / (cos (Lo2 - Lo1) sin L1 - cos L1 tan L2)
    >Hour angles and longitudes are 0 to 360 measured westwards (not a
    >convention that is accepted by all), so any E longitudes are negative. You
    >may note that as there is no mention of altitude (or distance) in this arc
    >tan formula, you can obtain the azimuth directly, without having to
    >previously calculate those quantities beforehand.
    >When making the above calculation, observe the sign of tan az, and if it's
    >negative, add 180 degrees to the result. If the hour angle was negative,
    >add another 180 degrees. The result should be the true azimuth measured 0
    >to 360 clockwise from North, with no ambiguities.
    >It's even simpler if your calculator or computer has rectangular-polar
    >conversion (sometimes labelled the atn2 function). In that case you put
    >into the y-coordinate the value (- sin hour angle) or (sin (L1 - L2)), and
    >into x goes cos hour angle*sin lat - cos lat*tan dec (or its equivalent).
    >The result emerges as a radius (which you ignore) and an angle, which is
    >the angle you require. All you have to do (it you're bothered to) is to add
    >360 if it is negative.
    >A footnote about calculations involving positions close to the antipode.
    >These are required rarely, but I can think of an example that came up on
    >another mailing list.
    >Muslims wish to face Mecca when they pray, but this must present problems
    >if they live on the other side of the world. There must be a
    >jump-discontinuity in the great-circle azimuth of Mecca at some spot at
    >Mecca's antipode, which one could consider as anti-Mecca. Near there, the
    >faithful would need to face away from that spot, at which one might erect
    >some sort of marker, if it was on land. Fortunately it isn't; it's
    >somewhere between the Tuamotus and the Gambier islands in the Pacific. But
    >imagine the problems that must face a pious Muslim crew member of an
    >inter-island vessel plying those waters, if he wishes to pray in the right
    >Note: everything in this posting assumes a spherical Earth. If its
    >ellipticity is taken into account, everything will get much more complex,
    >including the direction of Mecca.
    >George Huxtable.
    >Lu Abel said-
    > >Almost every celestial text provides the formulae for calculating great
    > >circle distance and direction between two points.  This is not surprising
    > >since the trigonometry is identical to that of sight reduction.  If one
    > >takes any standard celestial reduction method and substitutes the starting
    > >point for the AP and the ending point for the body GP, Zn gives the
    > >direction and 60*(90-Hc) gives the distance.
    > >
    > >A couple of years ago this august group was of great help in answering for
    > >me an obvious question I've never seen answered in any standard navigation
    > >text, namely how to calculate the L/Lo of intermediate points.
    > >
    > >Recently I've done a bit of further digging on great circles.  I learned
    > >that while the law of cosines formula typically used for Hc calculation is
    > >trigonometrically correct, it can produce incorrect answers on a computer
    > >or calculator when the starting and ending points are close because
    > >computers and calculators express numbers with a limited number of
    > >significant digits.  An alternate is the "Haversine" formula, called by
    > >that name because of its haversine-like terms  [recall that hav(x) =
    > >(1-cos(x))/2 = sin^2(x/2) where "^2" means squared ]
    > >
    > >The haversine distance formula goes as follows:
    > >
    > >L1, Lo1 are L/Lo of the starting point, L2/Lo2 are the L/Lo of the
    > destination
    > >
    > >DLat = L2 - L1
    > >DLo = Lo2 - Lo1
    > >
    > >A = sin^2(DLat/2) + cos(L2) * cos(L1) * sin^2(DLo/2)
    > >
    > >Distance  =  120 * arcsin (sqrt (A))
    > >
    > >
    > >For a great writeup on the haversine distance formula, see
    > >http://mathforum.org/dr.math/problems/neff.04.21.99.html (note no "www" at
    > >beginning)
    > >
    > >Having learned about this formula, I'm going to guess it's used in the
    > >majority of GPS sets, since they're often calculating small distances (like
    > >distance-to-go when approaching a waypoint)
    > >
    > >The standard reference on calculating great circles via the haversine
    > >formula seems to be R.W. Sinnott, "Virtues of the Haversine", Sky and
    > >Telescope, vol. 68, no. 2, 1984, p. 159.  I've seen it mentioned in several
    > >writeups of the haversine distance formula.
    > >
    > >I can't access a copy to see if Sinnott provided formulae beyond the
    > >distance formula, so here's my question for this group:
    > >
    > >Are there equivalent "haversine" formulae for initial direction and the
    > >L/Lo of intermediate points?
    > >
    > >Thanks
    > >
    > >Lu Abel
    >George Huxtable, 1 Sandy Lane, Southmoor, Abingdon, Oxon OX13 5HX, UK.
    >Tel. 01865 820222 or (int.) +44 1865 820222.
    Brian Whatcott
       Altus OK                      Eureka!

    Browse Files

    Drop Files


    What is NavList?

    Join NavList

    (please, no nicknames or handles)
    Do you want to receive all group messages by email?
    Yes No

    You can also join by posting. Your first on-topic post automatically makes you a member.

    Posting Code

    Enter the email address associated with your NavList messages. Your posting code will be emailed to you immediately.

    Email Settings

    Posting Code:

    Custom Index

    Start date: (yyyymm dd)
    End date: (yyyymm dd)

    Visit this site
    Visit this site
    Visit this site
    Visit this site
    Visit this site
    Visit this site