 NavList:

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

Message:αβγ
Message:abc
 Add Images & Files Posting Code: Name: Email:
Re: Sight Reduction Formula
From: Herbert Prinz
Date: 2003 Oct 22, 18:04 -0400

"Clampitt, Ralph" wrote:

> I am new to the list and have just recently started using a programmable
calculator for sight reductions.  I have found several useful formulae on
this list and in other references, however the one shown in the Nautical
Almanac for azimuth puzzles me.
> Can anyone explain the formula for azimuth from the Almanac?

With pleasure. But let me state one caveat at the outset. The mathematics
behind the formula is an objective fact. What I am going to say is valid,
unless I commit an error of logic. When I proceed to discuss the pragmatic
aspects, i.e. why the formula is preferable to another, we enter the
realm of personal opinion. Chances are that the reasoning of the authors of
the Nautical Almanac was similar to mine; but the true motivation for their
choice of a formula would have to be extracted from an HMNAO document
describing their algorithm. I am sure that such a document exists (at
least internally), whether it is readily available to the public, I do not know.

In the following, let

L = Latitude
D = Declination
T = Time angle (=local hour angle)
H = computed Altitude

The standard cosine formula for the azimuth, which follows trivially from the
cosine formula for the spherical triangle, is:

(I )    cos A = (sinD - sinL sinH) / (cosH cosL)

The formula on which the algorithm in the N.A. is based is

(II)    cos A = (sinD cosL - cosD sinL cosT) / cosH

The latter can be derived from the former through algebraic transformations.
It can also be demonstrated geometrically. (II) is an improvement over (I),
its advantage being the absence of the term cosL in the denominator. Formula
(II) is probably the best choice for the base of a reliable
and general algorithm that does not require extreme accuracy.

1. Equivalence of (I) and (II).

One way to "explain" a formula is to show that it is equivalent to another one
that we have already come to accept. This is what I will do first, because it
gives me a chance to point out the decisive difference between the formulae.

From Pythagoras:
cosL cosL  =  1 - sinL sinL

multiply with sin(D):
sinD cosL cosL  =  sinD - sinD sinL sinL

subtract cos(D)*sin(L)*cos(L)*cos(T):
sinD cosL cosL - cosD sinL cosL cosT = sinD - sinD sinL sinL - cosD sinL cosL cosT

collect terms in each member:
cosL * (sinD cosL - cosD sinL cosT) = sinD - sinL * (sinD sinL + cosD cosL cosT)

substitute sinH for sinD sinL + cosD cosL cosT in the second member:
cosL * (sinD cosL - cosD sinL cosT) = sinD - sinL sinH

divide by (cosL*cosH), assuming cosH cosL <> 0, thus H <> 90deg AND L <> 90deg !!!
(sinD cosL - cosD sinL cosT) / cosH = (sinD - sinL sinH) / (cosH cosL)

On the left side we now have the expression for the cosine of the azimuth from
(II), and on the right side that from (I). We have thus shown that one can be
derived from the other, for all altitudes and latitudes different from 90
deg.

2. Singularities. The essential difference between (I) and (II).

It is meaningless to speak of the azimuth of an object at 90 deg altitude,
i.e. an object directly in the zenith. There is also no way of consistently
assigning a value (such as 0 deg, or whatever) to the azimuth of such an
object by means of an arbitrary definition. Therefore, it is quite
appropriate that none of the above formulae (or any other one, for that
matter) work for H = 90 deg. This case must always be trapped and handled
specially, regardless of which formula is being used.

One might expect a similar problem for an observer directly at one of the
poles. But the situation is different. We immediately see that formula (I)
falters for a latitude of +/-90 deg, while (II) still returns some value. The
exact result it returns we find by substituting 0 for cosL, and
either +1 or -1 for sinL into (II):

cos A = (- cosD cosT) / cosH             .... North pole
cos A = (cosD cosT) / cosH               .... South pole

But an observer at the pole always finds H = Abs(D), hence cosH = cosD:

cos A = (- cosD cosT) / cosD             .... North pole
cos A = (cosD cosT) / cosD               .... South pole

which reduces to

cos A = - cosT             .... North pole
cos A = cosT               .... South pole

and further to

A = 180 - T             .... North pole
A = T                   .... South pole

Although it might at first seem meaningless to speak of azimuth when the
direction towards north is not defined, accepting the above relation allows
the algorithm in the N.A. to be applied smoothly to polar navigation without
modification. Note that in such an application the observer at the
pole is not the rare exception, but the rule, as it is standard practice to
identify the assumed position with the pole.

Besides being able to use the exact location of the pole as assumed position,
there may be an additional benefit for positions near it. Algorithms tend to
become numerically unstable in the vicinity of singularities if one does not
take great care in which order expressions are evaluated.

3. Geometrical derivation

The rather abstract derivation of formula (II) from (I)  may be pleasing to
some, others might be more satisfied to see a physical interpretation of it.
Before we are going to draw triangles, let me give a motivation for what I am
aiming for.

Looking at the term (sinD cosL - cosD sinL cosT) in formula (II), we recognize
some vague similarity with the form of the cosine formula for spher.
triangles. We could introduce an auxiliary variable X, such that

(A)    cosX = sinD cosL - cosD sinL cosT

and hope that we can press (A) into the shape of the cosine theorem. After
having done so, we can rewrite (II) as

(B)    cosA = cosX / cosH

which is the sine theorem for spherical triangles in disguise. A little kneading is called for:

(A)    cosX = sinD cosL - cosD sinL cosT
cosX = sinD cosL + cosD sinL cos(180-T)            cos(u) = -cos(180-u)
(A')   cosX = cos(90-D) cosL + sin(90-D) sinL cos(180-T)  sin(u) = cos(90-u)

(B)    cosA = cosX / cosH
sin(90-A) / 1 = sin(90-X) / sin(90-H)              cos(u) = sin(90-u)
(B')   sin(90-A) : sin(90) = sin(90-X) : sin(90-H)        rewriting as proportion

The latter is an application of the well known sine theorem.

So, first we are looking for a spherical triangle that is made up of  one side
equalling in length the latitude of the observer, another side equal to the
co-declination of the star, and the enclosed angle equal to the supplement of
the local hour angle. The third side of this triangle will
then be our auxiliary variable X. Once we have got this, we find a rectangular
spherical triangle with the side 90-H (= zenith distance) opposite the right
angle, and  the side 90-X  (=  the complement of the auxiliary variable)
opposite the angle which is the complement of the azimuth.

Now we proceed to construct the corresponding triangles. The first one can be
drawn on a globe as follows:

Mark P, the pole.
Mark O, the place of the observer at latitude L. Draw the meridian through O and P.
Mark S, the star at declination D. Draw the meridian through S and P.

Draw the prime vertical through the observer. Remember that the prime vertical
is a great circle perpendicular to the meridian.

Draw the pole of the prime vertical. Since this point is 90 degrees away from
any point on the prime vertical, it can be found by extending the meridian OP
beyond the pole by the length L. Mark this point V.

The triangle SPV is the first one that we are looking for. Obviously, PV
equals L, PS equals (90-D), and the angle SPV is (180-T). The arc VS then
represents the auxiliary variable X. It can be computed by means of formula
(A') given above.

To construct the second triangle,  intersect the great circle VS with the
prime vertical. Mark the intersection point as W. Since V is the pole of the
prime vertical, VW intersects the latter at a right angle. Furthermore, the
arc VW measures 90 degrees.

The triangle OWS is the one that has been sought. Arc OS is the zenith
distance of the star. It measures 90 - H degrees. This arc is opposite the
right angle OWS. Arc WS is equal to VW - X and therefore measures 90 - X
degrees. This arc is opposite the angle SOW, which is the angle between
the prime vertical and the direction of the star and hence the complement of
the azimuth. So, we know two sides and an angle that is not enclosed. The
triangle OWS can be resolved by means of (B') given above.

4. Numerical considerations.

The formula is part of a sight reduction algorithm for scientific calculators.
It can be used as a recipe for manual computation, but more likely is to be
programmed into a programmable calculator. (The same formula is also
recommended for the use in personal computers in HMNAO, "AstroNavPC,
Astro-Navigation Methods and Software for the PC", Willmann-Bell, 2000). Since
no particular hardware or software environment is addressed, the algorithm
has to be general enough in that it does not rely on the behaviour of any
particular machine or on the presence of any but the most basic
s/w  features.

The algorithm must also be robust enough to provide reliable results when
implemented by a competent engineer or navigator who is not necessarily a
skilled computer programmer. "Reliable" normally means reducing the branching
logic to an absolute minimum and correct handling of special
cases. (Note that this issue of "cases" is at the core of the vast literature
on ever "better" sight reduction methods in the 18th and 19th century.)

Finally, the algorithm must be designed so that any reasonable implementation
may provide the required accuracy of the results. The most accurate tables
that are commonly used in celestial navigation, Ho 229, provide an accuracy
of 0.1 degrees in azimuth. It is fair to assume that this is
sufficient. The cosines of two angles differing by 0.1 degrees differ from
each other at least in the 5th decimal. Every calculator can handle this. The
use of cosines is therefore perfectly adequate as far as accuracy is
concerned.

A more important aspect is numerical stability in limiting cases. Consider the expression

(sinD cosL - cosD sinL cosT) / cosH

The theoretical range of values of this expression is  -1 to +1. However, due
to rounding errors, the actual value that will be computed will almost always
differ from the correct value in the last significant digit of the internal
word. This poses a problem when the arc cosine function (the
inverse of the cosine) of the above expression is to be evaluated. The
slightest rounding error can bring the value for the expression just a tiny
bit outside of the permissible range where the arcus function is defined. The
consequences are fatal. In the best case an exception will be
raised. On the sillier calculators such as the HP48, the program will happily
continue in the domain of complex numbers, with highly questionable benefit
for the average navigator. This phenomenon has been well taken care of in the
algorithm that appears in the N.A.. The critical expression
is checked against the limits -1 and 1. If those are exceeded, the actual
value is replaced by the appropriate limit before the arcus function is
evaluated.

5. Tangent as alternative

George Huxtable recommends a formula given by Meeus employing the tangent
function as a "better" alternative to (II). He is correct in asserting that
this method provides better accuracy. But Meeus writes for the astronomer who
wants milliarcsecond precision. The navigator needs a tenth of a
degree. If the improved accuracy has to be bought at the cost of other
features, it might not be worth it.

George also points out that the tangent formula

(III)    Az = arctan ((-Sin HA) / (Cos Lat * Tan Dec - Sin Lat * Cos HA))

does not require the computed altitude as input. True again, but in the
context of  the sight reduction procedure we already have the altitude
anyway. On the other hand we don't have the tangent of the declination. So,
in the context in which our formula is proposed, this would be an extra
step, making the procedure as a whole longer.

Formula (III) uses the tangent of the declination. This function has a pole at
90 degrees. Admittedly, there is no star located exactly on either poles. But
one never knows how a function ends up being used, once it has been
programmed into a calculator. Just like HO. 229 can be used for
great circle sailing, so can (III). In fact, in my hp48, I use a common
subroutine for both purposes. A correct implementation of (III) must check
for D=90 and handle it as a special case. (II) does not suffer from this
problem.

A further problem arises if the denominator becomes zero or very small. The
first condition will produce a zero divide exception, the second may cause an
overflow. George warns us thus to check the denominator for zero. Since (II)
also has a denominator that can become zero, we have to
perform a similar check there too. But note the difference:

In (II) we check in order to catch a meaningless result that cannot be handled
(star in zenith). In (III), the check is performed in order to catch a
perfectly valid situation that the formula is not able to handle (body on
prime vertical). That means that in the latter case an additional
check for body in zenith has to be made eventually. The extra steps to handle
special cases make the algorithm cumbersome and diminish reliability. Even
the programmable ones amongst the calculators are much better suited to
evaluating formulae serially than for heavy program logic.

A similar argument applies to the selection of the quadrant. The cosine being
unique and invertible over the range from 0 to 180 degrees provides a more
fool proof solution than the tangent. Using the cosine, one needs to consider
only that the final result for the azimuth must be consistent
with the hemisphere in which the star is located. The azimuth of a body is
smaller than 180 deg if it's east of the observer; larger than 180 if it's
west. The solution via tangent requires an extra decision as to which
quadrant is to be used.

Functions like ATAN2, polar to rect conversion, etc. are not universal enough
to be relied on in a general procedure such as the one in the nautical
almanac. Some calculators don't have them, in some their implementation is
awkward (to say the least) as anybody with a hp48 will confirm.
Nevertheless, the HMNAO publication, op. cit., p.65, which is partly addressed
to computer programmers that have access to standard Fortran or C libraries
does actually mention cartesian to polar co-ordinate conversion as a
possibility. It is illuminating that their solution is NOT to simply
use the numerator and denominator of  formula (III), as suggested by George Huxtable.

Instead, they propose  to set:

x = cosL sinD - sinL cosD cosT        (for consistency, I changed LHA to T):
y = -cosD sinT

and subsequently convert (x,y) to polar co-ordinates. As one can see, these
are expressions that could be obtained from (III) after multiplication by
cosD. What looks like an unnecessary complication is actually a reduction of
computing steps, because cosD is already available from a
previous altitude computation, whereas tanD is not. But, more importantly, it
prevents the algorithm from faltering for D = 90 degrees.

Such are the subtleties to be considered when proposing a better formula. The
good people at the HMNAO have done their homework. In summa, I believe that
their choice of the cosine formula for sight reduction with a calculator is
the best one for the intended purpose, because of its
generality, simplicity, efficiency, reliability and sufficient accuracy.

Herbert Prinz Browse Files

Drop Files 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)