[comp.sys.handhelds] Celestial Navigation w/ HP-48SX

metcalf@akala.IFA.Hawaii.Edu (Tom Metcalf) (12/12/90)

Greetings.

I recently posted a program to do sight reductions for celestial navigation on
the HP-48SX calculator.  Version 2 is now ready for distribution.  

Apart from a few minor changes, the difference between the first post and
this version is that the program will now plot lines of position so that
you can see how good the fix is.

Please mail requests for the program to  metcalf@uhifa.ifa.hawaii.edu

Tom Metcalf
metcalf@uhifa.ifa.hawaii.edu

metcalf@akala.ifa.hawaii.edu (Tom Metcalf) (12/12/90)

Following is version 2 of the HP-48SX program I wrote to do sight
reductions for celestial navigation.  It computes a position fix from
observations of any number of celestial bodies using a least squares
fit to the altitude of the bodies as a function of time.  The routine
does all the standard corrections for dip, refraction, parallax etc. as
well as correcting for motion of the observer between sights.

The major change from version 1 is that the program will now plot the lines
of position, allowing you to *see* how good the fix is.

It does not compute the GHA/declination of celestial bodies, so a copy of
the nautical almanac is required to use these routines.  It will, however,
interpolate the GHA and declination from the hourly entries on the daily pages 
of the nautical almanac.

A detailed description of the mathematical basis of the algorithm is 
available upon request.

----------------------------------------------------------

 Instructions


There are several steps to go through to get a fix from a set of observations.
When prompted for input, key in the requested data, and press "ENTER".
All angular inputs must be in degrees and in the "hms" format: dddd.mmss, 
where dddd is degrees, mm is minutes and ss is seconds of arc.  All times 
must also be input in the "hms" format.
For example, 16.3427 is 16 degrees 34 minutes 27 seconds if an angle is 
being input and 14.0153 is 14:01:53 if a time is being input.

The output is the optimum latitude and longitude, both in "hms" format.
North latitude and west longitude are positive numbers, while south latitude
and east longitude are negative numbers.  For example,
157 deg 49 min 58 sec W, 21 deg 17 min 30 sec N would be output as

 LON:  157.4958
 LAT:   21.1730

1. If you want to start a new set of observations purge the variable
   "obs".  This variable stores all the observations, so, to start over, this
   variable must be removed.  You may want to rename it rather than remove it
   if it will be useful at a later time.

2. Run the program "setup".  This sets up the appropriate corrections
   and the GHA-declination interpolation for the observed body.  This 
   program must be run whenever a new body is observed or whenever the 
   observations have extended beyond the times given for the GHA/declination 
   interpolation (TIM1,TIM2) since the interpolation will become inaccurate.
   If you are observing more than one body, input all the observations for 
   each before proceeding to the next (temporal order does not matter).  
   The "setup" program asks for the following input:
   a) BODY: S is for the Sun, M is for the Moon, VM
      is for Venus or Mars, and anything else assumes a star or other planet.
      Note that alpha mode is automatically initiated.
   b) INDEX:  The index correction (degrees, in hms format) which is to be 
      *added* to the observed sextant altitude,  e.g. 1' should be
      input as 0.0100.
   c) SEMI-D: Semi-diameter in degrees, "hms" format: e.g. 16.2 min
      should be input as 0.1612 since 0.2 min = 12" (Sun only).
   d) HP: Parallax in degrees, "hms" format (Moon,Venus,Mars only).
   e) LIMB: For Upper limb enter "1", for lower limb enter "-1", and
      for disk center enter "0" (Sun/Moon only).
   f) HEIGHT: Height above water (in meters) at which the observations
      were taken (for the dip correction).
   g) PRESSURE TEMPERATURE: The atmospheric pressure in millibars and the 
      atmospheric temperature in Celsius.  These are used for the refraction
      correction: if you want to use standard conditions (usually good enough)
      simply hit ENTER without changing the displayed numbers.
   h) GHA1 DEC1 TIM1: The Greenwich Hour Angle and declination at time TIM1.
      The actual values for the observations will be interpolated linearly from 
      this value and the next.  TIM1 should be a whole hour near the 
      observation times.  All three numbers should be input on the appropriate 
      line *before* pressing "ENTER".  To move to the next item, use the 
      down-arrow key.  All three entries must be in "hms" format.
   i) GHA2 DEC2 TIM2: The second set of values for the linear interpolation.
      TIM2 should be a whole hour after TIM1 and in "hms" format
      (generally, TIM1 and TIM2 should be consecutive hours).  All observations
      must be between TIM1 and TIM2.  If this is not the case, the observations
      should be input in several groups, running "setup" between groups.
      If the GHA passes through zero between GHA1 and GHA2, 360 degrees should
      be added to GHA2.
   j) SPEED: Speed of vessel during the observations (knots). If the
      speed is zero, the program will terminate at this point.
   k) COURSE: True course of vessel during the observations (hms format).
   l) DR LAT LON:  Dead reckoning latitude and longitude to use in the
      correction of the observations for course and speed (hms format).
      Negative values indicate East longitude or South latitude.
   m) TIME OF FIX: The time to which the course and speed corrections 
      are made (hms format).  This will be the time at which the fix is valid.

3.  Enter the observations:
    a) Enter the time of the observation (hms format) and then the 
       uncorrected altitude (hms format) onto the stack.  
    b) Run the "correct" program to correct the observations for index, 
       dip, refraction, parallax, and course/speed.
    c) Run the "addob" program to add the observation to the "obs" variable.
    d) Repeat until all observations are input.
    e) If accurate dead reckoning information is available, it
       can be included in the fix by running the ADDDR program which includes
       the DR position in the "obs" variable.  Dead reckoning information is 
       only required when only two observations are available; otherwise, it 
       is optional.  The input should be in "hms" format.  Negative values 
       indicate East longitude or South latitude.  If the vessel is moving, 
       the dead reckoning position should be computed at the sime time used 
       in step 2(m) above.

Important note:  If observations of more than one body are input, "setup"
must be run before starting the input for each body.


4.  Get the fix by running the "solve" program.  This program can be run at 
    any time when there are at least 3 observations (including dead reckoning) 
    in "obs".  It does not affect "obs", so more observations can be input 
    after running "solve".  If "convergence error" appears or if the position 
    estimate is far from your dead reckoning position, there is probably an 
    error in the input data and it should be reentered (if the data is correct 
    and "convergence error" appears,  the position fix should not be trusted).
    Remember: the fix is, at best, only as good as the data you supply, and
    you should examine the results critically!

5.  If all observations are for a *single celestial body*, run the "error" 
    program to get an estimate of the position error (miles).  This program 
    assumes an error on the observations of one arc minute, and should be 
    multiplied by the actual sextant error in arc minutes if it is other than 
    one.

6.  Run the "plotp" program to plot the lines of position.  To run run this
    program the position fix from "solve" must be the first two numbers on
    the stack (longitude in level 2, latitude in level 1).  The "plotp" 
    program will then ask you for the scale of the plot in miles.  The
    program makes the scale of the plot as close as possible to this size.
    Hence, if you input 100 miles, the HP-48SX screen will display a region
    of the Earth's surface 100 miles on a side (latitude on the vertical 
    scale (North up), longitude on the horizontal scale (East right)).  The 
    default is 10 miles; if this is what you want, simply press return at 
    the prompt.  The program will work for very large scales (try 10000), but 
    the Earth's spherical surface is mapped onto a flat lat/lon grid and hence 
    the lines of position will appear somewhat distorted.  A cross-hair is 
    drawn at the position of the fix with line segments one mile long.  The 
    cross hair can be used to judge the scale of the plot but for large plots 
    it will disappear as the scale gets too large to resolve a one mile line 
    segment.  If no lines of position appear, you need to expand the plot by 
    using a large scale value.


Disclaimer:


This software is provided "as is" and is subject to change without
notice.  No warranty of any kind is made with regard to this software,
including, but not limited to, the implied warranties of merchantability 
and fitness for a particular purpose.  The author shall not be liable for 
any errors or for incidental or consequential damages in connection with 
the furnishing, performance, or use of this software.

Copyright 1990 by Thomas R. Metcalf.  Permission is granted to any
individual or institution to use, copy, or redistribute this software
so long as it is not sold for profit and provided that this copyright
notice and the above disclaimer are retained.

--------------------------CUT HERE--------------------------
%%HP: T(3)A(D)F(.);
DIR
  SOLVE
    \<< -22 SF 4 FIX
DEG 0 0 0 0 0 GSUM
a0 \->NUM 'A0' STO a1
\->NUM 'A1' STO EV1
\->NUM DUP '\Ga1' STO
EIGEN 'E1' STO EV3
\->NUM DUP '\Ga3' STO
EIGEN 'E3' STO EV2
\->NUM DUP '\Ga2' STO
EIGEN 'E2' STO R E1
DOT '\Gb1' STO R E2
DOT '\Gb2' STO R E3
DOT '\Gb3' STO 0
'NIT' STO 0 '\Gm' STO
      DO \Gm 'OLD'
STO ITER '\Gm' STO 1
'NIT' STO+
      UNTIL 'ABS((\Gm
-OLD)/\Gm)<.000001 OR
NIT>50'
      END
      IF 'NIT>50 OR
\Gm>\Ga1'
      THEN
"CONVERGENCE ERROR"
      END UVW OBJ\->
DROP OUT
    \>>
  ADDOB
    \<< \-> T A
      \<< T HMS\-> 'T'
STO A HMS\-> 'A' STO
OBS
        IFERR OBJ\->
        THEN T GHA1
GHA2 INTERP T DEC1
DEC2 INTERP A { 1 3
} \->ARRY SWAP STO
        ELSE OBJ\->
ROT 1 + ROT ROT
\->LIST T GHA1 GHA2
INTERP SWAP T DEC1
DEC2 INTERP SWAP A
SWAP \->ARRY 'OBS'
STO
        END
      \>>
    \>>
  CORRECT
    \<< DEG HMS\-> INDX
+ HGT \v/ .0293 * -
DUP DUP REFRACT
SWAP COS
      CASE BODY 'S'
SAME
        THEN
.002443 * SEMI
        END BODY
'M' SAME
        THEN HP *
HP .272476 *
        END BODY
'VM' SAME
        THEN HP * 0
        END 0 * 0
      END LU * +
SWAP - +
      IF 'SPD>0'
      THEN SWAP
HMS\-> DUP DUP GHA1
GHA2 INTERP SWAP
DEC1 DEC2 INTERP
SWAP DRLAT DRLON
AZIM DUP CSCORR ROT
SWAP - SWAP \->HMS
SWAP
      END \->HMS
    \>>
  SETUP
    \<< "BODY?" { ""
\Ga V } INPUT OBJ\->
'BODY' STO
"INDEX? (Deg)" { ""
V } INPUT OBJ\-> HMS\->
'INDX' STO
      IF BODY 'S'
SAME
      THEN
"SEMI-D? (Deg)" {
"" V } INPUT OBJ\->
HMS\-> 'SEMI' STO
      END
      IF BODY 'M'
SAME BODY 'VM' SAME
OR
      THEN
"HP? (Deg)" { "" V
} INPUT OBJ\-> HMS\->
'HP' STO
      END
      IF BODY 'M'
SAME BODY 'S' SAME
OR
      THEN
"LIMB (L/U/C=1/-1/0)?"
{ "" V } INPUT OBJ\->
'LU' STO
      END
"HEIGHT (m)?" { ""
V } INPUT OBJ\->
'HGT' STO
"ENTER for std cond"
{
":PRESS (mb): 1010
:TEMPER (C): 10"
-14 V } INPUT OBJ\->
'TMPTR' STO 'PRESS'
STO
"GHA1 DEC1 TIM1?" {
":GHA1:
:DEC1:
:TIM1:"
{ 1 0 } V } INPUT
OBJ\-> HMS\-> 'T1' STO
HMS\-> 'DEC1' STO
HMS\-> 'GHA1' STO
"GHA2 DEC2 TIM2" {
":GHA2:
:DEC2:
:TIM2:"
{ 1 0 } V } INPUT
OBJ\-> HMS\-> 'T2' STO
HMS\-> 'DEC2' STO
HMS\-> 'GHA2' STO
"SPEED? (Knots)" {
"" V } INPUT OBJ\->
'SPD' STO
      IF 'SPD\=/0'
      THEN
"COURSE? (True)" {
"" V } INPUT OBJ\->
HMS\-> 'CRS' STO
"DR LAT LON?" {
":LAT:
:LON:" { 1 0
} V } INPUT OBJ\->
HMS\-> 'DRLON' STO
HMS\-> 'DRLAT' STO
"TIME OF FIX?" { ""
V } INPUT OBJ\-> HMS\->
'TF' STO
      ELSE 0 'CRS'
STO 0 'DRLAT' STO 0
'DRLON' STO 0 'TF'
STO
      END
    \>>
  ADDDR
    \<< 0 \-> n
      \<< OBS OBJ\->
OBJ\-> DROP DROP 'n'
STO "dd.mmss" {
":DR_LAT: 
:DR_LON: "
{ 1 0 } V } INPUT
OBJ\-> HMS\-> SWAP HMS\->
90 n 1 + 3 2 \->LIST
\->ARRY 'OBS' STO
      \>>
    \>>
  ERROR
    \<< 0 0 0 0 0 0 0
0 \-> H1 H2 D1 D2 G1
G2 DT DH
      \<< OBS { 1 3 }
GET 'H1' STO OBS {
N 3 } GET 'H2' STO
OBS { 1 2 } GET
'D1' STO OBS { N 2
} GET 'D2' STO OBS
{ 1 1 } GET 'G1'
STO OBS { N 1 } GET
'G2' STO T2 T1 -
GHA2 GHA1 - / G2 G1
- * 'DT' STO H2 H1
- 'DH' STO 1 DT / N
\v/ / 57.3 H1 H2 + 2
/ COS * * 225 D1 D2
+ 2 / COS SQ * DH
DT / SQ - \v/ / "ERR"
\->TAG
      \>>
    \>>
  PLOTP
    \<< 2 DUPN HMS\->
'LAT' STO HMS\->
'LON' STO 0 0 0 0 0
0 0 0 0 0 0 0 0 \-> g
d a l n N sc sc\Gl
ssz d0 d1 ll lm
      \<<
"Scale? (Miles)" {
"10" -1 V } INPUT
OBJ\-> ABS 120 / DUP
'sc' STO LAT COS /
2.0469 * 180 MIN
NEG 'sc\Gl' STO ERASE
{ # 0h # 0h } PVIEW
LON sc\Gl + RANGE LAT
sc + 90 MIN DUP 3
ROLLD R\->C PMAX LON
sc\Gl - RANGE LAT sc
- -90 MAX DUP 3
ROLLD R\->C PMIN - 2
/ 'sc' STO OBS OBJ\->
OBJ\-> DROP2 DUP 'N'
STO 3 * DROPN 1 N
        FOR n DEPTH
'd0' STO OBS { n 1
} GET 'g' STO OBS {
n 2 } GET 'd' STO
OBS { n 3 } GET 'a'
STO
          IF 'LAT-
sc>d+90-a OR LAT+sc
<d-90+a'
          THEN
          ELSE LAT
sc + d 90 a - +
            IF DUP
90 >
            THEN
180 SWAP -
            END MIN
LAT sc - d 90 a - -
            IF DUP
-90 <
            THEN
180 + NEG
            END MAX
            IF LAT
d <
            THEN
SWAP
            END
DUP2 SWAP - DUP
SIGN
            IF DUP
0 ==
            THEN
DROP 1
            END
SWAP ABS 90 a -
PSCALE sc 32 / MAX
* 'ssz' STO DUP
'lm' STO SWAP DUP
'll' STO - ssz /
CEIL 0 SWAP
            FOR l g
d a l ssz * ll +
DUP lm
              IF '
ssz<0'
              THEN
SWAP
              END
              IF >
              THEN
DROP lm
              END
LOP DUP C\->R SWAP g
- NEG g + RANGE
SWAP R\->C DEPTH d0 -
ROLLD
            NEXT
DEPTH d0 - 2 / 2 +
'd1' STO
            WHILE
DEPTH d0 - DUP 1 >
            REPEAT
              IF d1
\=/
              THEN
OVER SWAP
              END
LIMIT LINE
            END
DEPTH d0 - DROPN
          END
        NEXT LAT
COS DUP LON .00833
ROT / - LAT R\->C
SWAP LON .00833 ROT
/ + LAT R\->C LINE
LON LAT .00833 -
R\->C LON LAT .00833
+ R\->C LINE
      \>> { } PVIEW
    \>>
  PSCALE
    \<< \-> s a
      \<<
        IF 's\=/0'
        THEN 'a/(30
+a/s)' \->NUM
        ELSE 0
        END
      \>>
    \>>
  LON 157.833138922
  LAT 21.2991627064
  LIMIT
    \<< 0 0 0 0 0 0 \->
g1 g2 d1 d2 d180 up
      \<< DUP2 C\->R
'd1' STO 'g1' STO
C\->R 'd2' STO 'g2'
STO
        IF 'ABS(g1-
g2)>180'
        THEN DROP2
LON 180
          IF 'g1>
LON'
          THEN +
          ELSE -
          END 'up'
STO 'd1+(up-g1)*(d1
-d2)/(g1-g2)' \->NUM
'd180' STO g2 d2
R\->C up 360
          IF 'up>
LON'
          THEN -
          ELSE +
          END d180
R\->C up d180 R\->C g1
d1 R\->C LINE
        END
      \>>
    \>>
  RANGE
    \<<
      WHILE DUP 180
LON + >
      REPEAT 360 -
      END
      WHILE DUP
-180 LON + <
      REPEAT 360 +
      END
    \>>
  NIT 3
  PPAR {
(158.01621892,21.2158293731)
(157.650058924,21.3824960397)
X 0 (0,0) FUNCTION
Y }
  LOP
    \<< \-> g d a l
      \<<
        IF 'ABS(l)\=/
90'
        THEN 'g+
ACOS((SIN(a)-SIN(l)
*SIN(d))/(COS(l)*
COS(d)))' \->NUM
        ELSE g
        END DUP IM
        IF 0 \=/
        THEN DROP g
        END
        IF 'ABS(l)>
90-ABS(d)+a'
        THEN 180 +
        END RANGE l
R\->C
      \>>
    \>>
  ITER
    \<< 0 0 \-> f fp
      \<< \Gb1 \Ga1 \Gm - /
SQ DUP 'f' STO+ 2 *
\Ga1 \Gm - / 'fp' STO+
\Gb2 \Ga2 \Gm - / SQ DUP
'f' STO+ 2 * \Ga2 \Gm -
/ 'fp' STO+ \Gb3 \Ga3 \Gm
- / SQ DUP 'f' STO+
2 * \Ga3 \Gm - / 'fp'
STO+ -1 'f' STO+ \Gm
f fp / -
      \>>
    \>>
  CST { OBS SOLVE
ADDOB CORRECT SETUP
ERROR TIME \->HMS
HMS\-> HMS+ HMS- }
  REFRACT
    \<< 0 \-> h rp
      \<< '1/TAN(h+
7.31/(h+4.4))' \->NUM
'rp' STO 'rp*((
PRESS-80)/930)/(1+
.00008*(rp+39)*(
TMPTR-10))' \->NUM 60
/
      \>>
    \>>
  PRESS 1010
  TMPTR 10
  a0 '-(G12*G23-G13
*G22)*G13+(G11*G23-
G12*G13)*G23-(G11*
G22-G12^2)*G33'
  a1 'G11*G22-G12^2
+G11*G33-G13^2+G22*
G33-G23^2'
  TF 0
  DRLON 0
  DRLAT 0
  CRS 0
  SPD 0
  CSCORR
    \<< \-> T
      \<< SPD T TF -
AZ CRS - COS 60 / *
*
      \>>
    \>>
  AZ 239.148905272
  AZIM
    \<< \-> D G L A
      \<< G A - 'A'
STO L COS D SIN * L
SIN D COS A COS * *
- A SIN D COS NEG *
R\->C ARG 'AZ' STO
        IF 'AZ<0'
        THEN 360
'AZ' STO+
        END
      \>>
    \>>
  EV3 '-2*\v/Q*COS((\Gh
+360)/3)+N/3'
  EV2 'N-\Ga1-\Ga3'
  EV1 '-2*\v/Q*COS(\Gh/
3)+N/3'
  OLD
-7.20287871696E-9
  \Gm
-7.20288321614E-9
  \Gb3 1.37838537944
  \Gb2 .575712622485
  \Gb1 -3.32596238E-7
  E3
[ .373606532854 -.192391048226 .90741602541 ]
  E2
[ -.147739742458 -.978108259972 -.146551015926 ]
  E1
[ -.915746213255 7.93089929862E-2 .393851439684 ]
  INTERP
    \<< \-> T V1 V2
      \<< V1 V2 V1 -
T2 T1 - / T T1 - *
+
      \>>
    \>>
  GSUM
    \<< \-> DS DC GS GC
HS
      \<< 0 'G11' STO
0 'G12' STO 0 'G13'
STO 0 'G22' STO 0
'G23' STO { 3 } 0
CON 'R' STO OBS
OBJ\-> OBJ\-> DROP DROP
'N' STO 1 N
        START SIN
'HS' STO DUP SIN
'DS' STO COS 'DC'
STO DUP SIN 'GS'
STO COS 'GC' STO DS
SQ 'G11' STO+ DS DC
GC * * 'G12' STO+
DS DC GS * * 'G13'
STO+ DC SQ GC SQ *
'G22' STO+ DC SQ GS
GC * * 'G23' STO+ R
OBJ\-> DROP DC GS HS
* * + ROT DS HS * +
ROT DC GC HS * * +
ROT { 3 } \->ARRY 'R'
STO
        NEXT N G11
G22 + - 'G33' STO
      \>>
    \>>
  OUT
    \<< \-> U V W
      \<<
        IF 'ABS(U)>
1'
        THEN U SIGN
'U' STO
        END U ASIN
V W R\->C ARG \->HMS
"LON" \->TAG SWAP
\->HMS "LAT" \->TAG
      \>>
    \>>
  UVW
    \<< \Gb1 \Ga1 \Gm - /
E1 * \Gb2 \Ga2 \Gm - / E2
* \Gb3 \Ga3 \Gm - / E3 *
+ +
    \>>
  EIGEN
    \<< \-> EV
      \<< 'G12*G23-
G13*G22+G13*EV'
\->NUM 'G13*G12-G11*
G23+G23*EV' \->NUM '
G11*G22-SQ(G12)-(
G11+G22)*EV+SQ(EV)'
\->NUM { 3 } \->ARRY
DUP ABS /
      \>>
    \>>
  \Ga2 .77929306844
  \Ga3 2.22070567227
  \Ga1 .000001259291
  \Gh 'ACOS(R1/Q^1.5)
'
  R1 'A0/2+N/3*(A1/
6-Q)'
  Q '(N/3)^2-A1/3'
  N 3
  A0 -.000002179305
  A1 1.73058431531
  G33 1.84527481377
  R
[ .429918452524 -.828298305696 1.16639758178 ]
  G23
-.275981896423
  G22 .827744288915
  G13 .769728325982
  G12
-.047009095401
  G11 .326980897319
  GHA2 1
  DEC2 1
  T2 1
  GHA1 1
  DEC1 1
  T1 1
  LU 0
  SEMI 1
  HP .986666666667
  HGT 1
  INDX 1
  BODY S
END
--------------------------CUT HERE--------------------------

metcalf@galileo.ifa.hawaii.edu (Tom Metcalf) (12/13/90)

There is an error in the instructions for using the Revision 2 celestial 
navigation program.  In section 2(e) on inputing which limb of the sun or 
moon was used, the upper and lower limb specifications are reversed.  The 
prompt in the program is correct, however.

Tom Metcalf
metcalf@uhifa.ifa.hawaii.edu