Geodesic2

Terrestrial Geodesic Distance(2) for the HP-41

Overview

0°) Method
1°) With Andoyer's formula of order 2
2°) With Andoyer's formula of order 3

a) Ellipsoid of Revolution
b) Triaxial Ellipsoid ( with parametric coordinates )
c) Triaxial Ellipsoid ( without parametric coordinates )

3°) Andoyer's Formula of order 2 + a few terms of order 3

a) Ellipsoid of Revolution
b) Triaxial Ellipsoid ( without parametric coordinates )
c) Triaxial Ellipsoid ( with parametric coordinates )

4°) With Vincenty's Formula
5°) Minimizing the Sum of 2 Distances
6°) Other Formulae for the Great Elliptic Distance

a) Formula of order 1 + Andoyer's formula of order 2
b) Formula of order 1 + Andoyer's formula of order 3
c) Formula of order 2 + Andoyer's formula of order 3

7°) Geodesic + Great Elliptic Distance ( Ellipsoid of Revolution )

-These programs evaluate the geodesic distance between 2 points M(L,B) & N(L',B') on a triaxial ellipsoid model of the Earth.
-The latitudes L and the longitudes B are supposed geodetic.
-The semi-axes are:

    a = 6378.172 km
b = 6378.102 km                 and the longitude of the equatorial major-axis = 14°92911 West
      c = 6356.752314 km

-So, the mean radius of the equator and the flattening correspond to WGS84.

0°) Method

-First, the distances of the 2 points M & N are computed on the great ellipses passing by M & N:

-For the triaxial model, it gives, say Δ
-For the WGS84 ellipsoid, it gives Δ'

-Then, the geodesic distance d is calculated on the WGS84 ellipsoid of revolution.

-Since the flattening of the Earth equator is very small, the geodesic distance D on the triaxial ellipsoid is very well approximated by the simple relation:

D ~ Δ + ( d - Δ' )

-In other words, the part of the formula giving the great elliptic arc distance on the ellipsoid of revolution in the WGS84 geodesic distance
is replaced by the great elliptic arc distance on the triaxial ellipsoid.

1°) With Andoyer's Formula of order 2

"TGD" determines the elements of the great ellipse (E) passing by M & N as follows:

With OM / OM = ( x , y , z ) & ON / ON = ( x' , y' ,z' ) ( the unit vectors defined by OM & ON ) we have

OP = ( Sin W ) / ( A² + B² + C² )^1/2 for a point P on the ellipsoid and with µ = ( OM , OP )

      A = [ x Sin ( W - µ ) + x' Sin µ ] / a
      B = [ y Sin ( W - µ ) + y' Sin µ ] / b           W = the angle ( OM , ON )
      C = [ z Sin ( W - µ ) + z' Sin µ ] / c

                         N P
                        /   /
                      / /
                    / / µ
               O//-----------------------------M

>>> The value of µ for the axes of the great ellipse (E) is found by equating to zero the derivative of OP with respect to µ

-After some calculations, it yields: ( Sin 2µ ) / ( Cos 2µ ) = P / Q with

P = -2 [ ( x Sin W ) ( x Cos W - x' ) / a² + ( y Sin W ) ( y Cos W - y' ) / b² + ( z Sin W ) ( z Cos W - z' ) / c² ]
Q = [ x² Sin² W - ( x Cos W - x' )² ] / a² + [ y² Sin² W - ( y Cos W - y' )² ] / b² + [ z² Sin² W - ( z Cos W - z' )² ] / c²

>>> However, it's simpler & faster to compute

    A² + B² + C² = T Sin² ( W - µ ) + S Sin² µ + 2 U Sin µ Sin ( W - µ )
   where S = x'² / a² + y'² / b² + z'² / c² ,   T = x² / a² + y² / b² + z² / c² ,   U = xx' / a² + yy' / b² + zz' / c²
    and   P = 2 U - 2 T Cos W , Q = T Sin W + ( 2 U Cos W - S - T Cos² W ) / Sin W

>>> But we can save a few bytes if we measure the angle µ from the bisector of the angle ( OM , ON )

The formulae become:       A² + B² + C² = T Sin² ( W/2 - µ ) + S Sin² ( W/2 + µ ) + 2 U Sin ( W/2 - µ ) Sin ( W/2 + µ )

and Tan 2µ = [ Sin W ( T - S ) ] / [ Cos W ( T + S ) - 2.U ]

>>> Then, the geodesic distances - on the great ellipses corresponding to WGS84 & the triaxial ellipsoid - between M & N ( say T' & T ) are computed by Andoyer's formula on order 2.

-We can transform the usual formulas employing parametric latitudes into formulas employing geocentric latitudes ( cf page 91 of reference [4] ). It yields:

D/a = W + (f/2) (-W+(Cos2F)(SinW)) +(f² /4) [W/4+(SinW)(CosW)/4-(Cos² (2F))(SinW)(CosW)/2]

a = major axis of the great ellipse , f = eccentricity
2F = b1+b2 where b = geocentric "latitudes"

( formulas of order 2 give a good accuracy to compute great elliptic distance for the Earth )

-Finally, the geodesic distance d on the WGS84 ellipsoid is calculated by 2nd order Andoyer's formula

and the geodesic distance D on the triaxial ellipsoid is approximated by D ~ Δ + ( d - Δ' )

Andoyer's Formula of order 2:

-With   H = ( L' - L )/2 ;   F = ( b + b' )/2 ; G = ( b' - b )/2

              S = sin² G cos² H + cos² F sin² H ;   µ = Arcsin S^1/2 ; R = ( S.(1-S) )^1/2

we define A = ( sin² G cos² F )/S + ( sin² F cos² G )/( 1-S )
      and      B = ( sin² G cos² F )/S - ( sin² F cos² G )/( 1-S )

-The geodesic distance on the biaxial ellipsoid is then:   d = 2a.µ ( 1 + eps )

where   eps =    (f/2) ( -A - 3B R/µ )
                       + (f²/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 }

   and   f = flattening = 1/298.25722

And we also have:

   d - Δ' = a (f²/4) µ [ (A²-A) ( 1 - µ Cot µ ) - B.(A-1) ( µ / Sin µ - (Sin µ) / µ ) ]      where µ , A , B are computed with parametric latitudes

( here, µ = 2 Arcsin S^1/2 )

-But with a formula of order 2, d - Δ' is computed by the same formula if µ , A , B are computed with geodetic latitudes !

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

01 LBL "TGD"
02 DEG
03 HR
04 STO 14
05 RDN
06 HR
07 14.929
08 STO 11
09 +
10 STO 13
11 RDN
12 HR
13 STO 12
14 X<>Y
15 HR
16 ST+ 11
17 6378.172
18 STO 01
19 .07
20 -
21 STO 02
22 21.349686
23 -
24 STO 03
25 RCL 11
26 RCL 12
27 XEQ 01
28 STO 05
29 RCL 09
30 STO 06
31 RCL 10
32 STO 07
33 RCL 13
34 ST- 11
35 RCL 14
36 XEQ 01
37 XEQ 00
38 STO 04
39 RCL 01
40 ST/ 05

41 ST/ 08
42 RCL 02
43 ST+ 01
44 ST/ 06
45 ST/ 09
46 RCL 03
47 ST/ 07
48 ST/ 10
49 RCL 08
50 X^2
51 RCL 09
52 X^2
53 RCL 10
54 X^2
55 +
56 +
57 STO 13
58 RCL 05
59 X^2
60 RCL 06
61 X^2
62 RCL 07
63 X^2
64 +
65 +
66 STO 00
67 -
68 RCL 04
69 ACOS
70 STO 02
71 SIN
72 STO 15
73 *
74 RCL 08
75 XEQ 00
76 ST+ X
77 STO 09
78 RCL 00
79 RCL 13
80 +

81 RCL 04
82 *
83 -
84 X#0?
85 /
86 ATAN
87 STO 08
88 2
89 ST/ 01
90 ST/ 11
91 /
92 STO 05
93 XEQ 02
94 STO 06
95 RCL 05
96 90
97 +
98 XEQ 02
99 STO 07
100 RCL 02
101 D-R
102 RCL 08
103 COS
104 STO 05
105 X^2
106 ST+ X
107 RCL 04
108 RCL 15
109 ST* 05
110 *
111 ST* Y
112 -
113 -
114 4
115 /
116 RCL 06
117 ST/ 15
118 RCL 07
119 ST- Y
120 ST+ X

121 /
122 STO 07
123 *
124 RCL 05
125 -
126 +
127 RCL 07
128 *
129 +
130 ST* 15
131 RCL 14
132 ENTER
133 GTO 03
134 LBL 00
135 RCL 05
136 *
137 RCL 09
138 RCL 06
139 *
140 +
141 RCL 10
142 RCL 07
143 *
144 +
145 RTN
146 LBL 01
147 1
148 P-R
149 X<>Y
150 RCL 03
151 X^2
152 *
153 STO 10
154 RDN
155 P-R
156 RCL 01
157 X^2
158 *
159 STO Z
160 X^2

161 X<>Y
162 RCL 02
163 X^2
164 *
165 STO 09
166 X^2
167 RCL 10
168 X^2
169 +
170 +
171 SQRT
172 ST/ 09
173 ST/ 10
174 /
175 STO 08
176 RTN
177 LBL 02
178 RCL 02
179 2
180 /
181 +
182 ENTER
183 SIN
184 ENTER
185 X^2
186 RCL 13
187 *
188 RCL 02
189 R^
190 -
191 SIN
192 ST* Z
193 X^2
194 RCL 00
195 *
196 +
197 X<>Y
198 RCL 09
199 *
200 +

201 SQRT
202 RTN
203 LBL 03
204 RCL 12
205 +
206 2
207 /
208 ST- Y
209 1
210 P-R
211 X^2
212 STO T
213 RDN
214 X^2
215 STO 08
216 SIGN
217 P-R
218 X^2
219 ST* 08
220 CLX
221 +
222 X^2
223 ST* T
224 -
225 RCL 11
226 COS
227 X^2
228 *
229 -
230 /
231 STO 02
232 RCL 08
233 1
234 LASTX
235 STO 08
236 -
237 /
238 ST- 02
239 +
240 X<> 08

241 SQRT
242 ASIN
243 ST+ X
244 D-R
245 STO 00
246 LASTX
247 1
248 P-R
249 RDN
250 /
251 STO 07
252 R^
253 *
254 RCL 08
255 ST* Y
256 -
257 RCL 07
258 ENTER
259 1/X
260 -
261 RCL 02
262 *
263 +
264 1
265 RCL 08
266 -
267 *
268 RCL 03
269 RCL 01
270 ST- Y
271 ST+ X
272 /
273 X^2
274 *
275 RCL 00
276 *
277 RCL 01
278 *
279 RCL 15
280 +
281 END

( 358 bytes / SIZE 016 )

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

Where D is the geodesic distance on the triaxial ellipsoid

Example: Calculate the geodesic distance between the Sydney Observatory ( Long = 151°12'17"8 E , Lat = 33°51'41"1 S )
and the Palomar Observatory ( Long = 116°51'50"4 W , Lat = 33°21'22"4 N )

     151.12178 ENTER^
    -33.51411   ENTER^
   -116.51504 ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65755 km                                                        ---Execution time = 23s---
                                            LASTX d = 12138.65875 km = great elliptic distance = R10

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

Notes:

-Longitudes are geodetic and positive East.

-When line 99 is executed ( STO 07 ), the semi axes of the great ellipse = R15/R06 & R15/R07.
-With the example above, 6378.151662 km & 6368.318461 km

-If the points are not nearly antipodal, the precision is of the order of one meter.

-With (L,b,L',b') = (70°,0°,270°,10°) "TGD" gives a distance of 17554.6159 instead of 17554.6140 km

-Here are a few results about the error of the results:

If | L' - L | < 120° rms error ~ 6 cm   and max error ~ 26 cm
If | L' - L | < 140° rms error ~ 7 cm   and max error ~ 38 cm
If | L' - L | < 150° rms error ~ 8 cm   and max error ~ 70 cm
If | L' - L | < 160° rms error ~ 11 cm and max error ~ 168 cm

-A better accuracy may be obtained by the programs listed in "Geodesics on a triaxial ellipsoid for the HP-41" but the program above is faster !

2°) With Andoyer's Formula of order 3

a) Ellipsoid of Revolution

-Sodano's formula is a 3rd order formula in the flattening to calculate the geodesic distance on an ellipsoid of revolution ( WGS84 ) - cf references [2] & [3]

-I've rewritten the formula as follows, after multiplying both sides by (1-f)

Dist = a { µ + f [ p Sin µ - (m/2) ( µ + Sin µ Cos µ )

+ f² [ p ((m-1)/2) µ² Csc µ + m ((1-m)/2) µ² Cot µ + ((m²-8.p²)/16) Sin µ Cos µ + (m²/16) µ - (m²/8) Sin µ Cos³ µ + (mp/2) Sin µ Cos² µ ]

          + f³ [ µ² Csc µ (p/4) (-2+5m-3m²) + µ³ Csc µ Cot µ (p/2) (1-4m+3m²) - µ² Cot µ (m/4) (-2+5m-3m²) + µ³ (m/6) (-1+2m-m²)
                + µ³ Cot² µ (m/2) (-1+3m-2m²) + Sinµ Cos µ (-8 p² + m² - 8 p²m - m³/2 )/16 + µ (8 p² + m² - 8 p²m - m³/2 )/16
                + µ³ Csc² µ (p²/2) (1-m) + Sin µ Cos³ µ (m²/16) (-2+m) + µ Cos² µ (m²/4) (1-m) + Sin µ Cos² µ (pm/2) (1+m) + p²m Sin³ µ Cos µ
                + µ Cos µ (3pm/4) (-1+m) + (pm²/2) Sin⁵ µ + Sin µ (p/2) (p²-m²) + (m³/12) Sin³ µ Cos³ µ - (2p³/3) Sin³ µ ] }

where µ = 2 Arc Sin [ sin² ((B'-B)/2) cos² ((L'-L)/2) + cos² ((B+B')/2) sin² ((L'-L)/2) ]^1/2 L , L' = longitudes & B , B' = parametric latitudes

Tan B = (1-f) Tan b where b = geodetic latitude

p = Sin B Sin B' m = 1 - [ Cos B Cos B' Sin (L'-L) ]² / Sin² µ , a = 6378.137 km , b = 6356.752314 km ( WGS84 )

-But Sodano's formula may be re-written as an Andoyer's formula.

( In reference [1] - page 87 - Paul D. Thomas gives the key to perform this transformation )

With L , L' = longitudes & p , p' = parametric latitudes , f = earth flattening

Tan p = (1-f) Tan b where b = geodetic latitude

H = ( L' - L )/2 ; F = ( p + p' )/2 ; G = ( p' - p )/2

S = sin² G cos² H + cos² F sin² H , µ = 2 Arc Sin S^1/2

we define A = ( sin² G cos² F )/S + ( sin² F cos² G )/( 1-S )
and B = ( sin² G cos² F )/S - ( sin² F cos² G )/( 1-S )

-After "some" calculations, it yields:

Dist = a { µ - (f/2) ( A.µ + B Sin µ )

+ (f²/4) [ A µ² Cot µ + B µ² / Sin µ + A² ( µ + Sin µ Cos µ - 4 µ² Cot µ ) / 4 - A.B µ² / Sin µ - ( B²/2 ) Sin µ Cos µ ]

                + (f³/8) [ A ( 2 µ² Cot µ - 2 µ³ Cot² µ - 8 µ³ / 6 ) + A² ( -5 µ² Cot µ + 5 µ³ Cot² µ + 8 µ³ / 3 + (1/2) Sin µ Cos µ + µ /2 )
                          + A³ ( 3 µ² Cot µ - 3 µ³ Cot² µ - 8 µ³ / 6 + (1/4) Sin µ Cos µ - µ / 4 )
                          + B ( 2 µ² / Sin µ - 2 µ³ Cos µ / Sin² µ ) + A.B ( -5 µ² / Sin µ + 6 µ³ Cos µ / Sin² µ + µ Cos µ )
                        + A² B ( 3 µ² / Sin µ - 4 µ³ Cos µ / Sin² µ - µ Cos µ + (1/2) Sin µ Cos² µ )
                        + B² ( - Sin µ Cos µ + µ³ / Sin² µ + µ ) + A B² ( (1/2) Sin µ Cos µ - µ³ / Sin² µ - µ )
                        + B³ ( (2/3) Sin³ µ - (1/2) Sin µ ) ] }

Data Registers: R00 = µ , R01 = A , R02 = B , R04 thru R09: temp
Flags: /
Subroutines: /

01 LBL "AND"
02 DEG
03 HR
04 STO 01
05 X<> T
06 HMS-
07 HR
08 2
09 /
10 SIN
11 X^2
12 STO 00
13 X<>Y
14 HR
15 1
16 STO 08
17 P-R
18 LASTX
19 298.25722
20 STO 04
21 ST+ 04
22 1/X
23 -
24 STO 05
25 /
26 R-P
27 X<> 01
28 RCL 08
29 P-R
30 RCL 05

31 /
32 R-P
33 RDN
34 +
35 2
36 /
37 ST- Y
38 RCL 08
39 P-R
40 X^2
41 STO 02
42 X<>Y
43 X^2
44 STO 01
45 RDN
46 X<>Y
47 RCL 08
48 P-R
49 X^2
50 ST* 01
51 RDN
52 X^2
53 ST* 02
54 STO Z
55 -
56 RCL 00
57 *
58 +
59 ST/ 02
60 RCL 02

61 RCL 01
62 RCL 08
63 R^
64 STO 01
65 -
66 /
67 ST- 02
68 +
69 X<> 01
70 SQRT
71 ASIN
72 ST+ X
73 D-R
74 STO 00
75 LASTX
76 RCL 08
77 P-R
78 STO 03
79 STO 06
80 RDN
81 STO 05
82 /
83 STO 07
84 ST/ Z
85 X^2
86 ST* 03
87 RCL 08
88 +
89 ENTER
90 R^

91 ST+ 08
92 ST- Z
93 2
94 /
95 STO 09
96 -
97 RCL 01
98 *
99 -
100 RCL 05
101 X^2
102 1.5
103 ST/ Y
104 FRC
105 -
106 RCL 07
107 /
108 RCL 02
109 *
110 +
111 RCL 02
112 *
113 RCL 06
114 X^2
115 RCL 07
116 ST- 03
117 ST+ X
118 /
119 RCL 03
120 4

121 ST/ 08
122 *
123 RCL 06
124 +
125 RCL 07
126 +
127 STO Z
128 -
129 RCL 01
130 *
131 +
132 RCL 03
133 ST+ X
134 STO 05
135 +
136 RCL 01
137 *
138 +
139 RCL 05
140 -
141 RCL 02
142 *
143 RCL 03
144 RCL 06
145 *
146 STO 03
147 RCL 00
148 ST+ X
149 X^2
150 3

151 ST* Z
152 /
153 STO Z
154 +
155 RCL 08
156 +
157 STO 05
158 RCL 01
159 *
160 RCL 05
161 ST+ X
162 -
163 RCL 03
164 ST+ Z
165 ST+ Z
166 +
167 RCL 01
168 *
169 +
170 RCL 01
171 *
172 -
173 RCL 04
174 /
175 RCL 08
176 RCL 06
177 RCL 07
178 *
179 STO Z
180 -
181 RCL 01

182 *
183 +
184 RCL 02
185 ST* 09
186 RCL 07
187 ST- 09
188 *
189 -
190 RCL 01
191 *
192 +
193 RCL 02
194 RCL 09
195 *
196 -
197 RCL 04
198 /
199 RCL 02
200 RCL 07
201 /
202 -
203 RCL 01
204 -
205 RCL 04
206 /
207 RCL 00
208 ST* Y
209 +
210 6378.137
211 *
212 END

( 268 bytes / SIZE 010 )

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

Example: Calculate the geodesic distance between the Sydney Observatory ( Long = 151°12'17"8 E , Lat = 33°51'41"1 S )
and the Palomar Observatory ( Long = 116°51'50"4 W , Lat = 33°21'22"4 N )

     151.12178 ENTER^
    -33.51411   ENTER^
   -116.51504 ENTER^
     33.21224   XEQ "AND"   >>>>   D = 12138.68432 km                                                      ---Execution time = 13s---

Notes:

-It is certainly possible to express the formula without calculating the parametric latitudes, but I've not found how.
-Without parametric latitudes, Andoyer's formula of order 2 gives a shorter & faster program,
but it's not at all sure with Andoyer's formula of order 3...

-We can also use an M-code routine to calculate the parametric latitudes:

( "PRM" doesn't check for alpha data )

0BD "M"
012 "R"
010 "P"
2A0 SETDEC
0F8 C=X
22D C=
060 1/C
2BE C=-C
02E C
0FA ->
0AE AB
001 C=
060 AB+1
128 L=C
0B8 C=Y
070 C
209 =
048 Tan(C)
138 C=L
13D C=
060 AB*C
10E A=C
0B8 C=Y
11E A=C MS
0AE A<>C
070 C
2A9 =
040 Atan(C)
0E8 X=C
078 C=Z
070 C
209 =
048 Tan(C)
138 C=L
13D C=
060 AB*C
10E A=C
078 C=Z
11E A=C MS
0AE A<>C
070 C
2A9 =
040 Atan(C)
068 Z=C
0A8 Y=C
3E0 RTN

STACK	INPUTS	OUTPUTS
Z	b	b_param
Y	b'	b_param
X	1/f	b'_param

b) Triaxial Ellipsoid ( With Parametric Coordinates )

-The following "TGD" employs Andoyer's formula of order 2 for the great elliptic distance and Andoyer's formula of order 3 to estimate d - Δ'

   (d - Δ') / a = (f²/4) µ [ (A²-A) ( 1 - µ Cot µ ) - B.(A-1) ( µ / Sin µ - (Sin µ) / µ ) ]       where µ , A , B are computed with parametric latitudes

+ (f³/8) [ A ( 2 µ² Cot µ - 2 µ³ Cot² µ - 8 µ³ / 6 ) + A² ( -5 µ² Cot µ + 5 µ³ Cot² µ + 8 µ³ / 3 + Sin µ Cos µ - µ )
                          + A³ ( 3 µ² Cot µ - 3 µ³ Cot² µ - 8 µ³ / 6 - Sin µ Cos µ + µ )
                          + B ( 2 µ² / Sin µ - 2 µ³ Cos µ / Sin² µ ) + A.B ( -5 µ² / Sin µ + 6 µ³ Cos µ / Sin² µ + µ Cos µ - 2 Sin µ )
                        + A² B ( 3 µ² / Sin µ - 4 µ³ Cos µ / Sin² µ - µ Cos µ + 2 Sin µ )
                        + B² ( - 2 Sin µ Cos µ + µ³ / Sin² µ + µ ) + A B² ( 2 Sin µ Cos µ - µ³ / Sin² µ - µ )

Data Registers: R00 thru R15: temp
Flags: /
Subroutines: /

01 LBL "TGD"
02 DEG
03 HR
04 STO 14
05 RDN
06 HR
07 14.92911
08 STO 11
09 +
10 STO 13
11 RDN
12 HR
13 STO 12
14 X<>Y
15 HR
16 ST+ 11
17 6378.172
18 STO 01
19 .07
20 -
21 STO 02
22 21.349686
23 -
24 STO 03
25 RCL 11
26 RCL 12
27 XEQ 01
28 STO 05
29 RCL 09
30 STO 06
31 RCL 10
32 STO 07
33 RCL 13
34 ST- 11
35 RCL 14
36 XEQ 01
37 XEQ 00
38 STO 04
39 RCL 01
40 ST/ 05
41 ST/ 08
42 RCL 02
43 ST+ 01
44 ST/ 06
45 ST/ 09
46 RCL 03
47 ST/ 07
48 ST/ 10
49 RCL 08
50 X^2
51 RCL 09
52 X^2
53 RCL 10

54 X^2
55 +
56 +
57 STO 13
58 RCL 05
59 X^2
60 RCL 06
61 X^2
62 RCL 07
63 X^2
64 +
65 +
66 STO 00
67 -
68 RCL 04
69 ACOS
70 STO 02
71 SIN
72 STO 15
73 *
74 RCL 08
75 XEQ 00
76 ST+ X
77 STO 09
78 RCL 00
79 RCL 13
80 +
81 RCL 04
82 *
83 -
84 X#0?
85 /
86 ATAN
87 STO 08
88 2
89 ST/ 01
90 ST/ 11
91 /
92 STO 05
93 XEQ 02
94 STO 06
95 RCL 05
96 90
97 +
98 XEQ 02
99 STO 07
100 RCL 02
101 D-R
102 RCL 08
103 COS
104 STO 05
105 X^2
106 ST+ X

107 RCL 04
108 RCL 15
109 ST* 05
110 *
111 ST* Y
112 -
113 -
114 4
115 /
116 RCL 06
117 ST/ 15
118 RCL 07
119 ST- Y
120 ST+ X
121 /
122 STO 07
123 *
124 RCL 05
125 -
126 +
127 RCL 07
128 *
129 +
130 ST* 15
131 GTO 03
132 LBL 00
133 RCL 05
134 *
135 RCL 09
136 RCL 06
137 *
138 +
139 RCL 10
140 RCL 07
141 *
142 +
143 RTN
144 LBL 01
145 1
146 P-R
147 X<>Y
148 RCL 03
149 X^2
150 *
151 STO 10
152 RDN
153 P-R
154 RCL 01
155 X^2
156 *
157 STO Z
158 X^2
159 X<>Y

160 RCL 02
161 X^2
162 *
163 STO 09
164 X^2
165 RCL 10
166 X^2
167 +
168 +
169 SQRT
170 ST/ 09
171 ST/ 10
172 /
173 STO 08
174 RTN
175 LBL 02
176 RCL 02
177 2
178 /
179 +
180 ENTER
181 SIN
182 ENTER
183 X^2
184 RCL 13
185 *
186 RCL 02
187 R^
188 -
189 SIN
190 ST* Z
191 X^2
192 RCL 00
193 *
194 +
195 X<>Y
196 RCL 09
197 *
198 +
199 SQRT
200 RTN
201 LBL 03
202 RCL 14
203 1
204 STO 08
205 STO 09
206 P-R
207 RCL 01
208 RCL 03
209 -
210 RCL 01
211 /
212 STO 04

213 SIGN
214 LASTX
215 -
216 STO 05
217 /
218 R-P
219 X<> 12
220 RCL 09
221 P-R
222 RCL 05
223 /
224 R-P
225 RDN
226 +
227 2
228 ST/ 04
229 /
230 ST- Y
231 RCL 09
232 P-R
233 X^2
234 STO T
235 RDN
236 X^2
237 STO 14
238 SIGN
239 P-R
240 X^2
241 ST* 14
242 CLX
243 +
244 X^2
245 ST* T
246 -
247 RCL 11
248 COS
249 X^2
250 *
251 -
252 /
253 STO 02
254 RCL 14
255 RCL 09
256 LASTX
257 STO 14
258 -
259 /
260 ST- 02
261 +
262 ST- 09
263 X<> 14
264 SQRT
265 ASIN

266 ST+ X
267 D-R
268 STO 00
269 LASTX
270 RCL 08
271 P-R
272 STO 03
273 STO 06
274 STO 13
275 RDN
276 STO 05
277 /
278 STO 07
279 ST/ Z
280 X^2
281 ST* 03
282 RCL 08
283 +
284 R^
285 ST- 08
286 ST+ X
287 -
288 RCL 02
289 *
290 2
291 RCL 07
292 ST- 03
293 ST/ Y
294 -
295 RCL 06
296 -
297 RCL 03
298 ST+ 03
299 ST* 13
300 4
301 *
302 -
303 RCL 09
304 ST* Z
305 *
306 RCL 03
307 -
308 RCL 14
309 *
310 -
311 RCL 03
312 -
313 RCL 02
314 *
315 RCL 13
316 RCL 00
317 ST* 01
318 ST+ X

319 X^2
320 3
321 ST* Z
322 /
323 STO 11
324 +
325 STO 12
326 RCL 08
327 -
328 RCL 14
329 *
330 RCL 12
331 ST+ X
332 -
333 RCL 13
334 +
335 RCL 08
336 +
337 RCL 14
338 *
339 RCL 13
340 ST+ X
341 +
342 RCL 11
343 +
344 RCL 14
345 *
346 -
347 RCL 04
348 *
349 RCL 06
350 RCL 07
351 *
352 RCL 14
353 ST* Y
354 -
355 RCL 07
356 ENTER
357 1/X
358 -
359 RCL 02
360 *
361 +
362 RCL 09
363 *
364 +
365 RCL 04
366 X^2
367 *
368 RCL 01
369 *
370 RCL 15
371 +
372 END

( 466 bytes / SIZE 016 )

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

Where D is the geodesic distance on the triaxial ellipsoid

Example: Calculate the geodesic distance between the Sydney Observatory ( Long = 151°12'17"8 E , Lat = 33°51'41"1 S )
and the Palomar Observatory ( Long = 116°51'50"4 W , Lat = 33°21'22"4 N )

     151.12178 ENTER^
    -33.51411   ENTER^
   -116.51504 ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65755 km                                                      ---Execution time = 28s---

LASTX = R10 = 12138.65875 km = great elliptic arc distance

Notes:

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

-When line 99 is executed ( STO 07 ), the semi axes of the great ellipse = R15/R06 & R15/R07.
-With the example above, 6378.151662 km & 6368.318461 km

-The root-mean-square error is about 2 cm if 80° < | L' - L | < 120° and the maximum error is about 13 cm
-The root-mean-square error is about 4 cm if 110° < | L' - L | < 150° and the maximum error is about 34 cm
-When the difference in the longitudes is 165°, the error can reach about 75 cm.

     c) Triaxial Ellipsoid ( Without Parametric Coordinates )

-The formula employed in the previous paragraph may be expressed without transforming geodetic coordinates into parametric coordinates.
-After some calculations, it yields:

   (d - Δ') / a = (f²/4) µ [ (A²-A) ( 1 - µ Cot µ ) - B.(A-1) ( µ / Sin µ - (Sin µ) / µ ) ]       where µ , A , B are computed with geodetic latitudes

+ (f³/8) [ A ( 4 µ - 2 µ² Cot µ - 2 µ³ Cot² µ - 8 µ³ / 6 )
       + A² ( 9 µ² Cot µ + 5 µ³ Cot² µ + 8 µ³ / 3 - Sin µ Cos µ - 13 µ )
                        + A³ ( -7 µ² Cot µ - 3 µ³ Cot² µ - 8 µ³ / 6 + Sin µ Cos µ + 9 µ )
                         + B ( -2 µ² / Sin µ - 2 µ³ Cos µ / Sin² µ + 4 Sin µ )
+ A.B ( 9 µ² / Sin µ + 6 µ³ Cos µ / Sin² µ - 3 µ Cos µ - 12 Sin µ )
                        + A² B ( -7 µ² / Sin µ - 4 µ³ Cos µ / Sin² µ + 3 µ Cos µ + 8 Sin µ )
                        + B² ( 2 Sin µ Cos µ + µ³ / Sin² µ - 3 µ )
+ A B² ( -2 Sin µ Cos µ - µ³ / Sin² µ + 3 µ )

Data Registers: R00 thru R15: temp
Flags: /
Subroutines: /

53 RCL 10
54 X^2
55 +
56 +
57 STO 13
58 RCL 05
59 X^2
60 RCL 06
61 X^2
62 RCL 07
63 X^2
64 +
65 +
66 STO 00
67 -
68 RCL 04
69 ACOS
70 STO 02
71 SIN
72 STO 15
73 *
74 RCL 08
75 XEQ 00
76 ST+ X
77 STO 09
78 RCL 00
79 RCL 13
80 +
81 RCL 04
82 *
83 -
84 X#0?
85 /
86 ATAN
87 STO 08
88 2
89 ST/ 01
90 ST/ 11
91 /
92 STO 05
93 XEQ 02
94 STO 06
95 RCL 05
96 90
97 +
98 XEQ 02
99 STO 07
100 RCL 02
101 D-R
102 RCL 08
103 COS
104 STO 05

105 X^2
106 ST+ X
107 RCL 04
108 RCL 15
109 ST* 05
110 *
111 ST* Y
112 -
113 -
114 4
115 /
116 RCL 06
117 ST/ 15
118 RCL 07
119 ST- Y
120 ST+ X
121 /
122 STO 07
123 *
124 RCL 05
125 -
126 +
127 RCL 07
128 *
129 +
130 ST* 15
131 RCL 14
132 ENTER
133 GTO 03
134 LBL 00
135 RCL 05
136 *
137 RCL 09
138 RCL 06
139 *
140 +
141 RCL 10
142 RCL 07
143 *
144 +
145 RTN
146 LBL 01
147 1
148 P-R
149 X<>Y
150 RCL 03
151 X^2
152 *
153 STO 10
154 RDN
155 P-R
156 RCL 01

157 X^2
158 *
159 STO Z
160 X^2
161 X<>Y
162 RCL 02
163 X^2
164 *
165 STO 09
166 X^2
167 RCL 10
168 X^2
169 +
170 +
171 SQRT
172 ST/ 09
173 ST/ 10
174 /
175 STO 08
176 RTN
177 LBL 02
178 RCL 02
179 2
180 /
181 +
182 ENTER
183 SIN
184 ENTER
185 X^2
186 RCL 13
187 *
188 RCL 02
189 R^
190 -
191 SIN
192 ST* Z
193 X^2
194 RCL 00
195 *
196 +
197 X<>Y
198 RCL 09
199 *
200 +
201 SQRT
202 RTN
203 LBL 03
204 RCL 12
205 +
206 2
207 /
208 ST- Y

209 1
210 STO 07
211 P-R
212 X^2
213 STO T
214 RDN
215 X^2
216 STO 08
217 SIGN
218 P-R
219 X^2
220 ST* 08
221 CLX
222 +
223 X^2
224 ST* T
225 -
226 RCL 11
227 COS
228 X^2
229 *
230 -
231 /
232 STO 02
233 RCL 08
234 RCL 07
235 LASTX
236 STO 08
237 -
238 /
239 ST- 02
240 +
241 ST- 07
242 X<> 08
243 SQRT
244 ASIN
245 ST+ X
246 D-R
247 STO 00
248 LASTX
249 1
250 P-R
251 STO 09
252 STO 12
253 STO 13
254 RDN
255 /
256 STO 06
257 ST/ 12
258 X^2
259 ST* 09
260 3

261 -
262 RCL 12
263 ST+ X
264 +
265 RCL 02
266 *
267 3
268 RCL 06
269 ST* Y
270 RCL 09
271 +
272 ST+ X
273 STO 10
274 STO 11
275 ST+ X
276 +
277 ST- 09
278 RCL 13
279 3
280 *
281 -
282 8
283 RCL 06
284 /
285 -
286 RCL 07
287 ST* Z
288 *
289 RCL 10
290 +
291 4
292 STO 14
293 RCL 06
294 /
295 ST- 11
296 -
297 RCL 08
298 *
299 +
300 RCL 11
301 -
302 RCL 02
303 *
304 RCL 09
305 RCL 13
306 ST* 10
307 *
308 RCL 12
309 +
310 9
311 +
312 RCL 00
313 ST+ X

314 X^2
315 3
316 /
317 ST- 14
318 -
319 RCL 07
320 *
321 RCL 10
322 -
323 RCL 14
324 ST- 10
325 +
326 RCL 08
327 *
328 RCL 10
329 +
330 RCL 08
331 *
332 -
333 RCL 01
334 RCL 03
335 -
336 RCL 01
337 ST+ X
338 /
339 STO 04
340 *
341 RCL 06
342 RCL 13
343 *
344 RCL 08
345 ST* Y
346 -
347 RCL 06
348 ENTER
349 1/X
350 -
351 RCL 02
352 *
353 +
354 RCL 07
355 *
356 +
357 RCL 04
358 X^2
359 *
360 RCL 00
361 *
362 RCL 01
363 *
364 RCL 15
365 +
366 END

( 458 bytes / SIZE 016 )

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

Where D is the geodesic distance on the triaxial ellipsoid

Example: Calculate the geodesic distance between the Sydney Observatory ( Long = 151°12'17"8 E , Lat = 33°51'41"1 S )
and the Palomar Observatory ( Long = 116°51'50"4 W , Lat = 33°21'22"4 N )

     151.12178 ENTER^
    -33.51411   ENTER^
   -116.51504 ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65755 km                                                      ---Execution time = 25s---

LASTX = R15 = 12138.65875 km = great elliptic arc distance

Note:

-When line 99 is executed ( STO 07 ), the semi axes of the great ellipse = R15/R06 & R15/R07.
-With the example above, 6378.151662 km & 6368.318461 km

3°) With Andoyer's Formula of order 2 + A Few Terms of Order 3

a) Ellipsoid of Revolution

-This version employs Andoyer's formula of order 2 ( where A , B , ... are computed with geodetic latitudes ) + a few terms of order 3

Formulae:

-With H = ( L' - L )/2 ; F = ( b + b' )/2 ; G = ( b' - b )/2

S = sin² G cos² H + cos² F sin² H ; µ = 2 Arcsin S^1/2

we define A = ( sin² G cos² F )/S + ( sin² F cos² G )/( 1-S )
and B = ( sin² G cos² F )/S - ( sin² F cos² G )/( 1-S )

Dist ~ a { µ - (f/2) ( A.µ + 3.B Sin µ )

+ (f²/4) [ A ( 4.µ + µ² Cot µ ) + B ( 6.Sin µ + µ² / Sin µ ) - A² ( 15.µ + 15 Sin µ Cos µ + 4 µ² Cot µ ) / 4 - A.B ( 6.Sin µ + µ² / Sin µ ) + ( B²/2 ) 15 Sin µ Cos µ ]

                + (f³/8) [ A ( 2 µ² Cot µ - 2 µ³ Cot² µ - 8 µ³ / 6 ) + A² ( -5 µ² Cot µ + 5 µ³ Cot² µ + 8 µ³ / 3 )
                          + A³ ( 3 µ² Cot µ - 3 µ³ Cot² µ - 8 µ³ / 6 )
                          + B ( 2 µ² / Sin µ - 2 µ³ Cos µ / Sin² µ ) + A.B ( -5 µ² / Sin µ + 6 µ³ Cos µ / Sin² µ )
                        + A² B ( 3 µ² / Sin µ - 4 µ³ Cos µ / Sin² µ )
                        + B² ( µ³ / Sin² µ ) + A B² ( - µ³ / Sin² µ )

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

01 LBL "TGD"
02 DEG
03 HR
04 X<> T
05 HMS-
06 HR
07 2
08 /
09 COS
10 X^2
11 RDN
12 HR
13 +
14 2
15 /
16 ST- Y
17 1
18 STO 07
19 P-R
20 X^2
21 STO 02
22 RDN
23 X^2
24 STO 01
25 SIGN

26 P-R
27 X^2
28 ST* 01
29 RDN
30 X^2
31 ST* 02
32 RCL Z
33 -
34 *
35 +
36 ST/ 02
37 RCL 02
38 RCL 01
39 RCL 07
40 R^
41 STO 01
42 -
43 /
44 ST- 02
45 +
46 ST- 07
47 X<> 01
48 SQRT
49 ASIN
50 ST+ X

51 D-R
52 STO 00
53 LASTX
54 1
55 P-R
56 STO 03
57 STO 05
58 RDN
59 /
60 STO 06
61 X^2
62 ST* 03
63 RCL 02
64 *
65 RCL 07
66 *
67 RCL 06
68 ST- 03
69 RCL 03
70 4
71 *
72 +
73 RCL 01
74 *
75 RCL 03

76 6
77 *
78 -
79 RCL 06
80 -
81 RCL 01
82 *
83 -
84 RCL 03
85 ST+ X
86 -
87 RCL 02
88 *
89 RCL 03
90 RCL 05
91 *
92 STO 03
93 3
94 *
95 RCL 00
96 X^2
97 .75
98 /
99 STO 04
100 +

101 RCL 01
102 *
103 RCL 04
104 ST+ X
105 -
106 RCL 03
107 ST+ 04
108 ST+ 04
109 5
110 *
111 -
112 RCL 01
113 *
114 RCL 04
115 +
116 RCL 01
117 *
118 -
119 596.51444
120 STO 04
121 /
122 6
123 RCL 06
124 ST/ Y
125 +

126 RCL 07
127 *
128 RCL 05
129 RCL 06
130 ST* 05
131 /
132 7.5
133 STO 03
134 *
135 ST+ 03
136 RCL 02
137 *
138 +
139 RCL 02
140 *
141 +
142 RCL 03
143 2
144 /
145 RCL 05
146 +
147 RCL 01
148 *
149 RCL 05

150 -
151 4
152 -
153 RCL 01
154 *
155 -
156 RCL 04
157 /
158 RCL 01
159 -
160 RCL 02
161 3
162 *
163 RCL 06
164 /
165 -
166 RCL 04
167 /
168 RCL 00
169 ST* Y
170 +
171 6378.137
172 *
173 END

( 220 bytes / SIZE 008 )

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

Where D is the geodesic distance on the WGS84 ellipsoid

Example: Calculate the geodesic distance between the Sydney Observatory ( Long = 151°12'17"8 E , Lat = 33°51'41"1 S )
and the Palomar Observatory ( Long = 116°51'50"4 W , Lat = 33°21'22"4 N )

     151.12178 ENTER^
    -33.51411   ENTER^
   -116.51504 ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.68425 km ---Execution time = 10s---

Notes:

-The flattening value ( x2: line 123 ) may be replaced by 596.514
-It will save 2 bytes without changing a lot the precision.

-The precision is the same as Andoyer's formula of order 2 for relatively long distances and almost the same as Andoyer's formula of order 3 for very long distances.
-For instance, with L = 0° b = 5° L' = 175° b' = 3° this program gives 18968.13264 km whereas Vincenty's formula gives 18968.13070 km
-So, an error of 1.94 m

-Here are a few maximum errors ( in meters ) for several L' - L values

-The second row = maximum error of Andoyer's formula of order 2 without any term of order 3
-Third row = maximum error with "TGD" above: Andoyer of order 2 + a few of order 3.
-Fourth row = maximum error with complete Andoyer's formula of order 3

-I found them with free42.

\| L'-L \|	90°	150°	155°	160°	165°	170°	171°	172°	173°	174°	175°	176°
And2	0.16	0.52	0.85	1.50	3.03	7.64	9.76	12.5	16.9	23.4	34.4	55.6
And2afew3	0.16	0.20	0.26	0.34	0.51	0.88	1.00	1.15	1.33	1.57	1.94	3.14
And3	~ 0	0.008	0.01	0.03	0.08	0.27	0.39	0.52	0.82	1.36	2.05	4.38

-Thus, if L'-L < 171° , the errors remain smaller than 1 meter ( except for nearly antipodal points ) !

-Surprisingly, if L' - L > 174°, the complete Andoyer's ( or Sodano's ) formula of order 3 may be less accurate than Andoyer's formula of order 2 + main terms of order 3 !

b) Triaxial Ellipsoid ( Without Parametric Coordinates )

-With H = ( L' - L )/2 ; F = ( b + b' )/2 ; G = ( b' - b )/2

S = sin² G cos² H + cos² F sin² H ; µ = 2 Arcsin S^1/2

we define A = ( sin² G cos² F )/S + ( sin² F cos² G )/( 1-S )
and B = ( sin² G cos² F )/S - ( sin² F cos² G )/( 1-S )

d - Δ' = a (f²/4) µ [ (A²-A) ( 1 - µ Cot µ ) - B.(A-1) ( µ / Sin µ - (Sin µ) / µ ) ]

                + a (f³/8) [ A ( 2 µ² Cot µ - 2 µ³ Cot² µ - 8 µ³ / 6 ) + A² ( -5 µ² Cot µ + 5 µ³ Cot² µ + 8 µ³ / 3 )
                         + A³ ( 3 µ² Cot µ - 3 µ³ Cot² µ - 8 µ³ / 6 )
                        + B ( 2 µ² / Sin µ - 2 µ³ Cos µ / Sin² µ ) + A.B ( -5 µ² / Sin µ + 6 µ³ Cos µ / Sin² µ )
                        + A² B ( 3 µ² / Sin µ - 4 µ³ Cos µ / Sin² µ )
                        + B² ( µ³ / Sin² µ ) + A B² ( - µ³ / Sin² µ )

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

49 RCL 08
50 X^2
51 RCL 09
52 X^2
53 RCL 10
54 X^2
55 +
56 +
57 STO 13
58 RCL 05
59 X^2
60 RCL 06
61 X^2
62 RCL 07
63 X^2
64 +
65 +
66 STO 00
67 -
68 RCL 04
69 ACOS
70 STO 02
71 SIN
72 STO 15
73 *
74 RCL 08
75 XEQ 00
76 ST+ X
77 STO 09
78 RCL 00
79 RCL 13
80 +
81 RCL 04
82 *
83 -
84 X#0?
85 /
86 ATAN
87 STO 08
88 2
89 ST/ 01
90 ST/ 11
91 /
92 STO 05
93 XEQ 02
94 STO 06
95 RCL 05
96 90

97 +
98 XEQ 02
99 STO 07
100 RCL 02
101 D-R
102 RCL 08
103 COS
104 STO 05
105 X^2
106 ST+ X
107 RCL 04
108 RCL 15
109 ST* 05
110 *
111 ST* Y
112 -
113 -
114 4
115 /
116 RCL 06
117 ST/ 15
118 RCL 07
119 ST- Y
120 ST+ X
121 /
122 STO 07
123 *
124 RCL 05
125 -
126 +
127 RCL 07
128 *
129 +
130 ST* 15
131 RCL 14
132 ENTER
133 GTO 03
134 LBL 00
135 RCL 05
136 *
137 RCL 09
138 RCL 06
139 *
140 +
141 RCL 10
142 RCL 07
143 *
144 +

145 RTN
146 LBL 01
147 1
148 P-R
149 X<>Y
150 RCL 03
151 X^2
152 *
153 STO 10
154 RDN
155 P-R
156 RCL 01
157 X^2
158 *
159 STO Z
160 X^2
161 X<>Y
162 RCL 02
163 X^2
164 *
165 STO 09
166 X^2
167 RCL 10
168 X^2
169 +
170 +
171 SQRT
172 ST/ 09
173 ST/ 10
174 /
175 STO 08
176 RTN
177 LBL 02
178 RCL 02
179 2
180 /
181 +
182 ENTER
183 SIN
184 ENTER
185 X^2
186 RCL 13
187 *
188 RCL 02
189 R^
190 -
191 SIN
192 ST* Z

193 X^2
194 RCL 00
195 *
196 +
197 X<>Y
198 RCL 09
199 *
200 +
201 SQRT
202 RTN
203 LBL 03
204 RCL 12
205 +
206 2
207 /
208 ST- Y
209 1
210 STO 07
211 P-R
212 X^2
213 STO T
214 RDN
215 X^2
216 STO 08
217 SIGN
218 P-R
219 X^2
220 ST* 08
221 CLX
222 +
223 X^2
224 ST* T
225 -
226 RCL 11
227 COS
228 X^2
229 *
230 -
231 /
232 STO 02
233 RCL 08
234 RCL 07
235 LASTX
236 STO 08
237 -
238 /
239 ST- 02
240 +

241 ST- 07
242 X<> 08
243 SQRT
244 ASIN
245 ST+ X
246 D-R
247 STO 00
248 LASTX
249 1
250 P-R
251 STO 04
252 STO 05
253 STO 09
254 RDN
255 /
256 ST* 05
257 STO 06
258 X^2
259 ST* 09
260 RCL 02
261 *
262 RCL 06
263 ST- 09
264 RCL 09
265 ST* 04
266 ST+ 09
267 4
268 *
269 +
270 RCL 07
271 ST* Z
272 *
273 RCL 09
274 +
275 RCL 08
276 *
277 +
278 RCL 09
279 -
280 RCL 02
281 *
282 RCL 04
283 RCL 00
284 ST+ X
285 X^2
286 3
287 ST* Z
288 /

289 STO Z
290 +
291 STO 09
292 RCL 08
293 *
294 RCL 09
295 ST+ X
296 -
297 RCL 04
298 ST+ Z
299 ST+ Z
300 +
301 RCL 08
302 ST* 05
303 *
304 +
305 RCL 08
306 ST- 05
307 *
308 -
309 RCL 01
310 RCL 03
311 -
312 RCL 01
313 ST+ X
314 /
315 STO 04
316 *
317 RCL 05
318 RCL 06
319 ENTER
320 1/X
321 -
322 RCL 02
323 *
324 +
325 RCL 07
326 *
327 +
328 RCL 04
329 X^2
330 *
331 RCL 00
332 *
333 RCL 01
334 *
335 RCL 15
336 +
337 END

( 427 bytes / SIZE 016 )

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

Where D is the geodesic distance on the triaxial ellipsoid

Example: Calculate the geodesic distance between the Sydney Observatory ( Long = 151°12'17"8 E , Lat = 33°51'41"1 S )
and the Palomar Observatory ( Long = 116°51'50"4 W , Lat = 33°21'22"4 N )

     151.12178 ENTER^
    -33.51411   ENTER^
   -116.51504 ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65755 km                                                        ---Execution time = 25s---
                                            LASTX d = 12138.65875 km = great elliptic distance = R10

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

Notes:

-When line 99 is executed ( STO 07 ), the semi axes of the great ellipse = R15/R06 & R15/R07.
-With the example above, 6378.151662 km & 6368.318461 km

-The precision is not so good as the programs listed in paragraphs 2b) & 2c)
-For example:

-If L = -40° & L' = 120° , the root mean square error = 15 cm and the maximum | error | = 55 cm
-If L = -95° & L' = 65° , rms = 31 cm and max | error | = 87 cm

c) Triaxial Ellipsoid ( With Parametric Coordinates )

-The method employed in paragraph above is not really rigorous !
-A better accuracy will be obtained if we use the same formulas but after transforming the geodetic coordinates into parametric coordinates to calculate µ , A , B ... etc ...

d - Δ' = a (f²/4) µ [ (A²-A) ( 1 - µ Cot µ ) - B.(A-1) ( µ / Sin µ - (Sin µ) / µ ) ]

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

51 RCL 09
52 X^2
53 RCL 10
54 X^2
55 +
56 +
57 STO 13
58 RCL 05
59 X^2
60 RCL 06
61 X^2
62 RCL 07
63 X^2
64 +
65 +
66 STO 00
67 -
68 RCL 04
69 ACOS
70 STO 02
71 SIN
72 STO 15
73 *
74 RCL 08
75 XEQ 00
76 ST+ X
77 STO 09
78 RCL 00
79 RCL 13
80 +
81 RCL 04
82 *
83 -
84 X#0?
85 /
86 ATAN
87 STO 08
88 2
89 ST/ 01
90 ST/ 11
91 /
92 STO 05
93 XEQ 02
94 STO 06
95 RCL 05
96 90
97 +
98 XEQ 02
99 STO 07
100 RCL 02

101 D-R
102 RCL 08
103 COS
104 STO 05
105 X^2
106 ST+ X
107 RCL 04
108 RCL 15
109 ST* 05
110 *
111 ST* Y
112 -
113 -
114 4
115 /
116 RCL 06
117 ST/ 15
118 RCL 07
119 ST- Y
120 ST+ X
121 /
122 STO 07
123 *
124 RCL 05
125 -
126 +
127 RCL 07
128 *
129 +
130 ST* 15
131 GTO 03
132 LBL 00
133 RCL 05
134 *
135 RCL 09
136 RCL 06
137 *
138 +
139 RCL 10
140 RCL 07
141 *
142 +
143 RTN
144 LBL 01
145 1
146 P-R
147 X<>Y
148 RCL 03
149 X^2
150 *

151 STO 10
152 RDN
153 P-R
154 RCL 01
155 X^2
156 *
157 STO Z
158 X^2
159 X<>Y
160 RCL 02
161 X^2
162 *
163 STO 09
164 X^2
165 RCL 10
166 X^2
167 +
168 +
169 SQRT
170 ST/ 09
171 ST/ 10
172 /
173 STO 08
174 RTN
175 LBL 02
176 RCL 02
177 2
178 /
179 +
180 ENTER
181 SIN
182 ENTER
183 X^2
184 RCL 13
185 *
186 RCL 02
187 R^
188 -
189 SIN
190 ST* Z
191 X^2
192 RCL 00
193 *
194 +
195 X<>Y
196 RCL 09
197 *
198 +
199 SQRT
200 RTN

201 LBL 03
202 RCL 14
203 1
204 STO 07
205 P-R
206 RCL 01
207 RCL 03
208 -
209 RCL 01
210 /
211 STO 04
212 SIGN
213 LASTX
214 -
215 STO 05
216 /
217 R-P
218 X<> 12
219 RCL 07
220 P-R
221 RCL 05
222 /
223 R-P
224 RDN
225 +
226 2
227 ST/ 04
228 /
229 ST- Y
230 RCL 07
231 P-R
232 X^2
233 STO T
234 RDN
235 X^2
236 STO 08
237 SIGN
238 P-R
239 X^2
240 ST* 08
241 CLX
242 +
243 X^2
244 ST* T
245 -
246 RCL 11
247 COS
248 X^2
249 *
250 -

251 /
252 STO 02
253 RCL 08
254 RCL 07
255 LASTX
256 STO 08
257 -
258 /
259 ST- 02
260 +
261 ST- 07
262 X<> 08
263 SQRT
264 ASIN
265 ST+ X
266 D-R
267 STO 00
268 LASTX
269 1
270 P-R
271 STO 05
272 STO 09
273 STO 13
274 RDN
275 /
276 ST* 05
277 STO 06
278 X^2
279 ST* 09
280 RCL 02
281 *
282 RCL 06
283 ST- 09
284 RCL 09
285 ST+ 09
286 ST* 13
287 4
288 *
289 +
290 RCL 07
291 ST* Z
292 *
293 RCL 09
294 +
295 RCL 08
296 *
297 +
298 RCL 09
299 -
300 RCL 02

301 *
302 RCL 13
303 RCL 00
304 ST* 01
305 ST+ X
306 X^2
307 3
308 ST* Z
309 /
310 STO Z
311 +
312 STO 09
313 RCL 08
314 *
315 RCL 09
316 ST+ X
317 -
318 RCL 13
319 ST+ Z
320 ST+ Z
321 +
322 RCL 08
323 ST* 05
324 *
325 +
326 RCL 08
327 ST- 05
328 *
329 -
330 RCL 04
331 *
332 RCL 05
333 RCL 06
334 ENTER
335 1/X
336 -
337 RCL 02
338 *
339 +
340 RCL 07
341 *
342 +
343 RCL 04
344 X^2
345 *
346 RCL 01
347 *
348 RCL 15
349 +
350 END

( 444 bytes / SIZE 016 )

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

Where D is the geodesic distance on the triaxial ellipsoid

Example: Calculate the geodesic distance between the Sydney Observatory ( Long = 151°12'17"8 E , Lat = 33°51'41"1 S )
and the Palomar Observatory ( Long = 116°51'50"4 W , Lat = 33°21'22"4 N )

     151.12178 ENTER^
    -33.51411   ENTER^
   -116.51504 ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65755 km                                                        ---Execution time = 27s---
                             LASTX = d = 12138.65875 km = great elliptic distance = R15

Notes:

-When line 99 is executed ( STO 07 ), the semi axes of the great ellipse = R15/R06 & R15/R07.
-With the example above, 6378.151662 km & 6368.318461 km

-The precision is the same as the programs listed in paragraphs 2b) & 2c) !!!
-For example:

-If L = -40° & L' = 120° , the root mean square error = 10 cm and the maximum | error | = 48 cm
-If L = -95° & L' = 65° , rms = 19 cm and max | error | = 55 cm

4°) With Vincenty's Formula

-Vincenty's formula is probably the best method for an ellipsoid of revolution.
-The following version of "TGD" - renamed "TGDV"- employs this method 3 times.

Data Registers: R00 thru R21: temp
Flags: /
Subroutines: /

-Line 47 is a three-byte GTO 14

01 LBL "TGDV"
02 DEG
03 HR
04 STO 21
05 RDN
06 HR
07 14.92911
08 STO 00
09 +
10 STO 20
11 RDN
12 HR
13 STO 19
14 X<>Y
15 HR
16 RCL 00
17 +
18 STO 18
19 6378.172
20 STO 01
21 .07
22 -
23 STO 02
24 6356.752314
25 STO 03
26 XEQ 01
27 STO 17
28 RCL 01
29 RCL 02
30 +
31 2
32 /
33 STO 01
34 STO 02
35 XEQ 01
36 ST- 17
37 RCL 01
38 STO 06
39 RCL 03
40 STO 07
41 RCL 21
42 RCL 19
43 RCL 20
44 RCL 18
45 XEQ 13
46 ENTER
47 GTO 14
48 LBL 01
49 RCL 18
50 RCL 19
51 XEQ 02
52 RCL 08
53 STO 05
54 RCL 09
55 STO 06
56 RCL 10
57 STO 07
58 RCL 20
59 RCL 21
60 XEQ 02
61 RCL 05
62 RCL 08
63 *
64 RCL 06

65 RCL 09
66 *
67 RCL 07
68 RCL 10
69 *
70 +
71 +
72 STO 00
73 ACOS
74 STO 04
75 RCL 01
76 ST/ 05
77 ST/ 08
78 RCL 02
79 ST/ 06
80 ST/ 09
81 RCL 03
82 ST/ 07
83 ST/ 10
84 RCL 05
85 X^2
86 RCL 06
87 X^2
88 RCL 07
89 X^2
90 +
91 +
92 STO 11
93 RCL 08
94 X^2
95 RCL 09
96 X^2
97 RCL 10
98 X^2
99 +
100 +
101 X<> 10
102 RCL 07
103 *
104 RCL 06
105 RCL 09
106 *
107 RCL 05
108 RCL 08
109 *
110 +
111 +
112 STO 09
113 RCL 00
114 *
115 ST+ X
116 RCL 10
117 -
118 RCL 11
119 RCL 00
120 X^2
121 *
122 -
123 RCL 11
124 RCL 04
125 SIN
126 STO 08
127 ST/ Z
128 *

129 +
130 RCL 09
131 RCL 11
132 RCL 00
133 *
134 -
135 ST+ X
136 X<>Y
137 X#0?
138 /
139 ATAN
140 2
141 /
142 90
143 STO 00
144 X<>Y
145 X<0?
146 +
147 STO 05
148 XEQ 03
149 STO 06
150 RCL 05
151 RCL 00
152 +
153 XEQ 03
154 STO 07
155 RCL 06
156 X>Y?
157 GTO 00
158 X<> 07
159 STO 06
160 RCL 00
161 ST+ 05
162 LBL 00
163 PI
164 R-D
165 RCL 05
166 X=Y?
167 CLX
168 STO 05
169 CHS
170 1
171 P-R
172 RCL 07
173 RCL 06
174 /
175 X^2
176 STO 09
177 *
178 R-P
179 CLX
180 RCL 04
181 RCL 05
182 -
183 1
184 P-R
185 RCL 09
186 *
187 R-P
188 CLX
189 ENTER
190 LBL 13
191 -
192 360

193 MOD
194 PI
195 R-D
196 STO 16
197 X>Y?
198 CLX
199 ST+ X
200 -
201 STO 11
202 STO 04
203 CLX
204 RCL 07
205 CHS
206 RCL 06
207 ST+ Y
208 /
209 STO 08
210 SIGN
211 LASTX
212 -
213 STO 15
214 P-R
215 LASTX
216 /
217 R-P
218 RDN
219 STO 09
220 CLX
221 RCL 15
222 P-R
223 LASTX
224 /
225 R-P
226 X<>Y
227 STO 10
228 LBL 12
229 RCL 16
230 RCL 04
231 X>Y?
232 ACOS
233 RCL 10
234 COS
235 STO 05
236 STO 12
237 RCL 04
238 SIN
239 *
240 STO 14
241 RCL 09
242 COS
243 ST* 05
244 ST* 14
245 RCL 10
246 SIN
247 STO 13
248 *
249 RCL 09
250 SIN
251 ST* 13
252 RCL 12
253 *
254 RCL 04
255 COS
256 ST* 05

257 *
258 -
259 R-P
260 RCL 05
261 RCL 13
262 +
263 R-P
264 X<> 14
265 X<>Y
266 STO 05
267 SIN
268 X#0?
269 /
270 ASIN
271 STO 12
272 COS
273 X^2
274 RCL 05
275 COS
276 RCL 13
277 ENTER
278 +
279 R^
280 X=0?
281 SIGN
282 /
283 -
284 STO 13
285 CLX
286 3
287 *
288 4
289 X<>Y
290 -
291 RCL 08
292 *
293 4
294 +
295 *
296 RCL 08
297 *
298 16
299 /
300 RCL 05
301 SIN
302 RCL 13
303 RCL Z
304 *
305 *
306 R-D
307 RCL 05
308 +
309 RCL 12
310 SIN
311 *
312 RCL 08
313 *
314 1
315 R^
316 -
317 *
318 RCL 11
319 +
320 ENTER

321 X<> 04
322 -
323 ABS
324 E-7
325 X<Y?
326 GTO 12
327 2
328 RCL 08
329 ST- Y
330 *
331 RCL 12
332 COS
333 RCL 15
334 /
335 X^2
336 *
337 12
338 RCL Y
339 5
340 *
341 -
342 *
343 64
344 -
345 *
346 256
347 ST- Y
348 /
349 STO 14
350 CLX
351 37
352 *
353 64
354 -
355 *
356 512
357 /
358 4
359 1/X
360 +
361 *
362 RCL 13
363 X^2
364 ST+ X
365 1
366 -
367 RCL 05
368 COS
369 4
370 /
371 *
372 *
373 RCL 13
374 +
375 *
376 RCL 05
377 SIN
378 *
379 RCL 05
380 D-R
381 -
382 RCL 14
383 *
384 RCL 15

385 *
386 RCL 06
387 *
388 RTN
389 LBL 02
390 1
391 P-R
392 X<>Y
393 RCL 03
394 X^2
395 *
396 STO 10
397 X^2
398 STO 04
399 RDN
400 P-R
401 RCL 01
402 X^2
403 *
404 STO 08
405 X^2
406 ST+ 04
407 X<>Y
408 RCL 02
409 X^2
410 *
411 STO 09
412 X^2
413 RCL 04
414 +
415 SQRT
416 ST/ 08
417 ST/ 09
418 ST/ 10
419 RTN
420 LBL 03
421 ENTER
422 SIN
423 STO 07
424 X^2
425 RCL 10
426 *
427 X<>Y
428 RCL 04
429 -
430 SIN
431 ST* 07
432 X^2
433 RCL 11
434 *
435 +
436 RCL 07
437 RCL 09
438 ST+ X
439 *
440 -
441 SQRT
442 RCL 08
443 X<>Y
444 /
445 RTN
446 LBL 14
447 RCL 17
448 +
449 END

( 564 bytes / SIZE 022 )

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

Where D is the geodesic distance on the triaxial ellipsoid
and d is the geodesic distance on the ellipsoid of revolution ( WGS 84 )

Example: Calculate the geodesic distance between the Sydney Observatory ( Long = 151°12'17"8 E , Lat = 33°51'41"1 S )
and the Palomar Observatory ( Long = 116°51'50"4 W , Lat = 33°21'22"4 N )

     151.12178 ENTER^
    -33.51411   ENTER^
   -116.51504 ENTER^
     33.21224   XEQ "TGDV"   >>>>   D = 12138.65757 km                                                      ---Execution time = 119s---
                                                X<>Y d = 12138.68433 km

Notes:

-This version gives the best accuarcy, provided the points are not ( almost ) antipodal.
-For nearly antipodal points, the routine stops and the HP41 displays DATA ERROR ( line 235 )

5°) Minimizing the Sum of 2 Distances

-No one of the programs above works well for nearly antipodal points.
-The following one uses the golden search to minimize the sum of 2 geodesic distances:

Let M(L,b) & N(L',b') , "D+D" minimizes the sum D = Dist(M,P) + Dist(P,N)

where P(L",b") with L" = the longitude of the point on a sphere in the middle of the geodesic MN on the sphere.

-The program starts with b" = -90° and +90° and employs the golden search to find the minimum D-value.

Data Registers: R00 thru R35: temp
Flags: /
Subroutine: "TGD" ( cf paragraph 2-b )

01 LBL "D+D"
02 STO 26
03 RDN
04 STO 25
05 RDN
06 STO 24
07 X<>Y
08 STO 23
09 R^
10 X<>Y
11 HMS-
12 HR
13 RCL 26
14 HR
15 COS
16 P-R
17 RCL 24
18 HR
19 COS
20 +
21 R-P
22 CLX

23 RCL 23
24 HR
25 +
26 HMS
27 STO 27
28 GTO 00
29 LBL 12
30 STO 35
31 RCL 27
32 X<>Y
33 RCL 25
34 RCL 26
35 XEQ "TGD"
36 STO 31
37 RCL 23
38 RCL 24
39 RCL 27
40 RCL 35
41 XEQ "TGD"
42 RCL 31
43 +
44 RTN

45 LBL 00
46 90
47 ENTER
48 STO 29
49 CHS
50 STO 28
51 -
52 5
53 SQRT
54 1
55 +
56 2
57 /
58 STO 34
59 /
60 RCL 28
61 HR
62 +
63 HMS
64 STO 32
65 XEQ 12
66 STO 33

67 LBL 10
68 RCL 32
69 RCL 28
70 HMS-
71 HR
72 RCL 34
73 /
74 RCL 28
75 HR
76 +
77 HMS
78 STO 30
79 XEQ 12
80 STO 31
81 VIEW 31
82 RCL 33
83 X<Y?
84 GTO 02
85 RCL 30
86 X<> 32
87 STO 29
88 RCL 31

89 STO 33
90 GTO 04
91 LBL 02
92 RCL 29
93 STO 28
94 RCL 30
95 STO 29
96 LBL 04
97 E-4
99 RCL 28
100 RCL 30
101 HMS-
102 HR
103 ABS
104 X>Y?
105 GTO 10
106 RCL 30
107 RCL 31
108 CLD
109 END

( 180 bytes / SIZE 036 )

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

Where D is the geodesic distance on the triaxial ellipsoid
and b" is the latitude of the "midpoint"

Example1: With (L,b) = (0°,0°) & (L',b') = (179.51°,0°)

        0         ENTER^ ENTER^
    179.51    ENTER^
        0         XEQ "D+D" >>>>   D = 20001.89956 km
                                        X<>Y b" ~ -75°37'

Example2: With (L,b) = (-40°,-40°) & (L',b') = (120°,50°)

   40 CHS ENTER^ ENTER^
      120       ENTER^
       50        XEQ "D+D" >>>>   D = 18095.49451 km
                                        X<>Y b" ~ 24°58'

Notes:

-The HP41 displays the successive approximations.

-The angle b" is not computed accurately: it's not necessary to get a good precision for the distance.
-You can press XEQ 10 to make another loop, or - for instance - replace lines 98-99 by a larger number.

-This program is very slow with a true HP41, so a good emulator is not superfluous !
-With free42, you get the results almost instantaneously...

-This program may also be used for very long geodesics when the points are not almost antipodal,
but when Andoyer's or Sodano's formulas give an error of a few meters.
-Since formulas of order 3 is very accurate when the distance is about 10000 km, the sum will remain very accurate too !

6°) Other Formulae for the Great Elliptic Distance

a) Formula of order 1 + Andoyer's Formula of order 2

-This variant employs Andoyer's formula of order 2 to compute the geodesic distance on WGS84 ellipsoid of revolution,
but the following formula of order 1 to calculate the distances Δ = MN on the great ellipses:

With u = OM / OM = ( x , y , z ) , v = ON / ON = ( x' , y' ,z' ) , W = angle ( OM , ON )

and f = (a-b) / a & f ' = (a-c) / c

Δ ~ a [ W - { ( f y² + f y'² + f ' z² + f ' z'² ) ( W/2 - (1/4) Sin 2.W ) + ( f y.y' + f ' z.z' ) ( Sin W - W Cos W ) } / Sin² W ]

-This formula is an approximation of the integral: Δ = §₀^W [ R² + ( dR/dµ )² ]^1/2 dµ

with R / Sin W = 1 / ( A² + B² + C² )^1/2 where

     A = [ x Sin ( W - µ ) + x' Sin µ ] / a
     B = [ y Sin ( W - µ ) + y' Sin µ ] / b
     C = [ z Sin ( W - µ ) + z' Sin µ ] / c

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

-Line 36 is a three-byte GTO 03

01 LBL "TGD"
02 DEG
03 HR
04 STO 04
05 RDN
06 HR
07 14.929
08 STO 01
09 +
10 STO 03
11 RDN
12 HR
13 STO 02
14 X<>Y
15 HR
16 ST+ 01
17 6378.172
18 STO 05
19 .07
20 -
21 STO 06
22 21.349686
23 -
24 STO 07
25 XEQ 01
26 STO 14
27 RCL 05
28 RCL 06
29 +
30 2
31 /
32 STO 05
33 STO 06
34 XEQ 01
35 ST- 14
36 GTO 03

37 LBL 01
38 RCL 01
39 RCL 02
40 XEQ 02
41 STO 08
42 RCL 11
43 STO 09
44 RCL 12
45 STO 10
46 RCL 03
47 RCL 04
48 XEQ 02
49 RCL 08
50 *
51 RCL 09
52 RCL 11
53 *
54 STO 13
55 +
56 RCL 10
57 RCL 12
58 *
59 STO 08
60 +
61 STO 00
62 RCL 10
63 X^2
64 RCL 12
65 X^2
66 +
67 RCL 05
68 RCL 07
69 -
70 ST* 08
71 STO 10
72 *

73 RCL 09
74 X^2
75 RCL 11
76 X^2
77 +
78 RCL 05
79 RCL 06
80 -
81 ST* 13
82 *
83 +
84 RCL 00
85 ENTER
86 ACOS
87 D-R
88 STO 11
89 LASTX
90 SIN
91 STO 12
92 /
93 ST* 00
94 -
95 *
96 2
97 /
98 RCL 00
99 RCL 08
100 RCL 13
101 +
102 ST* Y
103 -
104 +
105 RCL 12
106 /
107 RCL 11
108 RCL 05

109 ST/ 10
110 *
111 +
112 RTN
113 LBL 02
114 1
115 P-R
116 X<>Y
117 RCL 07
118 X^2
119 *
120 STO 12
121 X^2
122 X<> Z
123 X<>Y
124 P-R
125 RCL 05
126 X^2
127 *
128 STO 13
129 X^2
130 X<>Y
131 RCL 06
132 X^2
133 *
134 STO 11
135 X^2
136 +
137 +
138 SQRT
139 ST/ 11
140 ST/ 12
141 ST/ 13
142 RCL 13
143 RTN
144 LBL 03

145 RCL 04
146 RCL 02
147 RCL 03
148 RCL 01
149 -
150 2
151 ST/ 10
152 /
153 COS
154 X^2
155 RDN
156 +
157 2
158 /
159 ST- Y
160 1
161 P-R
162 X^2
163 STO 02
164 RDN
165 X^2
166 STO 01
167 SIGN
168 P-R
169 X^2
170 ST* 01
171 RDN
172 X^2
173 ST* 02
174 RCL Z
175 -
176 *
177 +
178 ST/ 02
179 STO 03
180 SQRT

181 ASIN
182 D-R
183 RCL 03
184 1
185 RCL 03
186 -
187 ST/ 01
188 ST- 03
189 *
190 SQRT
191 /
192 STO 00
193 CLX
194 RCL 02
195 RCL 01
196 ST- 02
197 +
198 STO 01
199 RCL 03
200 RCL 00
201 ST* 03
202 /
203 7.5
204 *
205 STO 04
206 LASTX
207 -
208 2
209 /
210 RCL 03
211 +
212 *
213 4
214 +
215 RCL 03
216 -

217 6
218 RCL 00
219 ST/ Y
220 +
221 RCL 02
222 *
223 STO Z
224 -
225 RCL 01
226 *
227 +
228 RCL 02
229 X^2
230 RCL 04
231 *
232 -
233 RCL 10
234 *
235 RCL 01
236 -
237 RCL 02
238 3
239 *
240 RCL 00
241 /
242 -
243 RCL 10
244 *
245 *
246 +
247 RCL 05
248 ST+ X
249 *
250 RCL 14
251 +
252 END

( 318 bytes / SIZE 015 )

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

Where D is the geodesic distance on the triaxial ellipsoid

Example: Calculate the geodesic distance between the Sydney Observatory ( Long = 151°12'17"8 E , Lat = 33°51'41"1 S )
and the Palomar Observatory ( Long = 116°51'50"4 W , Lat = 33°21'22"4 N )

     151.12178 ENTER^
    -33.51411   ENTER^
   -116.51504 ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65746 km                                                      ---Execution time = 25s---
                                 LASTX -     d = 12138.68425 km

where d is the geodesic distance on the ellipsoid of revolution ( WGS 84 )

-The exact results are D = 12138.65757 & d = 12138.68433 km ( rounded to 5 decimals )

Notes:

-The distance is computed by D ~ d + Δ₁ - Δ₂
-It might also be calculated by D ~ d . Δ₁/Δ₂ ( replace line 35 by ST/ 14 and line 251 by * )
-With example above, the results are the same.

-This program is shorter and slightly faster than the other versions.
-Though the precision is not the best, it's often accurate enough
-However, a formula of order 3 will give a better precision:

b) Formula of order 1 + Andoyer's Formula of order 3

-The same method is used for the great ellipse arc distances, but Andoyer's formula of order 3 is employed to compute the geodesic distance on WGS84 ellipsoid of revolution:

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

-Line 36 is a three-byte GTO 03

01 LBL "TGD"
02 DEG
03 HR
04 STO 14
05 RDN
06 HR
07 14.929
08 STO 11
09 +
10 STO 13
11 RDN
12 HR
13 STO 12
14 X<>Y
15 HR
16 ST+ 11
17 6378.172
18 STO 05
19 .07
20 -
21 STO 06
22 21.349686
23 -
24 STO 07
25 XEQ 01
26 STO 04
27 RCL 05
28 RCL 06
29 +
30 2
31 /
32 STO 05
33 STO 06
34 XEQ 01
35 ST- 04
36 GTO 03
37 LBL 01
38 RCL 11
39 RCL 12
40 XEQ 02
41 STO 08
42 RCL 02
43 STO 09
44 RCL 03
45 STO 10
46 RCL 13
47 RCL 14
48 XEQ 02
49 RCL 08

50 *
51 RCL 02
52 RCL 09
53 *
54 STO 01
55 +
56 RCL 03
57 RCL 10
58 *
59 STO 08
60 +
61 STO 00
62 RCL 03
63 X^2
64 RCL 10
65 X^2
66 +
67 RCL 05
68 RCL 07
69 -
70 ST* 08
71 STO 10
72 *
73 RCL 02
74 X^2
75 RCL 09
76 X^2
77 +
78 RCL 05
79 STO 09
80 ST/ 10
81 RCL 06
82 -
83 ST* 01
84 *
85 +
86 RCL 00
87 ENTER
88 ACOS
89 D-R
90 ST* 09
91 LASTX
92 SIN
93 STO 02
94 /
95 ST* 00
96 -
97 *
98 2

99 /
100 RCL 00
101 RCL 01
102 RCL 08
103 +
104 ST* Y
105 -
106 +
107 RCL 02
108 /
109 RCL 09
110 +
111 RTN
112 LBL 02
113 1
114 P-R
115 X<>Y
116 RCL 07
117 X^2
118 *
119 STO 03
120 RDN
121 P-R
122 RCL 05
123 X^2
124 *
125 STO 01
126 X^2
127 X<>Y
128 RCL 06
129 X^2
130 *
131 STO 02
132 X^2
133 RCL 03
134 X^2
135 +
136 +
137 SQRT
138 ST/ 01
139 ST/ 02
140 ST/ 03
141 RCL 01
142 RTN
143 LBL 03
144 RCL 12
145 1
146 STO 08
147 P-R

148 LASTX
149 RCL 10
150 -
151 STO 07
152 /
153 R-P
154 X<> 14
155 RCL 08
156 P-R
157 RCL 07
158 /
159 R-P
160 RDN
161 +
162 2
163 ST/ 10
164 /
165 ST- Y
166 RCL 08
167 P-R
168 X^2
169 STO 02
170 STO T
171 RDN
172 X^2
173 STO 01
174 SIGN
175 P-R
176 X^2
177 ST* 01
178 RDN
179 X^2
180 ST* 02
181 -
182 RCL 11
183 RCL 13
184 -
185 2
186 /
187 COS
188 X^2
189 *
190 -
191 ST/ 02
192 RCL 02
193 RCL 01
194 RCL 08
195 R^
196 STO 01

197 -
198 /
199 ST- 02
200 +
201 X<> 01
202 SQRT
203 ASIN
204 ST+ X
205 D-R
206 STO 00
207 LASTX
208 RCL 08
209 P-R
210 STO 03
211 STO 06
212 RDN
213 STO 11
214 /
215 STO 07
216 ST/ Z
217 X^2
218 ST* 03
219 RCL 08
220 +
221 ENTER
222 R^
223 ST+ 08
224 ST- Z
225 2
226 /
227 STO 09
228 -
229 RCL 01
230 *
231 -
232 RCL 11
233 X^2
234 1.5
235 ST/ Y
236 FRC
237 -
238 RCL 07
239 /
240 RCL 02
241 *
242 +
243 RCL 02
244 *
245 RCL 06
246 X^2

247 RCL 07
248 ST- 03
249 ST+ X
250 /
251 RCL 03
252 4
253 ST/ 08
254 *
255 RCL 06
256 +
257 RCL 07
258 +
259 STO Z
260 -
261 RCL 01
262 *
263 +
264 RCL 03
265 ST+ X
266 STO 11
267 +
268 RCL 01
269 *
270 +
271 RCL 11
272 -
273 RCL 02
274 *
275 RCL 03
276 RCL 06
277 *
278 STO 03
279 RCL 00
280 ST+ X
281 X^2
282 3
283 ST* Z
284 /
285 STO Z
286 +
287 RCL 08
288 +
289 STO 11
290 RCL 01
291 *
292 RCL 11
293 ST+ X
294 -
295 RCL 03
296 ST+ Z

297 ST+ Z
298 +
299 RCL 01
300 *
301 +
302 RCL 01
303 *
304 -
305 RCL 10
306 *
307 RCL 08
308 RCL 06
309 RCL 07
310 *
311 STO Z
312 -
313 RCL 01
314 *
315 +
316 RCL 02
317 ST* 09
318 RCL 07
319 ST- 09
320 *
321 -
322 RCL 01
323 *
324 +
325 RCL 02
326 RCL 09
327 *
328 -
329 RCL 10
330 *
331 RCL 02
332 RCL 07
333 /
334 -
335 RCL 01
336 -
337 RCL 10
338 *
339 RCL 00
340 ST* Y
341 +
342 RCL 05
343 *
344 RCL 04
345 +
346 END

( 429 bytes / SIZE 015 )

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

Where D is the geodesic distance on the triaxial ellipsoid

Example: Calculate the geodesic distance between the Sydney Observatory ( Long = 151°12'17"8 E , Lat = 33°51'41"1 S )
and the Palomar Observatory ( Long = 116°51'50"4 W , Lat = 33°21'22"4 N )

     151.12178 ENTER^
    -33.51411   ENTER^
   -116.51504 ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65753 km                                                      ---Execution time = 30s---
                                 LASTX -     d = 12138.68432 km

where d is the geodesic distance on the ellipsoid of revolution ( WGS 84 )

Notes:

-Surprisingly, the precision is almost as good as the results given by the program listed in §2°)b)
-For example, when the difference in the longitudes is 165°, the error can reach about 83 cm instead of 75 cm.

c) Formula of order 2 + Andoyer's Formula of order 3

-A formula of order 2 is used for the great ellipse arc distances and Andoyer's formula of order 3 is employed to compute the geodesic distance on WGS84 ellipsoid of revolution:

m = [ y² ( 2.f + 3.f² ) + z² ( 2.f ' + 3.f '² ) ] / ( 2 Sin² W )
With n = [ y'² ( 2.f + 3.f² ) + z'² ( 2.f ' + 3.f '² ) ] / ( 2 Sin² W ) [ In fact, 2.f + 3.f² & 2.f ' + 3.f '² may be replaced by (a²/b² - 1)/2 and (a²/c² - 1)/2 ]
p = [ y.y' ( 2.f + 3.f² ) + z.z' ( 2.f ' + 3.f '² ) ] / ( Sin² W )

Δ / (a.Sin W) ~ W / Sin W + ( p²/2 - m - n ) ( W/Sin W - Cos W ) / 2 + (p/2) ( W Cot W - 1 )
+ ( m² + n² ) [ (13/16) ( W/Sin W - Cos W ) + (1/8) Sin² W Cos W ) ]
+ 3 ( p²/2 ) [ (-1/4) W Sin W+ (3/8) ( W/Sin W - Cos W ) ]
+ m n [ (-7/4) W Sin W + (13/8) ( W/Sin W -Cos W ]
+ p ( m + n ) [ (3/4) Sin⁴ W + ( Cos W ) { (-13/8) ( W/Sin W - Cos W ) + (3/4) Sin² W Cos W } ]

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

-Line 36 is a three-byte GTO 03

01 LBL "TGD"
02 DEG
03 HR
04 STO 14
05 RDN
06 HR
07 14.92911
08 STO 11
09 +
10 STO 13
11 RDN
12 HR
13 STO 12
14 X<>Y
15 HR
16 ST+ 11
17 6378.172
18 STO 05
19 .07
20 -
21 STO 06
22 21.349686
23 -
24 STO 07
25 XEQ 01
26 STO 15
27 RCL 05
28 RCL 06
29 +
30 2
31 /
32 STO 05
33 STO 06
34 XEQ 01
35 ST- 15
36 GTO 03
37 LBL 01
38 RCL 11
39 RCL 12
40 XEQ 02
41 STO 08
42 RCL 02
43 STO 09
44 RCL 03
45 STO 10
46 RCL 13
47 RCL 14
48 XEQ 02
49 RCL 08
50 *
51 RCL 02
52 RCL 09
53 *
54 STO 01
55 +
56 RCL 03
57 RCL 10
58 *
59 STO 08
60 +
61 STO 00
62 RCL 09
63 X^2
64 RCL 05

65 RCL 06
66 /
67 X^2
68 1
69 -
70 ST* 01
71 STO 09
72 *
73 RCL 05
74 RCL 07
75 /
76 X^2
77 1
78 -
79 STO 04
80 ST* 08
81 RCL 10
82 X^2
83 *
84 +
85 STO 10
86 RCL 02
87 X^2
88 RCL 09
89 *
90 RCL 03
91 X^2
92 RCL 04
93 *
94 +
95 STO 02
96 RCL 01
97 RCL 08
98 +
99 RCL 00
100 ACOS
101 STO 08
102 SIN
103 STO 03
104 X^2
105 ST+ X
106 ST/ 02
107 ST/ 10
108 /
109 STO 09
110 RCL 08
111 D-R
112 STO 08
113 RCL 03
114 ST* 08
115 /
116 STO Y
117 RCL 00
118 -
119 2
120 /
121 STO 04
122 RCL 09
123 X^2
124 ST+ X
125 STO 01
126 RCL 02
127 -
128 RCL 10

129 -
130 *
131 +
132 X<>Y
133 RCL 00
134 *
135 RCL 09
136 ST* Y
137 -
138 +
139 RCL 04
140 3
141 *
142 RCL 08
143 -
144 RCL 01
145 *
146 .75
147 *
148 +
149 RCL 04
150 13
151 *
152 STO 01
153 RCL 08
154 7
155 *
156 -
157 RCL 02
158 *
159 RCL 10
160 *
161 4
162 /
163 +
164 RCL 01
165 RCL 00
166 RCL 03
167 X^2
168 *
169 ST+ 01
170 3
171 *
172 -
173 RCL 00
174 *
175 RCL 03
176 X^2
177 X^2
178 3
179 *
180 -
181 RCL 09
182 *
183 RCL 02
184 RCL 10
185 +
186 *
187 2
188 /
189 -
190 RCL 01
191 RCL 02
192 X^2

193 RCL 10
194 X^2
195 +
196 *
197 8
198 /
199 +
200 RCL 03
201 *
202 RCL 05
203 *
204 RTN
205 LBL 02
206 1
207 P-R
208 X<>Y
209 RCL 07
210 X^2
211 *
212 STO 03
213 RDN
214 P-R
215 RCL 05
216 X^2
217 *
218 STO 01
219 X^2
220 X<>Y
221 RCL 06
222 X^2
223 *
224 STO 02
225 X^2
226 RCL 03
227 X^2
228 +
229 +
230 SQRT
231 ST/ 01
232 ST/ 02
233 ST/ 03
234 RCL 01
235 RTN
236 LBL 03
237 RCL 05
238 RCL 07
239 -
240 RCL 05
241 /
242 STO 04
243 RCL 12
244 1
245 STO 08
246 P-R
247 LASTX
248 R^
249 -
250 STO 07
251 /
252 R-P
253 X<> 14
254 RCL 08
255 P-R
256 RCL 07

257 /
258 R-P
259 RDN
260 +
261 2
262 ST/ 04
263 /
264 ST- Y
265 RCL 08
266 P-R
267 X^2
268 STO 02
269 STO T
270 RDN
271 X^2
272 STO 01
273 SIGN
274 P-R
275 X^2
276 ST* 01
277 RDN
278 X^2
279 ST* 02
280 -
281 RCL 11
282 RCL 13
283 -
284 2
285 /
286 COS
287 X^2
288 *
289 -
290 ST/ 02
291 RCL 02
292 RCL 01
293 RCL 08
294 R^
295 STO 01
296 -
297 /
298 ST- 02
299 +
300 X<> 01
301 SQRT
302 ASIN
303 ST+ X
304 D-R
305 STO 00
306 LASTX
307 RCL 08
308 P-R
309 STO 03
310 STO 06
311 RDN
312 STO 10
313 /
314 STO 07
315 ST/ Z
316 X^2
317 ST* 03
318 RCL 08
319 +

320 ENTER
321 R^
322 ST+ 08
323 ST- Z
324 2
325 /
326 STO 09
327 -
328 RCL 01
329 *
330 -
331 RCL 10
332 X^2
333 1.5
334 ST/ Y
335 FRC
336 -
337 RCL 07
338 /
339 RCL 02
340 *
341 +
342 RCL 02
343 *
344 RCL 06
345 X^2
346 RCL 07
347 ST- 03
348 ST+ X
349 /
350 RCL 03
351 4
352 ST/ 08
353 *
354 RCL 06
355 +
356 RCL 07
357 +
358 STO Z
359 -
360 RCL 01
361 *
362 +
363 RCL 03
364 ST+ X
365 STO 10
366 +
367 RCL 01
368 *
369 +
370 RCL 10
371 -
372 RCL 02
373 *
374 RCL 03
375 RCL 06
376 *
377 STO 03
378 RCL 00
379 ST+ X
380 X^2
381 3
382 ST* Z

383 /
384 STO Z
385 +
386 RCL 08
387 +
388 STO 10
389 RCL 01
390 *
391 RCL 10
392 ST+ X
393 -
394 RCL 03
395 ST+ Z
396 ST+ Z
397 +
398 RCL 01
399 *
400 +
401 RCL 01
402 *
403 -
404 RCL 04
405 *
406 RCL 08
407 RCL 06
408 RCL 07
409 *
410 STO Z
411 -
412 RCL 01
413 *
414 +
415 RCL 02
416 ST* 09
417 RCL 07
418 ST- 09
419 *
420 -
421 RCL 01
422 *
423 +
424 RCL 02
425 RCL 09
426 *
427 -
428 RCL 04
429 *
430 RCL 02
431 RCL 07
432 /
433 -
434 RCL 01
435 -
436 RCL 04
437 *
438 RCL 00
439 ST* Y
440 +
441 RCL 05
442 *
443 RCL 15
444 +
445 END

( 537 bytes / SIZE 016 )

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

Where D is the geodesic distance on the triaxial ellipsoid

Example: Calculate the geodesic distance between the Sydney Observatory ( Long = 151°12'17"8 E , Lat = 33°51'41"1 S )
and the Palomar Observatory ( Long = 116°51'50"4 W , Lat = 33°21'22"4 N )

     151.12178 ENTER^
    -33.51411   ENTER^
   -116.51504 ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65756 km                                                      ---Execution time = 38s---
                                 LASTX -    d = 12138.68432 km

where d is the geodesic distance on the ellipsoid of revolution ( WGS 84 )

Note:

-Though the formula of order 2 doesn't give a very good precision for the great elliptic distances,
the differences between the great elliptic distances is accurately computed !

7°) Geodesic + Great Elliptic Distance

-With a triaxial ellipsoid, we can use programs listed above ( cf §1) §2b) §3b) §4) )
-But with an ellipsoid of revolution, a shoter & faster program may also be written.

-The routine hereunder employs Andoyer's formula of order 3 for the geodesic distance
and a formula of order 2 for the great elliptic distance ( cf "Great Elliptic Arc Length for the HP41" )

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

01 LBL "TGD+"
02 DEG
03 HR
04 STO 01
05 X<> T
06 HMS-
07 HR
08 2
09 /
10 SIN
11 X^2
12 STO 00
13 X<>Y
14 HR
15 1
16 STO 11
17 P-R
18 LASTX
19 298.25722
20 STO 04
21 ST+ 04
22 1/X
23 -
24 STO 05
25 /
26 R-P
27 X<> 01
28 RCL 11
29 P-R
30 RCL 05
31 /
32 R-P
33 RDN
34 +
35 2

36 /
37 ST- Y
38 RCL 11
39 P-R
40 X^2
41 STO 02
42 X<>Y
43 X^2
44 STO 01
45 RDN
46 X<>Y
47 RCL 11
48 P-R
49 X^2
50 ST* 01
51 RDN
52 X^2
53 ST* 02
54 STO Z
55 -
56 RCL 00
57 *
58 +
59 ST/ 02
60 RCL 02
61 RCL 01
62 RCL 11
63 R^
64 STO 01
65 -
66 /
67 ST- 02
68 +
69 X<> 01
70 SQRT

71 ASIN
72 ST+ X
73 D-R
74 STO 00
75 LASTX
76 RCL 11
77 P-R
78 STO 03
79 STO 06
80 STO 08
81 RDN
82 STO 05
83 /
84 STO 07
85 ST/ 08
86 X^2
87 ST* 03
88 RCL 11
89 +
90 STO Y
91 RCL 08
92 ST+ 11
93 ST- Z
94 2
95 /
96 STO 09
97 -
98 RCL 01
99 *
100 -
101 RCL 05
102 X^2
103 1.5
104 STO 10
105 ST/ Y

106 FRC
107 -
108 RCL 07
109 /
110 RCL 02
111 *
112 +
113 RCL 02
114 *
115 RCL 06
116 X^2
117 RCL 07
118 ST- 03
119 ST+ X
120 /
121 RCL 03
122 4
123 ST/ 11
124 *
125 RCL 06
126 +
127 RCL 07
128 +
129 STO 05
130 -
131 RCL 01
132 *
133 RCL 05
134 +
135 RCL 03
136 ST+ X
137 STO 05
138 +
139 RCL 01
140 *

141 +
142 RCL 05
143 -
144 RCL 02
145 *
146 RCL 03
147 RCL 06
148 *
149 STO 03
150 3
151 *
152 RCL 00
153 X^2
154 ST+ X
155 RCL 10
156 /
157 STO 10
158 +
159 RCL 11
160 +
161 STO 05
162 RCL 01
163 *
164 RCL 05
165 ST+ X
166 -
167 RCL 03
168 ST+ 10
169 ST+ 10
170 +
171 RCL 01
172 *
173 RCL 10
174 +
175 RCL 01
176 *

177 -
178 RCL 04
179 /
180 RCL 11
181 RCL 06
182 RCL 07
183 *
184 STO Z
185 -
186 RCL 01
187 *
188 +
189 RCL 02
190 ST* 09
191 RCL 07
192 ST- 09
193 *
194 -
195 RCL 01
196 *
197 +
198 RCL 02
199 RCL 09
200 *
201 -
202 RCL 04
203 /
204 RCL 02
205 ST* 06
206 RCL 07
207 /
208 STO 09
209 RCL 01
210 +
211 STO 10
212 -

213 RCL 08
214 3
215 -
216 4
217 /
218 RCL 01
219 *
220 RCL 09
221 -
222 RCL 01
223 ST* Y
224 +
225 RCL 06
226 2
227 /
228 RCL 09
229 ST* Y
230 -
231 -
232 RCL 04
233 /
234 RCL 10
235 -
236 X<>Y
237 RCL 04
238 ST/ Z
239 /
240 RCL 00
241 ST* Z
242 ST* Y
243 ST+ Z
244 +
245 6378.137
246 ST* Z
247 *
248 END

( 310 bytes / SIZE 012 )

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

Where D = geodesic distance & d = great elliptic distance

Example: Calculate the geodesic distance between the Sydney Observatory ( Long = 151°12'17"8 E , Lat = 33°51'41"1 S )
and the Palomar Observatory ( Long = 116°51'50"4 W , Lat = 33°21'22"4 N )

     151.12178 ENTER^
    -33.51411   ENTER^
   -116.51504 ENTER^
     33.21224   XEQ "TGD+"   >>>> D = 12138.68432 km                                                      ---Execution time = 14s---
                                 X<>Y    d = 12138.68551 km

Notes:

-With the great elliptic distance, the maximum error = 3 cm
-Of course, the precision would be better with the 3rd order formula ( cf "Great Elliptic Arc Length for the HP41" )

Conclusion:

-For a WGS84 ellipsoid, the more accurate results are obtained by Vincenty's formula ( cf "Terrestrial Geodesic Distance(1) for the HP41" paragraph 2 )
-Unfortunately, the iterations take a long time with an HP41.
-So "TGD" listed in §3-a) is often enough accurate.

-Likewise, with a triaxial ellipsoid model of the Earth, the program listed in §4 above is the most accurate, but very slow.
-So, "TGD" listed in §3-b) is satisfactory in many cases... it's up to you to judge !

Remark:

-The values of a & b that are used in all these programs are those given by WGS84.
-IERS2003 suggests other numbers: mean equatorial radius = 6378.13661 km & flattening = 1/298.256421 so, polar radius = 6356.751868 km

-For the triaxial ellipsoid, recent evaluations suggest: equatorial flattening = 1/91026
-So it yields::

a = 6378.171645 km
b = 6378.101575 km
c = 6356.751868 km

References:

[1] Paul D. Thomas - "Spheroidal Geodesics, Reference Systems & Local Geometry" - US Naval Oceanographic Office.
[2] Richard H. Rapp - "Geometric Geodesy, Part II - The Ohio State University
[3] Emanuel M. Sodano & Thelma A. Robinson - "Direct and Inverse Solutions of Geodesics" - Army Map Service, Technical Report n°7
[4] Paul D. Thomas - "MATHEMATICAL MODELS FOR NAVIGATION SYSTEMS" -

-Many thanks to Professor Richard H. Rapp who sent me reference [3] !

hp41programs

Terrestrial Geodesic Distance(2) for the HP-41