Geod-3axial

# Geodesics on a Triaxial Ellipsoid for the HP-41

Overview

1°)  Jacobi's Method

a)  Program#1
b)  Program#2
c)  Program(s)#3

2°)  Longitude & Latitude >>> Cartesian Coordinates

a)  Geodetic Coordinates
b)  Astronomical Coordinates

3°)  Cartesian Coordinates <> Jacobi Parameters ( u , v )
4°)  Approximate Methods

a)  Program#1
b)  Program#2
c)  Program#3
d)  Program#4
e)  For the Earth only ( 2 programs )

Latest Update:    A second program in §4e)

-The routines listed in paragraph 1 employ the rigorous method found by Jacobi in 1838.
-In paragraph 4, the results are only approximate but still very good for the Earth:
-With the last 3 programs, the accuracy is of the order of 20 cm, except for nearly antipodal points.

1°)  Jacobi's Method

a)  Program#1

-The geodesic distance between 2 points on a scalene - or triaxial - ellipsoid is much more difficult to calculate
than the geodesic distance on an ellpsoid of revolution.
-The program hereunder uses the formulae given by Jacobi in reference [1]

-Given an ellipsoid    x2 / a + y2 / b + z2 / c = 1  with  a > b > c   ( here, a , b , c  are the squares of the semi-axes ),

and 2 points M ( u1 , v1 ) & N ( u2 , v2 )  where  x y z  and  u v are related by:

x =  [ a / ( a - c ) ]1/2 Cos u  ( a - b Sin2 v - c Cos2 v )1/2
y =    b1/2  Sin u  Cos v
z =  [ c / ( a - c ) ]1/2 Sin v  ( a Sin2 u + b Cos2 u - c )1/2

-First, we have to determine the parameter  Beta  such that    ( § = Integral )

§u1u2  ( a Sin2 u + b Cos2 u )1/2  ( a Sin2 u + b Cos2 u - c ) -1/2 ( (a-b) Sin2 u + Beta ) -1/2  du

=  +/-  §v1v2  ( b Sin2 v + c Cos2 v )1/2  ( a - b Sin2 v - c Cos2 v ) -1/2 ( (b-c) Cos2 v - Beta ) -1/2  dv

-This is done by the routine "BETA"

-Then, the geodesic distance s between the 2 points may be computed by 2 formulae:

s =  §u1u2  [ P + R ( dv/du )2 ]1/2  du          (1)

s =  §v1v2  [ P ( du/dv )2 +R ]1/2  dv           (2)

where the geodesic parameters are given by

P = ( x / u )2 + ( y / u )2 + ( z / u )2
R = ( x / v )2 + ( y / v )2 + ( z / v )2

( the other geodesic parameter Q = ( x / u )( x / v ) + ( y / u )( y / v ) + ( z / u )( z / v ) = 0 here )

-It yields:

P = ( a / (a-c) ) Sin2 u ( a - b Sin2 v - c Cos2 v ) + b Cos2 u Cos2 v + ( c (a-b)2 / (a-c) ) Sin2 u Cos2 u Sin2 v ( a Sin2 u + b Cos2 u - c ) -1

R = ( c / (a-c) ) Cos2 v ( a Sin2 u + b Cos2 u - c ) + b Sin2 u Sin2 v + ( a (b-c)2 / (a-c) ) Sin2 v Cos2 v Cos2 u  ( a - b Sin2 v - c Cos2 v ) -1

-"DIST" calculate s with Gauss-Legendre 10-point formula.
-The intervall of integration is divided in n subintervals ( n must be stored in R00 )

-For each u-value, we have to calculate the corresponding v-value:
-So we must solve a non-linear equation of the form  §v1v f(v) dv - §u1u g(u) du = 0

-This is done by applying again Gauss-Lendre integration + Newton's method !
-So you can guess that a good emulator in turbo mode is not superfluous...

Data Registers:           •  R00 = n          ( Registers R00 thru R07 and R40 thru R49 are to be initialized before executing "BETA" & "DIST" )

•  R01 = a              •  R04 = u1          •  R06 = u2                      R09 to R39: temp
•  R02 = b              •  R05 = v1          •  R07 = v2
•  R03 = c                                   >>>  R08 = Beta

•  R40 = 0.2955242247             •  R45 = 0.6794095683
•  R41 = 0.1488743390             •  R46 = 0.1494513492
•  R42 = 0.2692667193             •  R47 = 0.8650633667            which are the coefficients for Gauss-Legendre 10-point formula
•  R43 = 0.4333953941             •  R48 = 0.06667134431
•  R44 = 0.2190863625             •  R49 = 0.9739065285

Flags:  F01  F02  F06  F07  F09  F10

CF 01  if   ( u1 -  u2 ).( v1 -  v2 ) > 0
SF 01  if   ( u1 -  u2 ).( v1 -  v2 ) < 0

CF 02  if you want to compute the geodesic distance with formula (1)
SF 02  if you want to compute the geodesic distance with formula (2)

Subroutines: /

 01  LBL "BETA"   02  DEG   03  CF 10   04  CF 07   05  49.039   06  STO 10    07  RDN   08  STO 21           09  X<>Y   10  STO 09   11  XEQ 00   12  STO 20    13  LBL 01   14  RCL 21   15  VIEW 21   16  XEQ 00   17  ENTER   18  ENTER   19  X<> 20    20  -   21  X#0?   22  /    23  RCL 09   24  RCL 21   25  STO 09   26  STO T   27  -   28  *   29  +   30  STO 21   31  X#Y?   32  GTO 01   33  STO 08   34  RTN   35  LBL "DIST"   36  LBL 12   37  CF 10   38  RCL 04   39  RCL 06   40  +   41  RCL 05   42  RCL 07   43  +   44  2   45  ST/ Z   46  /   47  FS? 02   48  X<>Y   49  STO 12   50  49.039   51  STO 10   52  RCL 04   53  FS? 02   54  RCL 05   55  STO 11   56  RCL 06   57  FS? 02   58  RCL 07   59   E-8   60  STO 22   61  SF 07   62  XEQ 14   63  ABS   64  D-R   65  CLD   66  RTN   67  GTO 12   68  LBL 00   69  STO 08   70  RCL 06   71  RCL 04   72  STO 11 73  CF 06    74  XEQ 14   75  STO 19    76  RCL 07    77  RCL 05    78  STO 11           79  SF 06   80  XEQ 14   81  FS? 01   82  CHS   83  ST- 19   84  RCL 19    85  RTN   86  GTO 00   87  LBL 13   88  CLX   89  STO 28    90  CLX   91  RCL 21   92  STO 23   93  -   94  RCL 00   95  STO 24   96  ST+ X   97  /   98  STO 25    99  LBL 06 100  ST+ 23 101  RCL 10 102  STO 26 103  LBL 07 104  RCL 23 105  STO 27 106  RCL 25 107  RCL IND 26 108  * 109  ST+ 27 110  - 111  XEQ 02 112  X<> 27 113  XEQ 02 114  RCL 27 115  + 116  DSE 26 117  RCL IND 26 118  * 119  ST+ 28 120  DSE 26 121  GTO 07 122  RCL 25 123  ST+ 23 124  DSE 24 125  GTO 06 126  ST* 28 127  RCL 28 128  RTN 129  LBL 14 130  CLX 131  STO 18 132  CLX 133  RCL 11 134  STO 13 135  - 136  RCL 00 137  STO 14 138  ST+ X 139  / 140  STO 15 141  LBL 08 142  ST+ 13 143  RCL 10 144  STO 16 145  LBL 09 146  RCL 13 147  STO 17  148  RCL 15         149  RCL IND 16 150  * 151  ST+ 17 152  - 153  FS? 07 154  XEQ 05 155  FC? 07 156  XEQ 02 157  X<> 17 158  FS? 07 159  XEQ 05 160  FC? 07 161  XEQ 02 162  RCL 17  163  + 164  DSE 16 165  RCL IND 16 166  * 167  ST+ 18 168  DSE 16 169  GTO 09 170  RCL 15  171  ST+ 13 172  DSE 14 173  GTO 08 174  ST* 18 175  RCL 18 176  RTN 177  LBL 02 178  FS? 06 179  GTO 03 180  SIN 181  X^2 182  RCL 01 183  RCL 02 184  - 185  * 186  STO Y 187  RCL 02 188  + 189  RCL X 190  RCL 03 191  - 192  / 193  X<>Y 194  RCL 08 195  + 196  / 197  X<0? 198  SF 09 199  FS? 10 200  ABS 201  SQRT 202  RTN 203  LBL 03 204  COS 205  X^2 206  RCL 02 207  RCL 03 208  - 209  * 210  ENTER 211  CHS 212  RCL 02 213  + 214  RCL 01 215  RCL Y 216  - 217  / 218  X<>Y 219  RCL 08         220  - 221  / 222  X<0? 223  SF 09 224  FS? 10 225  ABS 226  SQRT 227  RTN 228  LBL 05 229  STO 31  230  FC? 02 231  RCL 04  232  FS? 02 233  RCL 05  234  STO 21  235  CF 06 236  FS? 02 237  SF 06 238  XEQ 13 239  FS? 01 240  CHS 241  STO 20 242  RCL 05 243  FS? 02 244  RCL 04 245  STO 21 246  SF 10 247  RCL 12 248  STO 30 249  SF 06 250  FS? 02 251  CF 06 252  LBL 10 253  CF 09 254  RCL 30 255  FC? 02 256  XEQ 03 257  FS? 02 258  XEQ 02 259  STO 29 260  RCL 30 261  ENTER 262  XEQ 13 263  RCL 20 264  - 265  RCL 29 266  / 267  ST- 30 268  VIEW 30 269  RCL 30 270  X=0? 271  SIGN 272  / 273  ABS 274  RCL 22 275  X 30 282  STO 31 283  STO 12 284  1 285  P-R 286  X^2 287  STO 32 288  RCL 02 289  * 290  X<>Y 291  X^2 292  STO 33         293  RCL 01  294  * 295  + 296  STO 36  297  RCL 03 298  - 299  RCL 30 300  FC? 02 301  STO 12  302  COS 303  X^2 304  STO 34  305  * 306  RCL 03  307  * 308  RCL 01  309  RCL 03 310  - 311  / 312  RCL 30 313  SIN 314  X^2 315  STO 35 316  RCL 33 317  * 318  RCL 02 319  * 320  + 321  RCL 34 322  RCL 03 323  * 324  RCL 35 325  RCL 02 326  * 327  + 328  STO 37 329  RCL 01 330  X<>Y 331  - 332  RCL 34 333  RCL 35 334  * 335  X<>Y 336  / 337  RCL 32 338  * 339  RCL 01 340  * 341  RCL 02 342  RCL 03 343  - 344  X^2 345  * 346  RCL 01 347  RCL 03 348  - 349  / 350  + 351  STO 38 352  RCL 32 353  RCL 33 354  * 355  RCL 35 356  * 357  RCL 01 358  RCL 02 359  - 360  X^2 361  * 362  RCL 03         363  * 364  RCL 01  365  RCL 03 366  - 367  / 368  RCL 36 369  RCL 03  370  - 371  / 372  RCL 32  373  RCL 34  374  * 375  RCL 02  376  * 377  + 378  RCL 01  379  RCL 37 380  - 381  RCL 33 382  * 383  RCL 01 384  ST* Y 385  RCL 03 386  - 387  / 388  + 389  STO 39 390  RCL 36 391  RCL 01 392  RCL 37 393  - 394  * 395  RCL 02 396  RCL 03 397  - 398  RCL 34 399  * 400  RCL 08 401  - 402  * 403  RCL 37 404  RCL 36 405  RCL 03 406  - 407  * 408  RCL 01 409  RCL 02 410  - 411  RCL 33 412  * 413  RCL 08 414  + 415  * 416  FS? 02 417  X<>Y 418  / 419  RCL 39 420  RCL 38 421  FS? 02 422  X<>Y 423 RCL Z 424  * 425  + 426  SQRT 427  RTN 428  END

( 642 bytes / SIZE 050 )

 STACK BETA-INPUT BETA-OUTPUT DIST-OUTPUT Y beta1 / / X beta2 beta s

Where   beta1 & beta2 are 2 initial guesses,  beta the solution and s = geodesic distance.

Example1:   Ellipsoid  x2 / 41 + y2 / 37 + z2 / 35 = 1

u1 = +10°          u2 = +75°         ( all in decimal degrees )
v1 = -15°          v2 = +61°

41  STO 01      10  STO 04     75  STO 06
37  STO 02     -15  STO 05     61  STO 07
35  STO 03

-If you choose n = 4,    4  STO 00

CF 01  since ( u1 -  u2 ) & ( v1 -  v2 )  have the same sign

and if your initial guesses for beta are 0.1  and  0.2

0.1   ENTER^
0.2   XEQ "BETA"  >>>>   beta = 0.1986540209 = R08

-With n = 8  STO 00, the same guesses give  beta = 0.1986540199  which is stored in R08
-Let's keep this value for beta in R08.

-Then, simply press R/S or XEQ "DIST" to get the geodesic distance

-With  n = 2  STO 00  and  CF 02,  it yields  s = 8.594822585
-With  n = 4  STO 00          R/S        >>>>   s = 8.594822582

-You could want to check if using the 2nd formula confirms this value:

2  STO 00  SF 02  R/S  >>>>   s = 8.594822577
4  STO 00    R/S    >>>>     s = 8.594822577

-So the exact result is probably about  8.594822580

Example2:     On the Earth: a triaxial ellipsoid model with semi-axes  6378.172 km , 6378.102 km , 6356.752 km ,

and the equatorial major axis in the direction   14°93W  ,  165°07E

•  Store the squares of these semi-axes into R01-R02-R03

-Calculate the geodesic distance between the U.S. Naval Observatory at Washington (D.C.)    L1 = 77°03'56"0 W ;  b1 = 38°55'17"2 N
and Cape Town Observatory ( South Africa )  L2 = 18°28'35"7 E ;  b2 = 33°56'02"5 S

-With the programs listed in paragraph 2°)b) and 3°), we find the following Jacobi parameters:

u1 = -62°16090881        u2 = 33°42630683           store these 4 numbers     R04   R06        respectively
v1 = +38°84397819       v2 = -33°88879132                    into                      R05   R07

-Set flag F01    SF 01   since ( u1 -  u2 ).( v1 -  v2 ) < 0

-You will find, with  n = 4  STO 00 ,  that  beta = 136050.1048  and  s = 12709.56546  km

Notes:

-For n = 1, the execution time t for "DIST" is about 3 seconds with V41 in turbo mode
-For n = 2,  t = 8s
-For n = 4,  t = 26s

-Multiply these values by 600 to get the execution time with a real HP-41...

-In fact, getting Beta is seldom easy !
-Generally, you will have to calculate f(beta) for many beta-values before finding a small interval (b,b') for which  f(b) < 0 & f(b') > 0

-For that purpose, XEQ 00  with  CF 07  CF 10  and  R10 = 49.039  with your guess in X-register.

-Flag F10 is used to avoid a DATA ERROR line 198 or 223:
-When the Newton's method is applied - lines 252 to 276 - it may happen that the iterations lead to a negative radicand.
-However, the last iteration must not be the result of a wrong calculation:
-In this case, line 278 would produce a NONEXISTENT message.

-The method also works if  a = b > c  or if  a > b = c , but it's not the fastest way to solve the problem in this case...

Problem:

-This program cannot find the geodesic distance between Washington and Paris because the geodesic looks like this:

*
*                                               *
*                                                                     *
*                                                                         Paris
*
Washington

-In other words, the derivative  dv / du  changes its sign
-The program hereunder overcomes this limitation.

b)  Program#2

-So, we have to compute   §v1v0  f(v)  dv - §v0v2  f(v)  dv

with  f(v) = ( b Sin2 v + c Cos2 v )1/2  ( a - b Sin2 v - c Cos2 v ) -1/2 ( (b-c) Cos2 v - Beta ) -1/2

where the denominator  ( (b-c) Cos2 v - Beta )  vanishes between v1 & v2  i-e  ( (b-c) Cos2 v0 - Beta ) = 0

-In order to remove this singularity without producing another singularity for v = 0, I've made the following change of variable:

Cos w =   [ ( Cos2 v - K ) / ( 1 - K ) ] 1/2              with   K = beta / ( b - c )
Sin w =  Sin v /  ( 1 - K ) 1/2

-The integral becomes:    ( b-c ) -1/2 §w1w2  g(w)  dw

with  g(w) = [ b - ( b-c ).{ K + ( 1-K) Cos2 w } ]1/2  [ ( a-b ) - ( b-c ).{ K + ( 1-K) Cos2 w } ] -1/2  [ K + ( 1-K) Cos2 w } ] -1/2

-Flags F01-F02-F04 determine the form of the geodesic between the 2 points:

*                                                     *                                          *                            *
*                                                                *                             *             *                       *
*                                                                         *                    *                      *                       *                                *
*                                                                               *               *                                                       *                    *
*                                                                                    *           *                                                                   *

CF 01  CF 04                                               SF 01   CF 04         CF 01-CF 02-SF04                   SF 01-CF 02-SF 04
CF 02 or SF 02   <-- It's your choice -->     CF 02 or SF 02

-Thus, F02 must be cleared if F04 is set ( add FS? 04  CF 02 inside the listing below if you prefer )

-Of course, it's not always obvious to know the good case in advance.
-So, you'll probably have to make several tests before finding the right configuration.

Data Registers:           •  R00 = n          ( Registers R00 thru R07 and R42 thru R51 are to be initialized before executing "BETA" & "DIST" )

•  R01 = a              •  R04 = u1          •  R06 = u2                      R09 to R38: temp
•  R02 = b              •  R05 = v1          •  R07 = v2
•  R03 = c                                   >>>  R08 = Beta                    R39 = K  ,   R40 = w1  ,  R41 = w2

•  R42 = 0.2955242247             •  R47 = 0.6794095683
•  R43 = 0.1488743390             •  R48 = 0.1494513492
•  R44 = 0.2692667193             •  R49 = 0.8650633667            which are the coefficients for Gauss-Legendre 10-point formula
•  R45 = 0.4333953941             •  R50 = 0.06667134431
•  R46 = 0.2190863625             •  R51 = 0.9739065285

Flags:  F01  F02  F04  F06  F07  F09  F10

F01  &  F04  indicate the kind of geodesic ( see above )

CF 02  if you want to compute the geodesic distance with formula (1)
SF 02  if you want to compute the geodesic distance with formula (2)

Subroutines: /

 01  LBL "BETA"   02  DEG   03  CF 10   04  CF 07   05  51.041   06  STO 10           07  RDN   08  STO 21    09  X<>Y   10  STO 09   11  XEQ 00   12  STO 20    13  LBL 01   14  RCL 21    15  VIEW 21   16  XEQ 00   17  ENTER   18  ENTER   19  X<> 20    20  -   21  X#0?   22  /   23  RCL 09   24  RCL 21   25  STO 09   26  STO T   27  -   28  *   29  +   30  STO 21   31  X#Y?   32  GTO 01   33  STO 08   34  RTN   35  LBL "DIST"   36  LBL 12   37  CF 10   38  RCL 05   39  RCL 07   40  +   41  RCL 04   42  RCL 06   43  +   44  RCL 40   45  RCL 41   46  +   47  2   48  ST/ T   49  ST/ Z   50  /   51  FS? 02   52  X<> Z   53  STO 12   54  51.041   55  STO 10   56  RCL 04   57  FS? 02   58  RCL 05   59  STO 11   60  RCL 06   61  FS? 02   62  RCL 07   63   E-8   64  STO 22   65  SF 07   66  XEQ 14   67  D-R   68  ABS   69  CLD   70  RTN   71  GTO 12   72  LBL 00   73  STO 08   74  RCL 02   75  RCL 03   76  -   77  /   78  STO 39   79  RCL 06   80  RCL 04   81  STO 11   82  CF 06   83  XEQ 14   84  STO 19   85  RCL 07 86  RCL 05           87  FS? 02   88  GTO 11    89  X<>Y   90  1   91  P-R   92  X^2   93  RCL 39    94  -   95  SQRT   96  FS? 04   97  CHS   98  R-P   99  X<>Y 100  STO 41  101  RCL 05  102  1 103  P-R 104  X^2 105  RCL 39  106  - 107  SQRT 108  R-P 109  RDN 110  STO 40 111  LBL 11 112  STO 11 113  SF 06 114  XEQ 14 115  FS? 01 116  CHS 117  ST- 19 118  RCL 19 119  RTN 120  GTO 00 121  LBL 13 122  CLX 123  STO 28 124  CLX 125  RCL 21 126  STO 23 127  - 128  RCL 00 129  STO 24 130  ST+ X 131  / 132  STO 25 133  LBL 06 134  ST+ 23 135  RCL 10 136  STO 26 137  LBL 07 138  RCL 23 139  STO 27 140  RCL 25 141  RCL IND 26 142  * 143  ST+ 27 144  - 145  XEQ 02 146  X<> 27 147  XEQ 02 148  RCL 27 149  + 150  DSE 26 151  RCL IND 26 152  * 153  ST+ 28 154  DSE 26 155  GTO 07 156  RCL 25 157  ST+ 23 158  DSE 24 159  GTO 06 160  ST* 28 161  RCL 28 162  RTN 163  LBL 14 164  CLX 165  STO 18 166  CLX 167  RCL 11 168  STO 13 169  - 170  RCL 00 171  STO 14         172  ST+ X 173  / 174  STO 15  175  LBL 08 176  ST+ 13 177  RCL 10  178  STO 16  179  LBL 09 180  RCL 13 181  STO 17 182  RCL 15  183  RCL IND 16 184  * 185  ST+ 17 186  - 187  FS? 07 188  XEQ 05 189  FC? 07 190  XEQ 02 191  X<> 17 192  FS? 07 193  XEQ 05 194  FC? 07 195  XEQ 02 196  RCL 17  197  + 198  DSE 16 199  RCL IND 16 200  * 201  ST+ 18 202  DSE 16 203  GTO 09 204  RCL 15 205  ST+ 13 206  DSE 14 207  GTO 08 208  ST* 18 209  RCL 18 210  RTN 211  LBL 02 212  FS? 06 213  GTO 03 214  SIN 215  X^2 216  RCL 01 217  RCL 02 218  - 219  * 220  STO Y 221  RCL 02 222  + 223  RCL X 224  RCL 03 225  - 226  / 227  X<>Y 228  RCL 08 229  + 230  / 231  X<0? 232  SF 09 233  FS? 10 234  ABS 235  SQRT 236  RTN 237  LBL 03 238  FS? 02 239  GTO 04 240  COS 241  X^2 242  1 243  RCL 39 244  - 245  * 246  RCL 39 247  + 248  STO Y 249  RCL 03 250  RCL 02 251  - 252  * 253  RCL 02 254  + 255  ENTER 256  CHS 257  RCL 01         258  + 259  / 260  X<>Y 261  / 262  RCL 02  263  RCL 03  264  - 265  / 266  X<0? 267  SF 09 268  FS? 10 269  ABS 270  SQRT 271  RTN 272  LBL 04 273  COS 274  X^2 275  RCL 02 276  RCL 03 277  - 278  * 279  ENTER 280  CHS 281  RCL 02  282  + 283  RCL 01 284  RCL Y 285  - 286  / 287  X<>Y 288  RCL 08 289  - 290  / 291  X<0? 292  SF 09 293  FS? 10 294  ABS 295  SQRT 296  RTN 297  LBL 05 298  STO 31 299  FC? 02 300  RCL 04 301  FS? 02 302  RCL 05 303  STO 21 304  CF 06 305  FS? 02 306  SF 06 307  XEQ 13 308  FS? 01 309  CHS 310  STO 20 311  RCL 40 312  FS? 02 313  RCL 04 314  STO 21 315  SF 10 316  RCL 12 317  STO 30 318  SF 06 319  FS? 02 320  CF 06 321  LBL 10 322  CF 09 323  RCL 30 324  FC? 02 325  XEQ 03 326  FS? 02 327  XEQ 02 328  STO 29 329  RCL 30 330  ENTER 331  XEQ 13 332  RCL 20 333  - 334  RCL 29 335  / 336  ST- 30 337  VIEW 30 338  RCL 30 339  X=0? 340  SIGN 341  / 342  ABS 343  RCL 22         344  X 31 352  STO 30 353  FS? 02 354  GTO 05 355  SIN 356  X^2 357  1 358  RCL 39 359  - 360  * 361  SQRT 362  ASIN 363  STO 30  364  LBL 05 365  RCL 31  366  1 367  P-R 368  X^2 369  STO 32 370  RCL 02 371  * 372  X<>Y 373  X^2 374  STO 33 375  RCL 01 376  * 377  + 378   STO 36 379  RCL 03 380  - 381  RCL 30 382  COS 383  X^2 384  STO 34 385  * 386  RCL 03 387  * 388  RCL 01 389  RCL 03 390  - 391  / 392  RCL 30 393  SIN 394  X^2 395  STO 35 396  RCL 33 397  * 398  RCL 02 399  * 400  + 401  RCL 34 402  RCL 03 403  * 404  RCL 35 405  RCL 02 406  * 407  + 408  STO 37 409  RCL 01 410  X<>Y 411  - 412  RCL 34 413  RCL 35 414  * 415  X<>Y 416  / 417  RCL 32 418  * 419  RCL 01 420  * 421  RCL 02 422  RCL 03 423  - 424  X^2 425  * 426  RCL 01 427  RCL 03         428  - 429  / 430  + 431  STO 38  432  RCL 32  433  RCL 33  434  * 435  RCL 35 436  * 437  RCL 01 438  RCL 02 439  - 440  X^2 441  * 442  RCL 03 443  * 444  RCL 01  445  RCL 03 446  - 447  / 448  RCL 36  449  RCL 03 450  - 451 / 452  RCL 32 453  RCL 34 454  * 455  RCL 02 456  * 457  + 458  RCL 01 459  RCL 37 460  - 461  RCL 33 462  * 463  RCL 01 464  ST* Y 465  RCL 03 466  - 467  / 468  + 469  STO 35 470  RCL 36 471  RCL 01 472  RCL 37 473  - 474  * 475  RCL 02 476  RCL 03 477  - 478  RCL 34 479  * 480  RCL 08 481  - 482  * 483  RCL 37 484  RCL 36 485  RCL 03 486  - 487  * 488  RCL 01 489  RCL 02 490  - 491  RCL 33 492  * 493  RCL 08 494  + 495  * 496  FS? 02 497  X<>Y 498  / 499  RCL 35 500  RCL 38 501  FS? 02 502  X<>Y 503  RCL Z 504  * 505  + 506  SQRT 507  RTN 508  END

( 746 bytes / SIZE 052 )

 STACK BETA-INPUT BETA-OUTPUT DIST-OUTPUT Y beta1 / / X beta2 beta s

Where   beta1 & beta2 are 2 initial guesses,  beta the solution and s = geodesic distance.

Example1&2:  The same as in paragraph 1°) a)

-With the first example CF 01  CF 02 ( or SF 02 )  and CF 04 , you'll get the same results for beta & s
-With the second example  SF 01  CF 02 ( or SF 02 ) and CF 04, same results too.

Example3:    The same ellipsoidal model of the Earth:  semi-axes  6378.172 km , 6378.102 km , 6356.752 km ,

with the equatorial major axis in the direction   14°93W  ,  165°07E

-Calculate the geodesic distance between the U.S. Naval Observatory at Washington (D.C.)    L1 = 77°03'56"0 W ,  b1 = 38°55'17"2 N
and the "Observatoire de Paris"  L2 = 2°20'13"8 E ,  b2 = 48°50'11"2 N

-The routines in paragraph 2°) b) and 3°) return the Jacobi parameters

u1 = -62°16090881        u2 = +17°30207991           store these 4 numbers     R04   R06         respectively
v1 = +38°84397819       v2 =  +48°83869959                     into                      R05   R07

-With  n = 4  STO 00 ,  CF 01  CF 02  SF 04 ,  it yields   beta = 101348.0008  and  the geodesic distance  s = 6181.626001 km

Notes:

-With an oblate ellipsoid of revolution with a = b = 6378.137 km & c = 6356.752 km ,
it yields  12709.52608 km in example 2 and 6181.621921 km in example 3
which may be found more easily by Vincenty's formula.

-This routine is not perfect !
-The small number in R22 = E-8  may sometimes produce an infinite loop inside LBL 10
-Replace it by a less small number ( line 63 )  or

-Stop the routine
-Set flag 21 , R/S
-Place 0 in X-register when the program stops ( line 337 )
-Clear flag 21 and R/S again.

-Moreover, the change of variable w cannot be used if the 2 points are on the equator... except for nearly antipodal points !

-For example, if  L1 =  b1 = 0  and  L1 = 179°51' E ,  b2 = 0  you first find

u1 = 14°93015654        u2 = 194°7801551
v1 = 0                           v2 =  0

-With  n = 8  STO 00  it yields   beta = 16815.92370  &  the geodesic distance  s = 20001.90110 km

c)  Program(s)#3

-Here, the w-values are computed by the classical Runge-Kutta method of order 4 ( LBL 03 )
and the geodesic distance itself by the Newton-Côtes 7-point formula ( LBL 01 )
-We assume that the parameter "beta" has been already calculated by the first part of the program above,
but "FBET" may also be used ( with a root-finding program ) to solve f(beta) = 0

-Moreover, another version of "BETA" is also listed at the end of this paragraph ( see below )

-"DGD" uses the same change of argument of paragraph 1°)b) and flags F01 & F04 must be set or cleared as before:

*                                   *                                          *                            *
*                                               *                             *             *                       *
*                                                       *                    *                      *                       *                                *
*                                                            *               *                                                       *                    *
*                                                               *           *                                                                   *

CF 01-CF 04                             SF 01-CF 04              CF 01-SF04                                 SF 01-SF 04

-Flag F02 is unused here, so "DGD" cannot be used to apply the second formula for the distance.
-Here, we compute:

s =  §u1u2  [ P + R ( dv/du )2 ]1/2  du

Data Registers:           •  R00 = n          ( Registers R00 thru R08 are to be initialized before executing "DGD" - R00 to R07 for "FBET" )

•  R01 = a              •  R04 = u1          •  R06 = u2              R09 = s             R11 = u       R13 to R21: temp
•  R02 = b              •  R05 = v1          •  R07 = v2              R10 = h            R12 = w
•  R03 = c                                 >>>  •  R08 = Beta                                   R22 = K  ,   R23 = w2

Flags:  F00  F01 F04  F10  -  F01  &  F04  indicate the kind of geodesic ( see above )
Subroutines: /

 01  LBL "DGD"   02  CF 00   03  LBL 00   04  DEG   05  RCL 06    06  RCL 04    07  STO 11           08  -   09  RCL 00    10  STO 21   11  6   12  *   13  /   14  STO 10    15  RCL 08   16  RCL 02    17  RCL 03   18  -   19  /   20  STO 22   21  RCL 05   22  1   23  P-R   24  X^2   25  RCL Z   26  -   27  SQRT   28  R-P   29  X<>Y   30  STO 12   31  RCL 07   32  1   33  P-R   34  X^2   35  RCL 22   36  -   37  SQRT   38  FS? 04   39  CHS   40  R-P   41  X<>Y   42  STO 23   43  FS? 00   44  RTN   45  CLX   46  STO 09   47  SF 10   48  LBL 01   49  XEQ 02   50  41   51  FC?C 10   52  ST+ X   53  *   54  ST+ 09   55  XEQ 03   56  XEQ 02   57  216   58  * 59  ST+ 09   60  XEQ 03   61  XEQ 02   62  27   63  *   64  ST+ 09   65  XEQ 03   66  XEQ 02   67  272   68  *   69  ST+ 09   70  XEQ 03   71  XEQ 02   72  27   73  *   74  ST+ 09   75  XEQ 03   76  XEQ 02    77  216   78  *   79  ST+ 09   80  XEQ 03   81  DSE 21   82  GTO 01    83  XEQ 02   84  41   85  *   86  RCL 09           87  +   88  RCL 10    89  *   90  140   91  /   92  D-R   93  ABS   94  STO 09   95  RTN   96  STO 00    97  GTO 00   98  LBL 02   99  RCL 11 100  1 101  P-R 102  X^2 103  STO 13 104  RCL 02 105  * 106  X<>Y 107  X^2 108  STO 14 109  RCL 01 110  * 111  + 112  STO 17 113  RCL 03 114  - 115  STO 18 116  RCL 12 117  SIN 118  1 119  RCL 22 120  - 121  SQRT 122  * 123  STO 16 124  ST* 16 125  ASIN 126  COS 127  X^2 128  STO 15 129  RCL 03 130  * 131  RCL 16 132  RCL 02  133  * 134  + 135  STO 19 136  RCL 01 137  X<>Y 138  - 139  STO 20  140  RCL 01  141  RCL 01  142  RCL 03         143  - 144  / 145  RCL 14  146  * 147  * 148  RCL 02 149  RCL 13 150  * 151  RCL 15 152  * 153  + 154  RCL 01  155  RCL 02 156  - 157  X^2 158  RCL 01 159  RCL 03 160  ST* Z 161  - 162  / 163  RCL 14 164  * 165  RCL 13 166  * 167  RCL 16 168  * 169  RCL 18 170  / 171  + 172  RCL 02 173  RCL 03 174  - 175  X^2 176  RCL 01 177  ST* Y 178  RCL 03 179  - 180  / 181  RCL 15 182  * 183  RCL 16 184  * 185  RCL 13 186  * 187  RCL 20 188  ST* 17 189 / 190  RCL 02  191  RCL 14 192  * 193  RCL 16 194  * 195  + 196  RCL 01  197  RCL 03  198  ST- Y 199  X<>Y 200  / 201  RCL 15         202  * 203  RCL 18  204  * 205  + 206  RCL 01 207  RCL 02 208  - 209  RCL 14 210  * 211  RCL 08  212  + 213  RCL 19 214  * 215  ST* 18 216  CLX 217  RCL 02 218  RCL 03 219  - 220  RCL 15 221  * 222  RCL 08 223  - 224  RCL 17 225  * 226  RCL 18 227  / 228  * 229  + 230  SQRT 231  RTN 232  LBL 03 233  RCL 12 234  RCL 11 235  XEQ 04 236  STO 17 237  2 238  / 239  RCL 12 240  + 241  RCL 11 242  RCL 10 243  2 244  / 245  + 246  STO 18 247  XEQ 04  248  ST+ 17 249  ST+ 17 250  2 251  / 252  RCL 12  253  + 254  RCL 18  255  XEQ 04 256  ST+ 17 257  ST+ 17 258  RCL 12         259  + 260  RCL 10 261  RCL 11  262  + 263  STO 11 264  XEQ 04 265  RCL 17 266  + 267  6 268  / 269  ST+ 12 270  RTN 271  LBL 04 272  1 273  P-R 274  X^2 275  STO 13  276  RCL 02 277  * 278  X<>Y 279  X^2 280  STO 14 281  RCL 01 282  * 283  + 284  STO 16 285  RCL 03 286  - 287  ST/ 16 288  R^ 289  COS 290  X^2 291  1 292  RCL 22 293  STO T 294  - 295  * 296  + 297  STO 15 298  RCL 02 299  RCL 03 300  - 301  * 302  RCL 02 303  X<>Y 304  ST- Y 305  RCL 01  306  RCL 02 307  - 308  + 309  RCL 15  310  * 311  RCL 02 312  RCL 03         313  - 314  * 315  X<>Y 316  / 317  RCL 16  318  * 319  RCL 01 320  RCL 02 321  - 322  RCL 14 323  * 324  RCL 08 325  + 326  / 327  SQRT 328  RCL 10  329  * 330  FS? 01 331  CHS 332  RTN 333  LBL "FBET" 334  STO 08 335  SF 00 336  XEQ 00 337  6 338  ST* 21 339  LBL 05 340  XEQ 03 341  DSE 21 342  GTO 05 343  RCL 12 344  RCL 23 345  - 346  END

( 475 bytes / SIZE 024 )

 STACK DGD-INPUT DGD-OUTPUT FBET-INPUT FBET-OUTPUT Y / / / / X / s beta f(beta)

Where   s = geodesic distance    and      f(beta) = R12 - R23 = 0  if we achieve a perfect accuracy

Example:    The same ellipsoidal model of the Earth:  semi-axes  6378.172 km , 6378.102 km , 6356.752 km ,

with the equatorial major axis in the direction   14°93W  ,  165°07E

-Calculate the geodesic distance between the U.S. Naval Observatory at Washington (D.C.)    L1 = 77°03'56"0 W ,  b1 = 38°55'17"2 N
and the "Observatoire de Paris"  L2 = 2°20'13"8 E ,  b2 = 48°50'11"2 N

-The same Jacobi parameters

u1 = -62°16090881        u2 = +17°30207991           store these 4 numbers     R04   R06         respectively
v1 = +38°84397819       v2 =  +48°83869959                     into                      R05   R07

-If we already know that   beta = 101348.0008   STO 08

-With  n = 1  STO 00 ,  CF 01   SF 04

XEQ "DGD"  >>>>  s = 6181.533192 = R09  &  R12 = 108°0852527  instead of  R23 = 108°0854380

-To try another n-value, simply place it in X-register and  R/S

n = 2   R/S  >>>>   s = 6181.625172 = R09  &  R12 = 108°0854216           ---Execution time = 4mn36s---

n = 4   R/S  >>>>   s = 6181.625979 = R09  &  R12 = 108°0854368           ---Execution time = 8mn59s---

n = 8   R/S  >>>>   s = 6181.626012 = R09  &  R12 = 108°0854379

-So the results obtained with n = 8 or even  n = 4  are correct:  the errors are smaller than a few centimeters !

Notes:

-The geodesic distance is obtained much faster than with the first 2 programs !
-So, this program is better for a real HP-41... and V41 too !
-However, calculating beta will be easier with "BETA" listed in paragraph 1°)-b)

-Nevertheless, you can also use "FBET" to solve f(beta) = 0
-In this case, it's unuseful to compute the distance s ( flag F00 is set line 335 )

-For example  101348.0008  XEQ "FBET"  gives  -0.000001200   with  n = 4  in  R00
and   -0.000000100   with  n = 8  in  R00
-So, a root-finding program could call "FBET" as a subroutine to get the correct beta-value.

Variant:

-Instead of using the classical RK4, we can also use a Runge-Kutta method of order 6.
-For instance, replace lines 232 to 270 ( LBL 03 ) by

 232  LBL 03 233  RCL 12 234  RCL 11 235  XEQ 04 236  STO 17  237  3 238  / 239  RCL 12  240  + 241  RCL 10         242  3 243  / 244  RCL 11  245  + 246  STO 24 247  XEQ 04 248  STO 18 249  1.5 250  / 251  RCL 12 252  + 253  RCL 10  254  1.5 255  / 256  RCL 11 257  + 258  XEQ 04 259  STO 19 260  CHS 261  RCL 18 262  4 263  * 264  + 265  RCL 17  266  + 267  12 268  / 269  RCL 12  270  + 271  RCL 24         272  XEQ 04 273  STO 20  274  90 275  * 276  RCL 19  277  35 278  * 279  + 280  RCL 18  281  110 282  * 283  - 284  RCL 17 285  25 286  * 287  + 288  48 289  / 290  RCL 12 291  + 292  RCL 10 293  1.2 294  / 295  RCL 11  296  + 297  XEQ 04 298  STO 24  299  10 300  / 301  RCL 20         302  2 303  / 304  + 305  RCL 19  306  8 307  / 308  - 309  RCL 18  310   55 311  * 312  RCL 17 313  18 314  * 315  - 316  120 317   / 318  - 319  RCL 12 320  + 321  RCL 10 322  6 323  / 324  RCL 11  325  + 326  XEQ 04 327  STO 25  328  80 329  * 330  RCL 20         331  118 332  * 333  - 334  RCL 18  335  99 336  * 337  + 338  39 339  / 340  RCL 24  341  128 342  * 343  RCL 19 344  215 345  * 346  + 347  RCL 17 348  783 349  * 350  - 351  780 352  / 353  + 354  RCL 12  355  + 356  RCL 10  357  RCL 11 358  + 359  STO 11         360  XEQ 04 361  RCL 17 362 + 363  13 364  * 365  RCL 24  366  RCL 25 367  + 368  32 369  * 370  + 371  RCL 19 372  RCL 20  373  + 374  55 375  * 376   + 377  .5 378  % 379  ST+ 12 380  RTN

( 645 bytes / SIZE 026 )  for the whole program

-With the same example and  n = 1  STO 00 ,  CF 01   SF 04

XEQ "DGD"  >>>>  s = 6181.537185 = R09  &  R12 = 108°0854502  instead of  R23 = 108°0854380

n = 2   R/S  >>>>   s = 6181.625364 = R09  &  R12 = 108°0854381

n = 4   R/S  >>>>   s = 6181.625979 = R09  &  R12 = 108°0854380

n = 8   R/S  >>>>   s = 6181.626012 = R09  &  R12 = 108°0854380

-The convergence is obviously better for w but the precision is not really increased for the geodesic distance, except for small n-values.

-Here is also another variant of "BETA" that employs Ridder's method to solve f(beta) = 0
-You have to place in registers X & Y  2 initial guesses  b & b'  such that  f(b).f(b') < 0

Data Registers:           •  R00 = n          ( Registers R00 thru R07 and R26 thru R35 are to be initialized before executing "BETA" )

•  R01 = a              •  R04 = u1          •  R06 = u2               R10 to R25: temp
•  R02 = b              •  R05 = v1          •  R07 = v2
•  R03 = c                                   >>>  R08 = Beta             R09 = f(Beta)

•  R26 = 0.2955242247             •  R31 = 0.6794095683
•  R27 = 0.1488743390             •  R32 = 0.1494513492
•  R28 = 0.2692667193             •  R33 = 0.8650633667
•  R29 = 0.4333953941             •  R34 = 0.06667134431
•  R30 = 0.2190863625             •  R35 = 0.9739065285

Flags:  F01  F04  F10      F01 & F04 must be used as above.
Subroutines: /

 01  LBL "BETA"   02  DEG   03  35.025   04  STO 25    05  RDN   06  STO 11           07  X<>Y   08  STO 12   09  XEQ 00   10  STO 20    11   E-8   12  STO 24   13  RCL 11   14  XEQ 00   15  STO 09   16  LBL 01   17  VIEW 11   18  RCL 11   19  RCL 12   20  +   21  2   22  /   23  STO 21   24  XEQ 00   25  ABS   26  RCL 24   27  X>Y?   28  GTO 05   29  LASTX   30  STO 22   31  RCL 20   32  RCL 09   33  -   34  SIGN   35  *   36  RCL 22   37  X^2   38  RCL 09   39  RCL 20   40  * 41  -   42  SQRT   43  /   44  RCL 21    45  RCL 12           46  -   47  *   48  RCL 21    49  +   50  STO 23   51  XEQ 00   52  ABS   53  RCL 24    54  X>Y?   55  GTO 06   56  RCL 22   57  LASTX   58  *   59  X<0?   60  GTO 02   61  RCL 20   62  LASTX   63  *   64  X<0?   65  GTO 03   66  LASTX   67  X<> 09   68  STO 20   69  RCL 23   70  X<> 11   71  STO 12   72  GTO 04   73  LBL 02   74  RCL 22   75  STO 20   76  RCL 21   77  STO 12   78  LBL 03   79  LASTX   80  STO 09 81  RCL 23   82  STO 11   83  LBL 04   84  RCL 12    85  RCL 11           86  X#Y?   87  GTO 01   88  GTO 13   89  LBL 05   90  RCL 21    91  GTO 07   92  LBL 06   93  RCL 23    94  LBL 07   95  STO 11   96  LASTX   97  STO 09   98  LBL 13   99  RCL 09 100  RCL 11 101  STO 08 102  CLD 103  RTN 104  LBL 00 105  STO 08 106  RCL 02 107  RCL 03 108  - 109  / 110  STO 10 111  RCL 06 112  RCL 04 113  ENTER 114  CF 10 115  XEQ 14 116  STO 19 117  RCL 07 118  1 119  P-R 120  X^2 121  RCL 10         122  - 123  SQRT 124  FS? 04 125  CHS 126  R-P 127  X<>Y 128  RCL 05  129  1 130  P-R 131  X^2 132  RCL 10  133  - 134  SQRT 135  R-P 136  SF 10 137  XEQ 14 138  FS? 01 139  CHS 140  ST- 19 141  RCL 19 142  RTN 143  GTO 00 144  LBL 14 145  CLX 146  STO 18 147  RDN 148  STO 13 149  - 150  RCL 00 151  STO 14 152  ST+ X 153  / 154  STO 15 155  LBL 08 156  ST+ 13 157  RCL 25 158  STO 16 159  LBL 09 160  RCL 13 161  STO 17  162  RCL 15         163  RCL IND 16 164  * 165  ST+ 17 166  - 167  XEQ 02 168  X<> 17 169  XEQ 02 170  RCL 17  171  + 172  DSE 16 173  RCL IND 16 174  * 175  ST+ 18 176  DSE 16 177  GTO 09 178  RCL 15 179  ST+ 13 180  DSE 14 181  GTO 08 182  RCL 18 183  * 184  RTN 185  LBL 02 186  FS? 10 187  GTO 03 188  SIN 189  X^2 190  RCL 01 191  RCL 02 192  - 193  * 194  STO Y 195  RCL 02 196  + 197  RCL X 198  RCL 03 199  - 200  / 201  X<>Y 202  RCL 08         203  + 204  / 205  SQRT 206  RTN 207  LBL 03 208  COS 209  X^2 210  1 211  RCL 10  212  STO T 213  - 214  * 215  + 216  STO Y 217  RCL 03 218  RCL 02 219  - 220  * 221  RCL 02 222  + 223  ENTER 224  CHS 225  RCL 01 226  + 227  / 228  X<>Y 229  / 230  RCL 02 231  RCL 03 232  - 233  / 234  SQRT 235  RTN 236  END

( 335 bytes / SIZE 036 )

 STACK INPUTS OUTPUTS Y b f(beta) ~ 0 X b' beta

Example:   Again the same one

-With  n = 1  STO 00 ,  CF 01   SF 04   we find that  100000 < beta < 110000

100000  ENTER^
110000  XEQ "BETA"  >>>>   beta = 101348.0011    X<>Y    f(beta) = -4 E-9               ---Execution time = 6m42s---

-With  n = 2  STO 00

100000  ENTER^
110000  XEQ "BETA"  >>>>   beta = 101348.0010    X<>Y    f(beta) = -1 E-9

-With  n = 4  STO 00

100000  ENTER^
110000  XEQ "BETA"  >>>>   beta = 101348.0008    X<>Y    f(beta) = -4 E-9

Notes:

-"BETA" uses Ridder's method as root-finding routine ( lines 01 to 103 )
-Though not very fast, Gaussian integration is however less slow than using "FBET"
-When the program stops, register R08 = beta, ready for using "DGD" above.
-In order to find 2 values for which the results have opposite signs,
you can  XEQ 00  provided  R25 = 35.025 = control number of the Gauss-Legendre coefficients.

2°)  Longitude & Latitude  >>>>  Cartesian Coordinates

a) Geodetic Coordinates

-The formulae between the geodetic longitude l &  latitude b and the cartesian coordinates x , y , z
of a point on a triaxial ellipsoid are given in references [2] & [3]

x = w Cos b Cos l
y = w ( 1 - e2 ) Cos b Sin l                       In fact   l + 14°93  instead of l
z = w ( 1 - e' 2 ) Sin b

where    e2 = 1 - b2/a2    ,    e' 2 = 1 - c2/a2    ,     w = a / ( 1 - e' 2 Sin2b - e2 Cos2 b Sin2l ) 1/2

-This program is written for the Earth, assuming that the semi-axes are:

a = 6378.172 km     with  the equatorial major axis in the direction   14°93W  ,  165°07E
b = 6378.102 km     equatorial minor semi-axis
c = 6356.752 km     polar minor semi-axis

-Change these values if you know better approximations ( lines 08-14-17-19 ).
-The longitudes are measured positively Eastward from the meridian of Greenwich ( otherwise, replace line 09 + by a - ).

Data Registers: /
Flags: /
Subroutine:  "S-R"  Spherical-Rectangular conversion ( cf "Spherical Coordinates for the HP-41" paragraph 1°-a) )

 01  LBL "LB-XYZ"  02  DEG  03  HR  04  STO N              05  COS  06  X<>Y  07  HR  08  14.93  09  +  10  STO M   11  SIN 12  *  13  X^2  14  6356.752  15  X^2  16  X<>Y  17  6378.102   18  X^2  19  6378.172  20  STO P              21  X^2  22  ST- Y 23  ST- T  24  ST/ T  25  /  26  STO O              27  *  28  X<>Y  29  RCL N   30  SIN  31  X^2  32  *  33  + 34  PI  35  SIGN  36  ST+ O  37  ST+ Z  38  +  39  SQRT  40  RCL P              41  X<>Y  42  /  43  RCL M  44  RCL N 45  X<> Z  46  XEQ "S-R"  47  R^  48  ST* T  49  X<> O              50  ST* Z  51  RDN  52  CLA  53  END

( 110 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Z / z Y l ( ° ' " ) y X b ( ° ' " ) x

---Execution time = 6s---

Example:     The Observatoire de Paris l = 2°20'13"8 E ,  b = 48°50'11"2 N

2.20138   ENTER^
48.50112  XEQ "LB-XYZ" >>>>  x = 4016.614852 km
RDN   y = 1248.484248 km
RDN   z = 4778.596570 km

Note:

-The longitudes given in almanachs are more likely geocentric longitudes.
-The program below employs this definition.

b)  Astronomical Coordinates

-If the longitude l is geocentric i-e equals the angle between the meridian plane passing through the observer and the Prime Meridian plane,
and if the latitude b is the angle between the local vertical in the local meridian plane and the plane of the Equator,
the rectangular coordinates x , y , z of a point P on the triaxial ellipsoid may be found as follows:

r2 = b2 + ( a2-b2 ) / [ 1 + (a2/b2) Tan2 L ]                with   r = OQ   where  Q = point of intersection of the equator and the local meridian

R2 = c2 + ( r2-c2 ) / [ 1 + (c2/r2) Tan2 b ]                  and  L = l + 14°93

x = R Cos µ  Cos L
y = R Cos µ  Sin L                        where    R = OP  and  Tan µ = (c2/r2) Tan b
z = R Sin µ

Data Registers: /
Flags: /
Subroutine:  "S-R"  Spherical-Rectangular conversion ( cf "Spherical Coordinates for the HP-41" paragraph 1°-a) )

 01  LBL "LB-XYZ"  02  DEG  03  HR  04  STO N  05  X<>Y  06  HR  07  14.93  08  +  09  STO M             10  1  11  P-R  12  X^2  13  X<>Y 14  6378.172  15  STO O   16  6378.102  17  STO P  18  /  19  *  20  X^2  21  +  22  /  23  RCL O              24  X^2  25  RCL P   26  X^2 27  STO T  28  -  29  *  30  +  31  RCL N              32  1  33  P-R  34  X^2  35  X<>Y  36  X^2  37  6356.752  38  X^2  39  STO P 40  R^  41  STO O              42  /  43  *  44  +  45  /  46  RCL O   47  RCL P  48  -  49  *  50  RCL P   51  ST/ O  52  + 53  SQRT  54  RCL N              55  1  56  P-R  57  RCL O   58  *  59  R-P  60  X<> M   61  R^  62  XEQ "S-R"  63  CLA  64  END

( 121 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Z / z Y l ( ° ' " ) y X b ( ° ' " ) x

---Execution time = 7s---

Example:     Again the Observatoire de Paris l = 2°20'13"8 E ,  b = 48°50'11"2 N

2.20138   ENTER^
48.50112  XEQ "LB-XYZ" >>>>  x = 4016.607084 km
RDN   y = 1248.509238 km
RDN   z = 4778.596570 km

-We also get  R = OP = 6366.073591  in register T

Notes:

-Strictly speaking, the perpendicular to the local meridian plane is not perpendicular to the triaxial ellipsoid,
but the difference is negligible: both versions return the same z-values
-This would not be the case if the flattening were large.

-One could use  TAN  and  ATAN  instead of the polar <> rectangular conversion,
but there would be an OUT OF RANGE if the latitude = 90° and F24 is clear.
-Moreover, the HP-41 wrongly returns  ATAN ( -9.999999999 E99 ) = 90° instead of -90°
-So, using P-R & R-P is preferable

-I've used this program to get the coordinates in the examples of paragraph 1

-If the latitude b is strictly geodetic, i-e if it is the angle between the normal to the ellipsoid and the equatorial plane,
the geocentric longitude b' may be computed by the formula:

Tan b' = (c2/a2) (Tan b ) [ Cos2 L + (a4/b4) Sin2 L ] 1/2    where   L  is the geocentric longitude.

-This formula is used in paragraph 4°)d)

3°)  Cartesian Coordinates <> Jacobi Parameters ( u , v )

-The coordinates of a point on the ellipsoid      x2 / a + y2 / b + z2 / c = 1    ( a > b > c )   are seldom given by ( u , v )  where

x =  [ a / ( a - c ) ]1/2 Cos u  ( a - b Sin2 v - c Cos2 v )1/2
y =    b1/2  Sin u  Cos v
z =  [ c / ( a - c ) ]1/2 Sin v  ( a Sin2 u + b Cos2 u - c )1/2

-The 2 routines below perform the conversions

-In the conversion  ( x , y , z )  >>>  ( u , v ) , we have to solve the equation:

( b2 - b.c ) Sin4 v + [ b.y2 + b.c - b2 - a.y2 - b.(a-c)/c z2 ] Sin2 v + b.(a-c)/c  z2 = 0

-So, the 2nd part of the routine calls "P2" ( quadratic equation )
-The sign of v is determined by the sign of z ( lines 100-101-102 )

Data Registers:         R00 & R04 thru R08: unused       ( Registers R01 thru R03 are to be initialized before executing "UV-XYZ" & "XYZ-UV" )

•  R01 = a                                          R09 to R11: temp
•  R02 = b
•  R03 = c                                          x2 / a + y2 / b + z2 / c = 1    ( a > b > c )

Flag:  F00
Subroutine:  "P2"  ( quadratic equations ) - cf "Polynomials for the HP-41" paragraph 1°)-a)

 01  LBL "UV-XYZ"   02  1   03  P-R   04  STO 09    05  X^2   06  RCL 02    07  *   08  X<>Y   09  STO 10               10  X^2   11  RCL 01    12  *   13  +   14  RCL 03    15  ST- Y   16  *   17  SQRT   18  STO 11   19  SIGN   20  P-R   21  ST* 10   22  X^2   23  RCL 03   24  *   25  X<>Y   26  ST* 11   27  X^2 28  RCL 02    29  *   30  +   31  CHS   32  RCL 01    33  ST+ Y   34  *   35  SQRT   36  RCL 01               37  RCL 03    38  -   39  SQRT   40  ST/ 11   41  /   42  ST* 09   43  RCL 11    44  RCL 10   45  RCL 02   46  SQRT   47  *   48  RCL 09   49  RTN   50  LBL "XYZ-UV"   51  STO 09    52  RDN   53  STO 10   54  X<>Y 55  STO 11   56  RCL 02    57  RCL 03    58  -   59  RCL 02    60  *   61  ENTER^   62  CHS   63  RCL 02               64  RCL 01    65  -   66  RCL 10    67  X^2   68  *   69  +   70  RCL 01   71  RCL 03   72  ST- Y   73  /   74  RCL 02    75  *   76  RCL 11   77  X^2   78  *   79  ST- Y   80  R^   81  X=0? 82  GTO 01   83  RDN   84  XEQ "P2"   85  FS? 00   86  SF 41   87  X>Y?   88  X<>Y   89  X<0?   90  X<>Y   91  GTO 02   92  LBL 01   93  X<> Z   94  /   95  CHS   96  LBL 02   97  ENTER^   98  SQRT   99  ASIN 100  RCL 11             101  SIGN 102  * 103  X<>Y 104  RCL 03  105  RCL 02 106  - 107  * 108  RCL 03 109  - 110  RCL 01 111  ST+ Y 112  * 113  SQRT 114  RCL 01 115  RCL 03  116  - 117  SQRT 118  RCL 09  119  * 120  X<>Y 121  / 122  RCL 10             123  RCL 02  124  SQRT 125  / 126  R^ 127  COS 128  X=0? 129  SIGN 130  / 131  X<>Y 132  R-P 133  RDN 134  END

( 171 bytes / SIZE 012 )

 STACK INPUTS1 OUTPUTS1 INPUTS2 OUTPUTS2 Z / z z v Y v y y v X u x x u

Example:   The ellipsoid  x2 / 41 + y2 / 37 + z2 / 35 = 1              u = +10° ,   v = -15°

41  STO 01
37  STO 02
35  STO 03

-Set the HP-41 in DEG mode:  XEQ "DEG"

15  CHS  ENTER^
10  XEQ "UV-XYZ"  >>>>   x =  6.235047001
RDN    y =  1.020269420
RDN    z = -0.910302041

-With these values in registers  X , Y , Z  ( RDN RDN )

XEQ "XYZ-UV"  ( or simply R/S )  gives again  u = 10  RDN  v = -15 , with small roudoff-errors in the last decimals.

Notes:

-These programs work in all angular modes.
-"XYZ-UV" does not check that ( x , y , z ) is on the ellipsoid.

-Lines 80 to 83 and 91 to 96 are only useful if   b = c
-If the quadratic equation has 2 complex conjugate roots - which should not happen - the HP-41 displays NONEXISTENT ( line 86 )

-Registers R00 & R04 to R08 are unused because they contain important data for the main program of paragraph 1°)
-Otherwise, R09-R10-R11 could be replaced by R00-R04-R05.

4°)  Approximate Methods

a)  Program#1

-The center O of the ellipsoid and 2 points M & N define a plane ( at least if they are not aligned ) that intersects the ellipsoid.
and we can compute the distance of the arc of ellipse MN thus defined.
-This will not be - in general - the geodesic distance, but it will remain slightly larger if the flattenings are small.

-The following "TGD" uses this method for the Earth, assuming that the semi-axes are:

a = 6378.172 km     with  the equatorial major axis in the direction   14°93W  ,  165°07E
b = 6378.102 km     equatorial minor semi-axis
c = 6356.752 km     polar minor semi-axis                    ( same values as in paragraph 2 )

-A point P on the arc MN may be determined by the vectorial equality:

OP = x u + y v          where  u = OM / || OM ||   and  v = ON / || ON ||

-We have   x = R Sin ( W - µ ) / Sin W  and  y = R Sin µ / Sin W

with   W = the angle ( OM , ON )     µ = the angle ( OM , OP )    and    R = OP

-Writing that P(X,Y,Z)  is on the ellipsoid  X2/a2 + Y2/b2 + Z2/c2 = 1,  we find that

R / Sin W = 1 / ( A2 + B2 + C2 )1/2   with

A = [ x1 Sin ( W - µ ) + x2 Sin µ ] / a
B = [ y1 Sin ( W - µ ) + y2 Sin µ ] / b        where  u ( x1 , y1 , z1 )  &  v ( x2 , y2 , z2 )
C = [ z1 Sin ( W - µ ) + z2 Sin µ ] / c

-Then, we compute  dR / dµ  and the arc length is    §0W [ R2 + ( dR/dµ )2 ] 1/2

Data Registers:   R00 to R09 & R20 to R37: temp          ( Registers R10 thru R19 are to be initialized before executing "TGD" )

R03 = n = the number of subintervalls in which the Gaussian formula is applied.  n = 1 seems enough ( line 56 )

R20 = a , R21 = b , R22 = c ,  R23 to R28 contain the coordinates of the unit vectors u & v   , R29 = W

•  R10  = 0.2955242247                  •  R11 = 0.1488743390
•  R12  = 0.2692667193                  •  R13 = 0.4333953941
•  R14  = 0.2190863625                  •  R15 = 0.6794095683
•  R16  = 0.1494513492                 •  R17 = 0.8650633667
•  R18  = 0.06667134431                •  R19 = 0.9739065285

Flags: /
Subroutines:  "LB-XYZ" cf paragraph 2°)b) above
"GL10"  Gauss-Legendre 10-point formula  cf "Gaussian Integration for the HP-41" paragraph 1°)-b)

Note:

-If you have the M-Code routines  GX & GW no register must be initialized and you can replace
registers R20 , ......... , R38  by  R10 , .............. , R28 in the listing below

 01  LBL "TGD"   02  DEG   03  STO 08    04  RDN   05  STO 07    06  RDN   07  XEQ "LB-XYZ"   08  STO 23               09  X<>Y   10  STO 24    11  R-P   12  RCL Z   13  STO 25    14  R-P   15  ST/ 23   16  ST/ 24   17  ST/ 25   18  RCL 07   19  RCL 08   20  XEQ "LB-XYZ"   21  STO 26   22  X<>Y   23  STO 27   24  R-P   25  RCL Z   26  STO 28   27  R-P   28  ST/ 26   29  ST/ 27   30  ST/ 28   31  RCL 23   32  RCL 26 33  *   3 4 RCL 24    35  RCL 27    36  *   37  +   38  RCL 25    39  RCL 28               40  *   41  +   42  ACOS   43  STO 02    44  STO 29   45  6378.172   46  STO 20   47  .07   48  -   49  STO 21   50   6356.752   51  STO 22   52  CLX   53  STO 01   54  "S"   55  ASTO 00   56  SIGN   57  LBL 00   58  STO 03   59  XEQ "GL10"   60  RCL 29   61  SIN   62  *   63  D-R   64  RTN 65  GTO 00   66  LBL "S"   67  STO Y   68  1   69  P-R   70  STO 31               71  RDN   72  STO 30    73  RCL 26    74  *   75  RCL 29   76  RCL Z   77  -   78  1   79  P-R   80  STO 33   81  RDN   82  STO 32   83  RCL 23   84  *   85  +   86  RCL 20   87  /   88  STO 34   89  X^2   90  RCL 32   91  RCL 24   92  *   93  RCL 27   94  RCL 30   95  *   96  + 97  RCL 21   98  /   99  STO 35 100  X^2 101  + 102  RCL 25             103  RCL 32  104  * 105  RCL 28 106  RCL 30  107  * 108  + 109  RCL 22 110  / 111  STO 36 112  X^2 113  + 114  STO 37 115  RCL 26 116  RCL 31 117  * 118  RCL 23 119  RCL 33 120  * 121  - 122  RCL 34 123  * 124  RCL 20 125  / 126  RCL 27 127  RCL 31 128  * 129  RCL 24 130  RCL 33 131  * 132  - 133  RCL 35             134  * 135  RCL 21  136  / 137  + 138  RCL 28  139  RCL 31 140  * 141  RCL 25 142  RCL 33 143  * 144  - 145  RCL 36 146  * 147  RCL 22 148  / 149  + 150  X^2 151  RCL 37 152  3 153  Y^X 154  / 155  RCL 37 156  1/X 157  + 158  SQRT 159  RTN 160  END

( 279 bytes / SIZE 038 )

 STACK INPUTS OUTPUTS T l1 ( ° ' " ) / Z b1 ( ° ' " ) / Y l2 ( ° ' " ) / X b2 ( ° ' " ) D ( km )

Example:       Calculate the geodesic distance between the U.S. Naval Observatory at Washington (D.C.) and the "Observatoire de Paris"

l1 = 77°03'56"0 W ,  b1 = 38°55'17"2 N
l2 = 2°20'13"8 E  ,    b2 = 48°50'11"2 N

-77.0356    ENTER^
38.55172  ENTER^
2.20138   ENTER^
48.50112  XEQ "TGD"  >>>>   D = 6181.628624 km                                                       ---Execution time = 67s---

-To perform the calculation with another n-value, simply place it in X-register and  R/S
-For instance:  2  R/S  >>>>  D = 6181.628624 km , same result !

Notes:

-With this example, we've found in paragraph 1°)-b)  s = 6181.626001 km
-So, this simplified method gives an error < 3 meters in this case.
-For the geodesic distance between Washington & Cape Town Observatories, it yields  D = 12709.56721 km
whereas the exact value was  s = 12709.56545 km, an error < 2 meters for D

-However, with ( -40° , -40° ) ( +120° , +50° ) it yields 18095.54550 km
but the exact geodesic distance is 18095.48497 km
-So, the errors may reach at least 60 meters.

-They seem to remain relatively small ( except for nearly antipodal points )
but Andoyer's formula of order 2 ( for an ellipsoid of revolution ) remains preferable...

-If you apply this method to the first example at the top of this page ( modify the corresponding lines in "TGD" and "LB-XYZ" ),
you get, with n = 4 ,  D = 8.594880194  instead of  s = 8.594822580
-The difference is more important here since the flattenings are larger too.

-Gauss-Legendre 10-point formula may of course be replaced by another integration formula...

b)  Program#2

-To get a better accuracy, we can also divide the arc MN in several parts
-In the following program, 2 points P & Q are inserted between M & N
and we compute the sum of 3 arc lentgths  MP + PQ + QN

-P and Q are chosen to minimize this sum.

Ph                                Qh
M --------------------- P0 -------------------- Q0 ---------------------N
P -h                              Q -h

-If  W is the angle ( OM , ON ) , P0 & Q0 are in the plane ( OMN )  and  ( OM , OP0 ) = ( OP0 , OQ0 ) = ( OQ0 , ON ) = W / 3

-The points Ph , Qh  are above the plane ( OMN )  with  ( OP0 , OPh ) = ( OQ0 , OQh ) = h
-The points P -h , Q -h  are below the plane ( OMN )  with  ( OP0 , OP -h ) = ( OQ0 , OQ -h ) = h

where h is a small angle ( I've chosen h = W / 1000 , lines 64-65 )

----------------------------------------------------------------------------------------------------------------------------------------

-Let a point P ( on the ellipsoid ) with ( OM , OP' ) = µ , ( OP' , OP ) = h    ,    P' = orthogonal projection of P on the plane (OMN)

OP = x u + y v + z w       where  u = OM / || OM ||   ,  v = ON / || ON ||  ,  w = u x v          ( x = cross-product )

-We have, with R = OP

x = R Cos h Sin ( W - µ ) / Sin W
y = R Cos h Sin µ / Sin W
z = R Sin h / Sin W

-Writing that P(X,Y,Z)  is on the ellipsoid  X2/a2 + Y2/b2 + Z2/c2 = 1,  we find that

R / Sin W = 1 / ( A2 + B2 + C2 )1/2   with

A = [ x1 Cos h Sin ( W - µ ) + x2 Cos h Sin µ + ( y1z2 - y2z1 ) Sin h ] / a
B = [ y1 Cos h Sin ( W - µ ) + y2 Cos h Sin µ + ( z1x2 - z2x1 ) Sin h ] / b                       where  u ( x1 , y1 , z1 )  &  v ( x2 , y2 , z2 )
C = [ z1 Cos h Sin ( W - µ ) + z2 Cos h Sin µ + ( x1y2 - x2y1 ) Sin h ] / c

-----------------------------------------------------------------------------------------------------------------------------------------

-Several "distances" are calculated ( all the following "arcs" are elliptic arcs ).

f(0,0) = arc(MP0) + arc(P0Q0) + arc(Q0N)
f(h,0) = arc(MPh) + arc(PhQ0) + arc(Q0N)
f(-h,0) = arc(MP-h) + arc(P-hQ0) + arc(Q0N)
f(0,h) = arc(MP0) + arc(P0Qh) + arc(QhN)
f(0,-h) = arc(MP0) + arc(P0Q-h) + arc(Q-hN)
f(h,h) = arc(MPh) + arc(PhQh) + arc(QhN)

-f is approximated by a polynomial ( in 2 variables ) of degree 2
-Finally, the minimum is found after solving a 2x2 linear system.

Data Registers:   R00 = n = number of subintervals of integration

R01 = a             R05 = x1             R08 = x2           R21 = h                   R48 = f(h,0)         R51 = f(0,-h)
R02 = b             R06 = y1             R09 = y2           R22 = W                R49 = f(-h,0)        R52 = f(h,h)
R03 = c             R07 = z1             R10 = z2            R47 = f(0,0)           R50 = f(0,h)

>>>>    When the program stops,    R04 = s

Flag:  F00

CF 00 for the Earth:  n = 1 is sufficient and R01-R02-R03 are initialized   ( lines 08 to 17 )
SF 00 for another ellipsoid:  you'll have to initialize R00-R01-R02-R03

Subroutine:  "S-R"  Spherical-Rectangular conversion ( cf "Spherical Coordinates for the HP-41" paragraph 1°-a) )

-Line 259 is a three-byte GTO 10 ( not really important )
-LBL 14 is actually the routine "LB-XYZ" listed in paragraph 2°)b) - slightly modified.

 01  LBL "TGD2"   02  STO 10   03  RDN   04  STO 09    05  RDN   06  FS? 00   07  GTO 00   08  6378.172   09  STO 01   10  .07   11  -   12  STO 02           13  21.35   14  -   15  STO 03    16  SIGN   17  STO 00    18  RDN   19  LBL 00   20  XEQ 14   21  STO 05    22  X<>Y   23  STO 06   24  R-P   25  RCL Z   26  STO 07   27  R-P   28  ST/ 05   29  ST/ 06   30  ST/ 07   31  RCL 09   32  RCL 10   33  XEQ 14   34  STO 08   35  X<>Y   36  STO 09   37  R-P   38  RCL Z   39  STO 10   40  R-P   41  ST/ 08   42  ST/ 09   43  ST/ 10   44  LBL 10   45  RCL 05   46  STO 41   47  RCL 08   48  STO 44   49  *   50  RCL 06   51  STO 42   52  RCL 09   53  STO 45   54  *   55  +   56  RCL 07   57  STO 43   58  RCL 10   59  STO 46   60  *   61  +   62  ACOS   63  STO 22   64   E3   65  /   66  STO 21   67  RCL 22   68  3   69  /   70  XEQ 11   71  23.032009   72  REGMOVE   73  RCL 22   74  ST+ X   75  3   76  /   77  XEQ 11   78  RCL 38   79  STO 08   80  RCL 39    81  STO 09   82  RCL 40   83  STO 10   84  XEQ 12   85  STO 48   86  STO 52   87  RCL 32   88  STO 08   89  RCL 33   90  STO 09   91  RCL 34   92  STO 10   93  XEQ 12   94  STO 49   95  RCL 35   96  STO 08   97  RCL 36 98  STO 09   99  RCL 37 100  STO 10  101  XEQ 12 102  STO 47 103  STO 50 104  STO 51 105  RCL 26 106  STO 05 107  RCL 27 108  STO 06  109  RCL 28         110  STO 07 111  XEQ 12  112  ST+ 47 113  RCL 29  114  STO 05 115  RCL 30 116  STO 06 117  RCL 31  118  STO 07 119  XEQ 12 120  ST+ 50 121  RCL 23 122 STO 05 123  RCL 24 124  STO 06 125  RCL 25 126  STO 07 127  XEQ 12 128  ST+ 51 129  RCL 44 130  STO 08 131  RCL 45 132  STO 09 133  RCL 46 134  STO 10 135  XEQ 12 136  ST+ 51 137  RCL 29 138  STO 05 139  RCL 30 140  STO 06 141  RCL 31 142  STO 07 143  XEQ 12 144  ST+ 50 145  ST+ 52 146  RCL 26 147  STO 05 148  RCL 27 149  STO 06 150  RCL 28 151  STO 07 152  XEQ 12 153  ST+ 47 154  ST+ 48 155  ST+ 49 156  RCL 32 157  STO 08 158  RCL 33 159  STO 09 160  RCL 34 161  STO 10 162  XEQ 12 163  ST+ 49 164  RCL 38 165  STO 08 166  RCL 39 167  STO 09 168  RCL 40 169  STO 10 170  XEQ 12 171  ST+ 48 172  RCL 29 173  STO 05 174  RCL 30 175  STO 06 176  RCL 31 177  STO 07 178  XEQ 12 179  ST+ 52 180  RCL 48 181  RCL 49 182  + 183  RCL 47 184  ST+ X 185  - 186  STO 05 187  RCL 50 188  RCL 51 189  - 190  STO 06 191  * 192  RCL 52 193  RCL 48 194  - 195  RCL 47 196  + 197  RCL 50  198  - 199  STO 07 200  RCL 48 201  RCL 49 202  - 203  STO 09         204  * 205  - 206  STO 11  207  RCL 50  208  RCL 51  209  + 210  RCL 47 211  ST+ X 212  - 213  STO 08 214  RCL 09 215  * 216  RCL 06 217  RCL 07 218  * 219  - 220  RCL 07 221  X^2 222  RCL 05 223  RCL 08 224  * 225  - 226  ST+ X 227  ST/ 11 228  / 229  STO 10 230  RCL 05 231  * 232  RCL 09 233  + 234  RCL 10 235  * 236  RCL 08 237  RCL 11 238  * 239  RCL 06 240  + 241  RCL 11 242  * 243 + 244  2 245  / 246  RCL 10 247  RCL 11 248  * 249  RCL 07 250  * 251  + 252  RCL 47 253  + 254  STO 04 255  RTN 256  STO 00 257  41.005006 258  REGMOVE 259  GTO 10 260  LBL 11 261  STO 11 262  SIN 263  STO 13 264  RCL 08 265  * 266  RCL 22 267  RCL 11 268  - 269  SIN 270  STO 14 271  RCL 05 272  * 273  + 274  STO 26 275  RCL 21 276  COS 277  STO 15 278  * 279  RCL 06 280  RCL 10 281  * 282  RCL 07 283  RCL 09 284  * 285  - 286  RCL 21 287  SIN 288  STO 12 289  * 290  ST+ Z 291  - 292  STO 23 293  X<>Y 294  STO 29  295  RCL 06 296  RCL 14 297  * 298  RCL 09  299  RCL 13         300  * 301  + 302  STO 27 303  RCL 15  304  * 305  RCL 07  306  RCL 08 307  * 308  RCL 05 309  RCL 10  310 * 311  - 312  RCL 12 313  * 314  ST+ Z 315  - 316  STO 24 317  X<>Y 318  STO 30 319  RCL 07 320  RCL 14 321  * 322  RCL 10 323  RCL 13 324  * 325  + 326  STO 28 327  RCL 15 328  * 329  RCL 05 330  RCL 09 331  * 332  RCL 06 333  RCL 08 334  * 335  - 336  RCL 12 337  * 338  ST+ Z 339  - 340  STO 25 341  X<>Y 342  STO 31 343  RCL 30 344  R-P 345  RCL 29 346  R-P 347  ST/ 29 348  ST/ 30 349  ST/ 31 350  RCL 28 351  RCL 27 352  R-P 353  RCL 26 354  R-P 355  ST/ 26 356  ST/ 27 357  ST/ 28 358  RCL 25 359  RCL 24 360  R-P 361  RCL 23 362  R-P 363  ST/ 23 364  ST/ 24 365  ST/ 25 366  RTN 367  LBL 12 368  CLX 369  STO 04 370  RCL 05 371  RCL 08 372  * 373  RCL 06 374  RCL 09 375  * 376  + 377  RCL 07 378  RCL 10 379  * 380  + 381  ACOS 382  STO 11 383  RCL 00 384  STO M 385  ST+ X 386  / 387  STO 12 388  .6 389  SQRT 390  * 391  STO 16  392  LBL 01 393  RCL 12 394  RCL 16         395  - 396  XEQ 13 397  ST+ 04 398  RCL 12  399  XEQ 13 400  1.6 401  * 402  ST+ 04 403  RCL 12 404  RCL 16  405  + 406  XEQ 13 407  ST+ 04 408  RCL 11  409  RCL 00 410  / 411  ST+ 12 412  DSE M 413  GTO 01 414  RCL 04 415  * 416  RCL 11 417  SIN 418  * 419  3.6 420  / 421  D-R 422  STO 04 423  RTN 424  LBL 13 425  STO 13 426  1 427  P-R 428  STO 15 429  RDN 430  STO 14 431  RCL 08 432  * 433  RCL 11 434  RCL 13 435  - 436  1 437  P-R 438  STO 18 439  RDN 440  STO 17 441  RCL 05 442  * 443  + 444  RCL 01 445  / 446  STO 19 447  X^2 448  RCL 06 449  RCL 17 450  * 451  RCL 09 452  RCL 14 453  * 454  + 455  RCL 02 456  / 457  STO 20 458  X^2 459  + 460  RCL 07 461  RCL 17 462  * 463  RCL 10 464  RCL 14 465  * 466  + 467  RCL 03 468  / 469  STO 13 470  X^2 471  + 472  STO 14 473  RCL 08 474  RCL 15 475  * 476  RCL 05 477  RCL 18 478  * 479  - 480  RCL 19 481  * 482  RCL 01 483  / 484  RCL 09 485  RCL 15 486  * 487  RCL 06 488  RCL 18  489  * 490  - 491  RCL 20         492  * 493  RCL 02  494  / 495  + 496  RCL 10 497  RCL 15 498  * 499  RCL 07  500  RCL 18 501  * 502  - 503  RCL 13 504  * 505  RCL 03  506  / 507  + 508  X^2 509  RCL 14 510  3 511  Y^X 512  / 513  RCL 14 514   1/X 515  + 516  SQRT 517  RTN 518  LBL 14 519  DEG 520  HR 521  STO N 522  X<>Y 523  HR 524  14.93 525  FS? 00 526  CLX 527  + 528  STO M 529  1 530  P-R 531  X^2 532  X<>Y 533  RCL 01 534  RCL 02 535  / 536  * 537  X^2 538  + 539  / 540  RCL 01 541  X^2 542  RCL 02 543  X^2 544  STO T 545  - 546  * 547  + 548  RCL N 549  1 550  P-R 551  X^2 552  X<>Y 553  X^2 554  RCL 03 555  X^2 556  R^ 557  STO O 558  / 559  * 560  + 561  / 562  RCL O 563  RCL 03 564  X^2 565  - 566  * 567  RCL 03 568  X^2 569  ST/ O 570  + 571   SQRT 572   RCL N 573  1 574  P-R 575  RCL O 576  * 577  R-P 578  X<> M 579  R^ 580  XEQ "S-R" 581  CLA 582  END

( 835 bytes / SIZE 053 )

 STACK INPUTS OUTPUTS T l1 ( ° ' " ) / Z b1( ° ' " ) / Y l2 ( ° ' " ) / X b2 ( ° ' " ) D ( km )

Example1:       Calculate the geodesic distance between the U.S. Naval Observatory at Washington (D.C.) and the "Observatoire de Paris"

l1 = 77°03'56"0 W ,  b1 = 38°55'17"2 N
l2 = 2°20'13"8 E  ,    b2 = 48°50'11"2 N

CF 00

-77.0356    ENTER^
38.55172  ENTER^
2.20138   ENTER^
48.50112  XEQ "TGD2"  >>>>   D = 6181.626264 km                                             ---Execution time = 4mn04s---

-The error is about 26 centimeters.

Example2:         Calculate the geodesic distance between the U.S. Naval Observatory at Washington (D.C.)

L1 = 77°03'56"0 W ;  b1 = 38°55'17"2 N   and Cape Town Observatory ( South Africa )
L2 = 18°28'35"7 E ;  b2 = 33°56'02"5 S

-77.0356    ENTER^
38.55172   ENTER^
18.28357   ENTER^
-33.56025  XEQ "TGD2"  >>>>    D = 12709.56616 km

-The exact value was  s = 12709.56545 km, so, the error = 71 cm

Example3:    With ( -40° , -40° ) ( +120° , +50° ) it yields 18095.48823 km

-Exact geodesic distance = 18095.48497 km  -> Error = 3m26

-In this case, the error remains larger than 1 meter.

Notes:

-Gauss-Legendre 3-point formula is used to compute the integral.

-If you prefer the geodetic longitude, replace LBL 14 ( lines 518 to 581 ) by the routine listed in paragraph 2°)a)
with small modifications since the semi-axes a , b , c are in registers R01-R02-R03 here.

-If you have the cartesian coordinates  x , y , z  of the 2 points M & N, store these values in R05 to R10
then divide the 3 registers R05-R06-R07 by the norm of the vector OM and the 3 registers R08-R09-R10 by the norm of the vector ON
and press  XEQ 10

-With the first example at the top of this page, you will get something like  8.594838705  ( exact result = 8.594822580 )

-The sum of the elliptic arcs lengths arc(MP0) + arc(P0Q0) + arc(Q0N) is in register R47 and in L-register at the end.
-For instance, in example 3, you get  18095.54550  and  "TGD2"  has diminished this value by about  57 meters.

-This program should not be used for nearly antipodal points... though it may work in some cases.

-The union of such 3 arcs is not really differentiable at 2 points, so using a polynomial approximation remains doubtful.
-A better approach is given hereunder.

c)  Program#3

-This version employs the same angles µ & h to determine the position of a point P of the curve on the ellipsoid.
-But instead of adding several arc lengths in different planes, we search an approximation of the function  h(µ) under the form:

h(µ) = h1 Sin ( 180° µ / W ) + h2 Sin ( 360° µ / W )              ( in fact the first terms of a Fourier series )

-And the corresponding integral  s = §0W [ R2 Cos2 h  + R2 ( dh/dµ )2 + ( dR/dµ )2  ] 1/2

=  ( Sin W )    §0W [ ( Cos2 h ) / T + (1/T) ( dh/dµ )2 + ( dr/dµ )2  ] 1/2 dµ    is evaluated by a Gaussian quadrature.

where   T = ( A2 + B2 + C2 )   with the same notations as in paragraph above:

r = 1/sqrt(T)                  dr/dµ = r/µ + ( r/h ) ( dh/dµ )

r = 1 / ( A2 + B2 + C2 )1/2   with

A = [ x1 Cos h Sin ( W - µ ) + x2 Cos h Sin µ + ( y1z2 - y2z1 ) Sin h ] / a
B = [ y1 Cos h Sin ( W - µ ) + y2 Cos h Sin µ + ( z1x2 - z2x1 ) Sin h ] / b                       where  u ( x1 , y1 , z1 )  &  v ( x2 , y2 , z2 )
C = [ z1 Cos h Sin ( W - µ ) + z2 Cos h Sin µ + ( x1y2 - x2y1 ) Sin h ] / c

the partial derivatives

r/µ = ( -1/sqrt(T3) ) [ (A/a) { -x1 Cos h Cos ( W - µ ) + x2 Cos h Cos µ }
+ (B/b) { -y1 Cos h Cos ( W - µ ) + y2 Cos h Cos µ }
+ (C/c) { -z1 Cos h Cos ( W - µ ) + z2 Cos h Cos µ } ]
and

r/h = ( -1/sqrt(T3) ) [ (A/a) { -x1 Sin h Sin ( W - µ ) - x2 Sin h Sin µ + ( y1z2 - y2z1 ) Cos h }
+ (B/b) { -y1 Sin h Sin ( W - µ ) - y2 Sin h Sin µ + ( z1x2 - z2x1 ) Cos h }
+ (C/c) { -z1 Sin h Sin ( W - µ ) - z2 Sin h Sin µ + ( x1y2 - x2y1 ) Cos h } ]

-With  H = W / 1000  ( lines 89-90 ) , we calculate the following lengths

A = f(0,0)
B = f(H,0)
C = f(-H,0)                     f = length s corresponding to h1 = -H , 0 , +H ,  h2 = -H , 0 , +H
D = f(0,H)
E = f(0,-H)

-Then, the extremum is estimated by the formula:

s ~ A - (B-C)2 / [ 8.(B+C-2.A) ] - (D-E)2 / [ 8.(D+E-2.A) ]

Data Registers:   R00 = n = number of subintervals of integration        ( Registers R41 to R46 are to be initialized before executing "TGD3" )

R01 = a             R05 = x1             R08 = x2             R32 = W                  R35 = f(-H,0)
R02 = b             R06 = y1             R09 = y2             R33 = f(0,0)             R36 = f(0,H)
R03 = c             R07 = z1             R10 = z2              R34 = f(H,0)            R37 = f(0,-H)

>>>>    When the program stops,    R04 = s

•  R41 = 0.4679139346                •  R44 = 0.6612093865
•  R42 = 0.2386191861                •  R45 = 0.1713244924                   ( coefficients for Gauss-Legendre 6-point formula )
•  R43 = 0.3607615730                •  R46 = 0.9324695142

Flag:  F00

CF 00 for the Earth:  n = 1 is sufficient and R01-R02-R03 are initialized   ( lines 08 to 17 )
SF 00 for another ellipsoid:  you'll have to initialize R00-R01-R02-R03

Subroutine:  "S-R"  Spherical-Rectangular conversion ( cf "Spherical Coordinates for the HP-41" paragraph 1°-a) )

-Line 146 is a three-byte GTO 10 ( not really important )
-LBL 14 is actually the routine "LB-XYZ" listed in paragraph 2°)b) - slightly modified.
-You can also add several lines to initialize R41 to R46 , for instance after line 43.

 01  LBL "TGD3"   02  STO 10   03  RDN   04  STO 09    05  RDN   06  FS? 00   07  GTO 00   08  6378.172   09  STO 01   10  .07   11  -   12  STO 02           13  21.35   14  -   15  STO 03   16  SIGN   17  STO 00    18  RDN   19  LBL 00   20  XEQ 14   21  STO 05   22  X<>Y   23  STO 06   24  R-P   25  RCL Z   26  STO 07   27  R-P   28  ST/ 05   29  ST/ 06   30  ST/ 07   31  RCL 09   32  RCL 10   33  XEQ 14   34  STO 08   35  X<>Y   36  STO 09   37  R-P   38  RCL Z   39  STO 10   40  R-P   41  ST/ 08   42  ST/ 09   43  ST/ 10   44  LBL 10   45  RCL 06   46  RCL 10   47  *   48  RCL 07   49  RCL 09   50  *   51  -   52  STO 21   53  RCL 07   54  RCL 08   55  *   56  RCL 05   57  RCL 10   58  *   59  -   60  STO 22   61  RCL 05   62  RCL 09   63  *   64  RCL 06   65  RCL 08   66  *   67  -   68  STO 23   69  46.04   70  STO 38   71  RCL 05   72  RCL 08   73  * 74  RCL 06   75  RCL 09   76  *   77  +   78  RCL 07    79  RCL 10    80  *   81  +   82  ACOS   83  ENTER^   84  STO 32           85  SIN   86  D-R   87  STO 31    88  CLX   89   E3   90  /   91  STO 11   92  CLX   93  STO 12   94  XEQ 12   95  STO 34   96  RCL 11   97  CHS   98  STO 11   99  XEQ 12 100  STO 35 101  CLX 102  X<> 11 103  STO 12 104  XEQ 12 105  STO 37 106  RCL 12 107  CHS 108  STO 12 109  XEQ 12 110  STO 36 111  CLX 112  STO 12 113  XEQ 12 114  STO 33 115  RCL 34 116  RCL 35 117  - 118  X^2 119  RCL 34 120  RCL 35 121  + 122  R^ 123  ST+ X 124  - 125  / 126  RCL 36 127  RCL 37 128  - 129  X^2 130  RCL 36 131  R^ 132  ST+ X 133  - 134  RCL 37 135  + 136  / 137  + 138  8 139  / 140  CHS 141  RCL 33 142  + 143  STO 04 144  RTN 145  STO 00 146  GTO 10 147  LBL 12 148  CLX 149  STO 04 150  STO 40  151  RCL 32 152  RCL 00 153  STO 30 154  ST+ X 155  / 156  STO 39         157  LBL 01 158  ST+ 40 159  RCL 38  160  STO 47 161  LBL 02 162  RCL 40 163  STO 48 164  RCL 39 165  RCL IND 47 166  * 167  ST+ 48 168  - 169  XEQ 13 170  X<> 48 171  XEQ 13 172  RCL 48 173  + 174  DSE 47 175  RCL IND 47 176  * 177  ST+ 04 178  DSE 47 179  GTO 02 180  RCL 39 181  ST+ 40 182  DSE 30 183  GTO 01 184  RCL 04 185  * 186  RCL 31 187  * 188  RTN 189  LBL 13 190  STO 13 191  1 192  P-R 193  STO 15 194  RDN 195  STO 14 196  RCL 08 197  * 198  RCL 32 199  RCL 13 200  - 201  1 202  P-R 203  STO 17 204  RDN 205  STO 16 206  RCL 05 207  * 208  + 209  STO 18 210  PI 211  R-D 212  RCL 13 213  * 214  RCL 32 215  / 216  STO Z 217  RCL 11 218  P-R 219  STO 24 220  X<> T 221  ST+ X 222  RCL 12  223  P-R 224  ST+ X 225  ST+ 24 226  RDN 227  + 228  1 229  P-R 230  STO 26         231  ST* Z 232  RDN 233  STO 25 234  RCL 21 235  * 236  + 237  RCL 01 238  / 239  STO 27 240  X^2 241  RCL 06 242  RCL 16 243  * 244  RCL 09 245  RCL 14 246  * 247  + 248  STO 19 249  RCL 26 250  * 251  RCL 22 252  RCL 25 253  * 254  + 255  RCL 02 256  / 257  STO 28 258  X^2 259  + 260  RCL 07 261  RCL 16 262  * 263  RCL 10 264  RCL 14 265  * 266  + 267  STO 20 268  RCL 26 269  * 270  RCL 23 271  RCL 25 272  * 273  + 274  RCL 03 275  / 276  STO 29 277  X^2 278  + 279  X<> 18 280  RCL 25 281  * 282  RCL 21 283  RCL 26 284  * 285  - 286  RCL 27 287  RCL 01 288   / 289  STO 27 290  * 291  RCL 19 292  RCL 25 293  * 294  RCL 22 295  RCL 26  296  * 297  - 298  RCL 28 299  RCL 02 300  / 301  STO 28  302  * 303  + 304  RCL 20         305  RCL 25 306  * 307  RCL 23 308  RCL 26 309  * 310  - 311  RCL 29 312  RCL 03 313  / 314  STO 29 315  * 316  + 317  RCL 24 318  PI 319  * 320  RCL 32 321  / 322  STO 14 323  * 324  STO 13 325  RCL 05 326  RCL 17 327  * 328  RCL 08 329  RCL 15 330  * 331  - 332  RCL 27 333  * 334  RCL 06 335  RCL 17 336  * 337  RCL 09 338  RCL 15 339  * 340  - 341  RCL 28 342  * 343  + 344  RCL 07 345  RCL 17 346  * 347  RCL 10 348  RCL 15 349  * 350  - 351  RCL 29 352  * 353  + 354  RCL 26 355  * 356  RCL 13 357  + 358  RCL 18 359  / 360  X^2 361  RCL 14 362  X^2 363  + 364  RCL 26 365  X^2 366  + 367  RCL 18         368  / 369  SQRT 370  RTN 371  LBL 14 372  DEG 373  HR 374  STO N 375  X<>Y 376  HR 377  14.93 378  FS? 00 379  CLX 380  + 381  STO M 382  1 383  P-R 384  X^2 385  X<>Y 386  RCL 01  387  RCL 02 388  / 389  * 390  X^2 391  + 392  / 393  RCL 01 394  X^2 395  RCL 02 396  X^2 397  STO T 398  - 399  * 400  + 401  RCL N 402  1 403  P-R 404  X^2 405  X<>Y 406  X^2 407  RCL 03 408  X^2 409  R^ 410  STO O 411  / 412  * 413  + 414  / 415  RCL O 416  RCL 03 417  X^2 418  - 419  * 420  RCL 03 421  X^2 422  ST/ O 423  + 424  SQRT 425  RCL N 426  1 427  P-R 428  RCL O 429  * 430  R-P 431  X<> M 432  R^ 433  XEQ "S-R" 434  CLA 435  END

( 617 bytes / SIZE 049 )

 STACK INPUTS OUTPUTS T l1( ° ' " ) / Z b1 ( ° ' " ) / Y l2 ( ° ' " ) / X b2 ( ° ' " ) D ( km )

Example1:       Calculate the geodesic distance between the U.S. Naval Observatory at Washington (D.C.) and the "Observatoire de Paris"

l1 = 77°03'56"0 W ,  b1 = 38°55'17"2 N
l2 = 2°20'13"8 E  ,    b2 = 48°50'11"2 N

CF 00

-77.0356    ENTER^
38.55172  ENTER^
2.20138   ENTER^
48.50112  XEQ "TGD3"  >>>>   D = 6181.626034 km                                             ---Execution time = 5mn16s---

-The error is about 3 centimeters.

Example2:         Calculate the geodesic distance between the U.S. Naval Observatory at Washington (D.C.)

L1 = 77°03'56"0 W ;  b1 = 38°55'17"2 N   and Cape Town Observatory ( South Africa )
L2 = 18°28'35"7 E ;  b2 = 33°56'02"5 S

-77.0356    ENTER^
38.55172   ENTER^
18.28357   ENTER^
-33.56025  XEQ "TGD3"  >>>>    D = 12709.56552 km

-The exact value was  s = 12709.56545 km,  error = 7 cm

Example3:    With ( -40° , -40° ) ( +120° , +50° ) it yields 18095.48516 km

-Exact geodesic distance = 18095.48497 km  -> Error = 19 cm

-Now the errors are always much smaller than 1 meter ( except for nearly antipodal points ).

Notes:

-The elliptic arc length MN is in registers R33 and L

-Gauss-Legendre 6-point formula is used to compute the integral.
-  n = 1 seems sufficient for the Earth.
-With n = 1, Gauss-Legendre 5-point formula would add an extra error of a few decimeters for long geodesics.

-If you prefer the geodetic longitude, replace LBL 14 ( lines 371 to 434 ) by the routine listed in paragraph 2°)a)
with small modifications since the semi-axes a , b , c are already in registers R01-R02-R03.

-If you have the cartesian coordinates  x , y , z  of the 2 points M & N, store these values in R05 to R10
then divide the 3 registers R05-R06-R07 by the norm of the vector OM
and the 3 registers R08-R09-R10 by the norm of the vector ON, then press  XEQ 10

-You can also compute other arc lengths corresponding to different values of  h1 & h2
-Change lines 91 to 112 according to your choices and lines 114 to 143 to exploit these results.

Variants / Improvements:

•  With the geocentric longitudes that are used here, LBL 14 may be simplified:

-Replace lines 404 to 432 by

RCL 03        X<> T        X<> M
X^2              *                 1
ST* Z           R-P

and replace lines 21 to 43 with

STO 05        RCL 09        STO 09
RDN            RCL 10        X<>Y
STO 06        XEQ 14       STO 10
X<>Y          STO 08
STO 07        RDN

since the vectors are already unit vectors.

•  If the inputs are geocentric longitudes and latitudes, the simplifications are even more important:

-Simply replace lines  371 to 435 by

LBL 14      14.93           XEQ "S-R"
DEG          FS? 00         END
HR             CLX
X<>Y        +
HR             1

d)  Program#4

-In this variant, we search  h(µ)  under the form

h(µ) = h1 Sin ( 180° µ / W ) + h2 Sin ( 360° µ / W ) + h3 Sin ( 720° µ / W )

-Only 2 terms of the Fourier series if  CF 01 but 3 terms if  SF 01.
-So, with  SF 01  "TGD4" should return more accurate results.

Data Registers:   R00 = n = number of subintervals of integration        ( Registers R40 to R45 are to be initialized before executing "TGD4" )

R01 = a             R05 = x1             R08 = x2             R26 = W
R02 = b             R06 = y1             R09 = y2             R33 = f(0,0)             R35 to R39:  temp
R03 = c             R07 = z1             R10 = z2              R34 = f(H,0)

>>>>    When the program stops,    R04 = s

•  R40 = 0.4679139346                •  R43 = 0.6612093865
•  R41 = 0.2386191861                •  R44 = 0.1713244924                   ( coefficients for Gauss-Legendre 6-point formula )
•  R42 = 0.3607615730                •  R45 = 0.9324695142

Flags:  F00-F01-F03-F04

CF 00 for the Earth:  n = 1 is sufficient if CF 01 and R01-R02-R03 are initialized   ( lines 08 to 17 )
SF 00 for another ellipsoid:  you'll have to initialize R00-R01-R02-R03

CF 01  only 2 terms in the Fourier series: same results as  "TGD3"
SF 01  3 terms are found

SF 04   geocentric longitudes and latitudes

CF 04 & CF 03  geocentric longitudes &  geodetic latitudes
CF 04 & SF 03   geodetic longitudes and latitudes

Subroutine:  "S-R"  Spherical-Rectangular conversion ( cf "Spherical Coordinates for the HP-41" paragraph 1°-a) )

-Line 122 is a three-byte GTO 10 ( not really important )

 01  LBL "TGD4"   02  STO 08           03  RDN   04  STO 07    05  RDN   06  FS? 00   07  GTO 00   08  6378.172   09  STO 01   10  .07   11  -   12  STO 02   13  21.35   14  -   15  STO 03   16  SIGN   17  STO 00    18  RDN   19  LBL 00   20  XEQ 14   21  STO 05   22  RDN   23  STO 06   24  X<>Y   25  X<> 07   26  RCL 08   27  XEQ 14    28  STO 08   29  RDN    30  STO 09   31  X<>Y   32  STO 10   33  LBL 10    34  RCL 06   35  RCL 10   36  *   37  RCL 07   38  RCL 09   39  *   40  -    41  STO 11   42  RCL 07    43  RCL 08   44  *   45  RCL 05    46  RCL 10   47  *   48  -   49  STO 12   50  RCL 05   51  RCL 09   52  *    53  RCL 06   54  RCL 08   55  *   56  -   57  STO 13   58  45.039   59  STO 36   60  RCL 05   61  RCL 08   62  *   63  RCL 06   64  RCL 09   65  * 66  +   67  RCL 07           68  RCL 10    69  *   70  +   71  ACOS   72  STO 26   73  SIN    74  D-R   75  STO 27    76  CLX   77  STO 16    78  STO 17   79  STO 35   80  XEQ 12    81  STO 33    82  3    83  FC? 01   84  DSE X   85  STO 17   86  RCL 26   87   E3    88  /   89  STO 16   90  LBL 11    91  XEQ 12   92  STO 34   93  SIGN   94  CHS   95  ST* 16    96  XEQ 12   97  STO Y   98  RCL 34   99  ST+ Z 100  - 101  X^2 102  RCL 33 103  ENTER 104  + 105  R^ 106  - 107  X#0? 108  / 109  ST+ 35 110  DSE 17 111  GTO 11 112  RCL 35 113  8 114  ST/ Z 115  / 116  RCL 33 117  ST+ Z 118  + 119  STO 04 120  RTN 121  STO 00 122  GTO 10 123  LBL 12 124  CLX 125  STO 04 126  STO 28 127  RCL 26 128  RCL 00 129  STO 29 130  ST+ X 131  / 132  STO 30         133  LBL 01 134  ST+ 28 135  RCL 36  136  STO 31 137  LBL 02 138  RCL 28 139  STO 32 140  RCL 30 141  RCL IND 31 142  * 143  ST+ 32 144  - 145  XEQ 13 146  X<> 32 147  XEQ 13 148  RCL 32  149  + 150  DSE 31 151  RCL IND 31 152  * 153  ST+ 04 154  DSE 31 155  GTO 02 156  RCL 30 157  ST+ 28 158  DSE 29 159  GTO 01 160  RCL 04 161  * 162  RCL 27 163  * 164  RTN 165  LBL 13 166  STO 18 167  1 168  P-R 169  STO 20 170  RDN 171  STO 19 172  RCL 08 173  * 174  RCL 18 175  RCL 26 176  ST/ 18 177  - 178  1 179  P-R 180  STO 22 181  RDN 182  STO 21 183  RCL 05 184  * 185  - 186  STO 23 187  PI 188  R-D 189  RCL 18 190  * 191  RCL 17 192  * 193  RCL 16 194  P-R 195  RCL 17 196  * 197  PI 198  * 199  RCL 26         200  / 201  STO 18 202  CLX 203  SIGN 204  P-R 205  STO 15 206  ST* Z 207  RDN 208  STO 14 209  ST* 23 210  RCL 11 211  * 212  + 213  RCL 01  214  / 215  STO 37 216  X^2 217  RCL 09 218  RCL 19 219  * 220  RCL 06 221  RCL 21 222  * 223  - 224  STO 24 225  RCL 15 226  * 227  RCL 12 228  RCL 14 229  ST* 24 230  * 231  + 232  RCL 02 233  / 234  STO 38 235  X^2 236  + 237  RCL 10 238  RCL 19 239  * 240  RCL 07 241  RCL 21 242  * 243  - 244  STO 25 245  RCL 15 246  * 247  RCL 13 248  RCL 14 249  ST* 25 250  * 251  + 252  RCL 03 253  / 254  STO 39 255  X^2 256  + 257  RCL 01 258  ST/ 37 259  RCL 02 260  ST/ 38 261  RCL 03  262  ST/ 39 263  R^ 264  X<> 23 265  RCL 11         266  RCL 15 267  * 268  - 269  RCL 37 270  * 271  RCL 24 272  RCL 12 273  RCL 15 274  * 275  - 276  RCL 38  277  * 278  + 279  RCL 25 280  RCL 13 281  RCL 15 282  * 283  - 284  RCL 39 285  * 286  + 287  RCL 18 288  * 289  STO 14 290  RCL 05 291  RCL 22 292  * 293  RCL 08 294  RCL 20 295  * 296  - 297  RCL 37 298  * 299  RCL 06 300  RCL 22 301  * 302  RCL 09 303  RCL 20 304  * 305  - 306  RCL 38 307  * 308  + 309  RCL 07 310  RCL 22 311  * 312  RCL 10 313  RCL 20 314  * 315  - 316  RCL 39 317  * 318  + 319  RCL 15 320  * 321  RCL 14 322  + 323  RCL 23 324  / 325  X^2 326  RCL 18  327  X^2 328  + 329  RCL 15         330  X^2 331  + 332  RCL 23 333  / 334  SQRT 335  RTN 336  LBL 14 337  DEG 338  HR 339  STO N 340  X<>Y 341  HR 342  14.93 343  FS? 00 344  CLX 345  + 346  FS? 04 347  GTO 09 348  FC? 03 349  GTO 08 350  1 351  P-R 352  RCL 01  353  RCL 02 354  / 355  X^2 356  * 357  R-P 358  X<>Y 359  LBL 08 360  STO M 361  1 362  P-R 363  X^2 364  X<>Y 365  RCL 01 366  RCL 02 367  / 368  X^2 369  * 370  X^2 371  + 372  SQRT 373  RCL 03 374  RCL 01 375  / 376  X^2 377  * 378  RCL N 379  1 380  P-R 381  RCL Z 382  / 383  R-P 384  X<> M 385  LBL 09 386  1 387  XEQ "S-R" 388  CLA 389  END

( 554 bytes / SIZE 046 )

 STACK INPUTS OUTPUTS T l1 ( ° ' " ) / Z b1 ( ° ' " ) / Y l2 ( ° ' " ) D1 ( km ) X b2 ( ° ' " ) D2 or D3 ( km ) L / D0 ( km )

CF 01  ---Execution time = 4mn48s---       With  n = 1
SF 01   ---Execution time = 6mn41s---         ---------

Examples 1-2-3:

-The results are almost identical to those given by "TGD3" with or without F01 set

Example4:     On the ellipsoid    x2 / 41 + y2 / 37 + z2 / 35 = 1              the geocentric coordinates of the 2 points are:

l1 =  9°17'35"55660                  l2 = 63°20'07"3683
b1 = -8°11'55"79344                 b2 = 57°46'42"9935

41  SQRT  STO 01
37  SQRT  STO 02
35  SQRT  STO 03

-If you choose  n = 1     1  STO 00

SF 00   CF 01  SF 04

9.173555660  ENTER^
-8.115579344  ENTER^
63.20073683  ENTER^
57.46429935  XEQ "TGD4"  >>>>  D = 8.594824283

-To calculate D with a different n-value, say 2,  simply press  2  R/S   >>>>   8.594824312
-If  n = 4 ,  4  R/S  >>>>   8.594824330

-With  SF 01  ( more accurate )

1   R/S   >>>>   8.594823587
2   R/S   >>>>   8.594823575
4   R/S   >>>>   8.594823592

-In fact, this example is the 1st on this page and we have found the exact value: about  8.594822580 whence an error of  1 E-6

Notes:

-We could replace line 82 by 4 - instead of 3 - to get 4 terms of the Fourier series,
but the method to estimate the coefficients is not perfect and could lead to disappointing results.

-If you replace line 82 by 2 ( instead of 3 ) and if you delete lines 83-84, "TGD4" becomes equivalent to "TGD3",
apart from the fact that "TGD4" is shorter and faster than "TGD3".

-Lines 348 to 359 are useful if the longitudes are geodetic ( set flag F03 ). We have:

Tan Lgeocentric = (b2/a2)  Tan Lgeodetic

Variant:

-Instead of using flag F01, you can initialize a register - say R46 - with the number of terms
that you want to use in the Fourier series and replace lines 82 to 84 by  RCL 46

e)  For The Earth Only

-The integrand contains terms of the form:  u Sin ( W-µ ) + v Cos µ   ....
-They may be re-written  x Cos µ + x' Sin µ
-Moreover, the divisions by a , b , c can be performed at the beginning of the program i-e before the loops

-So, we save bytes and time ( this can be applied to "TGD3" & "TGD4" too )

-Performing the angular functions in RAD mode also seems preferable.

-For the Earth, there is no need to subdivise the interval of integration when Gauss-Legendre 6-point formula is used.
-Registers R00 thru R04 may also be re-used in the loops.
-All these simplifications are applied in "TGD5" hereunder.

Data Registers:   R00 to R37: temp      ( Registers R31 to R36 are to be initialized before executing "TGD5" )

>>>>    When the program stops,    R01 = D

•  R31 = 0.4679139346                •  R34 = 0.6612093865
•  R32 = 0.2386191861                •  R35 = 0.1713244924                   ( coefficients for Gauss-Legendre 6-point formula )
•  R33 = 0.3607615730                •  R36 = 0.9324695142

Flags:  F01-F03-F04

CF 01  2 terms in the Fourier series
SF 01  only 1 term is found

SF 04   geocentric longitudes and latitudes

CF 04 & CF 03  geocentric longitudes &  geodetic latitudes
CF 04 & SF 03   geodetic longitudes and latitudes

Subroutine:  "S-R"  Spherical-Rectangular conversion ( cf "Spherical Coordinates for the HP-41" paragraph 1°-a) )

 01  LBL "TGD5"   02  DEG   03  STO 08           04  RDN   05  STO 07   06  CLX   07  6378.172   08  STO 01   09  .07   10  -   11  STO 02   12  21.35   13  -   14  STO 03    15  CLX   16  14.93   17  STO 00   18  RDN   19  XEQ 05   20  STO 05   21  RDN   22  STO 06   23  X<>Y   24  X<> 07   25  RCL 08   26  XEQ 05   27  STO 08    28  RDN   29  STO 09   30  X<>Y    31  STO 10   32  RCL 06   33  *   34  RCL 07   35  RCL 09    36  *   37  -   38  STO 11   39  RCL 07   40  RCL 08    41  *   42  RCL 05   43  RCL 10    44  *   45  -   46  STO 12   47  RCL 05   48  RCL 09   49  *   50  RCL 06   51  RCL 08   52  *   53  -   54  STO 13   55  36.03   56  STO 37   57  RCL 05    58  RCL 08    59  *   60  RCL 06   61  RCL 09   62  *   63  +   64  RCL 07 65  RCL 10           66  *   67  +   68  STO 14   69  RAD   70  ACOS   71  STO 16   72  SIN   73  STO 15   74  LASTX   75  2   76  /   77  STO 17   78  RCL 05   79  RCL 14   80  *   81  ST- 08   82  RCL 01    83  ST/ 05   84  RCL 15   85  *   86  ST/ 08    87  ST/ 11   88  RCL 06   89  RCL 14   90  *   91  ST- 09   92  RCL 02   93  ST/ 06    94  RCL 15   95  *    96  ST/ 09   97  ST/ 12   98  RCL 07   99  RCL 14 100  * 101  ST- 10 102  RCL 03 103  ST/ 07 104  RCL 15 105  * 106  ST/ 10 107  ST/ 13 108  CLX 109  STO 18 110  STO 19 111  STO 20 112  XEQ 02 113  STO 21 114  2 115  FS? 01 116  SIGN 117  STO 18 118  RCL 16 119   E3 120  / 121  STO 19 122  LBL 01 123  XEQ 02 124  STO 22 125  SIGN 126  CHS 127  ST* 19 128  XEQ 02 129  STO Z 130  RCL 22         131  ST+ T 132  - 133  X^2 134  RCL 21 135  ST+ X 136  R^ 137  - 138  X#0? 139  / 140  ST+ 20 141  DSE 18 142  GTO 01 143  RCL 20 144  8 145  ST/ Z 146  / 147  RCL 21 148  ST+ Z 149  + 150  STO 01 151  DEG 152  RTN 153  LBL 02 154  CLX 155  STO 23 156  RCL 37 157  STO 24 158  LBL 03 159  RCL 17 160  ENTER^ 161  STO 25 162  RCL IND 24 163  * 164  ST- 25 165  + 166  XEQ 04 167  X<> 25 168  XEQ 04 169  RCL 25 170  + 171  DSE 24 172  RCL IND 24 173  * 174  ST+ 23 175  DSE 24 176  GTO 03 177  RCL 17 178  RCL 23 179  * 180  RTN 181  LBL 04 182  STO 14 183  RCL 18 184  RCL 16 185  / 186  STO 15 187  ST* 14 188  SIGN 189  P-R 190  STO 04 191  RCL 05 192  * 193  X<>Y 194  STO 00         195  RCL 08 196  * 197  + 198  STO 01 199  PI 200  RCL 14 201  * 202  RCL 19 203  P-R 204  RCL 15 205  * 206  PI 207  * 208  STO 14 209  CLX 210  SIGN 211  P-R 212  STO 15 213  ST* Z 214  RDN 215  STO 26 216  ST* 01 217  RCL 11 218  * 219  + 220  STO 27 221  X^2 222  RCL 09 223  RCL 00 224  * 225  RCL 06 226  RCL 04 227  * 228  + 229  STO 02 230  RCL 15 231  * 232  RCL 12 233  RCL 26 234  ST* 02 235  * 236  + 237  STO 28 238  X^2 239  + 240  RCL 10 241  RCL 00 242  * 243  RCL 07 244  RCL 04 245  * 246  + 247  STO 03 248  RCL 15 249  * 250  RCL 13 251  RCL 26 252  ST* 03 253  * 254  + 255  STO 29 256  X^2 257  + 258  STO 30         259  RCL 11  260  RCL 15 261  *  262  RCL 01 263  - 264  RCL 27 265  * 266  RCL 12 267  RCL 15 268  * 269  RCL 02 270  - 271  RCL 28 272  * 273  + 274  RCL 13 275  RCL 15 276  * 277  RCL 03 278  - 279  RCL 29 280  * 281  + 282  STO 01 283  RCL 04 284  RCL 08 285  * 286  RCL 00 287  RCL 05 288  * 289  - 290  RCL 27 291  * 292  RCL 04 293  RCL 09 294  * 295  RCL 00 296  RCL 06 297  * 298  - 299  RCL 28 300  * 301  + 302  RCL 04 303  RCL 10 304  * 305  RCL 00 306  RCL 07 307  * 308  - 309  RCL 29 310  * 311  + 312  RCL 15 313  * 314  RCL 01 315  RCL 14 316  * 317  + 318  RCL 30 319 / 320  X^2 321  RCL 14         322  X^2 323  + 324  RCL 15 325  X^2 326  + 327  RCL 30 328  / 329  SQRT 330  RTN 331  LBL 05 332  HR 333  STO 10 334  X<>Y 335  HR 336  RCL 00 337  + 338  FS? 04 339  GTO 07 340  FC? 03 341  GTO 06 342  1 343  P-R 344  RCL 01 345  RCL 02 346  / 347  X^2 348  * 349  R-P 350  X<>Y 351  LBL 06 352  STO 09 353  1 354  P-R 355  X^2 356  X<>Y 357  RCL 01 358  RCL 02 359  / 360  X^2 361  * 362  X^2 363  + 364  SQRT 365  RCL 03 366  RCL 01 367  / 368  X^2 369  * 370  RCL 10 371 1 372  P-R 373  RCL Z 374  / 375  R-P 376  X<> 09 377  LBL 07 378  1 379  XEQ "S-R" 380  END

( 512 bytes / SIZE 038 )

 STACK INPUTS OUTPUTS T l1 ( ° ' " ) / Z b1 ( ° ' " ) / Y l2 ( ° ' " ) D1 ( km ) X b2 ( ° ' " ) D2 or D1 ( km ) L / D0 ( km )

CF 01  ---Execution time = 3mn46s---
SF 01   ---Execution time = 2mn19s---

Example1:       Geodesic distance between the U.S. Naval Observatory at Washington (D.C.) and the "Observatoire de Paris"

l1 = 77°03'56"0 W ,  b1 = 38°55'17"2 N
l2 = 2°20'13"8 E  ,    b2 = 48°50'11"2 N

CF 01  CF 03  CF 04

-77.0356    ENTER^
38.55172  ENTER^
2.20138   ENTER^
48.50112  XEQ "TGD5"  >>>>   D2 = 6181.626035 km = R01
X<>Y  D1 = 6181.626036 km
LASTX   D0 = 6181.628623 km = R21

Example2:      Geodesic distance between the U.S. Naval Observatory at Washington (D.C.)

l1 = 77°03'56"0 W ,  b1 = 38°55'17"2 N   and Cape Town Observatory ( South Africa )
l2 = 18°28'35"7 E  ,  b2 = 33°56'02"5 S

CF 01  CF 03  CF 04

-77.0356    ENTER^
38.55172   ENTER^
18.28357   ENTER^
-33.56025  XEQ "TGD5"  >>>>    D2 = 12709.56553 km = R01
X<>Y   D1 = 12709.56694 km
LASTX   D0 = 12709.56721 km = R21

Example3:    With ( -40° , -40° ) ( +120° , +50° ) it yields:

•  If  CF 01  CF 03  CF 04      geocentric longitudes, geodetic latitudes

D2 = 18095.48517 km = R01
X<>Y   D1 = 18095.49054 km
LASTX   D0 = 18095.54547 km = R21

•  If  CF 01  SF 03  CF 04      geodetic longitudes and latitudes

D2 = 18095.49441 km = R01
X<>Y   D1 = 18095.49978 km
LASTX   D0 = 18095.55469 km = R21

•  If  CF 01  SF 04                  geocentric longitudes and latitudes

D2 = 18099.69164 km = R01
X<>Y   D1 = 18099.69700 km
LASTX   D0 = 18099.75156 km = R21

-If flag F01 is set, D2 is not computed ( faster but less accurate in general ).

Notes:

-The results should verify:  D2 <= D1 <= D0
-However, the formulae don't make the differences between maximum and minimum: they just return an approximate extremum with  D1 & D2
-So, you could add  X>0?  CLX  after line 139

-Lines 378-379 may be replaced by

X<>Y        RCL Z
1                X<>Y
P-R            P-R

thus avoiding the call to the subroutine "S-R" ...  which may be replaced by an M-Code routine S-R !

-The longest geodesic distance is  20003.98600 km  ( provided  a = 6378.172 km , b = 6378.102 km , c = 6356.752 km )
-It is twice the geodesic distance between ( -14°55'48" , 0 ) and ( -14°55'48" , 90° )

-According to other sources, a = 6378.171 km
-Change lines 07 to 14 if you know better evaluations of the semi-axes...

-Despite the title of this paragraph, "TGD5" may also be used for other celestial bodies, provided the flattenings remain small.
-Change lines 07 to 16 according to the characteristics of the triaxial ellipsoid.

Variant:

-The program may of course be simplified if you have an M-Code routine  DOT  ( dot product )
( cf for example "A few M-Code Routines for the HP-41" )

-For instance:

•   Replace lines 283 to 311 by

RCL 10        *                 RCL 06      RCL 04       -
RCL 04        -                 RCL 00       *                 DOT
*                  RCL 09       *                 RCL 05
RCL 07       RCL 04       -                  ST* 00
RCL 00       *                  RCL 08       X<> 00

•   Replace lines 259 to 281 by

RCL 13       -                   RCL 02       *
RCL 15       RCL 12        -                  RCL 01
*                 RCL 15        RCL 11        -
RCL 03       *                  RCL 15        DOT

•   Replace lines 220-237-255  with  STO M   STO N   STO O   respectively.

•   Add  CLA  after line 151

-The first lines of the program may be simplified further if you have also an M-Code routine  CROSS  ( cross product )

Remarks:

-The precision of these approximate methods seems quite acceptable, except for nearly antipodal points M & N.
-But in this case, we can always add an intermediate point P and find the minimum
of  MP + PN  by a golden section search - though it will be very slow with a real HP-41...
-Of course, to get very accurate results, Jacobi's method remains the best !

-This new version employs Gauss-Legendre 10-point integration formula

-It uses the following semi-axes:

a =     6378.172   km     with  the equatorial major axis in the direction   14°92911W
b =     6378.102   km     equatorial minor semi-axis
c = 6356.752314 km     polar minor semi-axis          ( like WGS84 )

Data Registers:   R00 to R49: temp
Flags:  /
Subroutines: /

-With an HP41, round the coefficients in lines 59-61-....-77 to 10 decimals.

 01 LBL "TGDX"  02 DEG  03 STO 08  04 RDN  05 STO 07  06 RDN  07 STO 06  08 X<>Y  09 STO 05  10 6378.172  11 STO 01  12 6378.102  13 STO 02  14 6356.752314  15 STO 03  16 14.92911  17 CHS  18 STO 00  19 RCL 05  20 RCL 06                  21 XEQ 05  22 STO 05  23 RDN  24 STO 06  25 X<>Y  26 X<> 07  27 RCL 08  28 XEQ 05  29 STO 08  30 RDN  31 STO 09  32 X<>Y  33 STO 10  34 RCL 06  35 *  36 RCL 07  37 RCL 09  38 *  39 -  40 STO 11  41 RCL 07  42 RCL 08  43 *  44 RCL 05  45 RCL 10        46 *  47 -  48 STO 12  49 RCL 05  50 RCL 09  51 *  52 RCL 06  53 RCL 08  54 *  55 -  56 STO 13  57 40.03  58 STO 49 59 0.295524224715  60 STO 31  61 0.148874338982  62 STO 32  63 0.26926671931  64 STO 33  65 0.433395394129  66 STO 34  67 0.219086362516  68 STO 35  69 0.679409568299  70 STO 36  71 0.149451349151  72 STO 37  73 0.865063366689  74 STO 38  75 0.0666713443087  76 STO 39  77 0.973906528517  78 STO 40                  79 RCL 05  80 RCL 08  81 *  82 RCL 06  83 RCL 09  84 *  85 +  86 RCL 07  87 RCL 10  88 *  89 +  90 STO 14  91 RAD  92 ACOS  93 STO 16  94 SIN  95 STO 15  96 LASTX  97 2  98 /  99 STO 17 100 RCL 05 101 RCL 14 102 * 103 ST- 08 104 RCL 01 105 ST/ 05 106 RCL 15 107 * 108 ST/ 08 109 ST/ 11 110 RCL 06 111 RCL 14 112 * 113 ST- 09 114 RCL 02 115 ST/ 06 116 RCL 15 117 * 118 ST/ 09 119 ST/ 12 120 RCL 07 121 RCL 14 122 * 123 ST- 10 124 RCL 03 125 ST/ 07 126 RCL 15 127 * 128 ST/ 10 129 ST/ 13 130 CLX 131 STO 18 132 STO 19 133 XEQ 02 134 STO 21 135 4 136 STO 18                 137 44 138 STO 20 139 RCL 16 140 1 E3 141 / 142 STO 19 143 GTO 01 144 LBL 02 145 CLX 146 STO 23 147 RCL 49 148 STO 24 149 LBL 03 150 RCL 17 151 ENTER 152 STO 25 153 RCL IND 24 154 * 155 ST- 25 156 + 157 XEQ 04 158 X<> 25 159 XEQ 04 160 RCL 25 161 + 162 DSE 24 163 RCL IND 24 164 * 165 ST+ 23 166 DSE 24 167 GTO 03 168 RCL 17 169 RCL 23 170 * 171 RTN 172 LBL 04 173 STO 14 174 RCL 18 175 RCL 16 176 / 177 STO 15 178 ST* 14 179 SIGN 180 P-R 181 STO 04 182 RCL 05 183 * 184 X<>Y 185 STO 48 186 RCL 08 187 * 188 + 189 STO 45 190 PI 191 RCL 14 192 * 193 RCL 19 194 P-R 195 RCL 15                 196 * 197 PI 198 * 199 STO 14 200 CLX 201 SIGN 202 P-R 203 STO 15 204 ST* Z 205 RDN 206 STO 26 207 ST* 45 208 RCL 11 209 * 210 + 211 STO 27 212 X^2 213 RCL 09 214 RCL 48 215 * 216 RCL 06 217 RCL 04 218 * 219 + 220 STO 46 221 RCL 15 222 * 223 RCL 12 224 RCL 26 225 ST* 46 226 * 227 + 228 STO 28 229 X^2 230 + 231 RCL 10 232 RCL 48 233 * 234 RCL 07 235 RCL 04 236 * 237 + 238 STO 47 239 RCL 15 240 * 241 RCL 13 242 RCL 26 243 ST* 47 244 * 245 + 246 STO 29 247 X^2 248 + 249 STO 30 250 RCL 11 251 RCL 15 252 * 253 RCL 45                 254 - 255 RCL 27 256 * 257 RCL 12 258 RCL 15 259 * 260 RCL 46 261 - 262 RCL 28 263 * 264 + 265 RCL 13 266 RCL 15 267 * 268 RCL 47 269 - 270 RCL 29 271 * 272 + 273 STO 45 274 RCL 04 275 RCL 08 276 * 277 RCL 05 278 RCL 48 279 * 280 - 281 RCL 27 282 * 283 RCL 04 284 RCL 09 285 * 286 RCL 06 287 RCL 48 288 * 289 - 290 RCL 28 291 * 292 + 293 RCL 04 294 RCL 10 295 * 296 RCL 07 297 RCL 48 298 * 299 - 300 RCL 29 301 * 302 + 303 RCL 15 304 * 305 RCL 14 306 RCL 45 307 * 308 + 309 RCL 30 310 / 311 X^2 312 RCL 14                 313 X^2 314 + 315 RCL 15 316 X^2 317 + 318 RCL 30 319 / 320 SQRT 321 RTN 322 LBL 05 323 HR 324 X<>Y 325 HR 326 RCL 00 327 - 328 1 329 X<>Y 330 RDN 331 P-R 332 R^ 333 X<>Y 334 P-R 335 RCL 01 336 X^2 337 * 338 X^2 339 STO 35 340 X<> L 341 X<>Y 342 RCL 02 343 X^2 344 * 345 X^2 346 ST+ 35 347 X<> L 348 R^ 349 RCL 03 350 X^2 351 * 352 X^2 353 ST+ 35 354 X<> L 355 X<> Z 356 RCL 35 357 SQRT 358 ST/ T 359 ST/ Z 360 / 361 RTN 362 LBL 07 363 1 364 X<>Y 365 RDN 366 P-R 367 R^ 368 X<>Y 369 P-R 370 RTN 371 LBL 01 372 XEQ 02 373 STO 22                 374 SIGN 375 CHS 376 ST* 19 377 XEQ 02 378 STO Z 379 RCL 22 380 ST+ T 381 - 382 X^2 383 RCL 21 384 ST+ X 385 R^ 386 - 387 8 388 * 389 X#0? 390 / 391 STO IND 20 392 DSE 20 393 DSE 18 394 GTO 01 395 DEG 396 RCL 21 397 ST+ 41 398 SIGN 399 RCL 41 400 ST+ 42 401 RCL 42 402 ST+ 43 403 RCL 43 404 ST+ 44 405 RCL 41 406 RCL 42 407 RCL 43 408 RCL 44 409 END

( SIZE 050 )

 STACK INPUTS OUTPUT T L ( ° ' " ) / Z b ( ° ' " ) / Y L' ( ° ' " ) / X b' ( ° ' " ) D ( km )

where  L , b , L' , b'  are geodetic coordinates

Example:      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   ENTER^
-116.51504  ENTER^
33.21224   XEQ "TGDX"   >>>>   D = 12138.657569 km

-The exact result ( computed with Jacobi's method ) is D = 12138.657551942....

Note:

-I've used free42 to compute these values.

References:

[1]  Jacobi - "De la ligne géodésique sur un ellipsoïde, et des différents usages d’une transformation analytique remarquable"
[2]  J. Feltens - Vector method to compute the Cartesian (X, Y, Z) to geodetic (f, l, h) transformation on a triaxial ellipsoid - J Geod (2009) 83:129–137
[3]  Marcin Ligas  - Cartesian to geodetic coordinates conversion on a triaxial ellipsoid - J Geod (2012) 86:249–256

-With many thanks to Charles Karney for the link that he sent me to get Jacobi's method:
http://lists.maptools.org/pipermail/proj/2012-October/006448.html

[4]  If you want to visualize the geodesics on a triaxial ellipsoid, go to this excellent page