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

       
    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