NavList:

   

Reply

"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

   

Reply