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 Distances ( Ellipsoid of Revolution )
8°) Geodesic + Great Elliptic Distances ( Triaxial Ellipsoid + Ellipsoid of Revolution )
9°) Best Program ?

Latest Update: §8°) & §9°)

-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 R-P
85 X<>Y
86 STO 08
87 2
88 ST/ 01
89 ST/ 11
90 /
91 STO 05
92 XEQ 02
93 STO 06
94 RCL 05
95 90
96 +
97 XEQ 02
98 STO 07
99 RCL 02
100 D-R
101 RCL 08
102 COS
103 STO 05
104 X^2
105 ST+ X
106 RCL 04
107 RCL 15
108 ST* 05
109 *
110 ST* Y
111 -
112 -
113 4
114 /
115 RCL 06
116 ST/ 15
117 RCL 07
118 ST- Y
119 ST+ X
120 /

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

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

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

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

( 357 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 98 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 R-P
85 X<>Y
86 STO 08
87 2
88 ST/ 01
89 ST/ 11
90 /
91 STO 05
92 XEQ 02
93 STO 06
94 RCL 05
95 90
96 +
97 XEQ 02
98 STO 07
99 RCL 02
100 D-R
101 RCL 08
102 COS
103 STO 05
104 X^2
105 ST+ X
106 RCL 04

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

160 X^2
161 *
162 STO 09
163 X^2
164 RCL 10
165 X^2
166 +
167 +
168 SQRT
169 ST/ 09
170 ST/ 10
171 /
172 STO 08
173 RTN
174 LBL 02
175 RCL 02
176 2
177 /
178 +
179 ENTER
180 SIN
181 ENTER
182 X^2
183 RCL 13
184 *
185 RCL 02
186 R^
187 -
188 SIN
189 ST* Z
190 X^2
191 RCL 00
192 *
193 +
194 X<>Y
195 RCL 09
196 *
197 +
198 SQRT
199 RTN
200 LBL 03
201 RCL 14
202 1
203 STO 08
204 STO 09
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 09
220 P-R
221 RCL 05
222 /
223 R-P
224 RDN
225 +
226 2
227 ST/ 04
228 /
229 ST- Y
230 RCL 09
231 P-R
232 X^2
233 STO T
234 RDN
235 X^2
236 STO 14
237 SIGN
238 P-R
239 X^2
240 ST* 14
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 14
254 RCL 09
255 LASTX
256 STO 14
257 -
258 /
259 ST- 02
260 +
261 ST- 09
262 X<> 14
263 SQRT
264 ASIN
265 ST+ X

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

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

( 465 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 98 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 µ )

-We may simplify this formula like this:

(d - Δ') / a = µ ( 1 - A ) { (f²/4) [ A ( µ Cot µ - 1 ) + B ( µ / Sin µ - (Sin µ) / µ ) ]
+ (f³/8) [ B² ( µ² / Sin² µ + 2 ( Sin µ Cos µ ) / µ - 3 )
+ B ( ( 4 - 8 A ) ( Sin µ ) / µ - 3 A Cos µ + ( 4 A - 2 ) µ² Cos µ / Sin² µ + ( 7 A - 2 ) µ / Sin µ )
+ A ( -A ( Sin µ Cos µ ) / µ + ( 7 A - 2 ) µ Cot µ - 8 ( 1 - A ) ( µ² / 6 ) + ( 3 A - 2 ) µ² Cot² µ + 4 - 9 A ) ] }

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 RCL 08
29 STO 05
30 RCL 09
31 STO 06
32 RCL 10
33 STO 07
34 RCL 13
35 ST- 11
36 RCL 14
37 XEQ 01
38 XEQ 00
39 STO 04
40 RCL 01
41 ST/ 05
42 ST/ 08
43 RCL 02
44 ST+ 01
45 ST/ 06
46 ST/ 09
47 RCL 03
48 ST/ 07
49 ST/ 10

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

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

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

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

246 P-R
247 STO 09
248 STO 12
249 STO 13
250 RDN
251 /
252 STO 06
253 ST/ 12
254 X^2
255 ST* 09
256 3
257 ST* 04
258 -
259 RCL 12
260 ST+ X
261 +
262 RCL 05
263 *
264 RCL 14
265 RCL 07
266 -
267 ST+ X
268 STO 10
269 ST+ X
270 RCL 06
271 /
272 -
273 RCL 09
274 RCL 10
275 *
276 +
277 RCL 04
278 ST+ 10
279 RCL 13
280 ST* 09
281 *
282 -
283 RCL 06
284 ST* 13
285 RCL 10
286 *
287 +
288 RCL 05
289 *
290 RCL 10
291 RCL 13
292 *
293 RCL 12
294 RCL 14
295 *

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

( 439 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.65754 km                                                      ---Execution time = 25s---

LASTX = R15 = 12138.65874 km = great elliptic arc distance

Note:

-When line 95 is executed, the semi axes of the great ellipse = R15/X & R15/Y.
-With the example above, 6378.151658 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 R-P
85 X<>Y
86 STO 08
87 2
88 ST/ 01
89 ST/ 11
90 /
91 STO 05
92 XEQ 02
93 STO 06
94 RCL 05
95 90
96 +

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

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

193 RCL 00
194 *
195 +
196 X<>Y
197 RCL 09
198 *
199 +
200 SQRT
201 RTN
202 LBL 03
203 RCL 12
204 +
205 2
206 /
207 ST- Y
208 1
209 STO 07
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 RCL 07
234 LASTX
235 STO 08
236 -
237 /
238 ST- 02
239 +
240 ST- 07

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

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

( 426 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 98 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 R-P
85 X<>Y
86 STO 08
87 2
88 ST/ 01
89 ST/ 11
90 /
91 STO 05
92 XEQ 02
93 STO 06
94 RCL 05
95 90
96 +
97 XEQ 02
98 STO 07
99 RCL 02
100 D-R

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

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

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

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

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

( 443 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 98 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 ENTER
138 SIGN
139 ST* Z
140 *
141 R-P
142 X<>Y
143 2
144 /
145 90
146 STO 00
147 X<>Y
148 X<0?
149 +
150 STO 05
151 XEQ 03
152 STO 06
153 RCL 05
154 RCL 00
155 +
156 XEQ 03
157 STO 07
158 RCL 06
159 X>Y?
160 GTO 00
161 X<> 07
162 STO 06
163 RCL 00
164 ST+ 05
165 LBL 00
166 PI
167 R-D
168 RCL 05
169 X=Y?
170 CLX
171 STO 05
172 CHS
173 1
174 P-R
175 RCL 07
176 RCL 06
177 /
178 X^2
179 STO 09
180 *
181 R-P
182 CLX
183 RCL 04
184 RCL 05
185 -
186 1
187 P-R
188 RCL 09
189 *
190 R-P
191 CLX
192 ENTER
193 LBL 13
194 -
195 360

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

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

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

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

( 568 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.
-With many 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) / a

Δ ~ 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 Distances ( Ellipsoid of Revolution )

-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" )

8°) Geodesic + Great Elliptic Distances ( Triaxial Ellipsoid + Ellipsoid of Revolution )

-This program computes the geodesic distance & the great elliptic arc distance of the triaxial ellipsoid ( semi-axes: a , b , c )
+ the geodesic distance & the great elliptic arc distance of the corresponding ellipsoid of revolution ( semi-axes (a+b)/2 , c )

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

-Line 47 is a three-byte GTO 04

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

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

131 *
132 -
133 ENTER
134 SIGN
135 ST* Z
136 *
137 R-P
138 CLX
139 2
140 /
141 STO 05
142 XEQ 02
143 STO 06
144 RCL 05
145 90
146 +
147 XEQ 02
148 SQRT
149 RCL 15
150 X<>Y
151 /
152 STO 07
153 RCL 15
154 RCL 06
155 SQRT
156 /
157 STO 06
158 RCL 05
159 RCL 04
160 ST+ 04
161 -
162 ST+ 04
163 1
164 P-R
165 RCL 07
166 RCL 06
167 /
168 X^2
169 STO 09
170 *
171 R-P
172 X<>Y
173 STO 14
174 RCL 04
175 1
176 P-R
177 RCL 09
178 *
179 R-P
180 CLX
181 LBL 03
182 2
183 /
184 SIN
185 X^2
186 STO 00
187 CLX
188 RCL 06
189 RCL 07
190 -
191 RCL 06
192 /
193 STO 07
194 CLX
195 SIGN

196 STO 11
197 P-R
198 LASTX
199 RCL 07
200 -
201 STO 08
202 /
203 R-P
204 X<> 14
205 RCL 11
206 P-R
207 RCL 08
208 /
209 R-P
210 RDN
211 +
212 2
213 ST/ 07
214 /
215 ST- Y
216 RCL 11
217 P-R
218 X^2
219 STO 05
220 X<>Y
221 X^2
222 STO 04
223 RDN
224 X<>Y
225 RCL 11
226 P-R
227 X^2
228 ST* 04
229 RDN
230 X^2
231 ST* 05
232 STO Z
233 -
234 RCL 00
235 *
236 +
237 ST/ 05
238 RCL 05
239 RCL 04
240 RCL 11
241 R^
242 STO 04
243 -
244 /
245 ST- 05
246 +
247 X<> 04
248 SQRT
249 ASIN
250 ST+ X
251 D-R
252 STO 00
253 LASTX
254 RCL 11
255 P-R
256 STO 13
257 STO 09
258 RDN
259 STO 08
260 /

261 STO 10
262 ST/ Z
263 X^2
264 ST* 13
265 RCL 11
266 +
267 ENTER
268 R^
269 ST+ 11
270 ST- Z
271 2
272 /
273 STO 12
274 -
275 RCL 04
276 *
277 -
278 RCL 08
279 X^2
280 1.5
281 ST/ Y
282 FRC
283 -
284 RCL 10
285 /
286 RCL 05
287 *
288 +
289 RCL 05
290 *
291 RCL 09
292 X^2
293 RCL 10
294 ST- 13
295 ST+ X
296 /
297 RCL 13
298 4
299 ST/ 11
300 *
301 RCL 09
302 +
303 RCL 10
304 +
305 STO Z
306 -
307 RCL 04
308 *
309 +
310 RCL 13
311 ST+ X
312 STO 08
313 +
314 RCL 04
315 *
316 +
317 RCL 08
318 -
319 RCL 05
320 *
321 RCL 09
322 RCL 13
323 *
324 STO 13
325 RCL 00

326 ST+ X
327 X^2
328 3
329 ST* Z
330 /
331 STO Z
332 +
333 RCL 11
334 +
335 STO 08
336 RCL 04
337 *
338 RCL 08
339 ST+ X
340 -
341 RCL 13
342 ST+ Z
343 ST+ Z
344 +
345 RCL 04
346 *
347 +
348 RCL 04
349 *
350 -
351 RCL 07
352 *
353 RCL 11
354 RCL 09
355 RCL 10
356 *
357 STO Z
358 -
359 RCL 04
360 *
361 +
362 RCL 05
363 ST* 12
364 RCL 10
365 ST- 12
366 *
367 -
368 RCL 04
369 *
370 +
371 RCL 05
372 RCL 12
373 *
374 -
375 RCL 07
376 *
377 RCL 05
378 RCL 10
379 /
380 -
381 RCL 04
382 -
383 RCL 07
384 *
385 RCL 00
386 ST* Y
387 +
388 RCL 06
389 *
390 RTN

391 LBL 01
392 1
393 P-R
394 X<>Y
395 RCL 03
396 X^2
397 *
398 STO 10
399 RDN
400 P-R
401 RCL 01
402 X^2
403 *
404 STO Z
405 X^2
406 X<>Y
407 RCL 02
408 X^2
409 *
410 STO 09
411 X^2
412 RCL 10
413 X^2
414 +
415 +
416 SQRT
417 ST/ 09
418 ST/ 10
419 /
420 STO 08
421 RTN
422 LBL 02
423 STO 12
424 RCL 04
425 +
426 SIN
427 ENTER
428 X^2
429 RCL 13
430 *
431 RCL 04
432 RCL 12
433 -
434 SIN
435 ST* Z
436 X^2
437 RCL 14
438 *
439 +
440 X<>Y
441 RCL 09
442 *
443 +
444 X<=0?
445 RCL 06
446 RTN
447 LBL 04
448 ENTER
449 STO 00
450 RCL 21
451 -
452 RCL 20
453 +
454 END

( 567 bytes / SIZE 022 )

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

Where D₃ = geodesic distance ( triaxial ellipsoid ) & d₃ = great elliptic distance ( triaxial ellipsoid ) = R20
D₂ = geodesic distance ( biaxial ellipsoid ) = R00 & d₂ = great elliptic distance ( biaxial ellipsoid ) = R21

Example: 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.65754 km                                                       ---Execution time = 72s---
                                     X<>Y    D₂ = 12138.68432 km = R00

and LASTX = d₃ = 12138.65876 km = R20 & RCL 21 = d₂ = 12138.68554 km

Notes:

-Andoyer's formula of order 3 is used 3 times in this program:
-Not very fast but less slow than with Vincenty's formula.

9°) Best Program ?

-We can make a shorter & faster program if we don't want to calculate the geodesic distance on the similar ellipsoid of revolution.
-Andoyer's formula of order 3 becomes less complex to compute the great elliptic arc length.
-In the program below, (d - Δ') is calculated like in paragraph 2°)c) ( All these formulas are of order 3 )

Data Registers:   R00 to 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 RCL 05
38 -
39 X^2
40 RCL 09
41 RCL 06
42 -
43 X^2
44 RCL 10
45 RCL 07
46 -
47 X^2
48 +
49 +
50 SQRT
51 2
52 /
53 ASIN
54 STO 04
55 RCL 01
56 ST/ 05
57 ST/ 08
58 RCL 02
59 ST+ 01
60 ST/ 06
61 ST/ 09

62 RCL 03
63 ST/ 07
64 ST/ 10
65 RCL 05
66 RCL 08
67 *
68 RCL 06
69 RCL 09
70 *
71 RCL 07
72 RCL 10
73 *
74 +
75 +
76 ST+ X
77 X<> 08
78 X^2
79 RCL 09
80 X^2
81 RCL 10
82 X^2
83 +
84 +
85 STO 10
86 RCL 05
87 X^2
88 RCL 06
89 X^2
90 RCL 07
91 X^2
92 +
93 +
94 STO 00
95 -
96 RCL 04
97 ST+ X
98 STO 07
99 1
100 STO 13
101 P-R
102 STO 02
103 RDN
104 STO 15
105 *
106 RCL 08
107 RCL 00
108 RCL 10
109 +
110 RCL 02
111 *
112 -
113 ENTER
114 SIGN
115 ST* Z
116 *
117 R-P
118 CLX
119 2
120 /
121 STO 05
122 XEQ 02

123 STO 06
124 RCL 05
125 90
126 +
127 XEQ 02
128 SQRT
129 X<> 05
130 RCL 04
131 -
132 ST+ 07
133 RCL 13
134 P-R
135 RCL 06
136 SQRT
137 STO 06
138 RCL 05
139 /
140 STO 09
141 *
142 R-P
143 X<>Y
144 ENTER
145 X<> 07
146 RCL 13
147 P-R
148 RCL 09
149 *
150 R-P
151 RDN
152 ST+ 07
153 -
154 ABS
155 STO 09
156 RCL 13
157 P-R
158 X<> 09
159 D-R
160 STO 00
161 RDN
162 STO 04
163 X^2
164 .75
165 /
166 RCL 13
167 -
168 RCL 04
169 *
170 RCL 07
171 COS
172 STO 07
173 *
174 RCL 04
175 RCL 09
176 *
177 STO 10
178 -
179 RCL 07
180 *
181 RCL 09
182 RCL 10
183 *

184 +
185 RCL 07
186 *
187 RCL 00
188 RCL 10
189 +
190 2
191 /
192 STO Z
193 +
194 RCL 05
195 RCL 06
196 -
197 RCL 05
198 ST+ X
199 /
200 STO T
201 *
202 +
203 RCL 07
204 ST* 04
205 X^2
206 RCL 10
207 *
208 -
209 *
210 2
211 /
212 RCL 04
213 -
214 RCL 00
215 -
216 *
217 RCL 00
218 +
219 RCL 06
220 /
221 ST* 15
222 GTO 03
223 LBL 01
224 1
225 P-R
226 X<>Y
227 RCL 03
228 X^2
229 *
230 STO 10
231 RDN
232 P-R
233 RCL 01
234 X^2
235 *
236 STO Z
237 X^2
238 X<>Y
239 RCL 02
240 X^2
241 *
242 STO 09
243 X^2
244 RCL 10

245 X^2
246 +
247 +
248 SQRT
249 ST/ 09
250 ST/ 10
251 /
252 STO 08
253 RTN
254 LBL 02
255 STO 09
256 RCL 04
257 +
258 SIN
259 ENTER
260 X^2
261 RCL 10
262 *
263 RCL 04
264 RCL 09
265 -
266 SIN
267 ST* Z
268 X^2
269 RCL 00
270 *
271 +
272 X<>Y
273 RCL 08
274 *
275 +
276 X<=0?
277 RCL 06
278 RTN
279 LBL 03
280 RCL 12
281 RCL 14
282 +
283 2
284 ST/ 01
285 ST/ 11
286 /
287 ST- 14
288 RCL 13
289 P-R
290 X^2
291 STO 05
292 X<>Y
293 X^2
294 X<> 14
295 RCL 13
296 P-R
297 X^2
298 ST* 14
299 RDN
300 X^2
301 ST* 05
302 STO Z
303 -
304 RCL 11
305 SIN

306 X^2
307 *
308 +
309 ST/ 05
310 RCL 05
311 RCL 14
312 RCL 13
313 R^
314 STO 14
315 -
316 /
317 ST- 05
318 +
319 STO 04
320 ST- 13
321 X<> 14
322 SQRT
323 ASIN
324 ST+ X
325 D-R
326 STO 00
327 LASTX
328 1
329 P-R
330 STO 07
331 STO 09
332 STO 12
333 RDN
334 /
335 STO 06
336 ST/ 12
337 X^2
338 ST* 09
339 3
340 ST* 04
341 -
342 RCL 12
343 ST+ X
344 +
345 RCL 05
346 *
347 RCL 14
348 RCL 13
349 -
350 ST+ X
351 STO 10
352 ST+ X
353 RCL 06
354 /
355 -
356 RCL 09
357 RCL 10
358 *
359 +
360 RCL 04
361 ST+ 10
362 RCL 07
363 ST* 09
364 *
365 -
366 RCL 06

367 ST* 07
368 RCL 10
369 *
370 +
371 RCL 05
372 *
373 RCL 07
374 RCL 10
375 *
376 RCL 12
377 RCL 14
378 *
379 -
380 RCL 00
381 ST+ X
382 X^2
383 3
384 /
385 RCL 13
386 ST* 00
387 ST- 10
388 ST- 10
389 *
390 -
391 RCL 04
392 2
393 -
394 RCL 09
395 *
396 +
397 RCL 10
398 -
399 RCL 14
400 ST* 07
401 ST- 07
402 *
403 +
404 RCL 01
405 RCL 03
406 -
407 RCL 01
408 ST* 00
409 ST+ X
410 /
411 STO Z
412 *
413 RCL 06
414 ENTER
415 1/X
416 -
417 RCL 05
418 *
419 +
420 RCL 07
421 +
422 *
423 *
424 RCL 00
425 *
426 RCL 15
427 +
428 END

( 529 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 = 30s---

   LASTX = R15 = 12138.65875 km = great elliptic arc distance

Note:

-This program also works well if D is small.

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