metcalf@akala.ifa.hawaii.edu (Tom Metcalf) (12/05/90)
I posted the following sight reduction program the other day with
a slight bug. Here is the corrected version. Sorry to waste bandwidth!
---------------------------------------------------------
Following is an HP-48SX program I wrote to do 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.
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) 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.
h) 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.
i) SPEED: Speed of vessel during the observations (knots). If the
speed is zero, the program will terminate at this point.
j) COURSE: True course of vessel during the observations (hms format).
k) 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.
l) 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 (for standard conditions), 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) Accurate dead reckoning information can be included by editing
the "obs" variable and inputing an additional row into the matrix with
the first element equal to the DR longitude (West positive), the second
element equal to the DR latitude (North positive), and the final element
equal to 90.00. The MatrixWriter application is particularly useful
for editing the "obs" variable.
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.
Important note: If observations of more than one body are input, "setup"
must be run before starting the input for each body.
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 DUP 4.4 +
7.31 SWAP / + TAN
.0167 SWAP / SWAP
COS
CASE BODY 'S'
SAME
THEN .0024
* SEMI
END BODY
'M' SAME
THEN HP *
HP .2724 *
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
"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
\>>
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
\>>
\>>
NIT 4
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 / -
\>>
\>>
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
-1.47280528459E-7
\Gm
-1.47296963855E-7
\Gb3 -13.5624809912
\Gb2
-1.50525051351E-2
\Gb1
-3.950284015E-7
E3
[ .188406852706 .980097318647 6.25468131232E-2 ]
E2
[ -1.40179729991E-2 6.63646826456E-2 -.997696960669 ]
E1
[ -.981991016574 .187096170084 2.62424562793E-2 ]
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 .0363326349
\Ga3 17.9636667352
\Ga1 .0000006299
\Gh 'ACOS(R1/Q^1.5)
'
R1 'A0/2+N/3*(A1/
6-Q)'
Q '(N/3)^2-A1/3'
N 18
A0
-4.110603687E-7
A1 .6526786832
G33 .1064412066
R
[ -2.55505296373 -13.2935502826 -.833272135648 ]
G23 1.09880240075
G22 17.255892215
G13 .212196428198
G12 3.31708381131
G11 .637666578414
GHA2
92.6916666667
DEC2
16.4916666667
T2 18
GHA1 77.65
DEC1
16.4916666667
T1 17
LU 1
SEMI
.266666666667
HP .986666666667
HGT 0
INDX 0
BODY T
END
--------------------------CUT HERE--------------------------