metcalf@akala.IFA.Hawaii.Edu (Tom Metcalf) (12/27/90)
Following is version 3.0 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 2 is that the program will now advance the lines
of position, allowing a running fix. I have also added some error checking,
simplified the input of stellar data, and fixed a slight error in the "PLOTP"
routine which cut off the bottom of the plot in a few cases.
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
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.
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. "OBS" is a matrix with each row
representing one observation in the format: GHA DEC ALT.
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; 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).
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 arc minutes
should be input as 0.1612 since 0.2 min = 12". Sun only.
d) HP: Parallax in degrees, "hms" 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. 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).
e) LIMB: Select lower limb (LL), upper limb (UL), or disk center (CENT)
from the menu keys. (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) If the body is a *star* the program will ask for GHA-ARIES 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. All inputs must be in hms format and all must be
input on the appropriate line *before* pressing "ENTER". To move to the
next item, use the down-arrow key. Skip to (j).
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 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
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. 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 terminate if the time and GHA do not increase.
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 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 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.
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
and should be the same for all observations.
3. Enter the observations:
a) Enter the time of the observation (hms format) and then the
uncorrected altitude (hms format) onto the stack. 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 2h and 2i. It should not pass though 00:00.
b) Run the "CORRECT" program to correct the observation for index,
dip, refraction, and parallax.
c) Run the "ADDOB" program to add the observation to the "OBS" variable.
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 the correction for
course/speed.
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 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(m) 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(m) for all bodies, unless the LOP's are advanced appropriately (see
step 7).
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 (very unlikely) 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 plot of the LOPs (step 6) 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".
5. If all observations are for a *single celestial body* (no DR) over a
relatively short period of time, 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 this
program the position fix from "SOLVE" must be left as the first two
numbers on the stack (longitude in level 2, latitude in level 1; this
defines the center of the plot). The "PLOTP" program will 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 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 cross-hair is drawn at the position of the fix 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 (*not* hms format) of the plot after pressing
<orange-shift GRAPH>.
7. 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 2m)
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(m) 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.
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.
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 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
IF '\Ga1==0 AND
\Gb1==0'
THEN
"AMBIGUOUS SOLUTION"
MESS KILL
END 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>30'
END
IF '\Gm>\Ga1 OR
NIT>30'
THEN Gx 'x' {
0 '-N' \Ga1 } ROOT
'\Gm' STO
END
IF '\Gm>\Ga1 OR \Gm
<-N'
THEN
"CONVERGENCE ERROR"
END UVW OBJ\->
DROP OUT CLLCD
"Update DR?" 2 DISP
DUP2 \->STR 4 DISP
\->STR 5 DISP ASK
IF 11.1 ==
THEN DUP2
HMS\-> 'DRLAT' STO
HMS\-> 'DRLON' STO
END
\>>
ADDOB
\<< 0 \-> TM A n
\<< TM HMS\->
'TM' STO
IF TM T1 <
TM T2 > BODY "T"
SAME NOT AND OR
THEN TM
\->HMS A
"Error:Bad Time"
MESS KILL
END A HMS\->
'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
\>>
\>>
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 - + \->HMS
\>>
SETUP
\<< CLLCD 2
FREEZE MBODY TMENU
"BODY?" PROMPT
'BODY' STO 0 MENU
"INDEX? (dd.mmss)"
{ "" V } INPUT OBJ\->
HMS\-> 'INDX' STO
IF BODY "S"
SAME
THEN
DO
"SEMI-D? (dd.mmss)"
{ "" V } INPUT OBJ\->
HMS\-> '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
"HP? (dd.mmss)" {
"" V } INPUT OBJ\->
HMS\-> '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
"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
IF BODY "T"
SAME
THEN "Star" {
":GHA\Gg:
:SHA:
:DEC:
:TIM: "
{ 1 0 } V } INPUT
OBJ\-> HMS\-> DUP 'T1'
STO 1 + 'T2' STO
HMS\-> DUP 'DEC1' STO
'DEC2' STO HMS\->
SWAP HMS\-> + DUP
'GHA1' STO
15.041067 + 'GHA2'
STO
ELSE
"Linear Interp 1" {
":GHA1:
:DEC1:
:TIM1:"
{ 1 0 } V } INPUT
OBJ\-> HMS\-> 'T1' STO
HMS\-> 'DEC1' STO
HMS\-> 'GHA1' STO
"Linear Interp 2" {
":GHA2:
:DEC2:
:TIM2:"
{ 1 0 } V } INPUT
OBJ\-> HMS\-> 'T2' STO
HMS\-> 'DEC2' STO
HMS\-> 'GHA2' STO
END
IF 'T1>T2'
THEN
"Error:T1>T2" MESS
KILL
END
IF 'GHA1>GHA2
'
THEN
"Error:GHA1>GHA2"
MESS KILL
END
"SPEED? (Knots)" {
"" V } INPUT OBJ\->
'SPD' STO
IF 'SPD\=/0'
THEN
"COURSE? (True)" {
"" V } INPUT OBJ\->
HMS\-> 180 RANGE
'CRS' STO DR
"TIME OF FIX?" { ""
V } INPUT OBJ\-> HMS\->
'TF' STO
ELSE 0 'CRS'
STO 0 'TF' STO
END
\>>
ADDDR
\<< 0 \-> n
\<< OBS
IFERR OBJ\->
THEN DROP 0
ELSE OBJ\->
DROP DROP
END 'n' STO
"dd.mmss" DRLAT
\->HMS "DR_LAT" \->TAG
\->STR "
" + DRLON
\->HMS "DR_LON" \->TAG
\->STR + { 1 0 } 'V'
3 \->LIST INPUT OBJ\->
HMS\-> SWAP HMS\-> 90 n
1 + 3 2 \->LIST \->ARRY
'OBS' STO
\>>
\>>
DR
\<<
"Dead Reckoning? (hms)"
DRLAT \->HMS "DR_LAT"
\->TAG \->STR "
" +
DRLON \->HMS "DR_LON"
\->TAG \->STR + { 1 0 }
'V' 3 \->LIST INPUT
OBJ\-> HMS\-> 'DRLON'
STO HMS\-> 'DRLAT'
STO
\>>
PLOTP
\<<
IF DEPTH 2 <
THEN
"LON/LAT NOT ON STACK"
MESS KILL
END 2 DUPN
HMS\-> 'LAT' STO HMS\->
'LON' STO 0 0 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
top bot
\<<
"Scale? (NMiles)" {
"9" -1 V } INPUT
OBJ\-> ABS
IF DUP 0 ==
THEN DROP
"SCALE\=/0 PLEASE"
MESS KILL
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 'LAT-
sc>d+90-a OR LAT+sc
<d-90+a'
THEN
ELSE top
d 90 a - +
IF DUP
90 >
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
\>>
ADV
\<< 0 0 0 0 0 0 \->
\Gh d \Gl l n n3
\<<
"Motion? (nmi,deg true)"
{
":DISTANCE:
:COURSE: "
{ 1 0 } V } INPUT
OBJ\-> HMS\-> 180 RANGE
'\Gh' STO 60 / 'd'
STO 2 FIX CLLCD
"Old DR: " DRLAT
\->HMS \->STR + " " +
DRLON \->HMS \->STR + 3
DISP OBS
IFERR OBJ\->
THEN DROP
ELSE OBJ\->
DROP SWAP DUP 'n'
STO * 'n3' STO 1 n
START 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 \->HMS \->STR +
" " + DRLON \->HMS
\->STR + 4 DISP 4 FIX
3 FREEZE
\>>
\>>
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
\>>
\>>
DRLAT
31.6262161969
DRLON
15.0182575167
NIT 3
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
\>> } }
Gx '(\Gb1/(\Ga1-x))^2
+(\Gb2/(\Ga2-x))^2+(\Gb3/
(\Ga3-x))^2-1'
x
-3.44048231887E-8
RMOVE
\<< 0 0 0 0 0 0 \->
\Gl l d \Gh d\Gl dl n\Gl nl
r n
\<< DRLON DRLAT
d \Gh CCMOVE DUP 'nl'
STO DRLAT - 'dl'
STO DUP 'n\Gl' STO
DRLON - 'd\Gl' STO l
COS \Gl COS * l COS \Gl
SIN * l SIN \->V3 'r'
STO 0 0 -1 \->V3 r d\Gl
SMOVE 'r' STO n\Gl 90
+ DUP COS SWAP SIN
0 \->V3 r 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
pl
\<< \Gh COS d * l
+
IF DUP ABS
90 \>=
THEN SIGN
90 * \Gl SWAP
ELSE
IF 'ABS(
COS(\Gh))<.01'
THEN '-
SIN(\Gh)*\.S(0,d,1/COS(
l+x*COS(\Gh)),x)'
\->NUM
ELSE 45 l
2 / + 'pl' STO '
-57.295779513*TAN(\Gh
)*LN(TAN(pl+d*COS(\Gh
)/2)/TAN(pl))' \->NUM
END \Gl +
SWAP
END
\>>
\>>
PSCALE
\<< \-> s a
\<<
IF 's\=/0'
THEN 'a/(
360+a/s)' \->NUM
ELSE 0
END
\>>
\>>
LON 15.0182575167
LAT 31.6262161969
IERR
9.78679938533E-2
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
\<< \-> \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)>
90-ABS(d)+a'
THEN 180 +
END LON
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 { SOLVE ADDOB
CORRECT SETUP ADV
ADDDR PLOTP 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
/
\>>
\>>
MESS
\<< 3 DISP 7
FREEZE 0 WAIT DROP
\>>
PPAR {
(25.0345573463,27.4595495302)
(5.0019576871,35.7928828636)
\Gm 0 (0,0) FUNCTION
Y }
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
CRS 0
SPD 0
EV3 '-2*\v/Q*COS((\Gh
+360)/3)+N/3'
EV2 'N-\Ga1-\Ga3'
EV1 '-2*\v/Q*COS(\Gh/
3)+N/3'
OLD
4.02732539364E-4
\Gm
4.02732539364E-4
\Gb3 .637635140172
\Gb2 -.684386502275
\Gb1 -.287515572887
E3
[ .788429842018 -.134721632922 .600190358001 ]
E2
[ .101065522502 -.934090506002 -.342433477858 ]
E1
[ -.606765312638 -.33064332367 .722849118343 ]
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
IF DUP 0 \=/
THEN /
ELSE DROP
END
\>>
\>>
\Ga2 .86588407145
\Ga3 1.46600912668
\Ga1 .66810680187
\Gh 'ACOS(R1/Q^1.5)
'
R1 'A0/2+N/3*(A1/
6-Q)'
Q '(N/3)^2-A1/3'
N 3
A0 -.84809073318
A1 2.8273476583
G33 .9787254307
R
[ .608017170664 .648440792514 .409228934633 ]
G23
-.001255484102
G22 .855154290614
G13 .370729035826
G12 -.10342306594
G11 1.16612027869
GHA2 1
DEC2 1
T2 1
GHA1 0
DEC1 0
T1 0
LU -1
SEMI 0
HP 0
HGT 0
INDX 0
BODY "S"
END
------------------------------ CUT HERE -----------------------