Welcome to the NavList Message Boards.

NavList:

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

Compose Your Message

Message:αβγ
Message:abc
Add Images & Files
    or...
       
    Reply
    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
    >destination."
    >
    >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
    >answer.
    >
    >=================================
    >
    >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
    >destination.
    >
    >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
    >direction.
    >
    >==================
    >
    >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---.u-net.com
    >George Huxtable, 1 Sandy Lane, Southmoor, Abingdon, Oxon OX13 5HX, UK.
    >Tel. 01865 820222 or (int.) +44 1865 820222.
    >------------------------------
    
    Brian Whatcott
       Altus OK                      Eureka!
    

       
    Reply
    Browse Files

    Drop Files

    NavList

    What is NavList?

    Join NavList

    Name:
    (please, no nicknames or handles)
    Email:
    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:

    Email Settings

    Posting Code:

    Custom Index

    Subject:
    Author:
    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