Following is version 4.3 of the HP48SX 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 3 is a much improved user interface (including a
choice of input/output formats) and some route planning programs (great circle,
rhumb line, etc). A bug which sometimes crashed the ADV program when the
speed from the INIT menu was set to zero has been corrected in version 4.2.
Version 4.3 includes better error trapping and a slight modification to
the PLOTP program to avoid the requirement that the fix be kept on the
stack. Also included in 4.3 is a new error estimate for all fixes.
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 (send a surface mail address).
Questions, comments, and suggestions are welcome
Tom Metcalf
metcalf@uhifa.ifa.hawaii.edu

Instructions
There are a number of programs included in the directory supplied, however,
there are even more variables etc used by the programs and these should not be
tampered with. To avoid confusion, press the custom menu key (CST) on the
calculator and only the programs (without all of the variables) will be in the
menu.
Note: these instructions are not a "tutorial" on celestial navigation.
A basic understanding of celestial navigation is assumed.
There are several steps to go through to get a fix from a set of observations.
1. First you need to set up several parameters used throughout the computation
of the fix. To do this, run the "INIT" program. You should examine the
displayed information and change anything which is incorrect. Select the
appropriate menu key to change the data:
Format: Selects HMS, HMT, or decimal format for angular input/output.
In HMS format angles are input in degreesminutesseconds format:
dd.mmss. In HMT format, angles are input/output in
degreesminutestenths format: dd.mmt. For example 45 degrees
16 minutes 12 seconds would be input/output as 45.1612 in HMS
format and 45.162 in HMT format. Note that times are
always input/output in HMS format regardless of the setting
in the INIT menu. For example, a time of 16:22:15 would
always be input as 16.2215. When you select HMT or HMS format,
all angular input must be in the selected format throughout the
sight reduction process.
Index: The index correction which is to be *added* to the observed
altitude, to be input in the current format (HMT or HMS).
For example, +1' should be input as 0.0100.
Height: The height of the observer's eye for use in the dip correction.
The default units are meters, however, any valid unit can be
specified. For example, to input a height of 10 feet, use
10_ft (use the orangeshift UNITS menu to avoid typing in the
"ft"). Note that the "ft" must be in lower case. If no
units are specified, or if the units are invalid, meters are
assumed.
C/S Course/Speed. This specifies the course steered and speed of
the vessel during the observations, used to correct for the
motion of the vessel. The course must be a true course in
degrees (HMS or HMT format). The default units for speed are
knots, but any valid unit of velocity may be used. Use the
arrow keys and delete key to specify both course and speed before
pressing "ENTER".
P/T Pressure/Temperature. The atmospheric pressure and temperature
used for the refraction correction: if you want to use standard
conditions (usually good enough) simply hit ENTER without
changing the displayed numbers. The default units for pressure
are millibars and the default units for temperature are Celsius.
Any valid units can be specified, however. Use the arrow keys
and delete key to specify both pressure and temperature before
pressing "ENTER".
EXIT Press EXIT when you're finished with the INIT menu.
Remember: all times must be input in the "hms" format!
2. 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. "OBS" is a matrix with each row
representing one observation in the format: GHA DEC ALT. It cannot be
accessed from the custom menu but, if present, it should be the first
variable in the VAR menu. Use the MatrixWriter application in the 48SX to
edit the "OBS" variable. In particular, if you want to delete an
observation, use the "ROW" key in the MatrixWriter application.
3. Run the program "SETUP". This sets up the appropriate corrections
and the GHAdeclination 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; see below) 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: Select the appropriate body by pressing a menu key. "planet"
means any planet other than Venus or Mars (i.e. Jupiter or Saturn).
edit the displayed number with the arrow and delete keys.
b) SEMID: Semidiameter in degrees, HMS or HMT format: e.g. 16.2 arc minutes
should be input as 0.1612 in HMS format (since 0.2 min = 12") and 0.162
in HMT format. This is requested for the Sun only. The default
value is computed for the current date in the calculator. If the
observations are not for the date currently in the calculator, the SD
may need to be edited before you press ENTER.
c) HP: Parallax in degrees, HMS or HMT format. Moon,Venus,Mars only.
For the Moon, HP is found on the daily pages of the Nautical Almanac.
For Venus/Mars, it is found in a short "parallax" table in the
explanation section of the Nautical Almanac. For the moon a typical
value might be 54.4 arc minutes which would be input as 0.5424 in
HMS format and 0.544 in HMT format. The correction is much smaller for
the planets: a typical value might be 0.3 arc minutes which would be
input as 0.0018 in HMS format (since 0.3 arc minutes is 18 arc seconds)
and 0.003 in HMT format.
d) LIMB: Select lower limb (LL), upper limb (UL), or disk center (CENT)
from the menu keys. (Sun/Moon only).
e) If the body is a *star* the program will ask for GHAARIES and
SHA, DEC for the star at time TIM (from the daily pages of the Nautical
Almanac). These should be at the whole hour (TIM) before the start of
the observations. GHA,DEC,SHA inputs must be in HMS or HMT format and the
time must be in HMS format. All must be input on the appropriate line
*before* pressing "ENTER". To move to the next item, use the downarrow
key. Since the stars move at a very regular rate, no interpolation
is necessary. Skip to (g).
e') 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 the whole hour before the
observation times. All three numbers should be input on the appropriate
line *before* pressing "ENTER". To move to the next item, use the
downarrow key. GHA and DEC must be in HMS or HMT format depending on
the selection in the INIT program (step 1), and the time must be in HMS
format.
f) 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 will be consecutive hours, however, this need not be the
case if you are willing to forego some accuracy). 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. Similarly, the time should not pass through zero:
if TIM1 is 23, TIM2 should be 24 not 0. As an error check, the program
will ask you to reinput the data if the time and GHA do not increase.
g) DR LAT LON: Dead reckoning latitude and longitude at the time of the
observations (not necessarily the time of the fix) to use in the
correction of the observations for course and speed (HMS or HMT format).
Negative values indicate East longitude or South latitude. The default
values displayed are the last known DR position. If these are correct,
simply hit ENTER. If they are incorrect, they can be edited with the
arrow and delete keys. These values need not be precise; use the best
information you have.
h) 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
and should be the same for all observations. Remember: time should not
pass through zero  add 24 hours if necessary.
4. Enter the observations:
a) Run the "ADDOB" program to add an observation to the "OBS" variable.
Enter the time of the observation (HMS format) and then the *uncorrected*
altitude. Care should be taken with the time: Since it is used in
the interpolation of the GHA and declination of the body, it must lie
within the range of times specified in steps 2e' and 2f. It should not
pass though 00:00. If the time is not within the range specified in the
"SETUP" program for the interpolation, the program will terminate
without adding the observation to the "OBS" variable. "ADDOB" does all
corrections (index, dip, refraction, parallax, semidiameter,
course/speed). The course/speed correction adjusts the GHA and
declination but not the altitude of the observation.
b) Repeat until all observations are input.
c) 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 just 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 same time used
in step 2(h) above.
Important note: If observations of more than one body are input, "SETUP"
must be run before starting the input for each body. Use the same time
in step 2(h) for all bodies, unless the LOP's are advanced appropriately (see
step 8).
5. 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 "root error" appears (theoretically impossible)
or if the position estimate is far from your dead reckoning position, there
may be an error in the input data and it should be reentered. If the data
is correct, but the error persists, the position fix should not be trusted;
note that the "ERROR" program (step 6) or the plot of the LOPs (step 7) can
be used to determine the validity of the fix. Remember: the fix is, at
best, only as good as the data you supply, and you should examine the
results critically!
"SOLVE" can update the DR position to keep a running fix; after the fix
is computed "SOLVE" will display the fix and ask if you want the DR
position updated. The "DR" program can also be used to update the DR
position "by hand".
The output is the optimum latitude and longitude, both in the current
format. North latitude and west longitude are positive numbers, while south
latitude and east longitude are negative numbers. For example,
157 deg 49 min 54 sec W, 21 deg 17 min 30 sec N would be output as
HMS: LON: 157.4954 HMT: LON: 157.499 Decimal: LON: 157.8317
LAT: 21.1730 LAT: 21.175 LAT: 21.2917
6. Run the "ERROR" program to get an estimate of the *random* position error
(miles). Systematic errors are not treated. This program gives the
formal errors on the least squares fit to the latitude/longitude and hence
gives a separate East/West and North/South error. The sextant error is
taken to be the difference between the computed and observed altitudes at
the position of the fix. Note that this difference is very small when only
two observations are specified since, by definition, both LOPs intersect
at the fix. Hence, this program gives meaningful results ONLY if more than
two observations (*not* including dead reckoning) are given. The computed
error is the standard "one sigma" error and hence you may want to multiply
by 2 or 3 in critical situations: assuming random errors with a normal
distribution, there is a 68.3% probability that the fix is within the
computed error range, a 95.4% probability that the fix is within twice the
computed range, and a 99.7% probability that the fix is within three times
the computed range. The "PLOTP" program (step 7) will give a visual, and
hence perhaps a better, indication of the reliability of the fix.
7. Run the "PLOTP" program to plot the lines of position.
The "PLOTP" program will ask you for the position of the center of the plot
and the scale of the plot in miles. The default position for the center
is the last fix from "SOLVE". If this is what you want simply press ENTER.
The program makes the scale of the plot as close as possible to the
specified size. Hence, if you input 100 miles, the HP48SX screen will
display a region of the Earth's surface about 100 miles on a side (latitude
on the vertical scale (North up), longitude on the horizontal scale (East
right)). The default is 9 miles; if this is what you want, simply press
ENTER 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 (really circles) of position will appear
somewhat distorted. The larger the scale the longer it takes to plot the
LOP's. A crosshair is drawn at the center position (the position
of the fix, if you used the default value) with line segments
one mile long (from tip to tip). The cross hair can be used to judge the
scale of the plot but for large plots it will disappear as it becomes
impossible to resolve a one mile line segment. If no lines of position
appear, you need to expand the plot by using a larger scale value. Press
"ATTN" when you are done with the plot. Note that the "coordinate" graph
utility and the arrow keys built into the calculator can be used to examine
the lat/lon coordinates (always decimal format) of the plot after pressing
.
"PLOTP" uses quite a bit of stack memory. If there is not enough
available, the program will terminate prematurely, though gracefully.
8. You can advance the LOP's contained in the "OBS" variable using the "ADV"
program. "ADV" assumes that the DR position is correct for the time implied
by the "OBS" data: if it is not, the DR position should be updated with the
"DR" program. "ADV" asks for a distance in nautical miles and a true course
describing the vessel's motion. All data in the "OBS" variable, as well as
the DR position, will be advanced along this track. After "ADV" is run,
another group of data may be added to "OBS" to update the running fix.
Note that the new data must be input with a "time of fix" (step 2h)
consistent with the *advanced* LOP's. For example, if you advance the LOP's
by the distance traveled in 6 hours, the time of fix in step 2(h) should
reflect this by having 6 hours added to the value used previously (so that
they all apply at the same time). Advancing lines of position by large
distances is not, in general, a good idea, since errors will be incurred
as your dead reckoning position becomes less certain. The "ADV" program
assumes the observer is moving along a constant course; if the course
varies, "ADV" must be run separately for each course steered.
If no "OBS" variable is present, "ADV" will update only the DR position
and hence is a useful tool for keeping track of your DR position, even when
you are not keeping a running fix. For completeness, "ADV" includes a
slight correction for the ellipticity of the Earth.
If the effects of a current need to be included, run "ADV" once for the
motion of the vessel through the water and once for the direction and
distance the current has moved the vessel since the last DR position.
If you make a mistake entering the data for "ADV", you can undo the damage
by running "ADV" again with the same course but the negative of the
distance.
Route Planning:
The program "SAIL" calculates the course and distance from one position
to another along a rhumb line or a great circle. Note that no correction
is made for the ellipticity of the Earth in this program. To use it run SAIL
and input the starting and ending lat/lon in HMS or HMT format. The default
starting position is the DR position and the default ending position is the
last known ending position. When asked for the type of course press the
appropriate menu key:
RHUMB: Returns the rhumb line course and distance
GC: Returns the great circle distance and initial course. Note that
the course along a great circle varies and the *initial* course
is given.
WAY: Computes waypoints along a great circle. You are asked for a "scale"
which defines the length of each computed rhumb line segment along
the great circle. Any valid units can be input: if no units are
specified or if the units are invalid, nautical miles are assumed.
The output is a list of waypoints. Each waypoint is itself a list in
the format {lat lon course} where "lat/lon" are the latitude and
longitude of the waypoint (HMS or HMT format) and "course" is the
rhumb line course to the next waypoint (HMS or HMT format). The
list of waypoints may be long and there are several ways to
examine it: (1) Use the editor or (2) execute OBJ> to break the list
into individual waypoints or (3) use the WVIEW program. Note that
the WVIEW program is the best way to view the data, but I wasn't
up to writing a data browser so I use Donnelly's DBR program
contained in his "Programmer's Toolkit". If you don't have the
"Programmer's Toolkit" installed, you won't be able to use WVIEW.
Also output is the total distance along the track and the additional
distance traveled due to approximating the great circle by a sequence
of rhumb lines (labeled ADDD).
VERT: Computes one vertex of the great circle connecting the starting and
ending points (point of highest latitude).
COMP: Computes a composite route which follows a great circle to a
maximum latitude and follows that parallel until another great
circle can be followed to the destination. The output is a list
of waypoints as in the WAY program. Total distance along the
track is also output. The program requests the limiting latitude
and the scale which gives the distance between waypoints along the
great circle routes. If the great circle from the starting
to the ending point never crosses the limiting latitude,
the message "GC is OK" is displayed. If the starting or ending
points are beyond the limiting latitude, the message "No sol'n"
is displayed.
EXIT: Exit the SAIL program when all desired computations have been
done.
Helpful Hints:
o Remember that South latitude, South declination and East longitude must be
negative numbers.
o If you need to compute a number for input into one of the programs (eg the
distance traveled in the ADV program) there are several ways to proceed:
a) Input an algebraic object, e.g. '(2+1/60)*6' (quotes required)
b) Compute before running the program and leave the number on the
stack. When the number is required use ^STK from the edit menu
to get the number off the stack (^STK displayed when the calculator
prompts you for info).
o Since the SOLVE program uses a least squares fit to the altitude of one or
more bodies, it is possible to get a fix from a relatively large number of
observations of a single body over a relatively short period of time (say
20 observations in 20 minutes).
o You should be able to exit all the programs gracefully by pressing the
"ATTN" (ON) key when the program is running.
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 direct, indirect, special, incidental or consequential
damages in connection with the furnishing, performance, or use of this software:
use it at your own risk.
Copyright 1991 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 AND AT END 
%%HP: T(3)A(D)F(.);
DIR
SOLVE
\<< DEPTH SAVES \>
dsv
\<<
IFERR FFIX
DEG 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
IF '\Ga1==0
AND \Gb1==0'
THEN
"AMBIGUOUS SOLUTION"
MESS 0 DOERR
END R E2
DOT '\Gb2' STO R E3
DOT '\Gb3' STO 'G\Gm'
'\Gm' { \GmST LBND UBND
} ROOT DROP
IF '\Gm>\Ga1
OR \GmSTR 4 DISP
\>STR 5 DISP ASK
IF 11.1
==
THEN DUP2
FMT\> 'DRLAT' STO
FMT\> 'DRLON' STO
END DUP2
FMT\> 'LAT' STO FMT\>
'LON' STO
THEN DEPTH
dsv  DROPN ABORT
END
\>> RESTS
\>>
ADDOB
\<< DEPTH SAVES
DEG RCLMENU 28 MENU
\> dsv om
\<<
IFERR
"Time/Altitude
(hh.mmss)/"
FMT +
":Time:
:H_s: " {
1 0 } 'V' 3 \>LIST
INPUT OBJ\> DTAG
SWAP DTAG SWAP 0 \>
TM A n
\<< TM HMS\>
'TM' STO
IF TM
T1 < TM T2 > BODY
"T" SAME NOT AND OR
THEN
"Error:Bad Time
Press ENTER"
MESS 0 DOERR
END A
CORRECT FMT\> 'A'
STO TM GHA1 GHA2
INTERP 180 RANGE TM
DEC1 DEC2 INTERP
IF 'SPD
\=/0'
THEN TF
TM  SPD * 60 / CRS
RMOVE SWAP 180
RANGE SWAP
END OBS
IFERR
OBJ\>
THEN 3
ROLLD A { 1 3 }
\>ARRY SWAP STO
ELSE
OBJ\> ROT 1 + DUP 3
* 'n' STO ROT ROT
\>LIST n ROLL n ROLL
ROT A SWAP \>ARRY
'OBS' STO
END
\>>
THEN DEPTH
dsv  DROPN ABORT
END om MENU
\>> RESTS
\>>
SETUP
\<< DEPTH RCLMENU
28 MENU \> dsv om
\<<
IFERR FFIX
CLLCD 2 FREEZE
MBODY TMENU "BODY?"
PROMPT 'BODY' STO 0
MENU
IF BODY
"S" SAME
THEN
DO
"SEMID? " FMT + SD
\>FMT \>STR 'V' 2
\>LIST INPUT OBJ\>
FMT\> 'SEMI' STO
IF '
SEMI>.55'
THEN
"TOO LARGE:PRESS ENTER"
MESS
END
UNTIL '
SEMI\<=.55'
END
END
IF BODY
"M" SAME BODY "VM"
SAME OR
THEN
DO
"HParallax? " FMT +
HP \>FMT \>STR 'V' 2
\>LIST INPUT OBJ\>
FMT\> 'HP' STO
IF '
HP>1.2'
THEN
"TOO LARGE:PRESS ENTER"
MESS
END
UNTIL '
HP<1.2'
END
END
IF BODY
"M" SAME BODY "S"
SAME OR
THEN
CLLCD 2 FREEZE
MLIMB TMENU "Limb?"
PROMPT 'LU' STO 0
MENU
END
DO
IF BODY
"T" SAME
THEN
"Star" ":GHA\Gg: " G\Gg
\>FMT \>STR +
"
:SHA:
:DEC:
" +
":TIM: " T\Gg \>HMS
\>STR + + { 1 0 }
'V' 3 \>LIST INPUT
OBJ\> HMS\> DUP 'T1'
STO DUP 'T\Gg' STO 1
+ 'T2' STO FMT\> DUP
'DEC1' STO 'DEC2'
STO FMT\> SWAP FMT\>
DUP 'G\Gg' STO + DUP
'GHA1' STO
15.041067 + 'GHA2'
STO
ELSE
"Linear Interp 1" {
":GHA1:
:DEC1:
:TIM1: "
{ 1 0 } V } INPUT
OBJ\> HMS\> 'T1' STO
FMT\> 'DEC1' STO
FMT\> 'GHA1' STO
"Linear Interp 2" {
":GHA2:
:DEC2:
:TIM2: "
{ 1 0 } V } INPUT
OBJ\> HMS\> 'T2' STO
FMT\> 'DEC2' STO
FMT\> 'GHA2' STO
END
IF 'T1\>=
T2'
THEN
"Err:T1\>=T2:Press ENTER"
MESS
END
IF '
GHA1>GHA2'
THEN
"GHA1>GHA2:Hit ENTER"
MESS
END
UNTIL 'T1
HMS \>STR 'V' 2
\>LIST INPUT OBJ\>
HMS\> 'TF' STO FFIX
END
THEN DEPTH
dsv  DROPN ABORT
END om MENU
\>>
\>>
INIT
\<< RCLMENU 28
MENU \> om
\<< FFIX { {
"INDEX" {
\<< 0 MENU
"INDEX? " FMT +
INDX \>FMT "INDEX"
\>TAG \>STR { 1 0 }
'V' 3 \>LIST INPUT
OBJ\> FMT\> 'INDX'
STO 0 CONT
\>> } } {
"HEIGHT" {
\<< 0 MENU
"HEIGHT? (m)" HGT
"HGT" \>TAG \>STR { 1
0 } 'V' 3 \>LIST
INPUT OBJ\> '1_m'
DOUNIT 'HGT' STO 0
CONT
\>> } } {
"C/S" {
\<< 0 MENU
"Motion? (True/Knots)"
":COURSE: " CRS
\>FMT \>STR +
"
:SPEED: " SPD
\>STR + + { 1 0 }
'V' 3 \>LIST INPUT
OBJ\> '1_knot'
DOUNIT 'SPD' STO
FMT\> 180 RANGE
'CRS' STO 0 CONT
\>> } } {
"P/T" {
\<< 0 MENU
"ENTER for std cond"
{
":PRESS (mb): 1010
:TEMPER (C): 10"
14 V } INPUT OBJ\>
'1_\^oC' DOUNIT
'TMPTR' STO '1_mbar
' DOUNIT 'PRESS'
STO 0 CONT
\>> } } {
"FORMAT" {
\<< 0 MENU
FFMT 1 +
IF DUP 3
==
THEN DROP
0
END
'FFMT' STO
CASE '
FFMT==2'
THEN
"(decimal)"
END '
FFMT==1'
THEN
"(dd.mmt)"
END '
FFMT==0'
THEN
"(dd.mmss)"
END
END 'FMT'
STO FFIX 0 CONT
\>> } } {
"EXIT" {
\<< 1 CONT
\>> } } }
TMENU
DO CLLCD
"INDEX " INDX \>FMT
\>STR + 2 DISP
"HEIGHT " HGT \>STR
"m" + + 3 DISP 1
FIX "MOTION " CRS
\>FMT \>STR + "T " +
SPD \>STR + "kn" + 4
DISP "P/T "
PRESS \>STR "mb " +
TMPTR \>STR + "C" +
+ 5 DISP FFIX
"FORMAT "
CASE '
FFMT==2'
THEN
"Decimal"
END '
FFMT==1'
THEN
"HMT"
END '
FFMT==0'
THEN
"HMS"
END "?"
END + 6
DISP 3 FREEZE HALT
0 MENU
UNTIL
END om MENU
\>>
\>>
ADDDR
\<< DEPTH SAVES 0
RCLMENU 28 MENU \>
dsv n om
\<<
IFERR OBS
IFERR
OBJ\>
THEN DROP
0
ELSE OBJ\>
DROP DROP
END 'n'
STO FMT DRLAT \>FMT
"DR_LAT" \>TAG \>STR
"
" + DRLON \>FMT
"DR_LON" \>TAG \>STR
+ { 1 0 } 'V' 3
\>LIST 28 MENU INPUT
0 MENU OBJ\> DTAG
FMT\> SWAP DTAG FMT\>
90 n 1 + 3 2 \>LIST
\>ARRY 'OBS' STO
THEN DEPTH
dsv  DROPN ABORT
END om MENU
\>> RESTS
\>>
DR
\<< DEPTH RCLMENU
28 MENU \> dsv om
\<<
IFERR FFIX
"Dead Reckoning?
"
FMT + DRLAT \>FMT
"DR_LAT" \>TAG \>STR
"
" + DRLON \>FMT
"DR_LON" \>TAG \>STR
+ { 1 0 } 'V' 3
\>LIST INPUT OBJ\>
FMT\> 'DRLON' STO
FMT\> 'DRLAT' STO
THEN DEPTH
dsv  DROPN ABORT
END om MENU
\>>
\>>
PLOTP
\<< DEPTH RCLMENU
28 MENU \> dsv smen
\<< SAVES DEG
IFERR
"Plot Center?
" FMT
+ LAT \>FMT "LAT"
\>TAG \>STR "
" + LON
\>FMT "LON" \>TAG
\>STR + { 1 0 } 'V'
3 \>LIST INPUT OBJ\>
FMT\> SWAP FMT\> 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 \> LON LAT g d
a l n N sc sc\Gl ssz
d0 d1 ll lm top bot
\<<
"Scale? (NMiles)" {
"9" 1 V } INPUT
OBJ\> ABS '1_nmi'
DOUNIT
IF DUP
0 ==
THEN
DROP
"SCALE\=/0 PLEASE"
MESS 0 DOERR
END 120
/ DUP 'sc' STO LAT
COS / 2.0469 * 180
MIN NEG 'sc\Gl' STO
ERASE { # 0h # 0h }
PVIEW LON sc\Gl + LON
RANGE LAT sc + 90
MIN DUP 'top' STO
DUP 3 ROLLD R\>C
PMAX LON sc\Gl  LON
RANGE LAT sc  90
MAX DUP 'bot' STO
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 '
LATsc>d+90a OR
LAT+sc
THEN 180 SWAP 
END MIN bot 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 +
LON 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
.0083333 ROT / 
LAT R\>C SWAP LON
.0083333 ROT / +
LAT R\>C LINE LON
LAT .0083333  R\>C
LON LAT .0083333 +
R\>C LINE
\>> { }
PVIEW
THEN DEPTH
dsv  DROPN
IF ERRN
# 0h \=/
THEN ERRM
CLLCD 4 DISP 2
FREEZE
END
END smen
MENU
\>> RESTS
\>>
ADV
\<< DEPTH SAVES
DEG RCLMENU 28 MENU
\> dsv om
\<<
IFERR 0 0 0
0 0 0 \> \Gh d \Gl l n
n3
\<<
"Motion? (nmi,deg true)"
{
":DISTANCE:
:COURSE: "
{ 1 0 } V } INPUT
OBJ\> FMT\> 180 RANGE
'\Gh' STO '1_nmi'
DOUNIT
IF 'SPD
\=/0'
THEN
DUP SPD / 'TF' STO+
END 60
/ 'd' STO 2 FIX
CLLCD "Old DR: "
DRLAT \>FMT \>STR +
" " + DRLON \>FMT
\>STR + 4 DISP OBS
IFERR
OBJ\>
THEN
DROP
ELSE
OBJ\> DROP SWAP DUP
'n' STO * 'n3' STO
1 n
FOR I
I 1 DISP 3 ROLLD
'l' STO '\Gl' STO \Gl l
d \Gh RMOVE SWAP 180
RANGE SWAP ROT n3
ROLLD n3 ROLLD n3
ROLLD
NEXT
{ n 3 } \>ARRY 'OBS'
STO
END
DRLON DRLAT d \Gh
CCMOVE 'DRLAT' STO
'DRLON' STO
"New DR: " DRLAT
\>FMT \>STR + " " +
DRLON \>FMT \>STR + 5
DISP FFIX 2 FREEZE
\>>
THEN DEPTH
dsv  DROPN ABORT
END om MENU
\>> RESTS
\>>
SAIL
\<< DEPTH SAVES
RCLMENU 28 MENU \>
dsv om
\<<
IFERR DEG 0
0 \> fr\Gl frl
\<<
"From? " FMT +
DRLAT \>FMT "Lat"
\>TAG \>STR "
" +
DRLON \>FMT "Lon"
\>TAG \>STR + { 1 0 }
'V' 3 \>LIST INPUT
OBJ\> FMT\> 'fr\Gl' STO
FMT\> 'frl' STO
"TO? " FMT + tol
\>FMT "Lat" \>TAG
\>STR "
" + to\Gl \>FMT
"Lon" \>TAG \>STR + {
1 0 } 'V' 3 \>LIST
INPUT OBJ\> FMT\>
'to\Gl' STO FMT\>
'tol' STO CLLCD 2
FREEZE { { "RHUMB"
\<< DEPTH
\> dsv
\<<
IFERR 0 MENU frl
fr\Gl tol to\Gl RHUMB
THEN DEPTH dsv 
DROPN ABORT
END 0 CONT
\>>
\>> } {
"GC"
\<< DEPTH
\> dsv
\<<
IFERR 0 MENU frl
fr\Gl tol to\Gl GC
THEN DEPTH dsv 
DROPN ABORT
END
\>> 0
CONT
\>> } {
"WAY"
\<< DEPTH
\> dsv
\<<
IFERR 0 MENU
"Scale? (nmi)" { ""
V } INPUT OBJ\> '1_
nmi' DOUNIT 60 /
frl fr\Gl tol to\Gl WAY
THEN DEPTH dsv 
DROPN ABORT
END
\>> 0
CONT
\>> } {
"VERT"
\<< DEPTH
\> dsv
\<<
IFERR 0 MENU frl
fr\Gl tol to\Gl VERTEX
THEN DEPTH dsv 
DROPN ABORT
END
\>> 0
CONT
\>> } {
"COMP"
\<< DEPTH
\> dsv
\<<
IFERR 0 MENU
"Composite" {
":Lat Limit:
:Scale: "
{ 1 0 } V } INPUT
OBJ\> '1_nmi' DOUNIT
60 / SWAP FMT\> frl
fr\Gl tol to\Gl COMP
THEN DEPTH dsv 
DROPN ABORT
END
\>> 0
CONT
\>> } {
"EXIT"
\<< 1
CONT
\>> } }
TMENU
DO
"Type?" PROMPT 0
MENU
UNTIL
END
\>>
THEN DEPTH
dsv  DROPN ABORT
END om MENU
\>> RESTS
\>>
WVIEW
\<< 2 FIX { }
SWAP { } 1 1 1 1
"Lat Lon Crs " FMT
+ 5 \>LIST DBR
IF 1 \=/
THEN DROP2
ELSE SWAP
DROP SWAP DUP ROT
GET
END FFIX
\>>
ERROR
\<< DEPTH \> dsv
\<<
IFERR DEG
OBS OBJ\> OBJ\> DROP
DROP 'N' STO { 2 2
} 0 CON 0 0 0 \> P
\GdH L \Gl
\<< 1 N
START
DERIVS '\GdH' STO '\Gl'
STO 'L' STO L SQ L
\Gl * DUP \Gl SQ { 2 2
} \>ARRY \GdH SQ / 'P'
STO+
NEXT P
INV OBJ\> DROP \v/ LAT
COS * '1_nmi' \>UNIT
"\177EW_Err" \>TAG SWAP
DROP SWAP DROP SWAP
\v/ '1_nmi' \>UNIT
"\177NS_Err" \>TAG
\>>
THEN DEPTH
dsv  DROPN ABORT
END
\>>
\>>
DRLAT
21.2922000006
DRLON
157.812100002
CORRECT
\<< DEG FMT\> 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  + \>FMT
\>>
DERIVS
\<< \> G D H
\<< 0 0 0 \> DC
HC LS
\<<
IF 'ABS(H
)\=/90'
THEN D
COS 'DC' STO H COS
'HC' STO LAT SIN
'LS' STO '(COS(LAT)
*SIN(D)LS*DC*COS(G
LON))/HC' \>NUM '
COS(LAT)*DC*SIN(G
LON)/HC' \>NUM 'ASIN
(SIN(D)*LS+DC*COS(
LAT)*COS(GLON))'
\>NUM H  ABS 60 *
ELSE 0 0
H
END
\>>
\>>
\>>
ABORT
\<< CLLCD
"Program Aborted
Press ENTER"
IF ERRN # 0h
\=/
THEN ERRM "
"
+ SWAP + ERR0
END MESS
\>>
RHUMB
\<< \> frl fr\Gl tol
to\Gl
\<< DEG to\Gl fr\Gl
RANGE 'to\Gl' STO 'LN
(TAN(45+tol/2)/TAN(
45+frl/2))' \>NUM '
\pi/180*(to\Glfr\Gl)'
\>NUM R\>C ARG 180
RANGE DUP \>FMT
"COURSE" \>TAG SWAP
IF DUP COS
ABS .0001 >
THEN COS
tol frl  SWAP /
ELSE to\Gl
fr\Gl  tol frl + 2 /
COS * SWAP SIN /
ABS
END 60 *
"DIST" \>TAG
\>>
\>>
GC
\<< \> frl fr\Gl tol
to\Gl
\<< DEG 'COS(
frl)*TAN(tol)SIN(
frl)*COS(to\Glfr\Gl)'
\>NUM 'SIN(fr\Glto\Gl)'
\>NUM R\>C ARG 180
RANGE \>FMT "COURSE"
\>TAG 'ACOS(SIN(frl)
*SIN(tol)+COS(frl)*
COS(tol)*COS(to\Gl
fr\Gl))' \>NUM 60 *
"DIST" \>TAG
\>>
\>>
COMP
\<< 0 0 0 0 0 0 0
0 \> scl ll frl fr\Gl
tol to\Gl vl v\Gl fc\Gl
tc\Gl n d d0 sn
\<< DEG frl fr\Gl
tol to\Gl VERTEX fr\Gl
RANGE 'v\Gl' STO 'vl'
STO to\Gl fr\Gl RANGE
'tc\Gl' STO
IF 'vl*SIGN
(ll)\<=ABS(ll)' 'ABS(
v\Gl(fr\Gl+tc\Gl)/2)>ABS
((fr\Gltc\Gl)/2)AND
ABS(vl)\=/90 AND ABS(
ll(frl+tol)/2)\>=ABS
((frltol)/2)' OR
THEN
"GC is OK: Hit ENTER"
MESS
ELSE DEPTH
'd0' STO to\Gl fr\Gl
RANGE fr\Gl
IF <
THEN 1
ELSE 1
END 'sn'
STO
IFERR ll
TAN INV DUP frl TAN
* ACOS NEG sn * fr\Gl
+ 0 RANGE 'fc\Gl' STO
tol TAN * ACOS sn *
to\Gl + 0 RANGE 'tc\Gl'
STO
THEN
DEPTH d0  DROPN
"No sol'n: Hit ENTER"
MESS
ELSE scl
frl fr\Gl ll fc\Gl WAY
DROP 'd' STO+ OBJ\>
'n' STO
IF 'RND
(fc\Gl,6)\=/RND(tc\Gl,6)'
THEN
OBJ\> SWAP DROP ll
fc\Gl ll tc\Gl RHUMB
'd' STO+ SWAP \>LIST
ELSE
DROP 1 'n' STO+
END scl
ll tc\Gl tol to\Gl WAY
DROP 'd' STO+ OBJ\>
n + \>LIST d "DIST"
\>TAG
END
END
\>>
\>>
VERTEX
\<< 0 \> frl fr\Gl
tol to\Gl C
\<< DEG frl fr\Gl
tol to\Gl GC DROP
FMT\> DUP 'C' STO
DUP SIN frl COS *
ABS ACOS frl 0 \>= 1
1 IFTE *
IF DUP 0 ==
THEN SWAP
DROP 0
ELSE DUP
ROT COS SWAP SIN /
ASIN NEG
IF 'C>180
'
THEN NEG
END fr\Gl +
IF 'ABS(
tol)>ABS(frl)AND
SIGN(tol)\=/SIGN(frl)
'
THEN 180
+ SWAP NEG SWAP
END 0
RANGE
END \>FMT
"V_Lon" \>TAG SWAP
\>FMT "V_Lat" \>TAG
SWAP
\>>
\>>
WAY
\<< \> scl frl fr\Gl
tol to\Gl
\<< DEG 0 frl
fr\Gl tol to\Gl GC SWAP
DROP 60 / frl fr\Gl
GETV DUP tol to\Gl
GETV CROSS DUP ABS
IF DUP 0 ==
THEN DROP2
IF 'RND(
frl,6)\=/RND(tol,6)OR
RND(fr\Gl,6)\=/RND(to\Gl,
6)'
THEN
"Ambiguous Sol'n" 3
DISP
END 0 fr\Gl
90  GETV
ELSE /
END NEG 0 0
\> d gcd r n d0 dsum
\<< DEPTH
'd0' STO
WHILE 'd<
gcd OR d==0'
REPEAT n
r d SMOVE V\> ASIN 3
ROLLD R\>C ARG 'd'
scl STO+
END tol
to\Gl gcd scl / FLOOR
2 + 'n' STO DUP2
"N/A" ROT \>FMT ROT
\>FMT ROT 3 \>LIST
DEPTH d0  ROLLD 1
n 1 
START 4
DUPN RHUMB 'dsum'
STO+ 3 ROLLD DROP2
3 ROLLD DUP2 5 ROLL
ROT \>FMT ROT \>FMT
ROT 3 \>LIST DEPTH
d0  ROLLD
NEXT
DROP2 n \>LIST dsum
DUP "DIST" \>TAG
SWAP gcd 60 *  '1_
nmi' \>UNIT "ADDD"
\>TAG
\>>
\>>
\>>
DOUNIT
\<< 55 CF
IFERR CONVERT
THEN DROP
END UVAL
\>>
SD
\<< 0 \> x
\<< DATE DUP
100 * FP 100 / 1.01
+ SWAP DDAYS 183 
183 / 'x' STO '(
15.762145+x*(
.02513+x*(1.15068+
x*(.02604+x*.62672
))))/60' \>NUM
\>>
\>>
RMOVE
\<< 0 0 0 0 \> \Gl l
d \Gh d\Gl dl n\Gl nl
\<< DRLON DRLAT
d \Gh CCMOVE DUP 'nl'
STO DRLAT  'dl'
STO DUP 'n\Gl' STO
DRLON  'd\Gl' STO l
\Gl d\Gl + GETV n\Gl 90 +
DUP COS SWAP SIN 0
\>V3 SWAP dl SMOVE
V\> ASIN 3 ROLLD R\>C
ARG SWAP
\>>
\>>
SMOVE
\<< \> n r d
\<< d COS r * n
n r DOT * 1 d COS 
* + r n CROSS d SIN
* +
\>>
\>>
CCMOVE
\<< 0 \> \Gl l d \Gh
l2
\<< l d \Gh MER l
+ DUP 'l2' STO
IF DUP ABS
90 \>=
THEN SIGN
90 * \Gl SWAP
ELSE
IF 'ABS(
COS(\Gh))<.0001'
THEN '
.998208257*d*SIN(\Gh
)/COS((l+l2)/2)*\v/(1
(ee*SIN((l+l2)/2))
^2)' \>NUM
ELSE l l2
\Gh DLo
END \Gl +
SWAP
END
\>>
\>>
MER
\<< \> l1 d \Gh
\<< '
.998208256722/(1ee
^2)*\.S(l1,l1+d*COS(\Gh
),(1(ee*SIN(l))^2)
^1.5,l)' \>NUM
\>>
\>>
DLo
\<< 0 0 \> l1 l2 \Gh
sl1 sl2
\<< l1 SIN
'sl1' STO l2 SIN
'sl2' STO '
57.2957795131*TAN(
\Gh)*(ATANH((sl2sl1)
/(1sl1*sl2))ee*
ATANH(ee*(sl2sl1)/
(1ee^2*sl2*sl1)))'
\>NUM
\>>
\>>
GETV
\<< \> l \Gl
\<< l COS \Gl COS
* l COS \Gl SIN * l
SIN \>V3
\>>
\>>
ee
8.18188106628E2
FMT "(dd.mmt)"
FFMT 1
FFIX
\<<
IF 'FFMT==1'
THEN 3 FIX
ELSE 4 FIX
END
\>>
FMT\>
\<<
CASE 'FFMT==1
'
THEN HMT\>
END 'FFMT==
0'
THEN HMS\>
END
END
\>>
\>FMT
\<<
CASE 'FFMT==1
'
THEN \>HMT
END 'FFMT==
0'
THEN \>HMS
END
END
\>>
\>HMT
\<< 4 RND DUP IP
SWAP FP .6 * +
\>>
HMT\>
\<< DUP IP SWAP
FP 1.66666667 * +
\>>
SVSTK {
# 81398000000FF0h
# 0h }
RESTS
\<< SVSTK STOF
FFIX
\>>
SAVES
\<< RCLF 'SVSTK'
STO 20 CF 21 CF
22 SF 55 CF 40
CF
\>>
\GmST
\<< 0 0 0 \> s2 s3
s4
\<< 2 SK 's2'
STO 3 SK 's3' STO 4
SK 's4' STO '(s3+\v/
(s3^23*s4*(s21)))
/(3*s4)' \>NUM RE
UBND MIN
\>>
\>>
UBND
\<< \Ga1 \Gb1 ABS 
\Ga2 \Gb2 ABS  \Ga3 \Gb3
ABS  MIN MIN
\>>
LBND
\<< \Ga1
1.73205080757 \Gb1
ABS *  \Ga2
1.73205080757 \Gb2
ABS *  \Ga3
1.73205080757 \Gb3
ABS *  MIN MIN
\>>
SK
\<< \> k
\<< '\Gb1^2/\Ga1^k+
\Gb2^2/\Ga2^k+\Gb3^2/\Ga3^k
' \>NUM
\>>
\>>
G\Gm
\<< \Gb1 \Ga1 \Gm  /
SQ \Gb2 \Ga2 \Gm  / SQ +
\Gb3 \Ga3 \Gm  / SQ + 1

\>>
ASK
\<< { "YES" "" ""
"" "" "NO" } TMENU
0
DO DROP 1
WAIT
UNTIL DUP {
11.1 16.1 } SWAP
POS DUP
IF NOT
THEN 880 .1
BEEP
END
END 0 MENU
\>>
MLIMB { { "LL"
\<< 1 CONT
\>> } "" { "UL"
\<< 1 CONT
\>> } "" { "CENT"
\<< 0 CONT
\>> } "" }
MBODY { { "SUN"
\<< "S" CONT
\>> } { "MOON"
\<< "M" CONT
\>> } { "VENUS"
\<< "VM" CONT
\>> } { "MARS"
\<< "VM" CONT
\>> } { "PLANET"
\<< "P" CONT
\>> } { "STAR"
\<< "T" CONT
\>> } }
PSCALE
\<< \> s a
\<<
IF 's\=/0'
THEN 'a/(
360+a/s)' \>NUM
ELSE 0
END
\>>
\>>
tol 36.9833333353
to\Gl 25.166666667
LON 157.812100002
LAT 21.2922000006
IERR
7.28357026453E2
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+(upg1)*(d1
d2)/(g1g2)' \>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
\<< \> \Gl
\<<
WHILE DUP
180 \Gl + >
REPEAT 360

END
WHILE DUP
180 \Gl + <
REPEAT 360
+
END
\>>
\>>
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)>
90ABS(d)+a'
THEN 180 +
END LON
RANGE l R\>C
\>>
\>>
CST { SOLVE ADDOB
SETUP INIT ADV
ADDDR DR PLOTP SAIL
WVIEW ERROR }
REFRACT
\<< 0 \> h rp
\<< '1/TAN(h+
7.31/(h+4.4))' \>NUM
'rp' STO 'rp*((
PRESS80)/930)/(1+
.00008*(rp+39)*(
TMPTR10))' \>NUM 60
/
\>>
\>>
MESS
\<< 3 DISP 7
FREEZE 0 WAIT DROP
\>>
PPAR {
(166.965199761,17.1250000005)
(148.658133575,25.4583333339)
X 0 (0,0) FUNCTION
Y }
T\Gg 6
G\Gg 231.103333334
PRESS 1010
TMPTR 10
a0 '(G12*G23G13
*G22)*G13+(G11*G23
G12*G13)*G23(G11*
G22G12^2)*G33'
a1 'G11*G22G12^2
+G11*G33G13^2+G22*
G33G23^2'
TF 162.361904762
CRS 315
SPD 7
EV3 '2*\v/Q*COS((\Gh
+360)/3)+N/3'
EV2 'N\Ga1\Ga3'
EV1 '2*\v/Q*COS(\Gh/
3)+N/3'
\Gm
5.11011480871E8
\Gb3 12.8516757606
\Gb2
5.92502509102E3
\Gb1
3.6381134513E7
E3
[ .390241406661 .920339419461 2.62106372201E2 ]
E2
[ 2.08258845121E3 2.93500423979E2 .999567025186 ]
E1
[ .920710196011 .390017549971 .013380795022 ]
INTERP
\<< \> T V1 V2
\<< V1 V2 V1 
T2 T1  / T T1  *
+
\>>
\>>
GSUM
\<< 0 0 0 0 0 \>
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
\<< OBJ\> DROP \> U
V W
\<<
IF 'ABS(U)>
1'
THEN U SIGN
'U' STO
END U ASIN
V W R\>C ARG \>FMT
"LON" \>TAG SWAP
\>FMT "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*G12G11*
G23+G23*EV' \>NUM '
G11*G22SQ(G12)(
G11+G22)*EV+SQ(EV)'
\>NUM { 3 } \>ARRY
DUP ABS
IF DUP 0 \=/
THEN /
ELSE DROP
END
\>>
\>>
\Ga2 .0156829589
\Ga3 19.984316546
\Ga1 .00000049507
\Gh 'ACOS(R1/Q^1.5)
'
R1 'A0/2+N/3*(A1/
6Q)'
Q '(N/3)^2A1/3'
N 20
A0
1.5585765436E7
A1 .313423116958
G33 .0293985547
G23 .48161523277
G22 16.9272222537
G13 .204441742038
G12 7.17745707216
G11 3.04337919163
R
[ 5.01524332669 11.8280778613 .330928157349 ]
GHA2
43.1883333337
DEC2
10.3033333339
T2 15
GHA1 28.686666668
DEC1
10.0483333334
T1 14
LU 1
SEMI .27166666721
HP .98666666864
HGT 2.7432
INDX 0
BODY "M"
END
 CUT HERE AND AT TOP 