hp41programs

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
 

-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<Y?
276  GTO 10
277  FS? 09
278  SF 41
279  RCL 31
280  FS? 02
281  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<Y?
345  GTO 10
346  FS? 09
347  SF 41
348  RCL 30 
349  STO 12 
350  FS? 02
351  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 !
 

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