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'2 )       ( v' = dv/du )

P = Cos2 v  ( a2 Sin2 u + b2 Cos2 u )
Q = (1/4) ( a2 - b2 ) Sin 2.u  Sin 2.v
R = Sin2 v  ( a2 Cos2 u + b2 Sin2 u ) + c2 Cos2 v

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

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

( dP/du =   P/uv' P/v   ... etc ... )

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

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

NC6:    §xx+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:    §xx+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'1  such that  v(u2) = v2

-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 = u1     •  R03 = u2     •  R05 = N  ( 6.N = number of steps )                   •  R06 = a    •  R07 = b    •  R08 = c
•  R02 = v1     •  R04 = v2

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

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'1 X x Dist

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

Example1:
a = 7   b = 6   c = 5   ,    u1 = 0 rd  ;  v1 = 0.1 rd  ;   u2 = 2 rd  ;  v2 = 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'1 =  0.539381209

-We have  R31 = v'2 = -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'1 =  0.436919794

and   R31 = v'2 = 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'1 = 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'1 =  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 = u1     •  R03 = u2     •  R05 = N  ( 6.N = number of steps )                   •  R06 = a    •  R07 = b    •  R08 = c
•  R02 = v1     •  R04 = v2

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

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'1 X x Dist

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

Example:
a = 7   b = 6   c = 5   ,    u1 = 0 rd  ;  v1 = 0.1 rd  ;   u2 = 2 rd  ;  v2 = 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 = 39s---      with V41 in TURBO mode
X<>Y   v'1 =  0.539381209

Notes:

-Thus, we have saved 3 seconds with V41 in turbo mode
-With a real HP41, probably about 30 minutes...

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 = u1     •  R03 = u2     •  R05 = N  ( 6.N = number of steps )                   •  R06 = a    •  R07 = b    •  R08 = c
•  R02 = v1     •  R04 = v2

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

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 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'1 X x Dist

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

Example:
a = 7   b = 6   c = 5   ,    u1 = 0 rd  ;  v1 = 0.1 rd  ;   u2 = 2 rd  ;  v2 = 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'1 =  0.539381183

-We have  R31 = v'2 = -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  u1 = u2  we have to express u as a function of v. The formulae become:

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

u''  ( P.R - Q2 ) / ( R + 2.Q u' + P u'2 )  =  - dQ/dv - u' dP/dv + (1/2) ( R/u + 2.u' Q/u + u'2 P/u ) + (1/2)  ( ( Q + P.u' ) / ( R + 2.Q u' + P u'2 ) ) ( dR/dv + 2.u' dQ/dv + u'2 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 = u1     •  R03 = u2     •  R05 = N  ( 6.N = number of steps )                   •  R06 = a    •  R07 = b    •  R08 = c
•  R02 = v1     •  R04 = v2

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

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 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 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'1 X x Dist

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

Example1:
a = 7   b = 6   c = 5   ,    u1 = 1 rd  ;  v1 = 0 rd  ;   u2 = 1 rd  ;  v2 = 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'1 = 0.059732533

-And   u'2 = 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'1 = 0.000001952

and  R31 = u'2 = -0.000002366

-With free42 and  N = 100 , it yields    D = 4540.56127748   &   u'1 = 0.00000195128  ,   u'2 = -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 = u1     •  R03 = u2     •  R05 = N  ( 6.N = number of steps )                   •  R06 = a    •  R07 = b    •  R08 = c
•  R02 = v1     •  R04 = v2

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

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 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 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'1 X x Dist

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

Example:
a = 7   b = 6   c = 5   ,    u1 = 1 rd  ;  v1 = 0 rd  ;   u2 = 1 rd  ;  v2 = 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'1 = 0.059732533

-And   u'2 = 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 = u1     •  R03 = u2     •  R05 = N  ( 6.N = number of steps )                   •  R06 = a    •  R07 = b    •  R08 = c
•  R02 = v1     •  R04 = v2

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

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 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 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'1 X x Dist

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

Example:
a = 7   b = 6   c = 5   ,    u1 = 1 rd  ;  v1 = 0 rd  ;   u2 = 1 rd  ;  v2 = 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

FIX 9

0    ENTER^
0.1   XEQ "DGD"   >>>>    D  = 5.375896774
X<>Y   u'1 = 0.059732531

-And   u'2 = R31 = -0.087357015

Remark:

-We can make similar modifications to use NC9 instead of NC6 for the integrals
and store a few numbers in data registers to slightly reduce execution time...