Geod-3axial

Geodesics on a Triaxial Ellipsoid(2) for the HP-41

Overview

0°) Geodetic Longitudes & Latitudes >>> U & V
1°) V = V(U)

a) Runge-Kutta of order 6
b) Runge-Kutta of order 4

2°) U = U(V)

a) Runge-Kutta of order 6
b) Runge-Kutta of order 4

-These programs employ Euler-Lagrange method to compute the geodesic distance between 2 points on a triaxial ellipsoid.
-It also works for biaxial ellipsoids.

-These routines use the coordinates u & v defined by:

x = a Cos u Cos v
y = b Sin u Cos v    where a , b , c are the semi-axes of the triaxial ellipsoid and x , y , z are the Cartesian coordinates
z = c Sin v

-We have to minimize the integral § L du where L = SQRT ( P + 2.Q v' + R v'² ) ( v' = dv/du )

P = Cos² v ( a² Sin² u + b² Cos² u )
    Q = (1/4) ( a² - b² ) Sin 2.u Sin 2.v
R = Sin² v ( a² Cos² u + b² Sin² u ) + c² Cos² v

-Euler-Lagrange equation: ¶L / ¶v = (d/du) ( ¶L / ¶v' )               ( where ¶ = partial derivative ) becomes:

v'' ( P.R - Q² ) / ( P + 2.Q v' + R v'² ) = - dQ/du - v' dR/du + (1/2) ( ¶P/¶v + 2.v' ¶Q/¶v + v'² ¶R/¶v ) + (1/2) ( ( Q + R.v' ) / ( P + 2.Q v' + R v'² ) ) ( dP/du + 2.v' dQ/du + v'² dR/du )

( dP/du = ¶P/¶u + v' ¶P/¶v ... etc ... )

-The second-order differential equation is solved by a Runge-Kutta formula and we use gardenhose method to find v'₁

-The integral is computed by Newton-Cotes6 ( or NC9 ) formula

NC6: §_x^x+6.h f(t).dt ~ (h/140) [41.f(x)+216.f(x+h)+27.f(x+2.h)+272.f(x+3.h)+27.f(x+4.h)+216.f(x+5.h)+41.f(x+6.h)]
NC9: §_x^x+9.h f(t).dt ~ (9.h/89600) [2857.f(x)+15741.f(x+h)+1080.f(x+2.h)+19344.f(x+3.h)+5778.f(x+4.h)+5778.f(x+5.h)+19344.f(x+6.h)+1080.f(x+7.h)+15741.f(x+8.h)+2857.f(x+9.h)]

-But first, let's transform the geodetic coordinates L & b into u & v

0°) Geodetic Longitudes & Latitudes >>> U & V

-"LB-UV" calculates u & v in radians

-This program also stores a = 6378.171645 km b = 6378.101575 km c = 6356.751868 km into R06 R07 R08

Data Registers: R11-R12-R13: temp

R06 = a R07 = b R08 = c

Flags: /
Subroutines: /

01 LBL "LB-UV"
02 DEG
03 HR
04 STO 12
05 COS
06 X<>Y
07 HR
08 14.92911
09 +
10 STO 11
11 SIN
12 *
13 X^2
14 6356.751868
15 STO 08

16 X^2
17 X<>Y
18 6378.101575
19 STO 07
20 X^2
21 6378.171645
22 STO 06
23 X^2
24 ST- Y
25 ST- T
26 ST/ T
27 /
28 STO 13
29 *
30 X<>Y

31 RCL 12
32 SIN
33 X^2
34 *
35 +
36 1
37 ST+ 13
38 ST+ Z
39 +
40 SQRT
41 RCL 06
42 X<>Y
43 /
44 RCL 11
45 RCL 12

46 X<> Z
47 X<>Y
48 RDN
49 P-R
50 R^
51 X<>Y
52 P-R
53 R^
54 ST* T
55 X<> 13
56 ST* Z
57 CLX
58 RCL 06
59 /

60 X<>Y
61 RCL 07
62 /
63 R^
64 RCL 08
65 /
66 X<> Z
67 RAD
68 R-P
69 X<>Y
70 RDN
71 R-P
72 X<> T
73 END

( 131 bytes / SIZE 014 )

STACK	INPUTS	OUTPUTS
Y	L ( ° ' " )	V ( rd )
X	b ( ° ' " )	U ( rd )

Example:    L = 141°23'45" & b = 64°12'34"

141.2345 ENTER^
64.1234 XEQ "LB-UV" >>>> U = 2.728389004 rd
X<>Y V = 1.119347728 rd

Notes:

-Change lines 08-14-18-21 if you want to use another Earth ellipsoid
-Line 08 is 14°92911 if the equatorial major axis in the direction   14°92911W

1°) V = V(U)

     a) Runge-Kutta of order 6

-We have to find v'₁ such that v(u₂) = v₂

-Coordinates u & v must be expressed in radians.
-"LB-UV" initializes registers R06-R07-R08

Data Registers: R00 = Geodesic Distance ( Registers R01 thru R08 are to be initialized before executing "DGD" )

    • R01 = u₁ • R03 = u₂ • R05 = N ( 6.N = number of steps )   • R06 = a • R07 = b • R08 = c
• R02 = v₁ • R04 = v₂

When the program stops, R26 = R27 = v'₁ R28 ~ 0 R29 = u R30 = v R31 = v' ( Check that R30 = v₂ )

Flags: /
Subroutines: /

-Lines 37 & 353 are three-byte GTOs

01 LBL "DGD"
02 RAD
03 STO 26
04 X<>Y
05 STO 27
06 RCL 03
07 RCL 01
08 -
09 RCL 05
10 6
11 *
12 /
13 STO 32
14 RCL 06
15 RCL 07
16 -
17 RCL 06
18 LASTX
19 +
20 *
21 STO 09
22 82
23 STO 51
24 216
25 STO 46
26 STO 50
27 27
28 STO 47
29 STO 49
30 272
31 STO 48
32 51.045
33 STO 25
34 RCL 27
35 XEQ 00
36 STO 28
37 GTO 04
38 LBL 00
39 STO 31
40 STO 34
41 RCL 05
42 6
43 *
44 STO 33
45 CLX
46 RCL 02
47 STO 30
48 RCL 01
49 STO 29
50 XEQ 02
51 RCL 15
52 SQRT
53 41
54 *
55 STO 00
56 DSE 25
57 GTO 03
58 LBL 01
59 RCL 31
60 STO 34
61 RCL 30
62 RCL 29
63 XEQ 02
64 RCL 15
65 SQRT
66 RCL IND 25
67 *
68 ST+ 00
69 DSE 25
70 GTO 03
71 CLX
72 6
73 ST+ 25
74 LBL 03
75 CLX
76 RCL 32
77 ST* 34
78 *
79 STO 39
80 3

81 /
82 RCL 31
83 +
84 STO 35
85 RCL 34
86 3
87 /
88 RCL 30
89 +
90 RCL 32
91 3
92 /
93 RCL 29
94 +
95 STO 38
96 XEQ 02
97 RCL 32
98 ST* 35
99 *
100 STO 40
101 1.5
102 /
103 RCL 31
104 +
105 STO 36
106 RCL 35
107 1.5
108 /
109 RCL 30
110 +
111 RCL 32
112 1.5
113 /
114 RCL 29
115 +
116 XEQ 02
117 RCL 32
118 ST* 36
119 *
120 STO 41
121 RCL 40
122 4
123 *
124 X<>Y
125 -
126 RCL 39
127 +
128 12
129 /
130 RCL 31
131 +
132 STO 37
133 RCL 35
134 4
135 *
136 RCL 36
137 -
138 RCL 34
139 +
140 12
141 /
142 RCL 30
143 +
144 RCL 38
145 XEQ 02
146 RCL 32
147 ST* 37
148 *
149 STO 42
150 18
151 *
152 RCL 41
153 7
154 *
155 +
156 RCL 40
157 22
158 *
159 -
160 RCL 39

161 5
162 *
163 +
164 9.6
165 /
166 RCL 31
167 +
168 STO 38
169 RCL 37
170 18
171 *
172 RCL 36
173 7
174 *
175 +
176 RCL 35
177 22
178 *
179 -
180 RCL 34
181 5
182 *
183 +
184 9.6
185 /
186 RCL 30
187 +
188 RCL 32
189 1.2
190 /
191 RCL 29
192 +
193 XEQ 02
194 RCL 32
195 ST* 38
196 *
197 STO 43
198 10
199 /
200 RCL 42
201 2
202 /
203 +
204 RCL 41
205 8
206 /
207 -
208 RCL 40
209 11
210 *
211 24
212 /
213 -
214 RCL 39
215 .15
216 *
217 +
218 RCL 31
219 +
220 STO 44
221 RCL 38
222 10
223 /
224 RCL 37
225 2
226 /
227 +
228 RCL 36
229 8
230 /
231 -
232 RCL 35
233 11
234 *
235 24
236 /
237 -
238 RCL 34
239 .15
240 *

241 +
242 RCL 30
243 +
244 RCL 32
245 6
246 /
247 RCL 29
248 +
249 XEQ 02
250 RCL 32
251 ST* 44
252 *
253 STO 45
254 80
255 *
256 RCL 40
257 99
258 *
259 +
260 RCL 42
261 118
262 *
263 -
264 20
265 *
266 RCL 43
267 ST+ 45
268 128
269 *
270 +
271 RCL 41
272 ST+ 42
273 215
274 *
275 +
276 RCL 39
277 783
278 *
279 -
280 780
281 /
282 RCL 31
283 +
284 RCL 35
285 99
286 *
287 RCL 44
288 80
289 *
290 +
291 RCL 37
292 118
293 *
294 -
295 20
296 *
297 RCL 38
298 ST+ 44
299 128
300 *
301 +
302 RCL 36
303 ST+ 37
304 215
305 *
306 +
307 RCL 34
308 783
309 *
310 -
311 780
312 /
313 RCL 30
314 +
315 RCL 29
316 RCL 32
317 +
318 XEQ 02
319 RCL 32
320 ST+ 29

321 ST* 12
322 *
323 RCL 39
324 +
325 13
326 *
327 RCL 45
328 32
329 ST* 44
330 *
331 +
332 RCL 42
333 55
334 ST* 37
335 *
336 +
337 .5
338 %
339 ST+ 31
340 RCL 34
341 RCL 12
342 +
343 13
344 *
345 RCL 44
346 +
347 RCL 37
348 +
349 .5
350 %
351 ST+ 30
352 DSE 33
353 GTO 01
354 RCL 31
355 RCL 30
356 RCL 29
357 XEQ 02
358 RCL 15
359 SQRT
360 41
361 *
362 ST+ 00
363 RCL 32
364 ABS
365 140
366 /
367 ST* 00
368 RCL 30
369 RCL 04
370 -
371 RTN
372 LBL 02
373 STO 10
374 X<> Z
375 STO 12
376 CLX
377 SIGN
378 P-R
379 STO 16
380 X<>Y
381 STO 17
382 RCL 10
383 LASTX
384 P-R
385 STO 13
386 X<>Y
387 STO 14
388 *
389 STO 11
390 *
391 *
392 RCL 09
393 *
394 STO 19
395 RCL 06
396 RCL 14
397 *
398 X^2
399 RCL 07
400 RCL 13

401 *
402 X^2
403 +
404 STO 15
405 RCL 16
406 X^2
407 *
408 STO 18
409 RCL 06
410 RCL 13
411 *
412 X^2
413 RCL 07
414 RCL 14
415 *
416 X^2
417 +
418 STO 10
419 RCL 17
420 X^2
421 *
422 RCL 08
423 RCL 16
424 *
425 X^2
426 +
427 STO 20
428 RCL 15
429 RCL 16
430 RCL 17
431 *
432 STO 24
433 ST+ X
434 *
435 CHS
436 STO 21
437 LASTX
438 RCL 10
439 RCL 08
440 X^2
441 -
442 *
443 STO 23
444 RCL 11
445 RCL 16
446 X^2
447 RCL 17
448 X^2
449 -
450 *
451 RCL 09
452 *
453 STO 22
454 RCL 12
455 *
456 RCL 13
457 X^2
458 RCL 14
459 X^2
460 -
461 RCL 24
462 *
463 RCL 09
464 *
465 +
466 STO 10
467 RCL 12
468 RCL 21
469 *
470 RCL 11
471 ST+ X
472 RCL 16
473 X^2
474 *
475 RCL 09
476 *
477 +
478 RCL 12
479 RCL 23

480 *
481 RCL 11
482 ST+ X
483 RCL 17
484 X^2
485 *
486 RCL 09
487 *
488 -
489 RCL 12
490 *
491 STO 24
492 RCL 10
493 ST+ 24
494 ST+ X
495 +
496 RCL 12
497 *
498 +
499 RCL 12
500 ST* 23
501 RCL 20
502 *
503 RCL 19
504 STO 11
505 +
506 ST+ 11
507 *
508 RCL 23
509 RCL 22
510 ST+ X
511 +
512 RCL 12
513 *
514 RCL 21
515 +
516 RCL 24
517 ST+ X
518 -
519 RCL 11
520 RCL 12
521 *
522 RCL 18
523 ST* 20
524 +
525 STO 15
526 *
527 +
528 RCL 20
529 RCL 19
530 X^2
531 -
532 ST+ X
533 /
534 RTN
535 LBL 04
536 VIEW 00
537 RCL 26
538 XEQ 00
539 ENTER
540 ENTER
541 X<> 28
542 -
543 X#0?
544 /
545 RCL 27
546 RCL 26
547 STO 27
548 STO T
549 -
550 *
551 +
552 STO 26
553 X#Y?
554 GTO 04
555 RCL 00
556 CLD
557 DEG
558 END

( 853 bytes / SIZE 052 )

STACK	INPUTS	OUTPUTS
Y	y	v'₁
X	x	Dist

Where x & y are 2 different initial guesses of v'₁

Example1: a = 7 b = 6 c = 5 , u₁ = 0 rd ; v₁ = 0.1 rd ; u₂ = 2 rd ; v₂ = 0.2 rd

7 STO 06 0 STO 01 2 STO 03
6 STO 07 0.1 STO 02 0.2 STO 04
5 STO 08

-If we choose N = 10 10 STO 05
-With initial guesses 0 & 1

0 ENTER^
1 XEQ "DGD" >>>> D = 12.91212647 ---Execution time = 42s--- with V41 in TURBO mode !
X<>Y v'₁ = 0.539381209

-We have R31 = v'₂ = -0.430961328

-With N = 20 , it yields the same distance: D = 12.91212647

Example2: On the Earth, with the values used in paragraph 0°)

Calculate the geodesic distance between the Sydney Observatory ( Long = 151°12'17"8 E , Lat = 33°51'41"1 S )
      and the Palomar Observatory ( Long = 116°51'50"4 W , Lat = 33°21'22"4 N )

     151.12178 ENTER^
    -33.51411   XEQ "LB-UV" >>>> 2.899588805 STO 01
   X<>Y -0.589438065 STO 02

   -116.51504 ENTER^
     33.21224   R/S >>>> -1.779101676 PI ENTER^ + + STO 03
X<>Y 0.580636810 STO 04

-With N = 10 STO 05 and initial guesses 0 1

0 ENTER^
1 XEQ "DGD" >>>> D = 12138.65673 km
X<>Y v'₁ = 0.436919794

and R31 = v'₂ = 0.451294935

-With N = 20 , it yields the same distance: D = 12138.65673 km

-With free42 and N = 100 or N = 200 , D = 12138.6567207 km

Example3: On the Earth, with the values used in paragraph 0°) L = 0° , b = 0° , L' = 179°51' , b = 0° ( nearly antipodal points )

     0 ENTER^
    0 XEQ "LB-UV" >>>> 0.260559389 STO 01
   X<>Y 0    STO 02

179.51 ENTER^
     0    R/S >>>> -2.883651233   PI + PI + STO 03
X<>Y 0    STO 04

-With N = 10 STO 05 and initial guesses    1 2

1 ENTER^
2 XEQ "DGD" >>>> D = 20001.79664 km
X<>Y v'₁ = 3.901657602

-With N = 20 , R00 remains around D = 20001.89623 km

-With free42 and N = 100 or N = 200 , D = 20001.8982845 km & v'₁ = 3.901001588

Notes:

-The successive D-approximations are displayed
-With a real HP41, execution time of the 1st example would be probably about 7 hours !
-We can modify this routine to reduce a little execution time if we replace line 553 X#Y? by;    - RND X#0? GTO 04 RCL 26
-The precision would be controlled by the display format.

-We can also use Newton-Cotes9 formula: replace

line 365 by 9 * 89600
line 360 by 2857
line 72 by 9
line 53 by 2857
line 42 by 9
lines 22 to 32 by 5714 STO 54 15741 STO 46 STO 53 1080 STO 47 STO 52 19344 STO 48 STO 51 5778 STO 49 STO 50 54.045
line 10 by 9

-Then, the number of steps is 9.N ( R05 = N )
-With N = 5, it yields almost the same geodesic distance D = 12.91212648 in the example1 above in 31 seconds with V41 ( turbo mode )
-So probably about 5h10m with a real HP41.

-Since digit entry is very slow, another way to make "DGD" slightly less slow is to store the numbers used by RK6 into data registers:

Data Registers: R00 = Geodesic Distance ( Registers R01 thru R08 are to be initialized before executing "DGD" )

• R01 = u₁ • R03 = u₂ • R05 = N ( 6.N = number of steps ) • R06 = a • R07 = b • R08 = c
• R02 = v₁ • R04 = v₂

When the program stops, R26 = R27 = v'₁ R28 ~ 0 R29 = u R30 = v R31 = v' ( Check that R30 = v₂ )

Flags: /
Subroutines: /

-Lines 117 & 413 are three-byte GTOs

01 LBL "DGD"
02 RAD
03 STO 26
04 X<>Y
05 STO 27
06 RCL 03
07 RCL 01
08 -
09 RCL 05
10 6
11 STO 65
12 *
13 /
14 STO 32
15 RCL 06
16 RCL 07
17 -
18 RCL 06
19 LASTX
20 +
21 *
22 STO 09
23 41
24 STO 76
25 ST+ X
26 STO 51
27 216
28 STO 46
29 STO 50
30 27
31 STO 47
32 STO 49
33 272
34 STO 48
35 51.045
36 STO 25
37 3
38 STO 52
39 2
40 STO 61
41 /
42 STO 53
43 12
44 STO 54
45 90
46 48
47 /
48 STO 55
49 35
50 LASTX
51 /
52 STO 56
53 110
54 LASTX
55 /
56 STO 57
57 25
58 LASTX
59 /
60 STO 58
61 40
62 LASTX
63 /
64 STO 59
65 10
66 STO 60
67 8
68 STO 62
69 11
70 24
71 /
72 STO 63
73 .15
74 STO 64
75 40
76 X^2
77 780
78 /
79 STO 66
80 128
81 LASTX
82 /
83 STO 67
84 2360
85 LASTX
86 /
87 STO 68
88 215

89 LASTX
90 /
91 STO 69
92 1980
93 LASTX
94 /
95 STO 70
96 783
97 LASTX
98 /
99 STO 71
100 13
101 200
102 /
103 STO 72
104 32
105 LASTX
106 /
107 STO 73
108 55
109 LASTX
110 /
111 STO 74
112 140
113 STO 75
114 RCL 27
115 XEQ 00
116 STO 28
117 GTO 04
118 LBL 00
119 STO 31
120 STO 34
121 RCL 05
122 RCL 65
123 *
124 STO 33
125 CLX
126 RCL 02
127 STO 30
128 RCL 01
129 STO 29
130 XEQ 02
131 RCL 15
132 SQRT
133 RCL 76
134 *
135 STO 00
136 DSE 25
137 GTO 03
138 LBL 01
139 RCL 31
140 STO 34
141 RCL 30
142 RCL 29
143 XEQ 02
144 RCL 15
145 SQRT
146 RCL IND 25
147 *
148 ST+ 00
149 DSE 25
150 GTO 03
151 CLX
152 RCL 65
153 ST+ 25
154 LBL 03
155 CLX
156 RCL 32
157 ST* 34
158 *
159 STO 39
160 RCL 52
161 /
162 RCL 31
163 +
164 STO 35
165 RCL 34
166 RCL 52
167 /
168 RCL 30
169 +
170 RCL 32
171 RCL 52
172 /
173 RCL 29
174 +
175 STO 38
176 XEQ 02

177 RCL 32
178 ST* 35
179 *
180 STO 40
181 RCL 53
182 /
183 RCL 31
184 +
185 STO 36
186 RCL 35
187 RCL 53
188 /
189 RCL 30
190 +
191 RCL 32
192 RCL 53
193 /
194 RCL 29
195 +
196 XEQ 02
197 RCL 32
198 ST* 36
199 *
200 STO 41
201 RCL 39
202 X<>Y
203 -
204 RCL 54
205 /
206 RCL 40
207 RCL 52
208 /
209 +
210 RCL 31
211 +
212 STO 37
213 RCL 34
214 RCL 36
215 -
216 RCL 54
217 /
218 RCL 35
219 RCL 52
220 /
221 +
222 RCL 30
223 +
224 RCL 38
225 XEQ 02
226 RCL 32
227 ST* 37
228 *
229 STO 42
230 RCL 55
231 *
232 RCL 41
233 RCL 56
234 *
235 +
236 RCL 40
237 RCL 57
238 *
239 -
240 RCL 39
241 RCL 58
242 *
243 +
244 RCL 31
245 +
246 STO 38
247 RCL 37
248 RCL 55
249 *
250 RCL 36
251 RCL 56
252 *
253 +
254 RCL 35
255 RCL 57
256 *
257 -
258 RCL 34
259 RCL 58
260 *
261 +
262 RCL 30
263 +
264 RCL 32

265 RCL 59
266 *
267 RCL 29
268 +
269 XEQ 02
270 RCL 32
271 ST* 38
272 *
273 STO 43
274 RCL 60
275 /
276 RCL 42
277 RCL 61
278 /
279 +
280 RCL 41
281 RCL 62
282 /
283 -
284 RCL 40
285 RCL 63
286 *
287 -
288 RCL 39
289 RCL 64
290 *
291 +
292 RCL 31
293 +
294 STO 44
295 RCL 38
296 RCL 60
297 /
298 RCL 37
299 RCL 61
300 /
301 +
302 RCL 36
303 RCL 62
304 /
305 -
306 RCL 35
307 RCL 63
308 *
309 -
310 RCL 34
311 RCL 64
312 *
313 +
314 RCL 30
315 +
316 RCL 32
317 RCL 65
318 /
319 RCL 29
320 +
321 XEQ 02
322 RCL 32
323 ST* 44
324 *
325 STO 45
326 RCL 66
327 *
328 RCL 40
329 RCL 70
330 *
331 +
332 RCL 42
333 RCL 68
334 *
335 -
336 RCL 43
337 ST+ 45
338 RCL 67
339 *
340 +
341 RCL 41
342 ST+ 42
343 RCL 69
344 *
345 +
346 RCL 39
347 RCL 71
348 *
349 -
350 RCL 31
351 +
352 RCL 35

353 RCL 70
354 *
355 RCL 44
356 RCL 66
357 *
358 +
359 RCL 37
360 RCL 68
361 *
362 -
363 RCL 38
364 ST+ 44
365 RCL 67
366 *
367 +
368 RCL 36
369 ST+ 37
370 RCL 69
371 *
372 +
373 RCL 34
374 RCL 71
375 *
376 -
377 RCL 30
378 +
379 RCL 29
380 RCL 32
381 +
382 XEQ 02
383 RCL 32
384 ST+ 29
385 ST* 12
386 *
387 RCL 39
388 +
389 RCL 72
390 *
391 RCL 45
392 RCL 73
393 ST* 44
394 *
395 +
396 RCL 42
397 RCL 74
398 ST* 37
399 *
400 +
401 ST+ 31
402 RCL 34
403 RCL 12
404 +
405 RCL 72
406 *
407 RCL 44
408 +
409 RCL 37
410 +
411 ST+ 30
412 DSE 33
413 GTO 01
414 RCL 31
415 RCL 30
416 RCL 29
417 XEQ 02
418 RCL 15
419 SQRT
420 RCL 76
421 *
422 ST+ 00
423 RCL 32
424 ABS
425 RCL 75
426 /
427 ST* 00
428 RCL 30
429 RCL 04
430 -
431 RTN
432 LBL 02
433 STO 10
434 X<> Z
435 STO 12
436 CLX
437 SIGN
438 P-R
439 STO 16
440 X<>Y

441 STO 17
442 RCL 10
443 LASTX
444 P-R
445 STO 13
446 X<>Y
447 STO 14
448 *
449 STO 11
450 *
451 *
452 RCL 09
453 *
454 STO 19
455 RCL 06
456 RCL 14
457 *
458 X^2
459 RCL 07
460 RCL 13
461 *
462 X^2
463 +
464 STO 15
465 RCL 16
466 X^2
467 *
468 STO 18
469 RCL 06
470 RCL 13
471 *
472 X^2
473 RCL 07
474 RCL 14
475 *
476 X^2
477 +
478 STO 10
479 RCL 17
480 X^2
481 *
482 RCL 08
483 RCL 16
484 *
485 X^2
486 +
487 STO 20
488 RCL 15
489 RCL 16
490 RCL 17
491 *
492 STO 24
493 ST+ X
494 *
495 CHS
496 STO 21
497 LASTX
498 RCL 10
499 RCL 08
500 X^2
501 -
502 *
503 STO 23
504 RCL 11
505 RCL 16
506 X^2
507 RCL 17
508 X^2
509 -
510 *
511 RCL 09
512 *
513 STO 22
514 RCL 12
515 *
516 RCL 13
517 X^2
518 RCL 14
519 X^2
520 -
521 RCL 24
522 *
523 RCL 09
524 *
525 +
526 STO 10
527 RCL 12
528 RCL 21
529 *

530 RCL 11
531 ST+ X
532 RCL 16
533 X^2
534 *
535 RCL 09
536 *
537 +
538 RCL 12
539 RCL 23
540 *
541 RCL 11
542 ST+ X
543 RCL 17
544 X^2
545 *
546 RCL 09
547 *
548 -
549 RCL 12
550 *
551 STO 24
552 RCL 10
553 ST+ 24
554 ST+ X
555 +
556 RCL 12
557 *
558 +
559 RCL 12
560 ST* 23
561 RCL 20
562 *
563 RCL 19
564 STO 11
565 +
566 ST+ 11
567 *
568 RCL 23
569 RCL 22
570 ST+ X
571 +
572 RCL 12
573 *
574 RCL 21
575 +
576 RCL 24
577 ST+ X
578 -
579 RCL 11
580 RCL 12
581 *
582 RCL 18
583 ST* 20
584 +
585 STO 15
586 *
587 +
588 RCL 20
589 RCL 19
590 X^2
591 -
592 ST+ X
593 /
594 RTN
595 LBL 04
596 VIEW 00
597 RCL 26
598 XEQ 00
599 ENTER
600 ENTER
601 X<> 28
602 -
603 X#0?
604 /
605 RCL 27
606 RCL 26
607 STO 27
608 STO T
609 -
610 *
611 +
612 STO 26
613 X#Y?
614 GTO 04
615 RCL 00
616 CLD
617 DEG
618 END

( 964 bytes / SIZE 077 )

STACK	INPUTS	OUTPUTS
Y	y	v'₁
X	x	Dist

01 LBL "DGD"
02 RAD
03 STO 26
04 X<>Y
05 STO 27
06 RCL 03
07 RCL 01
08 -
09 RCL 05
10 12
11 *
12 /
13 STO 32
14 RCL 06
15 RCL 07
16 -
17 RCL 06
18 LASTX
19 +
20 *
21 STO 09
22 82
23 STO 43
24 216
25 STO 38
26 STO 42
27 27
28 STO 39
29 STO 41
30 272
31 STO 40
32 43.037
33 STO 25
34 RCL 27
35 XEQ 00
36 STO 28
37 GTO 04
38 LBL 00
39 STO 31
40 RCL 05
41 6
42 *
43 STO 33
44 CLX
45 RCL 02
46 STO 30
47 RCL 01
48 STO 29
49 XEQ 02

50 RCL 15
51 SQRT
52 41
53 *
54 STO 00
55 DSE 25
56 GTO 03
57 LBL 01
58 RCL 31
59 RCL 30
60 RCL 29
61 XEQ 02
62 RCL 15
63 SQRT
64 RCL IND 25
65 *
66 ST+ 00
67 DSE 25
68 GTO 03
69 CLX
70 6
71 ST+ 25
72 LBL 03
73 CLX
74 RCL 32
75 ST+ 29
76 *
77 STO 35
78 RCL 31
79 +
80 STO 37
81 LASTX
82 RCL 32
83 *
84 STO 34
85 RCL 30
86 +
87 RCL 29
88 XEQ 02
89 RCL 32
90 ST* 37
91 *
92 ST+ 35
93 ST+ 35
94 RCL 31
95 +
96 ENTER
97 X<> 37
98 ST+ 34

99 ST+ 34
100 RCL 30
101 +
102 RCL 29
103 XEQ 02
104 RCL 32
105 ST+ 29
106 ST+ X
107 ST* 37
108 *
109 ST+ 35
110 RCL 31
111 +
112 ENTER
113 X<> 37
114 ST+ 34
115 RCL 30
116 +
117 RCL 29
118 XEQ 02
119 RCL 32
120 ST* 37
121 *
122 RCL 35
123 +
124 3
125 /
126 ST+ 31
127 RCL 37
128 RCL 34
129 +
130 3
131 /
132 ST+ 30
133 DSE 33
134 GTO 01
135 RCL 31
136 RCL 30
137 RCL 29
138 XEQ 02
139 RCL 15
140 SQRT
141 41
142 *
143 ST+ 00
144 RCL 32
145 ABS
146 70
147 /

148 ST* 00
149 RCL 30
150 RCL 04
151 -
152 RTN
153 LBL 02
154 STO 10
155 X<> Z
156 STO 12
157 CLX
158 SIGN
159 P-R
160 STO 16
161 X<>Y
162 STO 17
163 RCL 10
164 LASTX
165 P-R
166 STO 13
167 X<>Y
168 STO 14
169 *
170 STO 11
171 *
172 *
173 RCL 09
174 *
175 STO 19
176 RCL 06
177 RCL 14
178 *
179 X^2
180 RCL 07
181 RCL 13
182 *
183 X^2
184 +
185 STO 15
186 RCL 16
187 X^2
188 *
189 STO 18
190 RCL 06
191 RCL 13
192 *
193 X^2
194 RCL 07
195 RCL 14
196 *

197 X^2
198 +
199 STO 10
200 RCL 17
201 X^2
202 *
203 RCL 08
204 RCL 16
205 *
206 X^2
207 +
208 STO 20
209 RCL 15
210 RCL 16
211 RCL 17
212 *
213 STO 24
214 ST+ X
215 *
216 CHS
217 STO 21
218 LASTX
219 RCL 10
220 RCL 08
221 X^2
222 -
223 *
224 STO 23
225 RCL 11
226 RCL 16
227 X^2
228 RCL 17
229 X^2
230 -
231 *
232 RCL 09
233 *
234 STO 22
235 RCL 12
236 *
237 RCL 13
238 X^2
239 RCL 14
240 X^2
241 -
242 RCL 24
243 *
244 RCL 09
245 *

246 +
247 STO 10
248 RCL 12
249 RCL 21
250 *
251 RCL 11
252 ST+ X
253 RCL 16
254 X^2
255 *
256 RCL 09
257 *
258 +
259 RCL 12
260 RCL 23
261 *
262 RCL 11
263 ST+ X
264 RCL 17
265 X^2
266 *
267 RCL 09
268 *
269 -
270 RCL 12
271 *
272 STO 24
273 RCL 10
274 ST+ 24
275 ST+ X
276 +
277 RCL 12
278 *
279 +
280 RCL 12
281 ST* 23
282 RCL 20
283 *
284 RCL 19
285 STO 11
286 +
287 ST+ 11
288 *
289 RCL 23
290 RCL 22
291 ST+ X
292 +
293 RCL 12
294 *

295 RCL 21
296 +
297 RCL 24
298 ST+ X
299 -
300 RCL 11
301 RCL 12
302 *
303 RCL 18
304 ST* 20
305 +
306 STO 15
307 *
308 +
309 RCL 20
310 RCL 19
311 X^2
312 -
313 ST+ X
314 /
315 RTN
316 LBL 04
317 VIEW 00
318 RCL 26
319 XEQ 00
320 ENTER
321 ENTER
322 X<> 28
323 -
324 X#0?
325 /
326 RCL 27
327 RCL 26
328 STO 27
329 STO T
330 -
331 *
332 +
333 STO 26
334 -
335 RND
336 X#0?
337 GTO 04
338 RCL 26
339 RCL 00
340 CLD
341 DEG
342 END

( 508 bytes / SIZE 044 )

STACK	INPUTS	OUTPUTS
Y	y	v'₁
X	x	Dist

Where x & y are 2 different initial guesses of v'₁

Example: a = 7 b = 6 c = 5 , u₁ = 0 rd ; v₁ = 0.1 rd ; u₂ = 2 rd ; v₂ = 0.2 rd

7 STO 06 0 STO 01 2 STO 03
6 STO 07 0.1 STO 02 0.2 STO 04
5 STO 08

-If we choose N = 10 10 STO 05
-With initial guesses 0 & 1

FIX 9

0 ENTER^
1 XEQ "DGD" >>>> D = 12.91212648 ---Execution time = 17s--- with V41 in TURBO mode
X<>Y v'₁ = 0.539381183

-We have R31 = v'₂ = -0.430961303

Notes:

-Though v' are less accurate, the precision of the distance is almost the same !
-Execution time is 17 seconds instead of 42 seconds with V41 in turbo mode.
-With a real HP41, it would be about 2h50mn instead of 7 hours...

-We can also use Newton-Cotes9 formula: replace

line 146 by 9 * 44800
line 141 by 2857
line 70 by 9
line 52 by 2857
line 41 by 9
lines 22 to 32 by 5714 STO 46 15741 STO 38 STO 45 1080 STO 39 STO 44 19344 STO 40 STO 43 5778 STO 41 STO 42 46.037
line 10 by 18

-Then, the number of steps is 9.N ( R05 = N )
-With N = 5, it yields the same geodesic distance in the example above in about 12.8 seconds with V41 ( turbo mode )
-So probably about 2h09m with a real HP41.

2°) U = U(V)

     a) Runge-Kutta of order 6

-If u₁ = u₂ we have to express u as a function of v. The formulae become:

   § L dv where    L = SQRT ( R + 2.Q u' + P u'² ) ( u' = du/dv )

u'' ( P.R - Q² ) / ( R + 2.Q u' + P u'² ) = - dQ/dv - u' dP/dv + (1/2) ( ¶R/¶u + 2.u' ¶Q/¶u + u'² ¶P/¶u ) + (1/2) ( ( Q + P.u' ) / ( R + 2.Q u' + P u'² ) ) ( dR/dv + 2.u' dQ/dv + u'² dP/du )

( dP/dv = ¶P/¶v + u' ¶P/¶u ... etc ... )

Data Registers: R00 = +/-Geodesic Distance ( Registers R01 thru R08 are to be initialized before executing "DGD" )

    • R01 = u₁ • R03 = u₂ • R05 = N ( 6.N = number of steps )   • R06 = a • R07 = b • R08 = c
• R02 = v₁ • R04 = v₂

When the program stops, R26 = R27 = u'₁ R28 ~ 0 R29 = v R30 = u R31 = u'    ( Check that R30 = u₂ )

Flags: /
Subroutines: /

-Lines 37 & 353 are three-byte GTOs

01 LBL "DGD"
02 RAD
03 STO 26
04 X<>Y
05 STO 27
06 RCL 04
07 RCL 02
08 -
09 RCL 05
10 6
11 *
12 /
13 STO 32
14 RCL 06
15 RCL 07
16 -
17 RCL 06
18 LASTX
19 +
20 *
21 STO 09
22 82
23 STO 51
24 216
25 STO 46
26 STO 50
27 27
28 STO 47
29 STO 49
30 272
31 STO 48
32 51.045
33 STO 25
34 RCL 27
35 XEQ 00
36 STO 28
37 GTO 04
38 LBL 00
39 STO 31
40 STO 34
41 RCL 05
42 6
43 *
44 STO 33
45 CLX
46 RCL 01
47 STO 30
48 RCL 02
49 STO 29
50 XEQ 02
51 RCL 15
52 SQRT
53 41
54 *
55 STO 00
56 DSE 25
57 GTO 03
58 LBL 01
59 RCL 31
60 STO 34
61 RCL 30
62 RCL 29
63 XEQ 02
64 RCL 15
65 SQRT
66 RCL IND 25
67 *
68 ST+ 00
69 DSE 25
70 GTO 03
71 CLX
72 6
73 ST+ 25
74 LBL 03
75 CLX
76 RCL 32
77 ST* 34
78 *
79 STO 39
80 3

321 ST* 12
322 *
323 RCL 39
324 +
325 13
326 *
327 RCL 45
328 32
329 ST* 44
330 *
331 +
332 RCL 42
333 55
334 ST* 37
335 *
336 +
337 .5
338 %
339 ST+ 31
340 RCL 34
341 RCL 12
342 +
343 13
344 *
345 RCL 44
346 +
347 RCL 37
348 +
349 .5
350 %
351 ST+ 30
352 DSE 33
353 GTO 01
354 RCL 31
355 RCL 30
356 RCL 29
357 XEQ 02
358 RCL 15
359 SQRT
360 41
361 *
362 ST+ 00
363 RCL 32
364 ABS
365 140
366 /
367 ST* 00
368 RCL 30
369 RCL 03
370 -
371 RTN
372 LBL 02
373 STO 10
374 X<> Z
375 STO 12
376 CLX
377 SIGN
378 P-R
379 STO 13
380 X<>Y
381 STO 14
382 RCL 10
383 LASTX
384 P-R
385 STO 16
386 X<>Y
387 STO 17
388 *
389 STO 24
390 *
391 *
392 RCL 09
393 *
394 STO 19
395 RCL 06
396 RCL 14
397 *
398 X^2
399 RCL 07

400 RCL 13
401 *
402 X^2
403 +
404 STO 15
405 RCL 16
406 X^2
407 *
408 STO 18
409 RCL 06
410 RCL 13
411 *
412 X^2
413 RCL 07
414 RCL 14
415 *
416 X^2
417 +
418 STO 10
419 RCL 17
420 X^2
421 *
422 RCL 08
423 RCL 16
424 *
425 X^2
426 +
427 STO 20
428 RCL 13
429 RCL 14
430 *
431 STO 11
432 RCL 17
433 X^2
434 *
435 RCL 09
436 ST+ X
437 *
438 CHS
439 STO 21
440 RCL 13
441 X^2
442 RCL 14
443 X^2
444 -
445 RCL 09
446 *
447 RCL 24
448 *
449 STO 22
450 RCL 11
451 RCL 16
452 X^2
453 *
454 ST+ X
455 RCL 09
456 *
457 RCL 12
458 *
459 STO 23
460 RCL 15
461 RCL 24
462 *
463 ST+ X
464 -
465 X<> 16
466 X^2
467 RCL 17
468 X^2
469 -
470 RCL 09
471 *
472 RCL 11
473 *
474 RCL 12
475 RCL 22
476 *
477 +
478 X<> 24

479 ST+ X
480 RCL 10
481 RCL 08
482 X^2
483 -
484 *
485 RCL 12
486 ST* 16
487 RCL 21
488 *
489 +
490 RCL 16
491 RCL 24
492 ST+ X
493 +
494 RCL 12
495 *
496 +
497 RCL 12
498 RCL 18
499 *
500 RCL 19
501 STO 11
502 +
503 ST+ 11
504 *
505 RCL 23
506 RCL 22
507 ST+ X
508 +
509 RCL 12
510 *
511 RCL 21
512 +
513 RCL 16
514 RCL 24
515 +
516 ST+ X
517 -
518 RCL 11
519 RCL 12
520 *
521 RCL 20
522 ST* 18
523 +
524 STO 15
525 *
526 +
527 RCL 18
528 RCL 19
529 X^2
530 -
531 ST+ X
532 /
533 RTN
534 LBL 04
535 VIEW 00
536 RCL 26
537 XEQ 00
538 ENTER
539 ENTER
540 X<> 28
541 -
542 X#0?
543 /
544 RCL 27
545 RCL 26
546 STO 27
547 STO T
548 -
549 *
550 +
551 STO 26
552 X#Y?
553 GTO 04
554 RCL 00
555 CLD
556 DEG
557 END

( 854 bytes / SIZE 052 )

STACK	INPUTS	OUTPUTS
Y	y	u'₁
X	x	Dist

Where x & y are 2 different initial guesses of u'₁

Example1: a = 7 b = 6 c = 5 , u₁ = 1 rd ; v₁ = 0 rd ; u₂ = 1 rd ; v₂ = 1 rd

7 STO 06 1 STO 01 1 STO 03
6 STO 07 0 STO 02 1 STO 04
5 STO 08

-If we choose N = 10 10 STO 05
-With initial guesses 0 & 0.1

0    ENTER^
0.1 XEQ "DGD" >>>> D = 5.375896782
   X<>Y u'₁ = 0.059732533

-And u'₂ = R31 = -0.087357012

Example2: On the Earth, with the values used in paragraph 0°) L = 0° , b = 0° , L' = 0° , b = 41°

     0 ENTER^
    0 XEQ "LB-UV" >>>> 0.260559389 STO 01
   X<>Y 0    STO 02

0 ENTER^
     41 R/S >>>>    0.260559389 STO 03
X<>Y 0.713920143 STO 04

-With N = 10 STO 05 and initial guesses    0 0.1

0    ENTER^
0.1 XEQ "DGD" >>>> D = 4540.561276 km
   X<>Y u'₁ = 0.000001952

and R31 = u'₂ = -0.000002366

-With free42 and N = 100 , it yields D = 4540.56127748 & u'₁ = 0.00000195128 , u'₂ = -0.00000236680

Note:

-Like with the second program listed in paragraph 1°) we can store the numbers in data registers to make the routine a little faster:

Data Registers: R00 = +/-Geodesic Distance ( Registers R01 thru R08 are to be initialized before executing "DGD" )

    • R01 = u₁ • R03 = u₂ • R05 = N ( 6.N = number of steps )   • R06 = a • R07 = b • R08 = c
• R02 = v₁ • R04 = v₂

When the program stops, R26 = R27 = u'₁ R28 ~ 0 R29 = v R30 = u R31 = u' ( Check that R30 = u₂ )

Flags: /
Subroutines: /

-Lines 117 & 413 are three-byte GTOs

01 LBL "DGD"
02 RAD
03 STO 26
04 X<>Y
05 STO 27
06 RCL 04
07 RCL 02
08 -
09 RCL 05
10 6
11 STO 65
12 *
13 /
14 STO 32
15 RCL 06
16 RCL 07
17 -
18 RCL 06
19 LASTX
20 +
21 *
22 STO 09
23 41
24 STO 76
25 ST+ X
26 STO 51
27 216
28 STO 46
29 STO 50
30 27
31 STO 47
32 STO 49
33 272
34 STO 48
35 51.045
36 STO 25
37 3
38 STO 52
39 2
40 STO 61
41 /
42 STO 53
43 12
44 STO 54
45 90
46 48
47 /
48 STO 55
49 35
50 LASTX
51 /
52 STO 56
53 110
54 LASTX
55 /
56 STO 57
57 25
58 LASTX
59 /
60 STO 58
61 40
62 LASTX
63 /
64 STO 59
65 10
66 STO 60
67 8
68 STO 62
69 11
70 24
71 /
72 STO 63
73 .15
74 STO 64
75 40
76 X^2
77 780
78 /
79 STO 66
80 128
81 LASTX
82 /
83 STO 67
84 2360
85 LASTX
86 /
87 STO 68
88 215

89 LASTX
90 /
91 STO 69
92 1980
93 LASTX
94 /
95 STO 70
96 783
97 LASTX
98 /
99 STO 71
100 13
101 200
102 /
103 STO 72
104 32
105 LASTX
106 /
107 STO 73
108 55
109 LASTX
110 /
111 STO 74
112 140
113 STO 75
114 RCL 27
115 XEQ 00
116 STO 28
117 GTO 04
118 LBL 00
119 STO 31
120 STO 34
121 RCL 05
122 RCL 65
123 *
124 STO 33
125 CLX
126 RCL 01
127 STO 30
128 RCL 02
129 STO 29
130 XEQ 02
131 RCL 15
132 SQRT
133 RCL 76
134 *
135 STO 00
136 DSE 25
137 GTO 03
138 LBL 01
139 RCL 31
140 STO 34
141 RCL 30
142 RCL 29
143 XEQ 02
144 RCL 15
145 SQRT
146 RCL IND 25
147 *
148 ST+ 00
149 DSE 25
150 GTO 03
151 CLX
152 RCL 65
153 ST+ 25
154 LBL 03
155 CLX
156 RCL 32
157 ST* 34
158 *
159 STO 39
160 RCL 52
161 /
162 RCL 31
163 +
164 STO 35
165 RCL 34
166 RCL 52
167 /
168 RCL 30
169 +
170 RCL 32
171 RCL 52
172 /
173 RCL 29
174 +
175 STO 38
176 XEQ 02

353 RCL 70
354 *
355 RCL 44
356 RCL 66
357 *
358 +
359 RCL 37
360 RCL 68
361 *
362 -
363 RCL 38
364 ST+ 44
365 RCL 67
366 *
367 +
368 RCL 36
369 ST+ 37
370 RCL 69
371 *
372 +
373 RCL 34
374 RCL 71
375 *
376 -
377 RCL 30
378 +
379 RCL 29
380 RCL 32
381 +
382 XEQ 02
383 RCL 32
384 ST+ 29
385 ST* 12
386 *
387 RCL 39
388 +
389 RCL 72
390 *
391 RCL 45
392 RCL 73
393 ST* 44
394 *
395 +
396 RCL 42
397 RCL 74
398 ST* 37
399 *
400 +
401 ST+ 31
402 RCL 34
403 RCL 12
404 +
405 RCL 72
406 *
407 RCL 44
408 +
409 RCL 37
410 +
411 ST+ 30
412 DSE 33
413 GTO 01
414 RCL 31
415 RCL 30
416 RCL 29
417 XEQ 02
418 RCL 15
419 SQRT
420 RCL 76
421 *
422 ST+ 00
423 RCL 32
424 ABS
425 RCL 75
426 /
427 ST* 00
428 RCL 30
429 RCL 03
430 -
431 RTN
432 LBL 02
433 STO 10
434 X<> Z
435 STO 12
436 CLX
437 SIGN
438 P-R
439 STO 13
440 X<>Y

441 STO 14
442 RCL 10
443 LASTX
444 P-R
445 STO 16
446 X<>Y
447 STO 17
448 *
449 STO 24
450 *
451 *
452 RCL 09
453 *
454 STO 19
455 RCL 06
456 RCL 14
457 *
458 X^2
459 RCL 07
460 RCL 13
461 *
462 X^2
463 +
464 STO 15
465 RCL 16
466 X^2
467 *
468 STO 18
469 RCL 06
470 RCL 13
471 *
472 X^2
473 RCL 07
474 RCL 14
475 *
476 X^2
477 +
478 STO 10
479 RCL 17
480 X^2
481 *
482 RCL 08
483 RCL 16
484 *
485 X^2
486 +
487 STO 20
488 RCL 13
489 RCL 14
490 *
491 STO 11
492 RCL 17
493 X^2
494 *
495 RCL 09
496 ST+ X
497 *
498 CHS
499 STO 21
500 RCL 13
501 X^2
502 RCL 14
503 X^2
504 -
505 RCL 09
506 *
507 RCL 24
508 *
509 STO 22
510 RCL 11
511 RCL 16
512 X^2
513 *
514 ST+ X
515 RCL 09
516 *
517 RCL 12
518 *
519 STO 23
520 RCL 15
521 RCL 24
522 *
523 ST+ X
524 -
525 X<> 16
526 X^2
527 RCL 17
528 X^2

529 -
530 RCL 09
531 *
532 RCL 11
533 *
534 RCL 12
535 RCL 22
536 *
537 +
538 X<> 24
539 ST+ X
540 RCL 10
541 RCL 08
542 X^2
543 -
544 *
545 RCL 12
546 ST* 16
547 RCL 21
548 *
549 +
550 RCL 16
551 RCL 24
552 ST+ X
553 +
554 RCL 12
555 *
556 +
557 RCL 12
558 RCL 18
559 *
560 RCL 19
561 STO 11
562 +
563 ST+ 11
564 *
565 RCL 23
566 RCL 22
567 ST+ X
568 +
569 RCL 12
570 *
571 RCL 21
572 +
573 RCL 16
574 RCL 24
575 +
576 ST+ X
577 -
578 RCL 11
579 RCL 12
580 *
581 RCL 20
582 ST* 18
583 +
584 STO 15
585 *
586 +
587 RCL 18
588 RCL 19
589 X^2
590 -
591 ST+ X
592 /
593 RTN
594 LBL 04
595 VIEW 00
596 RCL 26
597 XEQ 00
598 ENTER
599 ENTER
600 X<> 28
601 -
602 X#0?
603 /
604 RCL 27
605 RCL 26
606 STO 27
607 STO T
608 -
609 *
610 +
611 STO 26
612 X#Y?
613 GTO 04
614 RCL 00
615 CLD
616 DEG
617 END

( 965 bytes / SIZE 077 )

STACK	INPUTS	OUTPUTS
Y	y	u'₁
X	x	Dist

Where x & y are 2 different initial guesses of u'₁

Example: a = 7 b = 6 c = 5 , u₁ = 1 rd ; v₁ = 0 rd ; u₂ = 1 rd ; v₂ = 1 rd

7 STO 06 1 STO 01 1 STO 03
6 STO 07 0 STO 02 1 STO 04
5 STO 08

-If we choose N = 10 10 STO 05
-With initial guesses 0 & 0.1

0    ENTER^
0.1 XEQ "DGD" >>>> D = 5.375896782
   X<>Y u'₁ = 0.059732533

-And u'₂ = R31 = -0.087357012

     b) Runge-Kutta of order 4

-We can save many bytes with 4th-order Runge-Kutta formula.
-And with a real HP41, the program will run faster.

Data Registers: R00 = Geodesic Distance ( Registers R01 thru R08 are to be initialized before executing "DGD" )

    • R01 = u₁ • R03 = u₂ • R05 = N ( 6.N = number of steps )   • R06 = a • R07 = b • R08 = c
• R02 = v₁ • R04 = v₂

When the program stops, R26 = R27 = u'₁ R28 ~ 0 R29 = v R30 = u R31 = u' ( Check that R30 = u₂ )

Flags: /
Subroutines: /

-Lines 37 & 134 are three-byte GTOs

01 LBL "DGD"
02 RAD
03 STO 26
04 X<>Y
05 STO 27
06 RCL 04
07 RCL 02
08 -
09 RCL 05
10 12
11 *
12 /
13 STO 32
14 RCL 06
15 RCL 07
16 -
17 RCL 06
18 LASTX
19 +
20 *
21 STO 09
22 82
23 STO 43
24 216
25 STO 38
26 STO 42
27 27
28 STO 39
29 STO 41
30 272
31 STO 40
32 43.037
33 STO 25
34 RCL 27
35 XEQ 00
36 STO 28
37 GTO 04
38 LBL 00
39 STO 31
40 RCL 05
41 6
42 *
43 STO 33
44 CLX
45 RCL 01
46 STO 30
47 RCL 02
48 STO 29
49 XEQ 02

148 ST* 00
149 RCL 30
150 RCL 03
151 -
152 RTN
153 LBL 02
154 STO 10
155 X<> Z
156 STO 12
157 CLX
158 SIGN
159 P-R
160 STO 13
161 X<>Y
162 STO 14
163 RCL 10
164 LASTX
165 P-R
166 STO 16
167 X<>Y
168 STO 17
169 *
170 STO 24
171 *
172 *
173 RCL 09
174 *
175 STO 19
176 RCL 06
177 RCL 14
178 *
179 X^2
180 RCL 07
181 RCL 13
182 *
183 X^2
184 +
185 STO 15
186 RCL 16
187 X^2
188 *
189 STO 18
190 RCL 06
191 RCL 13
192 *
193 X^2
194 RCL 07
195 RCL 14
196 *

197 X^2
198 +
199 STO 10
200 RCL 17
201 X^2
202 *
203 RCL 08
204 RCL 16
205 *
206 X^2
207 +
208 STO 20
209 RCL 13
210 RCL 14
211 *
212 STO 11
213 RCL 17
214 X^2
215 *
216 RCL 09
217 ST+ X
218 *
219 CHS
220 STO 21
221 RCL 13
222 X^2
223 RCL 14
224 X^2
225 -
226 RCL 09
227 *
228 RCL 24
229 *
230 STO 22
231 RCL 11
232 RCL 16
233 X^2
234 *
235 ST+ X
236 RCL 09
237 *
238 RCL 12
239 *
240 STO 23
241 RCL 15
242 RCL 24
243 *
244 ST+ X

245 -
246 X<> 16
247 X^2
248 RCL 17
249 X^2
250 -
251 RCL 09
252 *
253 RCL 11
254 *
255 RCL 12
256 RCL 22
257 *
258 +
259 X<> 24
260 ST+ X
261 RCL 10
262 RCL 08
263 X^2
264 -
265 *
266 RCL 12
267 ST* 16
268 RCL 21
269 *
270 +
271 RCL 16
272 RCL 24
273 ST+ 16
274 ST+ X
275 +
276 RCL 12
277 *
278 +
279 RCL 12
280 RCL 18
281 *
282 RCL 19
283 STO 11
284 +
285 ST+ 11
286 *
287 RCL 23
288 RCL 22
289 ST+ X
290 +
291 RCL 12
292 *

293 RCL 21
294 +
295 RCL 16
296 ST+ X
297 -
298 RCL 11
299 RCL 12
300 *
301 RCL 20
302 ST* 18
303 +
304 STO 15
305 *
306 +
307 RCL 18
308 RCL 19
309 X^2
310 -
311 ST+ X
312 /
313 RTN
314 LBL 04
315 VIEW 00
316 RCL 26
317 XEQ 00
318 ENTER
319 ENTER
320 X<> 28
321 -
322 X#0?
323 /
324 RCL 27
325 RCL 26
326 STO 27
327 STO T
328 -
329 *
330 +
331 STO 26
332 -
333 RND
334 X#0?
335 GTO 04
336 RCL 26
337 RCL 00
338 CLD
339 DEG
340 END

( 508 bytes / SIZE 044 )

STACK	INPUTS	OUTPUTS
Y	y	u'₁
X	x	Dist

hp41programs

Geodesics on a Triaxial Ellipsoid(2) for the HP-41