The Museum of HP Calculators

Terrestrial Geodesic Distance for the HP-41C

This program is Copyright © 1999-2007 by Jean-Marc Baillard and may be freely used by permission.

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


Overview
 

   1°) Andoyer's Formulae

     a) Geodesic distance only

      a-1) Program#1 ( First order accuracy )
      a-2) Program#2 ( Second order accuracy )  ( new )

     b) Geodesic distance and azimuths

      b-1) Program#1 ( new )
      b-2) Program#2 ( new )

   2°) Vincenty's Formulae

     a) Geodesic distance and azimuths  ( improved )
     b) Forward Geodesic Problem  ( new )

   3°) Numerical Integration
   4°) A Series Expansion
   5°) Euclidean Distance

    a) Geocentric Coordinates ( new )
    b) 3D-Distance  ( new )
 

-The programs listed in paragraphs 1 to 4 solve geodesic problems  between 2 points.
-These points are defined by their geographic coordinates:  ( L1 , b1 )  ( L2 , b2 )   where   Li = Longitudes ,  bi = latitudes.
-The geodesic distance D between these 2 points P1 & P2 ( on the Earth's surface , at mean sea level ) is expressed in kilometers.
-The azimuths  are expressed in  ° ' "
-An ellipsoidal model of the Earth is used with

  a = equatorial radius = 6378.137 km     ( DE403 )
  f  = Earth's flattening = 1/298.257       ( IERS 1992 )

-If you only want to compute the geodesic distance, the longitudes can be measured positively westward or positively eastward from the meridian of Greenwich.

-Otherwise, the longitudes are measured positively Eastwards and the azimuths are reckoned clockwise from North.

-In paragraph 5, the observers' heights are taken into account to compute the Euclidean 3D-distance.
 

1°) Andoyer's Formulae
 

     a) Geodesic Distance only
 

       a-1) Program#1 ( First order accuracy )
 

-Let                    L = ( L2 - L1 )/2  ;   F = ( b1 + b2 )/2  ;  G = ( b2 - b1 )/2
-We calculate     S = sin2 G  cos2 L  +  cos2 F  sin2 L  ;   µ = Arcsin S1/2  ;  R = ( S.(1-S) )1/2

-Then,                D = a.{ 2µ + f.[ (3R-µ)/(1-S)  sin2 F  cos2 G  - (3R+µ)/S  sin2 G  cos2 F ] }
 

Data Registers: /
Flags: /
Subroutines: /

-Synthetic register M may be replaced by any data register like R00 ...
 

01  LBL "TGD"
02  DEG
03  HR
04  X<> T
05  HMS-
06  HR
07  2
08  /
09  SIN
10  X^2
11  RDN
12  HR
13  +
14  2
15  /
16  ST- Y
17  COS
18  X^2
19  ST* T
20  RDN
21  SIN
22  X^2
23  STO M
24  ST* Y
25  -
26  -
27  ENTER^
28  SQRT
29  ASIN
30  D-R
31  RCL M
32  R^
33  ST* M
34  +
35  RCL M
36  -
37  1
38  ST- Y
39  R^
40  ST/ M
41  -
42  ST/ Y
43  LASTX
44  *
45  SQRT
46  3
47  *
48  R^
49  ST+ T
50  -
51  ST* Y
52  R^
53  +
54  0
55  X<> M
56  *
57  +
58  298.257
59  /
60  -
61  6378.137
62  *
63  END

( 98 bytes / SIZE 000 )
 
 
      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.)    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

   77.0356    ENTER^
   38.55172  ENTER^
   -2.20138   ENTER^
   48.50112  XEQ "TGD"  >>>>   D = 6181.6156 km   ( execution time = 5 seconds )

-In this example, the error is about  6 meters.
-If the points are not ( nearly ) antipodal, the relative error is less than  10 -5 and the maximum error is of the order of 60 meters in absolute value:
 with ( 0° , -40° ) ( 120° , 40° )   the error = 66 meters
-However, errors can reach several kilometers for nearly antipodal points.
-If the 2 points are on the equator:  ( L1 , 0 ) ( L2 , 0 ) , Andoyer's formula gives a perfect accuracy provided  |  L2 - L1 |  <  179°23'47"377
 

        a-2) Program#2 ( Second order accuracy )
 

-This program uses a formula which is second-order in the flattening f:
 

-With   L = ( L2 - L1 )/2  ;   F = ( b1 + b2 )/2  ;  G = ( b2 - b1 )/2
           S = sin2 G  cos2 L  +  cos2 F  sin2 L  ;   µ = Arcsin S1/2  ;  R = ( S.(1-S) )1/2

 we define  A = ( sin2 G  cos2 F )/S + ( sin2 F  cos2 G )/( 1-S )
         and   B = ( sin2 G  cos2 F )/S - ( sin2 F  cos2 G )/( 1-S )

-The geodesic distance is then   D = 2a.µ ( 1 + eps )

  where   eps =    (f/2) ( -A - 3B R/µ )
                       + (f2/4) { { - ( µ/R + 6R/µ ).B + [ -15/4 + ( 2S - 1 ) ( µ/R + (15/4) R/µ ) ] A + 4 - ( 2S - 1 ) µ/R } A
                                      - [ (15/2) ( 2S - 1 ).B  R/µ - ( µ/R + 6R/µ ) ] B }
 

Data Registers:   R00 = µ/R    R01 = A    R02 = B    R03 = 2S-1    R04 = 2/f = 596.514
Flags: /
Subroutines: /
 

  01  LBL "TGD1"
  02  DEG
  03  HR
  04  X<> T
  05  HMS-
  06  HR
  07  2
  08  /
  09  SIN
  10  X^2
  11  RDN
  12  HR
  13  +
  14  2
  15  /
  16  ST- Y
  17  1
  18  P-R
  19  X^2
  20  STO 02
  21  X<>Y
  22  X^2
  23  STO 01
  24  RDN
  25  X<>Y
  26  1
  27  P-R
  28  X^2
  29  ST* 01
  30  RDN
  31  X^2
  32  ST* 02
  33  STO T
  34  -
  35  *
  36  +
  37  ST/ 02
  38  STO 03
  39  SQRT
  40  ASIN
  41  D-R
  42  STO 00
  43  RCL 03
  44  1
  45  RCL 03
  46  -
  47  ST/ 01
  48  ST- 03
  49  *
  50  SQRT
  51  ST/ 00
  52  CLX
  53  RCL 02
  54  RCL 01
  55  ST- 02
  56  +
  57  STO 01
  58  3.75
  59  STO 04
  60  RCL 00
  61  ST/ Y
  62  +
  63  RCL 03
  64  *
  65  RCL 04
  66  -
  67  *
  68  6
  69  RCL 00
  70  ST/ Y
  71  +
  72  STO 04
  73  RCL 02
  74  *
  75  -
  76  4
  77  +
  78  RCL 00
  79  RCL 03
  80  *
  81  -
  82  RCL 01
  83  *
  84  RCL 02
  85  RCL 03
  86  *
  87  7.5
  88  *
  89  RCL 00
  90  /
  91  RCL 04
  92  -
  93  RCL 02
  94  *
  95  -
  96  596.514
  97  STO 04
  98  /
  99  RCL 01
100  -
101  RCL 02
102  3
103  *
104  RCL 00
105  /
106  -
107  RCL 04
108  /
109  *
110  +
111  12756.274
112  *
113  END

( 153 bytes / SIZE 005 )
 
 
      STACK        INPUTS      OUTPUTS
           T       L1 ( ° ' " )             /
           Z       b1 ( ° ' " )             /
           Y       L2 ( ° ' " )             /
           X       b2 ( ° ' " )        D ( km )

Example:

   77.0356    ENTER^
   38.55172  ENTER^
   -2.20138   ENTER^
   48.50112  XEQ "TGD1"  >>>>   D = 6181.621809 km   ( execution time = 7.6 seconds )

-The error is only 15 mm!
-The accuracy is very good ( less than 1 meter ), provided D is not too large:
-Otherwise, the errors may be greater: for example, ( 0° , 0° ) ( 179° , 1° )  give 19860.3438 km whereas the correct result is 19860.5092 km
 

     b) Geodesic Distance and Azimuths
 

-With the same notations as above, the forward azimuth Az1 in the direction from P1 towards P2 is obtained by the following formulae:    Az1 = Az1sph + dAz

      where  dAz = (f/2).( cos2 F - sin2 G ).[ ( 1+µ/R ).( sin G cos F )/S + ( 1-µ/R ).( sin F cos G )/(1-S) ].( sin 2L )

       and    Tan Az1sph = cos b2 sin 2L / ( cos b1 sin b2 - sin b1 cos b2 cos 2L )     >>>  the R-P function returns the angle in the proper quadrant

  Az1sph is the forward azimuth in a spherical model of the Earth,
  dAz is a correction that takes the Earth's flattening f into account.

-The back azimuth ( from P2 to P1 ) is obtained by replacing L by -L and exchanging b1 & b2 in the formulas above.
 

         b-1) Program#1
 

-This routine computes D with the 1st-order formula used by "TGD"
 

Data Registers:   R00 = 2L , µ , µ/R         R03 = ( sin F cos G )/(1-S)               R06 = (f/2).( cos2 F - sin2 G ).( sin 2L )
                              R01 = Az1sph               R04 =  ( sin G cos F )/S , temp          R07 = 1 - S
                              R02 = Az2sph               R05 = S
Flags: /
Subroutines: /
 

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

( 173 bytes / SIZE 008 )
 
 
      STACK        INPUTS      OUTPUTS
           T       L1 ( ° ' " )      Az2 ( ° ' " )
           Z       b1 ( ° ' " )      Az2 ( ° ' " )
           Y       L2 ( ° ' " )      Az1 ( ° ' " )
           X       b2 ( ° ' " )        D ( km )

Example:

   -77.0356    ENTER^
    38.55172  ENTER^
     2.20138   ENTER^
    48.50112  XEQ "TGD+"  >>>>    D   =   6181.6156 km                              ( execution time = 13 seconds )
                                            RDN    Az1 =  51°47'36"76 = forward azimuth      Az1-error = -0"05
                                            RDN    Az2 = -68°09'59"26 = back azimuth          Az2-error = -0"29

Notes:

-Like "TGD" , "TGD+" gives accurate results if the 2 points are not nearly antipodal.
-Az is obtained with an accuracy of the order of a few arcseconds ( the error is smaller than 1" in this example ).
-If the distance doesn't exceed 17000 km, the maximum errors I've seen is 66 meters for D and 20" for the azimuths.

-Otherwise, the results may still be accurate:       with ( 0° , -60° ) ( 170° , 70° )   D-error = 13m ,    Az1-error = 2" ,       Az2-error = 3"
  but in other examples, it may be disappointing:  with ( 0° , 0° )    ( 179° , 1° )     D-error = 3303m , Az1-error = 49'06" , Az2-error = 49'09"
 
 

         b-2) Program#2
 

-This program computes D with the 2nd-order formula used by "TGD1"
 

Data Registers:   R00 = 2L , µ , µ/R         R03 = ( sin F cos G )/(1-S)               R06 = (f/2).( cos2 F - sin2 G ).( sin 2L )          R08 = A
                              R01 = Az1sph               R04 =  ( sin G cos F )/S , temp          R07 = 1 - S , temp , 2/f = 596.514                 R09 = B
                              R02 = Az2sph               R05 = S , 2S-1
Flags: /
Subroutines: /
 

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

( 230 bytes / SIZE 010 )
 
 
      STACK        INPUTS      OUTPUTS
           T       L1 ( ° ' " )      Az2 ( ° ' " )
           Z       b1 ( ° ' " )      Az2 ( ° ' " )
           Y       L2 ( ° ' " )      Az1 ( ° ' " )
           X       b2 ( ° ' " )        D ( km )

Example:

   -77.0356    ENTER^
    38.55172  ENTER^
     2.20138   ENTER^
    48.50112  XEQ "TGD1+"  >>>>    D   =   6181.621809 km                              ( execution time = 15 seconds )
                                              RDN    Az1 =  51°47'36"76 = forward azimuth
                                              RDN    Az2 = -68°09'59"26 = back azimuth

Notes:

-Andoyer's formulas have one advantage: they are non-iterative!
-In most cases, their accuracy is satisfactory, all the more that the Earth's irregularities cannot be taken into account.
-If, however, a better precision is required, this can be achieved via Vincenty's formulas below.
 

2°) Vincenty's Formulae
 

     a) Geodesic Distance and Azimuths
 

-Let   L =  L2 - L1  = L0   ;   U = Arctan((1-f).tan b1)  ;   V = Arctan((1-f).tan b2)

-We compute recursively:      Ln+1 = L + (1-C).f.(sin µ).{ c + C.sin c [ cos 2d + C.cos c ( -1 + 2.cos2 2d ) ] }

   where      sin c = [ ( cos V sin Ln )2 + ( cos U sin V - sin U cos V cos Ln )2 ]1/2
                 cos c = sin U sin V + cos U cos V cos Ln
                       µ = Arcsin [ ( cos U cos V sin Ln )/sin c ]  ;  µ = azimuth of the geodesic at the equator
                      C = (f/16).(cos2 µ ).[ 4 + f.(4-3.cos2 µ ) ]
              cos 2d = cos c - 2.sin U sin V / cos2 µ

   until   Ln+1 = Ln = lambda

-Then,     D = a.A.(1-f).(c-D)

    with     A = 1 + (u/16384).{ 4096 + u.[-768 + u.(320-175.u)] }      ;       B = (u/1024).{ 256 + u.[ -128 + u.(74-47.u)] }
               D = (B.sin c).{ cos 2d + (B/4). [ ( cos c ).( -1 + 2.cos2 2d ) - (B/6).(cos 2d).(-3+4.sin2 c ).(-3+4.cos2 2d )] }
                u =  f(2-f)/(1-f)2 . cos2 µ

-Actually, the HP-41 works with 10 digits only. So several terms may be omitted without reducing the accuracy,  namely the terms in  u4 and  B3
-The azimuths  Az1 &  Az2 are finally calculated by the formulas:

       Tan Az1 = [ cos V sin (lambda) ] / [ cos U sin V - sin U cos V cos (lambda) ]
       Tan Az2 = [ -cos U sin (lambda) ] / [ sin U cos V - cos U sin V cos (lambda) ]
 

Data Registers:  R00 = f  ;  R01 = U  ;  R02 = V  ;  R04 = lambda  ;  R06 = µ = Azequator  ;   R03 , R05 and R07 to R09: temp
Flags:  /
Subroutines: /
 

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

( 289 bytes / SIZE 010 )
 
 
      STACK        INPUTS      OUTPUTS
           T       L1 ( ° ' " )     Az2 ( ° ' " )
           Z       b1 ( ° ' " )     Az2 ( ° ' " )
           Y       L2 ( ° ' " )     Az1 ( ° ' " )
           X       b2 ( ° ' " )        D ( km )

 where   Az1  =  forward azimth  ( from P1 to P2 )
   and     Az2  =  back azimuth  ( from P2 to P1 )

Example:

   -77.0356    ENTER^
    38.55172   ENTER^
      2.20138   ENTER^
   48.50112  XEQ "TGD2"  >>>>   the successive Ln-values are displayed  and  ( after 47 seconds )

                                           D   =  6181.621787 km
                             RDN    Az1 =  51°47'36"8132
                             RDN    Az2 = -68°09'58"9656

-We also have  R06 = µ = azimuth at the equator = 31.74571879°

Notes:

-The accuracy is excellent, provided the 2 points are not nearly antipodal.
-Otherwise, the algorithm doesn't converge! It happens when | Ln | exceeds 180° ( in register R04 )
-If you don't want to compute the azimuths, replace lines 193 to 229 by  RCL 09  *
 

     b) Forward Geodesic Problem
 

-Given the coordinates ( L1 , b1 ) of the first point  P1 , the forward azimuth Az1 and the length D of the geodesic,
 "FWD" computes the longitude and the latitude ( L2 , b2 ) of the second point  P2

-The formulas are given in the reference [3]  ( the terms in  u4 and  B3 have been neglected )
 

Data Registers:   R00 thru R13
Flags: /
Subroutines: /
 

  01  LBL "FWD"
  02  DEG
  03  R-D
  04  STO 01
  05  RDN
  06  HR
  07  STO 02
  08  RDN
  09  HR
  10  X<>Y
  11  HR
  12  STO 03
  13  CLX
  14  SIGN
  15  P-R
  16  LASTX
  17  298.257
  18  1/X
  19  STO 00
  20  -
  21  STO 09
  22  ST/ 01
  23  /
  24  R-P
  25  X<>Y
  26  STO 04
  27  1
  28  P-R
  29  STO 13
  30  RCL 02
  31  COS
  32  *
  33  R-P
  34  X<>Y
  35  STO 05
  36  6378.137
  37  ST/ 01
  38  RCL 13
  39  RCL 02
  40  SIN
  41  *
  42  STO 06
  43  ASIN
  44  2
  45  RCL 00
  46  ST- Y
  47  *
  48  X<>Y
  49  COS
  50  STO 07
  51  RCL 09
  52  /
  53  X^2
  54  *
  55  12
  56  RCL Y
  57  5
  58  *
  59  X<>Y
  60  -
  61  *
  62  64
  63  +
  64  *
  65  256
  66  ST+ Y
  67  /
  68  ST/ 01
  69  CLX
  70  37
  71  *
  72  64
  73  -
  74  *
  75  512
  76  /
  77  4
  78  1/X
  79  +
  80  *
  81  STO 08
  82  RCL 01
  83  STO 10
  84  LBL 01
  85  VIEW 10
  86  RCL 10
  87  RCL 05
  88  ST+ X
  89  +
  90  COS
  91  STO 11
  92  X^2
  93  ST+ X
  94  RCL 10
  95  COS
  96  ST* Y
  97  -
  98  STO 12
  99  RCL 08
100  *
101  4
102  /
103  RCL 11
104  +
105  RCL 10
106  SIN
107  *
108  RCL 08
109  *
110  R-D
111  RCL 01
112  +
113  ENTER^
114  X<> 10
115  -
116  ABS
117   E-7
118  X<Y?
119  GTO 01
120  RCL 04
121  SIN
122  STO 04
123  RCL 10
124  COS
125  STO 01
126  *
127  RCL 13
128  ST* 01
129  RCL 10
130  SIN
131  STO 05
132  *
133  RCL 02
134  COS
135  STO 13
136  *
137  +
138  RCL 04
139  RCL 05
140  *
141  RCL 01
142  RCL 13
143  *
144  -
145  X^2                     One byte can be saved if you replace lines 145 to 149 by   RCL 06   R-P   X<>Y   CLX  but it is a little slower...
146  RCL 06
147  X^2
148  +
149  SQRT
150  RCL 09
151  *
152  R-P
153  X<>Y
154  HMS
155  RCL 02
156  RCL 05
157  P-R
158  RCL 04
159  *
160  RCL 01
161  X<>Y
162  -
163  R-P
164  CLX
165  RCL 07
166  X^2
167  STO 07
168  3
169  *
170  4
171  -
172  RCL 00
173  ST* 06
174  ST* 07
175  *
176  4
177  -
178  RCL 07
179  *
180  16
181  /
182  STO 08
183  ST* 05
184  RCL 12
185  *
186  RCL 11
187  -
188  RCL 05
189  *
190  R-D
191  RCL 10
192  +
193  RCL 06
194  *
195  ST- Y
196  RCL 08
197  *
198  -
199  RCL 03
200  +
201  360                     Lines 201 to 208 are only useful to return a longitude between   -180° and +180°
202  MOD
203  PI
204  R-D
205  X>Y?
206  CLX
207  ST+ X
208  -
209  HMS
210  CLD
211  END

( 263 bytes / SIZE 014 )
 
 
      STACK        INPUTS      OUTPUTS
           T       L1 ( ° ' " )       b2 ( ° ' " )
           Z       b1 ( ° ' " )       b2 ( ° ' " )
           Y     Az1 ( ° ' " )       b2 ( ° ' " )
           X        D ( km )       L2 ( ° ' " )

 where   Az1  is the forward azimuth in the direction from the first point P1 towards the second point P2

Example:     L1 = 10°30'     b1 = 49°41'     Az1 = 12°24'     D = 16000 km       Calculate  L2 & b2

   10.30   ENTER^
   49.41   ENTER^
   12.24   ENTER^
  16000   XEQ "FWD"  >>>  the successive sigma-approximations are displayed and eventually, it yields:

               L2 = -177°03'07"987     X<>Y    b2 =  -14°06'40"7154       ( execution time = 25 seconds )

-This method works for antipodal points too and - more generally - for any length D !
 

3°) Numerical Integration
 

-The rigorous formulas are:

       D/a = s = §r1r2  r [ ( 1-e2.r2 )/( 1 - r2 )/ ( r2 - k2 ) ] 1/2 dr                                                             ( § = Integral )

     where  e = f.(2-f) = the eccentricity of the ellipsoid,
                (a.r)  is the radius of the parallel of latitude b ,   r = ( cos b ).( 1 - e2.sin2 b ) -1/2

  and  k  is the solution of    | L1 - L2 | =  §r1r2  (k/r) [ ( 1-e2.r2 )/( 1 - r2 )/ ( r2 - k2 ) ] 1/2 dr

-Actually,  a given geodesic cuts all the parallels at an angle phi such that   r.cos(phi) = k = Constant
-In order to remove the singularities of these integrals, I've made the following change of argument:

      sin 2µ = (2/( 1 - k2 )).[ ( 1 - r2 ).( r2 - k2 ) ] 1/2    or   sin µ = +/- [ ( 1 - r2 )/( 1 - k2) ] 1/2   ;    cos µ = +/- [ ( r2 - k2 )/( 1 - k2 ) ] 1/2

-And the integrals become:

            s = §µ1µ2  ( 1 - e2.( cos2 µ + k2.sin2 µ ) )1/2 dµ            ( E )

    | L1 - L2 | =  §µ1µ2  ( k/( cos2 µ + k2.sin2 µ ) ).( 1 - e2.( cos2 µ + k2.sin2 µ ) )1/2 dµ

-However, this last integral is difficult if the geodesic passes near the Pole, so I've split it into 2 parts which eventually yields:

   | L1 - L2 | =  [ Arctan ( k.tan µ ) ]µ1µ2 -   §µ1µ2  k.e2/[ 1 + ( 1 - e2.( cos2 µ + k2.sin2 µ ) )1/2 ] dµ     ( E' )

-"TGD3"  solves equation (E') by the secant method, thus giving k , µ1 , µ2
-And Gaussian integration produces  s   through equation (E).
 

Data Registers:  R00 thru R08 are used by "GL3"    R00 = "S" ;  R01 = µ1 ;  R02 = µ2
                             R03 = 2  for solving (E')      ( line 64 )              ( number of subintervals over which
                             R03 = 4  for evaluating  s    ( line 191 )               the Gauss-Legendre formula is applied )

                          These values produce ( almost ) 10 significant figures for the Earth.
                          They should be increased if you want to use this program for another planet with a greater flattening.

                             R09 = -e2  ;  R10 = r1 ( or -r1 if µ2 > 90° )  ;  R11 = sign(b1).(1-r1)1/2  ;  R12 = r2 ( calling r2 < r1 ) ;  R13 = (1-r2)1/2
                             R14 & R15 = the successive k-values  ;  R16 = f(k) ;  R17 = | L1 - L2 |
Flag:  F10
Subroutine:  "GL3"  Gauss-Legendre 3-point formula ( cf "Numerical Integration for the HP-41" )
 

  01  LBL "TGD3"
  02  DEG
  03  HR
  04  X<> T
  05  HMS-
  06  HR
  07  ABS
  08  PI
  09  R-D
  10  X>Y?
  11  CLX
  12  ST+ X
  13  -
  14  ABS
  15  STO 17
  16  RDN
  17  HR
  18  ST* Z
  19  ABS
  20  X<>Y
  21  ABS
  22  X>Y?
  23  X<>Y
  24  RCL Z
  25  SIGN
  26  *
  27  COS
  28  LASTX
  29  SIN
  30  STO 11
  31  X^2
  32  .006694385       ( e2  with  f = 1/298.257 )
  33  CHS
  34  STO 09
  35  *
  36  1
  37  +
  38  SQRT
  39  ST/ 11
  40  /
  41  STO 10
  42  X<>Y
  43  COS
  44  LASTX
  45  SIN
  46  STO 13
  47  X^2
  48  RCL 09
  49  *
  50  1
  51  +
  52  SQRT
  53  ST/ 13
  54  /
  55  STO 12
  56  STO 14
  57  STO 15
  58  SIGN
  59  RCL 09
  60  +
  61  SQRT
  62  ST* 11
  63  ST* 13
  64  2
  65  STO 03
  66  "S"
  67  ASTO 00
  68  RCL 13
  69  X#0?
  70  GTO 00
  71  179.3964936       if  | L1 - L2 | < 179°23'47"377 ,  s = | L1 - L2 | radians
  72  RCL 17
  73  X<=Y?
  74  GTO 04
  75  LBL 00
  76  SF 10
  77  XEQ 01
  78  STO 16
  79  SIGN
  80  ST* 10
  81  .9                         or another number between 0 and 1 ,  I don't know which value is statistically the best!
  82  ST* 15
  83  GTO 03
  84  LBL 01
  85  XEQ 02
  86  RCL 15
  87  *
  88  RCL 09
  89  *
  90  RCL 02
  91  1
  92  P-R
  93  RCL 15
  94  ST* Z
  95  RDN
  96  R-P
  97  RDN
  98  +
  99  RCL 01
100  1
101  P-R
102  RCL 15
103  ST* Z
104  RDN
105  R-P
106  RDN
107  -
108  RCL 17
109  -
110  RTN
111  LBL 02
112  RCL 13
113  RCL 10
114  SIGN
115  RCL 12
116  X^2
117  RCL 15
118  X^2
119  -
120  SQRT
121  *
122  R-P
123  X<>Y
124  STO 02
125  RCL 11
126  RCL 10
127  X^2
128  RCL 15
129  X^2
130  -
131  SQRT
132  R-P
133  X<>Y
134  STO 01
135  XEQ "GL3"
136  RTN
137  LBL "S"
138  1
139  P-R
140  X^2
141  X<>Y
142  RCL 15
143  *
144  X^2
145  +
146  RCL 09
147  *
148  1
149  +
150  SQRT
151  FC? 10
152  RTN
153  1
154  +
155  1/X
156  RTN
157  LBL 03
158  CLX
159  SIGN
160   E-10
161  -
162  RCL 15
163  ABS
164  X>Y?
165  X<>Y
166  STO 15
167  RCL 12
168  X<Y?
169  STO 15
170  VIEW 15
171  XEQ 01
172  ENTER^
173  ENTER^
174  X<> 16
175  -
176  X#0?
177  /
178  RCL 14
179  RCL 15
180  STO 14
181  -
182  *
183  ST+ 15
184  ABS
185   E-9
186  X<Y?
187  GTO 03
188  VIEW 15
189  XEQ 01
190  STO 16
191  4
192  STO 03
193  CF 10
194  XEQ 02
195  RCL 15             lines 195 to 198 may look strange, but they are useful to correct
196  RCL 16             the inaccuracy in the µ-values, especially if  k is close to 1.
197  *
198  -
199  LBL 04
200  D-R
201  6378.137
202  *
203  CLA
204  CLD
205  END

( 295 bytes / SIZE 018 )
 
 
      STACK        INPUTS      OUTPUTS
           T       L1 ( ° ' " )             /
           Z       b1 ( ° ' " )             /
           Y       L2 ( ° ' " )             /
           X       b2 ( ° ' " )        D ( km )

Example1:

   77.0356    ENTER^
   38.55172  ENTER^
   -2.20138   ENTER^
   48.50112  XEQ "TGD3"  >>>>   the successive k-values are displayed  and eventually,    D = 6181.621797 km   ( in 2mn45s ) , the last digit should be a 4.

-We also have in register R15    k = 0.612158197

Example2:      L1 = b1 = 0  ;   L2 = 179°51'   b2 = 0

        0      ENTER^  ENTER^
   179.51  ENTER^
        0      XEQ "TGD3"  >>>>   D = 20001.85464 km ( in 2mn10s ) , the last digit should be a 3.
                                  and  R15 = k = 0.248743230

-Unlike  TGD & TGD2 , TGD3 produces exact results for nearly antipodal points too ... But it is much slower!

-If the 2 points are on the equator and  | L1 - L2 |  <  179°23'47"377  we have  k = 1  whereas  k = 0  if  | L1 - L2 | = 180°
-The maximum geodesic distance occurs with ( 0 ; 0 ) ( 180 ; 0 )  and  Dmax = 6378.137 * 3.136328278 = 20003.93143 km
 

4°) A Series Expansion
 

-"TGD4" employs the same integrals as "TGD3", but they are evaluated by a series expansion: the program runs faster and no subroutine is used.
-The first neglected term is proportinal to e6 :
 

Formulae:

   | L1 - L2 | =  [ Arctan ( k.tan µ ) ]µ1µ2 - k.(e2/2).( µ2 - µ1 ) - k.(e4/8).[ ( µ2 - µ1 ) - ( 1 - k2).( µ2 - µ1 )/2 + (1/4).( 1 - k2).( Sin 2µ2 - Sin 2µ1 ) ]

                s = ( µ2 - µ1 ) - (e2/2).[ ( µ2 - µ1 ) - ( 1 - k2).( µ2 - µ1 )/2 + ( 1 - k2).( Sin 2µ2 - Sin 2µ1 )/4 ]
                                    - (e4/8).[ ( µ2 - µ1 ) - ( 1 - k2).( µ2 - µ1 ) + ( 1 - k2).( Sin 2µ2 - Sin 2µ1 )/2 + (3/8).( 1 - k2)2.( µ2 - µ1 )
                                                  - ( 1 - k2)2.( Sin 2µ2 - Sin 2µ1 )/4 + ( 1 - k2)2.( Sin 4µ2 - Sin 4µ1 )/32 ]
 
 

Data Registers:  R00 = | L1 - L2 | ;  R01 = µ1 ;  R02 = µ2 ;  R03 & R04 = the successive k-values ; R05 = f(k)
                             R06 = r1 ( or -r1 if µ2 > 90° )  ;  R07 = sign(b1).(1-r1)1/2  ;  R08 = r2 ( calling r2 < r1 ) ;  R09 = (1-r2)1/2
                             R10 = -e2 ;  R11 = Sin 2µ2 - Sin 2µ1 ;  R12 = µ2 - µ1 ;  R13 = 1 - k2  ;  R14 = temp ;  R15 = 10 -9
Flags:  /
Subroutines:  /
 

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

( 293 bytes / SIZE 016 )
 
 
      STACK        INPUTS      OUTPUTS
           T       L1 ( ° ' " )             /
           Z       b1 ( ° ' " )             /
           Y       L2 ( ° ' " )             /
           X       b2 ( ° ' " )        D ( km )

Example1:

   77.0356    ENTER^
   38.55172  ENTER^
   -2.20138   ENTER^
   48.50112  XEQ "TGD4"  >>>>   the successive k-values are displayed  and eventually,    D = 6181.621792 km   ( in 64s ) , the last digit should be a 4.

-We also have in register R04    k = 0.612158197

Example2:      L1 = b1 = 0  ;   L2 = 179°51'   b2 = 0

        0      ENTER^  ENTER^
   179.51  ENTER^
        0      XEQ "TGD4"  >>>>   D = 20001.85474 km ( in 40s ) , error = 11cm.
                                  and  R04 = k = 0.248743896

-The maximum error is about 13 centimeters.
 

5°) Euclidean Distance
 

     a) Geocentric Coordinates
 

-If the latitude b and the height h of an observer O are known, the following routine calculates
 the geocentric latitude b' and the distance rho = OC  where C is the center of the Earth.
 

Formulae:      (rho)  sin b' = a(1-f) sin u + h sin b       where   tan u = (1-f) tan b              a = 6378.137 km
                       (rho) cos b' =  a cos u + h cos b                                                                f = 1/298.257
 
 

Data Registers: /
Flags: /
Subroutines: /
 

01  LBL "GEOC"
02  DEG
03  HR
04  RCL Y
05  STO M                synthetic register M can be replaced by R00
06  CLX
07  SIGN
08  P-R
09  ST* M
10  X<>Y
11  ST* Z
12  298.257
13  1/X
14  SIGN
15  ST- L
16  X<> L
17  CHS
18  ST* Y
19  X<> Z
20  R-P
21  CLX
22  6378.137
23  P-R
24  ST+ M
25  RDN
26  *
27  +
28  0
29  X<> M
30  R-P
31  X<>Y
32  HMS
33  END

( 65 bytes / SIZE 000 )
 
 
      STACK        INPUTS      OUTPUTS
           Y         h (km)        rho (km)
           X        b ( ° ' " )        b' ( ° ' " )

Example:   Find the geocentric coordinates of the Mount Palomar Observatory     b = 33°21'22"4     h = 1706 m = 1.706 km

   1.706      ENTER^
 33.21224  XEQ "GEOC"  >>>>    b'  = 33°10'47"12
                                         X<>Y   rho = 6373.4156 km

Notes:

-If you prefer h & rho in meters, replace line 22 by 6378137
-In several applications, (rho) sin b' & (rho) cos b' are more useful than rho or b':  Simply delete lines 30-31-32
  (rho) sin b'  will be inY-register,
  (rho) cos b' will be in X-register.
 
 

     b) 3D-Distance
 

-"EDST" computes the geocentric rectangular coordinates of 2 points and their Euclidean 3D-distance
  knowing their longitudes, latitudes and heights.
 

Data Registers:               R00 is unused                                         ( Registers R01 thru R06 are to be initialized before executing "EDST" )

                                      •  R01 = L1 ( ° ' " )   •  R04 = L2 ( ° ' " )        R07 = x1     R10 = x2
                                      •  R02 = b1 ( ° ' " )   •  R05 = b2 ( ° ' " )        R08 = y1     R11 = y2         xi  yi  zi   in km
                                      •  R03 = h1 ( km )    •  R06 = h2  ( km )        R09 = z1     R12 = z2
Flags: /
Subroutine:   "GEOC"
 

01  LBL "EDST"
02  RCL 03
03  RCL 02
04  XEQ "GEOC"
05  HR
06  X<>Y
07  P-R
08  RCL 01
09  HR
10  X<>Y
11  P-R
12  STO 07
13  RDN
14  STO 08
15  X<>Y
16  STO 09
17  RCL 06
18  RCL 05
19  XEQ "GEOC"
20  HR
21  X<>Y
22  P-R
23  RCL 04
24  HR
25  X<>Y
26  P-R
27  STO 10
28  RCL 07
29  -
30  X^2
31  X<>Y
32  STO 11
33  RCL 08
34  -
35  X^2
36  +
37  X<>Y
38  STO 12
39  RCL 09
40  -
41  X^2
42  +
43  SQRT
44  END

( 63 bytes / SIZE 013 )
 
 
      STACK        INPUT      OUTPUT
           X             /        D (km)

Example:     Find the euclidean distance between the Mount Palomar Observatory:   L1 = -116°51'50"4  ,  b1 = 33°21'22"4  ,  h1 = 1706 m = 1.706 km
                     and the "Observatoire du Pic du Midi"   L2 = 0°08'32"4  ,  b2 = 42°56'12"0  ,  h2 = 2861 m = 2.861 km

    -116.51504    STO 01        0.08324      STO 04
       33.21224    STO 02       42.56120     STO 05
         1.706        STO 03         2.861         STO 06

  XEQ "EDST"  >>>>   D = 8585.5760 km   ( execution time = 11 seconds )

-And we have the geocentric coordinates:       R07 = x1 = -2410.4237                  R10 = x2 = 4678.8290        ( all in kilometers )
                                                                      R08 = y1 = -4758.6127                  R11 = y2 =     11.6231
                                                                      R09 = z1 =   3487.9636                  R12 = z2 = 4324.3023

-In comparison, the geodesic distance is  9410.6525 km  and the "loxodromic distance" is  10284.7558 km
 

Note:

-Delete lines 20-21-22 and 05-06-07 if you have deleted lines 30-31-32 in the "GEOC" listing.
-This reduces the execution time to 8.6 seconds.
 

References:

[1]  Jean Meeus - "Astronomical Algorithms" - Willmann-Bell -  ISBN 0-943396-35-2
[2]  Jacqueline Lelong-Ferrand - "Geometrie Differentielle" - Masson - ( in French )
[3]  Vincenty's formula
 
 

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