The Museum of HP Calculators

Calculating a Fix for the HP-41

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Overview

-After measuring the altitudes of n = 2 , 3 ( or more ) stars,
the 2 following programs compute your geographical position: Longitude & Latitude
-The second routine is longer but faster than the first one.

Program #1

Data Registers: R00 = "T" R01 to R05 are used by "EXY" R06 = 4n+10 R07 to R10: temp

                                R11 = ST₁ -RA₁                     R15 = ST₂ -RA₂                R19 = ST₃ -RA₃           .......................   ST = Sidereal Time
                                R12 = Cos (decl₁)                   R16 = Cos (decl₂)              R20 = Cos (decl₃)        .......................   RA = Right-Ascension
                                R13 = Sin (decl₁)                    R17 = Sin (decl₂)                R21 = Sin (decl₃)          .......................   decl = declination
                                R14 = Sin (alt₁)                       R18 = Sin (alt₂)                  R22 = Sin (alt₃)             .......................    alt = altitude

-When the program stops, R01 = your longitude & R02 = your latitude ( in degrees and decimals )

Flags:              F22 F29
Subroutines: "MST"   Mean Sidereal Time at Greenwich ( cf "Rising-Transit-Setting for the HP-41" )
                        "H0-H" ( cf "Atmospheric Refraction for the HP-41" )
                        "EXY"   ( cf "Extrema for the HP-41" , 2 dimensional problems )

01 LBL "FIX"
02 1
03 STO 04
04 11.9
05 STO 05
06 FIX 0
07 CF 29
08 DEG
09 LBL 01
10 "DATE^TIME"
11 ARCL 04
12 "~?"                        ( append ? )
13 FIX 4
14 CF 22
15 PROMPT
16 FC?C 22
17 GTO 03
18 FIX 0
19 XEQ "MST"
20 LASTX
21 STO IND 05
22 "RA^DECL"
23 ARCL 04
24 "~?"                        ( append ? )
25 FIX 4
26 PROMPT
27 FIX 0
28 X<>Y
29 HR
30 ST- IND 05
31 15
32 ST* IND 05
33 ISG 05
34 X<> Z
35 HR
36 1
37 P-R
38 STO IND 05
39 ISG 05
40 X<>Y
41 STO IND 05
42 ISG 05
43 "ALT"
44 ARCL 04
45 "~?"                        ( append ? )
46 PROMPT
47 XEQ "H0-H"          In the example below, I've used the refraction of a standard atmosphere,
48 HR                         other choices may of course be better according to the atmospueric conditions.
49 SIN
50 STO IND 05
51 ISG 04
52 CLX
53 RCL 05
54 INT
55 STO 06
56 ISG 05
57 GTO 01
58 LBL "T"                 ( or another global name, the same as line 97 )
59 STO 07
60 CLX
61 SIGN
62 P-R
63 STO 08
64 X<>Y
65 STO 09
66 RCL 06
67   E3
68 /
69 11
70 +
71 STO 10
72 CLX
73 LBL 02
74 RCL 07
75 RCL IND 10
76 +
77 COS
78 RCL 08
79 *
80 ISG 10
81 RCL IND 10
82 *
83 RCL 09
84 ISG 10
85 RCL IND 10
86 *
87 +
88 ISG 10
89 RCL IND 10
90 -
91 X^2
92 +
93 ISG 10
94 GTO 02
95 RTN
96 LBL 03
97 "T"
98 ASTO 00
99 "STEP^LON^LAT"
100 PROMPT
101 X<>Y
102 CLA
103 XEQ "EXY"
104 X<>Y
105 HMS
106 X<>Y
107 HMS
108 SF 29
109 CLD
110 END

( 215 bytes / SIZE 4n+21 )

      STACK        INPUTS      OUTPUTS

           Z             /        min(f)

           Y             / latitude ( ° ' " )

           X             / longitude ( ° ' " )

Example: On 2006/04/04 you measure the altitude of n = 3 stars:

   at 6h00m (UT) alpha Lyr       ( RA = 18h37m09s , Decl = 38°47'21" )      Altitude = 29°41'59"
   at 6h01m (UT) zeta UMaj      ( RA = 13h24m10s , Decl = 54°53'34" )      Altitude = 60°29'56"
   at 6h02m (UT) gamma UMin ( RA = 15h20m43s , Decl = 71°48'43" )      Altitude = 42°26'03"

-SIZE 023 or greater, then:

XEQ "FIX" the HP-41 displays "DATE^TIME1?" ( enter the different data and press R/S )

2006.0404 ENTER^
       6.00        R/S           -------------------- "RA^DECL1?"
     18.3709 ENTER^
     38.4721    R/S           -------------------- "ALT1?"
     29.4159    R/S           -------------------- "DATE^TIME2?"

2006.0404 ENTER^
       6.01        R/S           -------------------   "RA^DECL2?"
     13.2410 ENTER^
     54.5334    R/S           -------------------   "ALT2?"
     60.2956    R/S           -------------------   "DATE^TIME3?"

2006.0404 ENTER^
       6.02        R/S           -------------------   "RA^DECL3?"
     15.2043 ENTER^
     71.4843    R/S           -------------------   "ALT3?"
     42.2603    R/S           -------------------   "DATE^TIME4?"                      ( since we have only 3 stars, press R/S without any digit entry )

                      R/S           -------------------   "STEP^LON^LAT"                  ( enter a stepsize, say 2° and your estimated position, for instance -70° ; +20° )
           2       ENTER^
        -70      ENTER^
         20         R/S           the HP-41 executes "EXY" and displays the successive longitude-values, expressed in degrees and decimals
                                       and 8mn49s later:

     >>>   your longitude = X-register = -74°40'30"                 So, you are on the Atlantic Ocean, near the Bahamas!
     >>>   your latitude    = Y-register = +25°49'41"
                                         Z-register = f_min ~ 2.4 10^-11 = R03   this value must be a very small number.

-Longitudes are reckoned positively Eastwards from the meridian of Greenwich.

-If the process seems to diverge, stop the program, GTO "FIX" , XEQ 03 and key in new estimated step/longitude/latitude.
-Too bad initial guesses may lead to a relative minimum and a wrong position but in this case, Z-output is not very small.
-However, if n = 2, the system has 2 solutions. Therefore, measuring 3 altitudes is safer.

-This program solves the system of non-linear equations:

f_i(Long,Lat) = -sin(alt_i) + sin(Lat) sin(decl_i) + cos(Lat) cos(decl_i) cos(ST_i - RA_i + Long) = 0 i = 2 , 3 , ....

by a least-squares method: It minimizes SUM_i (f_i)² using the simple ( and slow ) EXY-algorithm
-The next program solves the system we get after equating to zero the partial derivatives.

Program #2

Data Registers: R00 = "T" R01 to R11 are used by "SXY" R12 = 4n+22 R13 to R22: temp

                                R23 = ST₁ -RA₁                     R27 = ST₂ -RA₂                R31 = ST₃ -RA₃           .......................   ST = Sidereal Time
                                R24 = Cos (decl₁)                   R28 = Cos (decl₂)              R32 = Cos (decl₃)        .......................   RA = Right-Ascension
                                R25 = Sin (decl₁)                    R29 = Sin (decl₂)                R33 = Sin (decl₃)          .......................   decl = declination
                                R26 = Sin (alt₁)                       R30 = Sin (alt₂)                  R34 = Sin (alt₃)             .......................    alt = altitude

-When the program stops, R01 = your longitude & R02 = your latitude ( in degrees and decimals )

Flags:              F22 F29
Subroutines: "MST"   Mean Sidereal Time at Greenwich ( cf "Rising-Transit-Setting for the HP-41" )
                        "H0-H" ( cf "Atmospheric Refraction for the HP-41" )
                        "SXY"   2x2 non-linear systems ( cf "Linear & Non-Linear Systems for the HP-41" )

01 LBL "FIX2"
02 1
03 STO 10
04 23.9
05 STO 11
06 FIX 0
07 CF 29
08 DEG
09 LBL 01
10 "DATE^TIME"
11 ARCL 10
12 "~?"                        ( append ? )
13 FIX 4
14 CF 22
15 PROMPT
16 FC?C 22
17 GTO 03
18 FIX 0
19 XEQ "MST"
20 LASTX
21 STO IND 11
22 "RA^DECL"
23 ARCL 10
24 "~?"                        ( append ? )
25 FIX 4
26 PROMPT
27 FIX 0
28 X<>Y
29 HR
30 ST- IND 11
31 15
32 ST* IND 11
33 ISG 11
34 X<> Z
35 HR
36 1
37 P-R
38 STO IND 11
39 ISG 11
40 X<>Y
41 STO IND 11
42 ISG 11
43 "ALT"
44 ARCL 10
45 "~?"                        ( append ? )
46 PROMPT
47 XEQ "H0-H"
48 HR
49 SIN
50 STO IND 11
51 ISG 10
52 CLX
53 RCL 11
54 INT
55 STO 12
56 ISG 11
57 GTO 01
58 LBL "T"                 ( or another global name, the same as line 126 )
59 STO 13
60 CLX
61 SIGN
62 P-R
63 STO 14
64 X<>Y
65 STO 15
66 RCL 12
67   E3
68 /
69 23
70 +
71 STO 18
72 CLX
73 STO 16
74 STO 17
75 STO 19
76 LBL 02
77 RCL IND 18
78 RCL 13
79 +
80 STO 20
81 COS
82 ISG 18
83 RCL IND 18
84 STO 21
85 *
86 RCL 14
87 *
88 RCL 15
89 ISG 18
90 RCL IND 18
91 STO 22
92 *
93 +
94 ISG 18
95 RCL IND 18
96 -
97 STO 11
98 X^2
99 ST+ 19
100 LASTX
101 RCL 20
102 SIN
103 *
104 RCL 21
105 *
106 ST+ 16
107 LASTX
108 RCL 20
109 COS
110 *
111 RCL 15
112 *
113 RCL 14
114 RCL 22
115 *
116 -
117 RCL 11
118 *
119 ST+ 17
120 ISG 18
121 GTO 02
122 RCL 17
123 RCL 16
124 RTN
125 LBL 03
126 "T"
127 ASTO 00
128 "LON^LAT GUESS1?"
129 PROMPT
130 STO 02
131 X<>Y
132 STO 01
133 "LON^LAT GUESS2?"
134 PROMPT
135 STO 04
136 X<>Y
137 STO 03
138 FIX 9
139 CLA
140 XEQ "SXY"
141 RCL 19
142 X<> Z
143 HMS
144 X<>Y
145 HMS
146 FIX 4
147 SF 29
148 CLD
149 END

( 296 bytes /SIZE 023+4n )

      STACK        INPUTS      OUTPUTS

           Z             /        min(f)

           Y             / latitude ( ° ' " )

           X             / longitude ( ° ' " )

Example: The same example as above: On 2006/04/04,

-SIZE 035 or greater, then:

XEQ "FIX" the HP-41 displays "DATE^TIME1?" ( enter the different data and press R/S )

                      R/S           -------------------   "LON^LAT GUESS1?"             ( enter your estimated position, for instance -70 ; +20 )
        -70      ENTER^
          20        R/S           -------------------   "LON^LAT GUESS2?"             ( enter a second estimated position, say -80 ; +30 )
        -80      ENTER^
         30         R/S           the HP-41 executes "SXY" and displays the successive longitude-values, expressed in degrees and decimals
                                      ( each iteration requires about 35 seconds if n = 3 ) and eventually,

     >>>   your longitude = X-register = -74°40'30"
     >>>   your latitude    = Y-register = +25°49'41"
                                         Z-register = f_min ~ 2.3 10^-11 = R19   this value must be a very small number.

-Longitudes are reckoned positively Eastwards from the meridian of Greenwich.

-If the process seems to diverge, stop the program, GTO "FIX2" , XEQ 03 and key in new estimated longitudes/latitudes.
-Too bad initial guesses may lead to a relative minimum and a wrong position but in this case, Z-output is not very small.
-However, if n = 2, the system has 2 solutions. Therefore 3 measured altitudes are safer.
-"SXY" may also stop prematurely if the jacobian = 0 ( always check Z-register or R19 )
-"SXY" tries to produce very accurate results, but here, 8 decimals are superfluous.
So, execution time can be reduced if you modify the final test in the "SXY" listing.

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall

STACK	INPUTS	OUTPUTS
Z	/	min(f)
Y	/	latitude ( ° ' " )
X	/	longitude ( ° ' " )