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

-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