Template

Terrestrial Geodesic Distance(4) for the HP-41

Overview

1°) Geodetic Longitudes & Latiudes

a) Program#1
b) Program#2
c) Program#3
d) Program#4

2°) Geocentric Longitudes & Geodetic Latitudes

a) Program#1
b) Program#2

3°) Geocentric Longitudes & Parametric Latitudes

4°) Geocentric Longitudes & Latitudes

a) Program(s)#1
b) Program(s)#2
c) Program#3
d) Program#4
e) Program#5 ( with great elliptic arc length )
f) My favorite Programs

f1) Program#1
f2) Program(s)#2

5°) Transformation of Coordinates

a) Geocentric Coordinates <> Geodetic Coordinates
b) WGS84 coordinates -> Triaxial Ellipsoid Geocentric Coordinates

Latest Update: paragraph 4°)a)

-These programs employ a little generalization of Andoyer's formulae to a triaxial ellipsoid model of the Earth.

-The semi-axes are:

a = 6378.172 km
b = 6378.102 km           the longitudes are measured from the equatorial major-axis = 14°92911 West
c = 6356.752314 km

-I've recently found relatively simple formulae to calculate the geodesic distance on a triaxial ellipsoid model of the Earth.

-With geocentric coordinates, the relative error is smaller than with geodetic coordinates.
-So my favorite programs are in paragraph 4°)f)

-If we compare the precision of this formula with the non-iterative formulae for an ellipsoid of revolution,
it is more accurate than Andoyer's formula of order 2 but less accurate than Sodano's formula of order 3:
-So, I think it is interesting but it is not perfect.

-Of course, if the points are ~ antipodal and if the difference between the longitudes is almost 180°, errors may be large.

-However, if we use the program listed in §4°e) as a subroutine of the program listed in "Terrestrial Geodesic Distance(2)" (§5) to minimize the sum of 2 distances,
we will have a good accuracy, even if the 2 points are nearly antipodal because the error are very small if the difference between the longitudes is <= 90°
-But it is very slow with an HP41: free42 or DM42 are much faster.

1°) Geodedic Longitudes & Latitudes

a) Program#1

-This is one of the shortest programs... but not the best formula !

   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 ( 1 - A ) (8/6) µ³ + ( B + A Cos µ ) ( µ³ / Sin² µ ) ( B + ( 3 A - 2 ) Cos µ ) ] ( 1 - A )

+   (f'/2) [ -A' µ - 3.B' Sin µ ]
+ ( f f'/4 ) ( -6 AB' - 5 A'B - AA'/4 - BB'/2 ) µ } f/2 ~ 1/595.543 & f'/2 ~ 1/182233 & a = 6378.172

-The angle µ is computed with:

Sin² µ = Sin² ( λ - λ' ) Cos² β' + [ Cos β Sin β' - Cos β' Sin β Cos ( λ - λ' ) ]²
   Cos µ = Sin β Sin β' + Cos β Cos β' Cos ( λ - λ' )

-And A & B and A' & B' are calculated with:

A = [ Sin² β + Sin² β' - 2 Sin β Sin β' Cos µ ] / Sin² µ
B = [ ( Sin² β + Sin² β' ) Cos µ - 2 Sin β Sin β' ] / Sin² µ

   A' = ( cos² β Sin² λ + cos² β' Sin² λ' - 2 cos β Sin λ cos β' Sin λ' Cos µ ) / Sin² µ
   B' = [ ( cos² β Sin² λ + cos² β' Sin² λ' ) Cos µ - 2 cos β Sin λ cos β' Sin λ' ] / Sin² µ

Data Registers: R00 thru R09: temp R00 = µ R01 = -A/4 R02 = -B R08 = -A' R09 = -B'
Flags: /
Subroutines: /

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

28 X<> 07
29 ST- 05
30 SIN
31 *
32 ST* 08
33 X^2
34 +
35 STO 03
36 X<>Y
37 STO 01
38 ST+ 01
39 X^2
40 RCL 07
41 ST* 01
42 ST* 02
43 X^2
44 +
45 X<> 05
46 RCL 06
47 P-R
48 ST* 00
49 R^
50 ST* 07
51 *
52 RCL 02
53 -
54 X^2

55 X<>Y
56 X^2
57 +
58 ST/ 01
59 ST/ 03
60 ST/ 05
61 ST/ 08
62 SQRT
63 RCL 08
64 RCL 00
65 RCL 07
66 +
67 STO 09
68 ST* 08
69 RCL 03
70 ST- 08
71 *
72 -
73 X<> 09
74 STO 03
75 R-P
76 X<>Y
77 D-R
78 STO 00
79 R^
80 /
81 RCL 01

82 RCL 03
83 ST* 01
84 RCL 05
85 ST- 01
86 *
87 -
88 STO 02
89 RCL 01
90 ST+ 04
91 RCL 03
92 *
93 +
94 STO 07
95 *
96 STO 05
97 *
98 RCL 03
99 RCL 04
100 ST+ X
101 *
102 RCL 07
103 +
104 *
105 RCL 00
106 ST+ X
107 X^2
108 RCL 01

109 *
110 RCL 04
111 *
112 3
113 /
114 +
115 595.543
116 STO 07
117 /
118 RCL 05
119 3
120 R^
121 ST/ 03
122 /
123 STO 06
124 RCL 02
125 *
126 STO 05
127 ST+ X
128 +
129 -
130 RCL 04
131 *
132 RCL 01
133 ST+ 05
134 ENTER
135 X^2

136 STO 04
137 RCL 02
138 X^2
139 ST+ X
140 -
141 RCL 03
142 *
143 RCL 04
144 +
145 ST+ Y
146 4
147 ST* Z
148 ST/ 01
149 /
150 -
151 -
152 RCL 07
153 /
154 RCL 05
155 +
156 RCL 08
157 5
158 *
159 RCL 09
160 2
161 /
162 +

163 RCL 02
164 *
165 RCL 09
166 ST* 06
167 24
168 *
169 RCL 08
170 ST+ 06
171 +
172 RCL 01
173 *
174 +
175 RCL 07
176 ST/ Z
177 /
178 RCL 06
179 -
180 182233
181 /
182 -
183 RCL 00
184 ST* Y
185 +
186 6378.172
187 *
188 END

( 248 bytes / SIZE 010 )

STACK	INPUTS	OUTPUTS
T	λ ( ° ' " )	/
Z	β ( ° ' " )	/
Y	λ' ( ° ' " )	/
X	β' ( ° ' " )	D ( km )

Where D is the geodesic distance on the triaxial ellipsoid

Example1: λ = +166°08'02"2 , β = -33°51'41"1
λ' = -101°56'06"0 , β' = +33°21'22"4

     166.08022 ENTER^
    -33.51411   ENTER^
   -101.5606    ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65759 km                                                      ---Execution time = 11.3s---

-Exact result = 12138.657552.....

Example2:    λ = -24°    , β = -24°
                      λ' = +144° ,  β' = +13°

24 ENTER^ ENTER^
144 ENTER^
13 XEQ "TGD" >>>> D = 18263.94639 km

-Exact D = 18263.946445....
-In this example, the great elliptic arc length = 18264.171677...

Example3:    λ = 0°    , β = 5°
                      λ' = +175° ,  β' = 3°

0 ENTER^
5 ENTER^
175 ENTER^
3 XEQ "TGD" >>>> D = 18968.18007 km

-Exact D = 18968.180062....

Notes:

-If | λ - λ' | = 160° , the maximum relative error is about 3.78 10^-8 and on average, the relative error is about 0.91 10^-8
-If | λ - λ' | = 120° , the maximum relative error is also ~ 3.78 10^-8 and on average, the relative error is about 0.73 10^-8

-Here are a few maximum errors ( in centimeters ):

| λ - λ' | = 120°

( +30°, +88° ) ( +150° , -88° ) M = +28.7 cm

| λ - λ' | = 160°

( -65°, +43° ) ( +95° , -27° ) M = +37.9 cm
( -22°, +50° ) ( +138° , -55° ) m = -37.8 cm

| λ - λ' | = 164°

( -62°, +48° ) ( +102°, -39° ) M = +51.2 cm
( +8° , +4° ) ( +172° , +4° ) m = -61.1 cm

| λ - λ' | = 168°

( -84°, +46° ) ( +84° , -38° ) M = +77.7 cm
( +6° , +3° ) ( +174° , +4° ) m = -105 cm

-If A is close to 1 ( when the geodesics are almost perpendicular to the equator ),
  the precision is often better without the term with f.f'/4 (!)

b) Program#2

-We can find a formula with   ( f f'/4 ) µ ( 1-A ) to get a better precision if A is close to 1 ( when the geodesics are almost perpendicular to the equator ):

   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 ( 1 - A ) µ³ + ( B + A Cos µ ) ( µ³ / Sin² µ ) ( B + ( 3 A - 2 ) Cos µ ) ] ( 1 - A )

+   (f'/2) [ -A' µ - 3.B' Sin µ ]
+ ( f f'/4 ) µ² ( 1-A ) ( -6 AB' - 4 A'B - 2 AA' ) } f/2 ~ 1/595.543 & f'/2 ~ 1/182233 & a = 6378.172

-The angle µ is computed with:

Sin² µ = Sin² ( λ - λ' ) Cos² β' + [ Cos β Sin β' - Cos β' Sin β Cos ( λ - λ' ) ]²
   Cos µ = Sin β Sin β' + Cos β Cos β' Cos ( λ - λ' )

-And A & B and A' & B' are calculated with:

A = [ Sin² β + Sin² β' - 2 Sin β Sin β' Cos µ ] / Sin² µ
B = [ ( Sin² β + Sin² β' ) Cos µ - 2 Sin β Sin β' ] / Sin² µ

   A' = ( cos² β Sin² λ + cos² β' Sin² λ' - 2 cos β Sin λ cos β' Sin λ' Cos µ ) / Sin² µ
   B' = [ ( cos² β Sin² λ + cos² β' Sin² λ' ) Cos µ - 2 cos β Sin λ cos β' Sin λ' ] / Sin² µ

Data Registers: R00 thru R10: temp R00 = µ R01 = -A R02 = A' B R08 = -A' R09 = -B'
Flags: /
Subroutines: /

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

27 X<>Y
28 X<> 07
29 ST- 05
30 SIN
31 *
32 ST* 08
33 X^2
34 +
35 STO 03
36 X<>Y
37 STO 01
38 ST+ 01
39 X^2
40 RCL 07
41 ST* 01
42 ST* 02
43 X^2
44 +
45 X<> 05
46 RCL 06
47 P-R
48 ST* 00
49 R^
50 ST* 07
51 *
52 RCL 02

53 -
54 X^2
55 X<>Y
56 X^2
57 +
58 ST/ 01
59 ST/ 03
60 ST/ 05
61 ST/ 08
62 SQRT
63 RCL 08
64 RCL 00
65 RCL 07
66 +
67 STO 09
68 ST* 08
69 RCL 03
70 ST- 08
71 *
72 -
73 X<> 09
74 STO 03
75 R-P
76 X<>Y
77 D-R
78 STO 00

79 STO 10
80 R^
81 /
82 RCL 01
83 RCL 03
84 ST* 01
85 RCL 05
86 ST- 01
87 *
88 -
89 STO 02
90 RCL 01
91 ST+ 04
92 RCL 03
93 *
94 +
95 STO 07
96 *
97 STO 05
98 *
99 RCL 03
100 RCL 04
101 ST* 10
102 ST+ X
103 *
104 RCL 07

105 +
106 *
107 RCL 00
108 RCL 01
109 *
110 RCL 10
111 *
112 +
113 595.543
114 STO 07
115 /
116 RCL 05
117 3
118 R^
119 ST/ 03
120 /
121 STO 06
122 RCL 02
123 *
124 STO 05
125 ST+ X
126 +
127 -
128 RCL 04
129 *
130 RCL 01

131 ST+ 05
132 ENTER
133 X^2
134 STO 04
135 RCL 02
136 X^2
137 ST+ X
138 -
139 RCL 03
140 *
141 RCL 04
142 +
143 ST+ Y
144 4
145 ST* Z
146 /
147 -
148 -
149 RCL 07
150 /
151 RCL 05
152 +
153 RCL 09
154 ST* 06
155 3
156 *

157 RCL 08
158 ST+ 06
159 ST* 02
160 +
161 RCL 01
162 *
163 RCL 02
164 ST+ X
165 +
166 ST+ X
167 RCL 10
168 *
169 RCL 07
170 ST/ Z
171 /
172 RCL 06
173 -
174 182233
175 /
176 -
177 RCL 00
178 ST* Y
179 +
180 6378.172
181 *
182 END

( 243 bytes / SIZE 011 )

STACK	INPUTS	OUTPUTS
T	λ ( ° ' " )	/
Z	β ( ° ' " )	/
Y	λ' ( ° ' " )	/
X	β' ( ° ' " )	D ( km )

Where D is the geodesic distance on the triaxial ellipsoid

Example1: λ = +166°08'02"2 , β = -33°51'41"1
λ' = -101°56'06"0 , β' = +33°21'22"4

     166.08022 ENTER^
    -33.51411   ENTER^
   -101.5606    ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65757 km                                                      ---Execution time = 10.9s---

-Exact result = 12138.657552.....

Example2:    λ = -24°    , β = -24°
                      λ' = +144° ,  β' = +13°

24 ENTER^ ENTER^
144 ENTER^
13 XEQ "TGD" >>>> D = 18263.94632 km

-Exact D = 18263.946445....
-In this example, the great elliptic arc length = 18264.171677...

Example3:    λ = 0°    , β = 5°
                      λ' = +175° ,  β' = 3°

0 ENTER^
5 ENTER^
175 ENTER^
3 XEQ "TGD" >>>> D = 18968.17996 km

-Exact D = 18968.180062....

Notes:

-If | λ - λ' | = 160° , the maximum relative error is about 3.8 10^-8 and on average, the relative error is about 10^-8

-Here are a few maximum errors ( in centimeters ):

| λ - λ' | = 120°

( +30° , +53° ) ( +150° , -53° ) m = -23.6 cm
( +30°, +19° ) ( +150° , +19° ) M = +22.3 cm

| λ - λ' | = 160°

( -80°, +34° ) ( +80° , -34° ) m = -46.2 cm
( -80°, +54° ) ( +80° , -63° ) M = +43.6 cm

| λ - λ' | = 168°

( -84°, +58° ) ( +84° , -52° ) M = +84.6 cm
( +6° , +4° ) ( +174° , +3° ) m = -79.1 cm

-The relative precision is a little smaller
-The absolute precision is better if | λ - λ' | = 120° & | λ - λ' | = 168° but not if | λ - λ' | = 160°

c) Program#3

-We can also use the mean radius of the equator ( 6378.137 instead of 6378.172 )

-Formula:

   Dist ~ a' µ { 1 + (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 ( 1 - A ) (8/6) µ + ( B + A Cos µ ) ( µ² / Sin² µ ) ( B + ( 3 A - 2 ) Cos µ ) ) ( 1 - A )
+ B² (3A-3) + B ( - 3A(1-A) Cos µ + (1-A) (7A-2) µ/Sin µ ) + A(1-A) (7A-2) µ Cot µ
+ ( 1-A ) ( A ( 7A-2 )sqrt(7) + ( 11B² - 8A² ) (Sin µ) (Cos µ)/µ ) ]

+   (f'/2) [ -A" - 3.B" (Sin µ)/µ ]
+ ( f f'/4 ) ( -3 AB" - 3 A"B - AA"/2 ) µ } f/2 ~ 1/596.51444 & f'/2 ~ 1/364465 & a' = 6378.137 km

       Sin² µ = Sin² ( λ - λ' ) Cos² β' + [ Cos β Sin β' - Cos β' Sin β Cos ( λ - λ' ) ]²
      Cos µ = Sin β Sin β' + Cos β Cos β' Cos ( λ - λ' )

A" = [ cos² β ( Sin² λ - Cos² λ ) + cos² β' ( Sin² λ' - Cos² λ' ) + 2 cos β cos β' Cos ( λ + λ' ) Cos µ ] / Sin² µ
B" = [ ( cos² β ( Sin² λ - Cos² λ ) + cos² β' ( Sin² λ' - Cos² λ' ) ) Cos µ + 2 cos β cos β' Cos ( λ + λ' ) ] / Sin² µ

A = [ Sin² β + Sin² β' - 2 Sin β Sin β' Cos µ ] / Sin² µ
B = [ ( Sin² β + Sin² β' ) Cos µ - 2 Sin β Sin β' ] / Sin² µ

-In fact, we can express -A" = 2 ( 1 - A' ) - A and -B" = -2 B' - B with:

   A' = ( cos² β Sin² λ + cos² β' Sin² λ' - 2 cos β Sin λ cos β' Sin λ' Cos µ ) / Sin² µ
   B' = [ ( cos² β Sin² λ + cos² β' Sin² λ' ) Cos µ - 2 cos β Sin λ cos β' Sin λ' ] / Sin² µ

-We save 0.5 second in execution time ( same numbers of bytes )

Data Registers: R00 thru R12: temp R01 = -A R02 = -B R08 = -A" R09 = -3.B"
Flags: /
Subroutines: /

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

63 RCL 08
64 RCL 00
65 RCL 07
66 +
67 STO 09
68 ST* 08
69 RCL 03
70 ST- 08
71 *
72 -
73 ST+ X
74 X<> 09
75 STO 03
76 R-P
77 RCL 08
78 +
79 ST+ X
80 STO 08
81 X<>Y
82 D-R
83 STO 00
84 R^
85 /
86 STO 06
87 RCL 01
88 RCL 03
89 ST* 01
90 RCL 05
91 ST- 01
92 *
93 -

94 STO 02
95 ST+ 09
96 RCL 01
97 ST+ 04
98 ST+ 08
99 RCL 03
100 *
101 +
102 STO 05
103 LASTX
104 RCL 03
105 +
106 ST+ X
107 +
108 *
109 RCL 01
110 7
111 *
112 2
113 +
114 STO 07
115 +
116 *
117 RCL 02
118 3
119 ST* 09
120 *
121 STO 10
122 STO 12
123 -
124 RCL 05

125 *
126 RCL 07
127 7
128 SQRT
129 *
130 RCL 00
131 ST+ X
132 X^2
133 RCL 04
134 *
135 3
136 /
137 +
138 RCL 01
139 *
140 +
141 RCL 10
142 RCL 02
143 ST* Y
144 X^2
145 STO 07
146 RCL 01
147 X^2
148 STO 11
149 -
150 ST+ 07
151 8
152 *
153 +
154 RCL 03
155 RCL 06

156 ST* 05
157 ST/ 10
158 /
159 ST* 07
160 *
161 +
162 596.51444
163 STO 03
164 /
165 RCL 05
166 -
167 RCL 10
168 ST+ X
169 -
170 RCL 04
171 *
172 RCL 01
173 ST+ 10
174 RCL 07
175 RCL 11
176 -
177 ST- Y
178 4
179 ST* Z
180 /
181 +
182 -
183 RCL 03
184 /
185 RCL 10
186 +

187 RCL 09
188 RCL 08
189 ST* 12
190 2
191 /
192 +
193 RCL 01
194 *
195 RCL 12
196 +
197 RCL 00
198 *
199 RCL 03
200 ST/ Z
201 /
202 RCL 08
203 -
204 RCL 09
205 RCL 06
206 /
207 -
208 364465
209 /
210 -
211 RCL 00
212 ST* Y
213 +
214 6378.137
215 *
216 END

( 283 bytes / SIZE 013 )

STACK	INPUTS	OUTPUTS
T	λ ( ° ' " )	/
Z	β ( ° ' " )	/
Y	λ' ( ° ' " )	/
X	β' ( ° ' " )	D ( km )

Where D is the geodesic distance on the triaxial ellipsoid

Example: λ = +166°08'02"2 , β = -33°51'41"1
λ' = -101°56'06"0 , β' = +33°21'22"4

     166.08022 ENTER^
    -33.51411   ENTER^
   -101.5606    ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65757 km                                                      ---Execution time = 12.4 s---

-Exact D = 12138.657552....

Notes:

-The precision becomes:

| λ - λ' | = 120°

( +30°, +70° ) ( +150° , +70° ) m = -17.9 cm
( +30°, +28° ) ( +150° , +27° ) M = +15.2 cm

| λ - λ' | = 160°

( -80°, +67° ) ( +80° , -72° ) M = +23.6 cm
( +10°, +66° ) ( +170° , -71° ) m = -21.8 cm

| λ - λ' | = 168°

( -39°, +42° ) ( +129° , -38° ) M = +52.5 cm
( +6° , +65° ) ( +174° , -61° ) m = -50.9 cm

| λ - λ' | = 170°

( -85°, +48° ) ( +85° , -43° ) M = +70.8 cm
( +5° , +57° ) ( +175° , -61° ) m = -69.5 cm

-Though the absolute error is smaller, the relative error may be larger:
-If | λ - λ' | = 160° , the maximum relative error is about 7.4 10^-8 & on average, the relative error is ~ 0.94 10^-8

d) Program#4

-Here is a formula with   ( f f'/4 ) µ ( 1-A ) to get a better precision if A is close to 1 ( when the geodesics are almost perpendicular to the equator ):

-Formula:

   Dist ~ a' µ { 1 + (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 ( 1 - A ) (8/6) µ + ( B + A Cos µ ) ( µ² / Sin² µ ) ( B + ( 3 A - 2 ) Cos µ ) ) ( 1 - A )
+ B² (3A-3) + B ( - 3A(1-A) Cos µ + (1-A) (7A-2) µ/Sin µ ) + A(1-A) (7A-2) µ Cot µ
+ ( 1-A ) ( A ( 7A-2 )sqrt(7) + ( 11B² - 8A² ) (Sin µ) (Cos µ)/µ ) ]

+   (f'/2) [ -A" - 3.B" (Sin µ)/µ ]
+ ( f f'/4 ) ( -8 AB" - 5 A"B - 6 AA" ) µ ( 1-A ) } f/2 ~ 1/596.51444 & f'/2 ~ 1/364465 & a' = 6378.137 km

       Sin² µ = Sin² ( λ - λ' ) Cos² β' + [ Cos β Sin β' - Cos β' Sin β Cos ( λ - λ' ) ]²
      Cos µ = Sin β Sin β' + Cos β Cos β' Cos ( λ - λ' )

A = [ Sin² β + Sin² β' - 2 Sin β Sin β' Cos µ ] / Sin² µ
B = [ ( Sin² β + Sin² β' ) Cos µ - 2 Sin β Sin β' ] / Sin² µ

-A" = 2 ( 1 - A' ) - A and -B" = -2 B' - B with:

   A' = ( cos² β Sin² λ + cos² β' Sin² λ' - 2 cos β Sin λ cos β' Sin λ' Cos µ ) / Sin² µ
   B' = [ ( cos² β Sin² λ + cos² β' Sin² λ' ) Cos µ - 2 cos β Sin λ cos β' Sin λ' ] / Sin² µ

Data Registers: R00 thru R12: temp R01 = -A R02 = A".B R08 = -A" R09 = -B"
Flags: /
Subroutines: /

33 X^2
34 +
35 STO 03
36 X<>Y
37 STO 01
38 ST+ 01
39 X^2
40 RCL 07
41 ST* 01
42 ST* 02
43 X^2
44 +
45 X<> 05
46 RCL 06
47 P-R
48 ST* 00
49 R^
50 ST* 07
51 *
52 RCL 02
53 -
54 X^2
55 X<>Y
56 X^2
57 +
58 ST/ 01
59 ST/ 03
60 ST/ 05
61 ST/ 08
62 SQRT
63 RCL 08
64 RCL 00

65 RCL 07
66 +
67 STO 09
68 ST* 08
69 RCL 03
70 ST- 08
71 *
72 -
73 ST+ X
74 X<> 09
75 STO 03
76 R-P
77 RCL 08
78 +
79 ST+ X
80 STO 08
81 X<>Y
82 D-R
83 STO 00
84 R^
85 /
86 STO 06
87 RCL 01
88 RCL 03
89 ST* 01
90 RCL 05
91 ST- 01
92 *
93 -
94 STO 02
95 ST+ 09
96 RCL 01

97 ST+ 04
98 ST+ 08
99 RCL 03
100 *
101 +
102 STO 05
103 LASTX
104 RCL 03
105 +
106 ST+ X
107 +
108 *
109 RCL 01
110 7
111 *
112 2
113 +
114 STO 07
115 +
116 *
117 RCL 02
118 3
119 *
120 STO 10
121 -
122 RCL 05
123 *
124 RCL 07
125 7
126 SQRT
127 *
128 RCL 00

129 ST+ X
130 X^2
131 RCL 04
132 *
133 3
134 /
135 +
136 RCL 01
137 *
138 +
139 RCL 10
140 RCL 02
141 ST* Y
142 X^2
143 STO 07
144 RCL 01
145 X^2
146 STO 11
147 -
148 ST+ 07
149 8
150 *
151 +
152 RCL 03
153 RCL 06
154 ST* 05
155 ST/ 10
156 /
157 ST* 07
158 *
159 +
160 596.51444

161 STO 03
162 /
163 RCL 05
164 -
165 RCL 10
166 ST+ X
167 -
168 RCL 04
169 *
170 RCL 01
171 ST+ 10
172 RCL 07
173 RCL 11
174 -
175 ST- Y
176 4
177 ST* Z
178 /
179 +
180 -
181 RCL 03
182 /
183 RCL 10
184 +
185 RCL 08
186 ST* 02
187 6
188 *
189 RCL 09
190 8
191 *

192 +
193 RCL 01
194 *
195 RCL 02
196 5
197 *
198 +
199 RCL 00
200 *
201 RCL 04
202 *
203 RCL 03
204 ST/ Z
205 /
206 RCL 08
207 -
208 RCL 09
209 3
210 *
211 RCL 06
212 /
213 -
214 364465
215 /
216 -
217 RCL 00
218 ST* Y
219 +
220 6378.137
221 *
222 END

( 288 bytes / SIZE 012 )

STACK	INPUTS	OUTPUTS
T	λ ( ° ' " )	/
Z	β ( ° ' " )	/
Y	λ' ( ° ' " )	/
X	β' ( ° ' " )	D ( km )

Where D is the geodesic distance on the triaxial ellipsoid

Example: λ = +166°08'02"2 , β = -33°51'41"1
λ' = -101°56'06"0 , β' = +33°21'22"4

     166.08022 ENTER^
    -33.51411   ENTER^
   -101.5606    ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65755 km                                                      ---Execution time = 12.7 s---

-Exact D = 12138.657552....

Notes:

-Unfortunately, if | λ - λ' | = 160° , the maximum relative error remains about 7.2 10^-8
-And the absolute errors become:

| λ - λ' | = 160°

( -80°, +66° ) ( +80° , -59° ) M = +41.0 cm
( +10°, +66° ) ( +170° , -59° ) m = -40.4 cm

| λ - λ' | = 162°

( -82°, +63° ) ( +82° , -57° ) M = +52.9 cm
( +8° , +64° ) ( +172° , -58° ) m = -52.0 cm

2°) Geocentric Longitudes & Geodetic Latitudes

a) Program#1

-This program takes the geocentric longitudes and geodetic latitudes on WGS84 ( ellipsoid of revolution ) ( cf JPL )
and computes the geodesic distance on the triaxial ellipsoid.

-Formula:

   Dist ~ a' µ { 1 + (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 ( 1 - A ) (8/6) µ + ( B + A Cos µ ) ( µ² / Sin² µ ) ( B + ( 3 A - 2 ) Cos µ ) ) ( 1 - A )
+ B² (3A-3) + B ( - 3A(1-A) Cos µ + (1-A) (7A-2) µ/Sin µ ) + A(1-A) (7A-2) µ Cot µ
+ ( 1-A ) ( A ( 7A-2 )sqrt(7) + ( 11B² - 8A² ) (Sin µ) (Cos µ)/µ ) ]

+   (f'/2) [ -A" + B" (Sin µ)/µ ]
+ ( f f'/4 ) ( -4 AB" - 3 A"B - 3.5 AA" ) } f/2 ~ 1/596.51444 & f'/2 ~ 1/364465 & a' = 6378.137 km

       Sin² µ = Sin² ( λ - λ' ) Cos² β' + [ Cos β Sin β' - Cos β' Sin β Cos ( λ - λ' ) ]²
      Cos µ = Sin β Sin β' + Cos β Cos β' Cos ( λ - λ' )

   -A" = 2 ( 1 - A' ) - A
   -B" = -2 B' - B with:

       A' = ( cos² β Sin² λ + cos² β' Sin² λ' - 2 cos β Sin λ cos β' Sin λ' Cos µ ) / Sin² µ
   B' = [ ( cos² β Sin² λ + cos² β' Sin² λ' ) Cos µ - 2 cos β Sin λ cos β' Sin λ' ] / Sin² µ

A = [ Sin² β + Sin² β' - 2 Sin β Sin β' Cos µ ] / Sin² µ
B = [ ( Sin² β + Sin² β' ) Cos µ - 2 Sin β Sin β' ] / Sin² µ

Data Registers: R00 thru R09: temp R01 = A.A" R02 = -B R08 = -A" R09 = -B"
Flags: /
Subroutines: /

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

94 STO 02
95 ST+ 09
96 RCL 01
97 ST+ 04
98 ST+ 08
99 RCL 03
100 *
101 +
102 STO 05
103 LASTX
104 RCL 03
105 +
106 ST+ X
107 +
108 *
109 RCL 01
110 7
111 *
112 2
113 +
114 STO 07
115 +
116 *
117 RCL 02
118 3
119 *
120 STO 10
121 -
122 RCL 05
123 *
124 RCL 07

125 7
126 SQRT
127 *
128 RCL 00
129 ST+ X
130 X^2
131 RCL 04
132 *
133 3
134 /
135 +
136 RCL 01
137 *
138 +
139 RCL 10
140 RCL 02
141 ST* Y
142 X^2
143 STO 07
144 RCL 01
145 X^2
146 STO 11
147 -
148 ST+ 07
149 8
150 *
151 +
152 RCL 03
153 RCL 06
154 ST* 05
155 ST/ 10

156 /
157 ST* 07
158 *
159 +
160 596.51444
161 STO 03
162 /
163 RCL 05
164 -
165 RCL 10
166 ST+ X
167 -
168 RCL 04
169 *
170 RCL 01
171 ST+ 10
172 RCL 07
173 RCL 11
174 -
175 ST- Y
176 4
177 ST* Z
178 /
179 +
180 -
181 RCL 03
182 /
183 RCL 10
184 +
185 RCL 01
186 RCL 09

187 *
188 4
189 *
190 RCL 02
191 RCL 08
192 ST* 01
193 *
194 3
195 *
196 +
197 RCL 01
198 3.5
199 *
200 +
201 RCL 03
202 ST/ Z
203 /
204 RCL 08
205 -
206 RCL 09
207 RCL 06
208 /
209 +
210 364465
211 /
212 -
213 RCL 00
214 ST* Y
215 +
216 6378.137
217 *
218 END

( 286 bytes / SIZE 012 )

STACK	INPUTS	OUTPUTS
T	λ ( ° ' " )	/
Z	β ( ° ' " )	/
Y	λ' ( ° ' " )	/
X	β' ( ° ' " )	D ( km )

Where D is the geodesic distance on the triaxial ellipsoid

Example: λ = +166°08'02"2 , β = -33°51'41"1
λ' = -101°56'06"0 , β' = +33°21'22"4

     166.08022 ENTER^
    -33.51411   ENTER^
   -101.5606    ENTER^
     33.21224   XEQ "TGB"   >>>>   D = 12138.70169 km                                                      ---Execution time = 12.5s---

-Exact D = 12138.701763....

Notes:

-The precision is:

-If | λ - λ' | = 168° , the maximum error is about 73 cm.
-If | λ - λ' | = 160° , the maximum error is about 32 cm
-If | λ - λ' | = 120° , the maximum error is about 20 cm.

-If | λ - λ' | = 160° , the maximum relative error is about 5.2 10^-8

>>> If we replace ( f f'/4 ) ( -4 AB' - 3 A'B - 3.5 AA' ) by ( f f'/4 ) µ ( - AB' + 0.75 A'B - 2 AA' ) the precision becomes:

-If | λ - λ' | = 168° , the maximum error is about 76 cm.
-If | λ - λ' | = 160° , the maximum error is about 36 cm
-If | λ - λ' | = 120° , the maximum error is about 15 cm.

-If | λ - λ' | = 160° , the maximum relative error is about 3.3 10^-8

( Replace lines 185 to 200 by RCL 01 RCL 09 * RCL 02 RCL 08 ST* 01 * .75 * - RCL 01 ST+ X + RCL 00 *

b) Program#2

-Formula:

   Dist ~ a' µ { 1 + (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 ( 1 - A ) (8/6) µ + ( B + A Cos µ ) ( µ² / Sin² µ ) ( B + ( 3 A - 2 ) Cos µ ) ) ( 1 - A )
+ B² (3A-3) + B ( - 3A(1-A) Cos µ + (1-A) (7A-2) µ/Sin µ ) + A(1-A) (7A-2) µ Cot µ
+ ( 1-A ) ( A ( 7A-2 )sqrt(7) + ( 11B² - 8A² ) (Sin µ) (Cos µ)/µ ) ]

+   (f'/2) [ -A" + B" (Sin µ)/µ ]
+ ( f f'/4 ) µ ( 1 - A ) ( -4 AB" - 3 A"B - 2 AA" ) } f/2 ~ 1/596.51444 & f'/2 ~ 1/364465 & a' = 6378.137 km

       Sin² µ = Sin² ( λ - λ' ) Cos² β' + [ Cos β Sin β' - Cos β' Sin β Cos ( λ - λ' ) ]²
      Cos µ = Sin β Sin β' + Cos β Cos β' Cos ( λ - λ' )

   -A" = 2 ( 1 - A' ) - A
   -B" = -2 B' - B with:

       A' = ( cos² β Sin² λ + cos² β' Sin² λ' - 2 cos β Sin λ cos β' Sin λ' Cos µ ) / Sin² µ
   B' = [ ( cos² β Sin² λ + cos² β' Sin² λ' ) Cos µ - 2 cos β Sin λ cos β' Sin λ' ] / Sin² µ

A = [ Sin² β + Sin² β' - 2 Sin β Sin β' Cos µ ] / Sin² µ
B = [ ( Sin² β + Sin² β' ) Cos µ - 2 Sin β Sin β' ] / Sin² µ

Data Registers: R00 thru R09: temp R01 = -A R02 = A".B R08 = -A" R09 = -B"
Flags: /
Subroutines: /

187 RCL 08
188 ST* 02
189 +
190 ST+ X
191 RCL 01
192 *
193 RCL 02
194 3
195 *
196 +
197 RCL 00
198 *
199 RCL 04
200 *
201 RCL 03
202 ST/ Z
203 /
204 RCL 08
205 -
206 RCL 09
207 RCL 06
208 /
209 +
210 364465
211 /
212 -
213 RCL 00
214 ST* Y
215 +
216 6378.137
217 *
218 END

( 286 bytes / SIZE 012 )

STACK	INPUTS	OUTPUTS
T	λ ( ° ' " )	/
Z	β ( ° ' " )	/
Y	λ' ( ° ' " )	/
X	β' ( ° ' " )	D ( km )

Where D is the geodesic distance on the triaxial ellipsoid

Example: λ = +166°08'02"2 , β = -33°51'41"1
λ' = -101°56'06"0 , β' = +33°21'22"4

     166.08022 ENTER^
    -33.51411   ENTER^
   -101.5606    ENTER^
     33.21224   XEQ "TGB"   >>>>   D = 12138.70168 km                                                      ---Execution time = 12.4s---

-Exact D = 12138.701763....

Notes:

-If | λ - λ' | = 160° , the maximum relative error is about 3.8 10^-8 and the root-mean-square relative error is ~ 10^-8

-Here are a few maximum errors ( in centimeters ):

| λ - λ' | = 120°

( -60°, +37° ) ( +60° , -37° ) m = -22.1 cm
( +30°, +36° ) ( +150° , -22° ) M = +14.2 cm

| λ - λ' | = 160°

( -80°, +53° ) ( +80° , -62° ) M = +38.2 cm
( +10°, +63° ) ( +170° , -54° ) m = -38.7 cm

| λ - λ' | = 164°

( -82°, +49° ) ( +82° , -57° ) M = +52.6 cm
( +8° , +55° ) ( +172° , -62° ) m = -52.3 cm

| λ - λ' | = 168°

( -84°, +48° ) ( +84° , -56° ) M = +78.1 cm
( +6° , +52° ) ( +174° , -58° ) m = -76.9 cm

3°) Geocentric Longitudes & Parametric Latitudes

-This program takes the geocentric longitudes & latitudes and computes the parametric latitudes ( WGS84 )
-Then it uses the Andoyer's formula of order 3 and 2 terms with f'/2 & f.f'/4

>>> If you want to use the geodetic latitudes, simply replace lines 23 & 31 ( * ) by /

Formula:

   Dist / ( a' µ ) ~ 1 - (f/2) ( A + B ( Sin µ ) / µ )
+ (f²/4) [ ( 1 - A ) ( µ / Sin µ ) ( B + A Cos µ ) + (1/2) [ ( A²/2 ) ( 1 + (Sin µ Cos µ) / µ ) - B² (Sin µ Cos µ) / µ ) ) ]
+ (f³/8) { ( 1 - A ) { B² + A [ B Cos µ - ( 1 - A ) (8/6) µ² ] + ( B + A Cos µ ) ( µ / Sin µ ) [ - ( 3 A - 2 ) + ( µ / Sin µ ) ( ( 3 A - 2 ) Cos µ + B ) ] }
+ [ B³ ( (2/3) Sin² µ - 1/2 ) + ( A/2 - 1 ) B² Cos µ ] (Sin µ) / µ ) - (A²/2) [ ( A/2 - 1 ) ( 1 + (Sin µ Cos µ) / µ ) - B (Sin µ Cos² µ) / µ ] }

+   (f'/2) [ -A"+ B" (Sin µ)/µ ]
+ ( f f'/4 ) µ (1-A) [ µ ( - AB" - A"B - 2 AA" ) - BB" ] f/2 ~ 1/596.51444 & f'/2 ~ 1/364465 & a' = 6378.137 km

     Sin² µ = Sin² ( λ - λ' ) Cos² β' + [ Cos β Sin β' - Cos β' Sin β Cos ( λ - λ' ) ]²
      Cos µ = Sin β Sin β' + Cos β Cos β' Cos ( λ - λ' )

A = [ Sin² β + Sin² β' - 2 Sin β Sin β' Cos µ ] / Sin² µ
B = [ ( Sin² β + Sin² β' ) Cos µ - 2 Sin β Sin β' ] / Sin² µ

   -A" = 2 ( 1 - A' ) - A
   -B" = -2 B' - B with:

   A' = ( cos² β Sin² λ + cos² β' Sin² λ' - 2 cos β Sin λ cos β' Sin λ' Cos µ ) / Sin² µ
   B' = [ ( cos² β Sin² λ + cos² β' Sin² λ' ) Cos µ - 2 cos β Sin λ cos β' Sin λ' ] / Sin² µ

Data Registers: R00 thru R09: temp R00 = µ R01 = A R02 = -B.B" R11 = -A" R12 = -B"
Flags: /
Subroutines: /

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

38 STO 00
39 STO 05
40 RCL 01
41 SIN
42 *
43 STO 08
44 ST+ 08
45 X^2
46 RCL 09
47 RCL 04
48 P-R
49 STO 06
50 X<>Y
51 X<> 07
52 ST- 01
53 SIN
54 *
55 ST* 08
56 X^2
57 +
58 STO 03
59 X<>Y
60 STO 02
61 ST+ 02
62 X^2
63 RCL 07
64 ST* 02
65 ST* 05
66 X^2
67 +
68 X<> 01
69 RCL 06
70 P-R
71 ST* 00
72 R^
73 ST* 07
74 *

75 RCL 05
76 -
77 X^2
78 X<>Y
79 X^2
80 +
81 ST/ 01
82 ST/ 02
83 ST/ 03
84 ST/ 08
85 SQRT
86 STO 09
87 RCL 08
88 RCL 00
89 RCL 07
90 +
91 STO 12
92 ST* 08
93 RCL 03
94 ST- 08
95 *
96 -
97 ST+ X
98 X<> 12
99 STO 03
100 R-P
101 RCL 08
102 +
103 ST+ X
104 STO 11
105 X<>Y
106 D-R
107 STO 00
108 R^
109 /
110 STO 06
111 RCL 01

112 RCL 02
113 RCL 03
114 ST* 01
115 *
116 -
117 ST- 04
118 ST- 11
119 X<> 01
120 RCL 02
121 -
122 STO 02
123 ST- 12
124 RCL 01
125 RCL 03
126 *
127 +
128 *
129 STO 08
130 RCL 01
131 3
132 *
133 2
134 -
135 STO 07
136 RCL 03
137 *
138 RCL 02
139 +
140 R^
141 *
142 RCL 07
143 -
144 *
145 RCL 02
146 RCL 03
147 *
148 RCL 00

149 ST+ X
150 X^2
151 RCL 04
152 *
153 3
154 /
155 -
156 RCL 01
157 *
158 +
159 RCL 02
160 X^2
161 STO 07
162 +
163 RCL 04
164 ST* 08
165 *
166 RCL 09
167 X^2
168 1.5
169 ST/ Y
170 FRC
171 -
172 RCL 02
173 *
174 RCL 01
175 2
176 ST- Y
177 /
178 STO 09
179 RCL 03
180 *
181 +
182 RCL 07
183 *
184 RCL 06
185 /

186 +
187 RCL 03
188 RCL 06
189 /
190 ST* 03
191 ST* 07
192 1
193 +
194 STO 05
195 RCL 09
196 *
197 RCL 02
198 RCL 03
199 *
200 -
201 RCL 01
202 X^2
203 2
204 /
205 ST* 05
206 *
207 -
208 RCL 10
209 /
210 RCL 05
211 RCL 07
212 -
213 2
214 /
215 +
216 RCL 08
217 +
218 RCL 10
219 /
220 RCL 02
221 RCL 06
222 /

223 RCL 01
224 +
225 -
226 RCL 02
227 RCL 11
228 ST* Y
229 ST+ X
230 RCL 12
231 ST* 02
232 +
233 RCL 01
234 *
235 +
236 RCL 00
237 ST* 04
238 *
239 RCL 02
240 +
241 RCL 04
242 *
243 RCL 10
244 ST/ Z
245 /
246 RCL 11
247 +
248 RCL 12
249 RCL 06
250 /
251 -
252 364465
253 /
254 +
255 RCL 00
256 ST* Y
257 +
258 6378.137
259 *
260 END

( 329 bytes / SIZE 013 )

STACK	INPUTS	OUTPUTS
T	λ ( ° ' " )	/
Z	β ( ° ' " )	/
Y	λ' ( ° ' " )	/
X	β' ( ° ' " )	D ( km )

Where D is the geodesic distance on the triaxial ellipsoid

Example1: λ = +166°08'02"2 , β = -33°41'00"9 ( with geocentric latitudes )
λ' = -101°56'06"0 , β' = +33°10'47"0

     166.08022 ENTER^
    -33.41009 ENTER^
   -101.56060    ENTER^
     33.10470 XEQ "TGCB"   >>>>   D = 12138.70246 km                                                      ---Execution time = 16s---

Example2:     λ = +166°08'02"2 , β = -33°51'41"1 ( with geodetic latitudes )
                   λ' = -101°56'06"0 ,  β' = +33°21'22"4

-After replacing lines 23 & 31 with /

     166.08022 ENTER^
    -33.51411 ENTER^
   -101.56060    ENTER^
     33.21224 XEQ "TGCB"   >>>>   D = 12138.70173 km

Notes:

-If | λ - λ' | = 160° , the maximum relative error is about 2.22 10^-8 and the root-mean-square relative error is ~ 0.73 10^-8
-If | λ - λ' | = 120° , the root-mean-square relative error is ~ 0.36 10^-8

-Here are a few maximum errors ( in centimeters ):

| λ - λ' | = 0°

( 0°, 26° ) ( 0° , -26° ) M = +6.8 cm
( 90° , 26° ) ( 90° , -26° ) m = -6.8 cm

| λ - λ' | = 120°

( -60°, +34° ) ( +60° , -34° ) m = -16.5 cm
( +30°, +34° ) ( +150° , -34° ) M = +16.5 cm

| λ - λ' | = 160°

( -65°, +58° ) ( +95° , -52° ) M = +41.3 cm
( +10°, +62° ) ( +170° , -54° ) m = -40.6 cm

| λ - λ' | = 168°

( -69°, +52° ) ( +99° , -47° ) M = +80.6 cm
( +6° , +5° ) ( +174° , +5° ) m = -84.3 cm

4°) Geocentric Longitudes & Latitudes

a) Program(s)#1

Formula:

   Dist/(a µ) ~ 1 + (f/2) ( -A + B (Sin µ)/µ )
                + (f²/4) [ A ( -4 + µ Cot µ ) + B ( 2 (Sin µ)/µ + µ / Sin µ ) + (A²/4) ( 17 + (Sin µ) (Cos µ)/µ - 4 µ Cot µ ) - A.B ( 2 (Sin µ)/µ + µ / Sin µ ) - ( B²/2 ) (Sin µ) (Cos µ)/µ ]
                  + (f³/8) { ( 1 - A ) [ ( -A ( 1 - A ) (4/3) µ² + A ( 21.A - 13 ) + ( B + A Cos µ ) [ ( µ / Sin µ ) ( 6 - 13.A) + 5.B + ( µ² / Sin² µ ) ( B + ( 3 A - 2 ) Cos µ ) ) ]
      + ( 3.A² - 9.B² ) ( Sin µ ) ( Cos µ ) / µ ] + B³ [ ( Sin µ ) / µ + 2 ( Sin³ µ ) / µ ] }

+   (f'/2) [ -A'+ B' (Sin µ)/µ ]
+ ( f f'/4 ) (1-A) (1-A') µ² ( -5 A B' - 5 A' B + 6 A A' + 5 B B ' ) f/2 ~ 1/595.543 & f'/2 ~ 1/182233 & a = 6378.172 km

     Sin² µ = Sin² ( λ - λ' ) Cos² β' + [ Cos β Sin β' - Cos β' Sin β Cos ( λ - λ' ) ]²
      Cos µ = Sin β Sin β' + Cos β Cos β' Cos ( λ - λ' )

A = [ Sin² β + Sin² β' - 2 Sin β Sin β' Cos µ ] / Sin² µ
B = [ ( Sin² β + Sin² β' ) Cos µ - 2 Sin β Sin β' ] / Sin² µ

   A' = ( cos² β Sin² λ + cos² β' Sin² λ' - 2 cos β Sin λ cos β' Sin λ' Cos µ ) / Sin² µ
   B' = [ ( cos² β Sin² λ + cos² β' Sin² λ' ) Cos µ - 2 cos β Sin λ cos β' Sin λ' ] / Sin² µ

Data Registers: R00 thru R09: temp R00 = µ R01 = A.A' R02 = -B R08 = -A' R09 = -B'
Flags: /
Subroutines: /

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

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

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

91 +
92 STO 05
93 LASTX
94 RCL 03
95 +
96 ST+ X
97 +
98 *
99 RCL 01
100 13
101 *
102 STO T
103 -
104 6
105 -
106 *
107 RCL 02
108 STO 11
109 5
110 *
111 +
112 RCL 05
113 *
114 +
115 RCL 01
116 X^2
117 STO 07
118 21
119 *
120 +

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

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

181 RCL 09
182 -
183 *
184 5
185 *
186 RCL 01
187 +
188 RCL 08
189 RCL 03
190 ST* Y
191 +
192 *
193 RCL 10
194 ST/ Z
195 /
196 RCL 08
197 +
198 RCL 09
199 RCL 06
200 /
201 -
202 182233
203 /
204 +
205 RCL 00
206 ST* Y
207 +
208 6378.172
209 *
210 END

( 282 bytes / SIZE 013 )

STACK	INPUTS	OUTPUTS
T	λ ( deg )	/
Z	β ( deg )	/
Y	λ' ( deg )	/
X	β' ( deg )	D ( km )

Where D is the geodesic distance on the triaxial ellipsoid

Example: λ = -24° , β = -10°
                λ' = +144° ,  β' = +41°

-24 ENTER^
-10 ENTER^
144 ENTER^
41 XEQ "TGC" >>>> D = 16362.07202 km                                                       ---Execution time = 12s---

-Exact D = 16362.072084...

Notes:

-Here are a few maximum errors ( in centimeters ):

| λ - λ' | = 120°

( -30°, +32° ) ( +90°, -48° ) m = -18.7 cm
( -45°, +62° ) ( +75° , -34° ) M = +18.1 cm

| λ - λ' | = 140°

( -25°, +39° ) ( +115°, -52° ) m = -28.0 cm
( -55°, +60° ) ( +85° , -43° ) M = +26.9 cm

| λ - λ' | = 160°

( -65°, +57° ) ( +95° , -49° ) M = +52.5 cm
( -20°, +42° ) ( +140° , -50° ) m = -57.5 cm

| λ - λ' | = 168°

( -69°, +55° ) ( +99° , -50° ) M = +85.2 cm
( -9°, +38° ) ( +159° , -44° ) m = -107.5 cm

-Here is another program slightly more accurate:

Formula:

   Dist/(a µ) ~ 1 + (f/2) ( -A + B (Sin µ)/µ )
                + (f²/4) [ A ( -4 + µ Cot µ ) + B ( 2 (Sin µ)/µ + µ / Sin µ ) + (A²/4) ( 17 + (Sin µ) (Cos µ)/µ - 4 µ Cot µ ) - A.B ( 2 (Sin µ)/µ + µ / Sin µ ) - ( B²/2 ) (Sin µ) (Cos µ)/µ ]
                  + (f³/8) { ( 1 - A ) [ ( -A ( 1 - A ) (4/3) µ² + A ( 21.A - 13 ) + ( B + A Cos µ ) [ ( µ / Sin µ ) ( 6 - 13.A) + 5.B + ( µ² / Sin² µ ) ( B + ( 3 A - 2 ) Cos µ ) ) ]
      + ( 3.A² - 9.B² ) ( Sin µ ) ( Cos µ ) / µ ] + B³ [ ( Sin µ ) / µ + 2 ( Sin³ µ ) / µ ] }

+   (f'/2) [ -A'+ B' (Sin µ)/µ ]
+ ( f f'/4 ) (1-A) (1-A') µ² [ ( - 3 A B' - 3 A' B + 4 A A' + 3 B B ' ) + µ ( - A B' - A' B + A A' + B B ' ) ]

f/2 ~ 1/595.543 & f'/2 ~ 1/182233 & a = 6378.172 km

     Sin² µ = Sin² ( λ - λ' ) Cos² β' + [ Cos β Sin β' - Cos β' Sin β Cos ( λ - λ' ) ]²
      Cos µ = Sin β Sin β' + Cos β Cos β' Cos ( λ - λ' )

A = [ Sin² β + Sin² β' - 2 Sin β Sin β' Cos µ ] / Sin² µ
B = [ ( Sin² β + Sin² β' ) Cos µ - 2 Sin β Sin β' ] / Sin² µ

   A' = ( cos² β Sin² λ + cos² β' Sin² λ' - 2 cos β Sin λ cos β' Sin λ' Cos µ ) / Sin² µ
   B' = [ ( cos² β Sin² λ + cos² β' Sin² λ' ) Cos µ - 2 cos β Sin λ cos β' Sin λ' ] / Sin² µ

Data Registers: R00 thru R11: temp R00 = µ R01 = A.A' R02 = -B R08 = -A' R09 = -B'
Flags: /
Subroutines: /

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

91 +
92 STO 05
93 LASTX
94 RCL 03
95 +
96 ST+ X
97 +
98 *
99 RCL 01
100 13
101 *
102 STO T
103 -
104 6
105 -
106 *
107 RCL 02
108 5
109 *
110 +
111 RCL 05
112 *
113 +
114 RCL 01
115 X^2
116 STO 07
117 21
118 *
119 +
120 RCL 07

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

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

183 -
184 *
185 RCL 00
186 3
187 +
188 *
189 RCL 01
190 +
191 RCL 08
192 RCL 03
193 ST* Y
194 +
195 *
196 RCL 10
197 ST/ Z
198 /
199 RCL 08
200 +
201 RCL 09
202 RCL 06
203 /
204 -
205 182233
206 /
207 +
208 RCL 00
209 ST* Y
210 +
211 6378.172
212 *
213 END

( 284 bytes / SIZE 012 )

STACK	INPUTS	OUTPUTS
T	λ ( deg )	/
Z	β ( deg )	/
Y	λ' ( deg )	/
X	β' ( deg )	D ( km )

Where D is the geodesic distance on the triaxial ellipsoid

Example: λ = -24° , β = -10°
λ' = +144° , β' = +41°

-24 ENTER^
-10 ENTER^
144 ENTER^
41 XEQ "TGC" >>>> D = 16362.07202 km

-Exact D = 16362.072084...

Notes:

-Here are a few maximum errors ( in centimeters ):

| λ - λ' | = 120°

( -30°, +32° ) ( +90°, -48° ) m = -18.6 cm
( -45°, +62° ) ( +75° , -34° ) M = +18.4 cm

| λ - λ' | = 140°

( -25°, +40° ) ( +115°, -53° ) m = -27.0 cm
( -55°, +60° ) ( +85° , -43° ) M = +27.8 cm

| λ - λ' | = 160°

( -65°, +57° ) ( +95° , -49° ) M = +53.9 cm
( -20°, +46° ) ( +140° , -53° ) m = -54.2 cm

| λ - λ' | = 168°

( -69°, +55° ) ( +99° , -50° ) M = +87.0 cm
( -24°, +40° ) ( +144° , -45° ) m = -101.6 cm

-The maximum errors are smaller in the following programs.

b) Program(s)#2

Formula:

   Dist ~ a' µ { 1 + (f/2) ( -A + B (Sin µ)/µ )
               + (f²/4) [ A ( -4 + µ Cot µ ) + B ( 2 (Sin µ)/µ + µ / Sin µ ) + (A²/4) ( 17 + (Sin µ) (Cos µ)/µ - 4 µ Cot µ ) - A.B ( 2 (Sin µ)/µ + µ / Sin µ ) - ( B²/2 ) (Sin µ) (Cos µ)/µ ]
               + (f³/8) ( 1 - A ) [ ( -A ( 1 - A ) (4/3) µ² + ( B + A Cos µ ) ( µ² / Sin² µ ) ( B + ( 3 A - 2 ) Cos µ ) )
+   ( 7.6 B² + 7.6 A.B Cos µ ) + (µ/Sin µ) ( B (5.6 - 12.2 A) + A (Cos µ) (7.6-15.2 A) ) ]

+   (f'/2) [ -A"+ B" (Sin µ)/µ ]
+ ( f f'/4 ) µ (1-A) [ µ ( -2 AB" - 4 AA" ) + 2 BB" ] } f/2 ~ 1/596.51444 & f'/2 ~ 1/364465 & a' = 6378.137 km

     Sin² µ = Sin² ( λ - λ' ) Cos² β' + [ Cos β Sin β' - Cos β' Sin β Cos ( λ - λ' ) ]²
      Cos µ = Sin β Sin β' + Cos β Cos β' Cos ( λ - λ' )

A = [ Sin² β + Sin² β' - 2 Sin β Sin β' Cos µ ] / Sin² µ
B = [ ( Sin² β + Sin² β' ) Cos µ - 2 Sin β Sin β' ] / Sin² µ

   -A" = 2 ( 1 - A' ) - A
   -B" = -2 B' - B with:

   A' = ( cos² β Sin² λ + cos² β' Sin² λ' - 2 cos β Sin λ cos β' Sin λ' Cos µ ) / Sin² µ
   B' = [ ( cos² β Sin² λ + cos² β' Sin² λ' ) Cos µ - 2 cos β Sin λ cos β' Sin λ' ] / Sin² µ

Data Registers: R00 thru R09: temp R00 = µ R01 = -A R02 = -B R08 = -A" R09 = -B"
Flags: /
Subroutines: /

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

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

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

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

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

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

179 +
180 RCL 01
181 *
182 RCL 00
183 ST* 04
184 *
185 -
186 ST+ X
187 RCL 04
188 *
189 RCL 07
190 ST/ Z
191 /
192 RCL 08
193 +
194 RCL 09
195 RCL 06
196 /
197 -
198 364465
199 /
200 +
201 RCL 00
202 ST* Y
203 +
204 6378.137
205 *
206 END

( 277 bytes / SIZE 010 )

STACK	INPUTS	OUTPUTS
T	λ ( ° ' " )	/
Z	β ( ° ' " )	/
Y	λ' ( ° ' " )	/
X	β' ( ° ' " )	D ( km )

Where D is the geodesic distance on the triaxial ellipsoid

Example1: λ = +166°08'02"2 , β = -33°41'00"9
λ' = -101°56'06"0 , β' = +33°10'47"0

     166.08022 ENTER^
    -33.41009 ENTER^
   -101.56060    ENTER^
     33.10470 XEQ "TGC"   >>>>   D = 12138.70252 km                                                      ---Execution time = 11.9s---

-Exact D = 12138.702497.... ( with free42, TGC gives 12138.702511... )

Example2: λ = -16° , β = -16°
λ' = +157° , β' = +13°

     -16 ENTER^ ENTER^
   157 ENTER^
     13   R/S >>>>   D = 19205.13660 km

-Exact D = 19205.13653....
-In this example, the great elliptic arc length = 19205.463988...

Notes:

-If | λ - λ' | = 160° , the maximum relative error is about 2.28 10^-8 and the root-mean-square relative error is ~ 0.63 10^-8
-If | λ - λ' | = 120° , the root-mean-square relative error is ~ 0.47 10^-8

-Here are a few maximum errors ( in centimeters ):

| λ - λ' | = 0°

( 0° , 51° ) ( 0° , -51° ) m = -11.5 cm
( 90° , 90° ) ( 90° , 43° ) M = +4.3 cm

| λ - λ' | = 120°

( +15°, +69° ) ( +135° , -46° ) m = -15.7 cm
( +0°, +30° ) ( +120° , -9° ) M = +19.2 cm

| λ - λ' | = 160°

( -65°, +67° ) ( +95° , -62° ) M = +43.2 cm
( 10°, +53° ) ( +170° , -61° ) m = -41.6 cm

| λ - λ' | = 168°

( -69°, +61° ) ( +99° , -57° ) M = +83.1 cm
( -9° , +45° ) ( +159° , -50° ) m = -80.1 cm

Remarks:

-If we replace lines 201 to 206 with:

RCL 00 ST* Y ST* Z ST+ Z + 6378.137 ST* Z * END

Y-output = Geodesic distance on the ellipsoid of revolution ( WGS84 )

-For instance:

     166.08022 ENTER^
    -33.41009 ENTER^
   -101.56060    ENTER^
     33.10470 XEQ "TGC"   >>>>   D = 12138.70252 km
X<>Y d = 12138.68509 km ( WGS84 )

( cf "Terrestrial Geodesic Distance(3)" paragraph 1°)d) )

-With the ellipsoid model of the Earth, f.f'/4 ~ f³/8 , so we can save 10 bytes and ~ 0.4 second with the following program

Data Registers: R00 thru R09: temp R00 = µ R01 = -A R02 = -B R08 = -A" R09 = -B" (Sin µ)/µ
Flags: /
Subroutines: /

29 ST- 05
30 SIN
31 *
32 ST* 08
33 X^2
34 +
35 STO 03
36 X<>Y
37 STO 01
38 ST+ 01
39 X^2
40 RCL 07
41 ST* 01
42 ST* 02
43 X^2
44 +
45 X<> 05
46 RCL 06
47 P-R
48 ST* 00
49 R^
50 ST* 07
51 *
52 RCL 02
53 -
54 X^2
55 X<>Y
56 X^2

57 +
58 ST/ 01
59 ST/ 03
60 ST/ 05
61 ST/ 08
62 SQRT
63 RCL 08
64 RCL 00
65 RCL 07
66 +
67 STO 09
68 ST* 08
69 RCL 03
70 ST- 08
71 *
72 -
73 ST+ X
74 X<> 09
75 STO 03
76 R-P
77 RCL 08
78 +
79 ST+ X
80 STO 08
81 X<>Y
82 D-R
83 STO 00
84 R^

85 /
86 STO 06
87 RCL 01
88 RCL 03
89 ST* 01
90 RCL 05
91 ST- 01
92 *
93 -
94 STO 02
95 ST+ 09
96 RCL 01
97 ST+ 04
98 ST+ 08
99 RCL 03
100 *
101 +
102 STO 07
103 *
104 STO 05
105 RCL 03
106 RCL 04
107 ST+ X
108 *
109 RCL 07
110 +
111 *
112 RCL 07

113 7.6
114 *
115 STO 07
116 RCL 02
117 ST* 07
118 ST- Y
119 -
120 ST- Y
121 ST+ X
122 LASTX
123 +
124 RCL 01
125 *
126 -
127 *
128 RCL 07
129 +
130 RCL 02
131 RCL 04
132 3
133 /
134 RCL 08
135 -
136 ST+ X
137 RCL 09
138 ST* Z
139 -
140 RCL 00
141 *

142 RCL 01
143 *
144 +
145 ST+ X
146 RCL 00
147 *
148 +
149 596.51444
150 STO 07
151 /
152 RCL 05
153 RCL 02
154 RCL 06
155 ST/ 03
156 ST/ 09
157 /
158 STO 05
159 ST+ X
160 +
161 -
162 RCL 04
163 *
164 RCL 01
165 ST* 04
166 ST- 05
167 X^2
168 RCL X
169 RCL 02
170 X^2

171 ST+ X
172 -
173 RCL 03
174 *
175 +
176 RCL 04
177 4
178 ST/ Z
179 *
180 +
181 +
182 RCL 07
183 /
184 RCL 05
185 -
186 RCL 08
187 RCL 09
188 -
189 611
190 /
191 +
192 RCL 07
193 /
194 RCL 00
195 ST* Y
196 +
197 6378.137
198 *
199 END

( 267 bytes / SIZE 010 )

STACK	INPUTS	OUTPUTS
T	λ ( ° ' " )	/
Z	β ( ° ' " )	/
Y	λ' ( ° ' " )	/
X	β' ( ° ' " )	D ( km )

Where D is the geodesic distance on the triaxial ellipsoid

Example: λ = +166°08'02"2 , β = -33°41'00"9
λ' = -101°56'06"0 , β' = +33°10'47"0

     166.08022 ENTER^
    -33.41009 ENTER^
   -101.56060    ENTER^
     33.10470 XEQ "TGC"   >>>>   D = 12138.70252 km                                                      ---Execution time = 11.5s---

-Exact D = 12138.702497....

c) Program#3

Formula:

   Dist ~ a' µ { 1 + (f/2) ( -A + B (Sin µ)/µ )
               + (f²/4) [ A ( -4 + µ Cot µ ) + B ( 2 (Sin µ)/µ + µ / Sin µ ) + (A²/4) ( 17 + (Sin µ) (Cos µ)/µ - 4 µ Cot µ ) - A.B ( 2 (Sin µ)/µ + µ / Sin µ ) - ( B²/2 ) (Sin µ) (Cos µ)/µ ]
               + (f³/8) ( 1 - A ) [ ( -A ( 1 - A ) (4/3) µ² + A ( 24 A - 12 ) + ( B + A Cos µ ) [ 8 B + ( µ / Sin µ ) ( 6 - 12 A ) + ( µ² / Sin² µ ) ( B + ( 3 A - 2 ) Cos µ ) ) ] ]

+   (f'/2) [ -A"+ B" (Sin µ)/µ ]
+ ( f f'/4 ) µ (1-A) [ µ ( -2.AB" - 4 AA" ) + 2 BB" ] } f/2 ~ 1/596.51444 & f'/2 ~ 1/364465 & a' = 6378.137 km

     Sin² µ = Sin² ( λ - λ' ) Cos² β' + [ Cos β Sin β' - Cos β' Sin β Cos ( λ - λ' ) ]²
      Cos µ = Sin β Sin β' + Cos β Cos β' Cos ( λ - λ' )

A = [ Sin² β + Sin² β' - 2 Sin β Sin β' Cos µ ] / Sin² µ
B = [ ( Sin² β + Sin² β' ) Cos µ - 2 Sin β Sin β' ] / Sin² µ

   -A" = 2 ( 1 - A' ) - A
   -B" = -2 B' - B with:

   A' = ( cos² β Sin² λ + cos² β' Sin² λ' - 2 cos β Sin λ cos β' Sin λ' Cos µ ) / Sin² µ
   B' = [ ( cos² β Sin² λ + cos² β' Sin² λ' ) Cos µ - 2 cos β Sin λ cos β' Sin λ' ] / Sin² µ

Data Registers: R00 thru R09: temp R00 = µ R01 = -A R02 = -B R08 = -A" R09 = -B"
Flags: /
Subroutines: /

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

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

59 ST/ 03
60 ST/ 05
61 ST/ 08
62 SQRT
63 RCL 08
64 RCL 00
65 RCL 09
66 +
67 ST* 08
68 STO 09
69 RCL 03
70 ST- 08
71 *
72 -
73 ST+ X
74 X<> 09
75 STO 03
76 R-P
77 RCL 08
78 +
79 ST+ X
80 STO 08
81 X<>Y
82 D-R
83 STO 00
84 R^
85 /
86 STO 06
87 RCL 01

88 RCL 03
89 ST* 01
90 RCL 05
91 ST- 01
92 *
93 -
94 STO 02
95 ST+ 09
96 RCL 01
97 ST+ 04
98 ST+ 08
99 RCL 03
100 *
101 +
102 STO 05
103 LASTX
104 RCL 03
105 +
106 ST+ X
107 +
108 *
109 RCL 01
110 ST+ X
111 STO 07
112 6
113 ST* Y
114 +
115 ST* 07
116 -

117 *
118 RCL 02
119 8
120 *
121 +
122 RCL 05
123 *
124 ST+ 07
125 RCL 01
126 4
127 *
128 RCL 00
129 X^2
130 RCL 04
131 3
132 /
133 *
134 *
135 RCL 07
136 +
137 596.51444
138 STO 07
139 /
140 +
141 RCL 05
142 RCL 02
143 RCL 06
144 ST* Z
145 ST/ 03

146 /
147 STO 05
148 ST+ X
149 +
150 -
151 RCL 04
152 *
153 RCL 01
154 ST- 05
155 X^2
156 RCL X
157 RCL 02
158 X^2
159 ST+ X
160 -
161 RCL 03
162 *
163 +
164 4
165 /
166 +
167 RCL 07
168 /
169 RCL 05
170 -
171 RCL 02
172 RCL 08
173 ST+ X
174 RCL 09

175 ST* Z
176 +
177 RCL 01
178 *
179 RCL 00
180 ST* 04
181 *
182 -
183 ST+ X
184 RCL 04
185 *
186 RCL 07
187 ST/ Z
188 /
189 RCL 08
190 +
191 RCL 09
192 RCL 06
193 /
194 -
195 364465
196 /
197 +
198 RCL 00
199 ST* Y
200 +
201 6378.137
202 *
203 END

( 272 bytes / SIZE 010 )

STACK	INPUTS	OUTPUTS
T	λ ( ° ' " )	/
Z	β ( ° ' " )	/
Y	λ' ( ° ' " )	/
X	β' ( ° ' " )	D ( km )

Where D is the geodesic distance on the triaxial ellipsoid

Example: λ = +166°08'02"2 , β = -33°41'00"9
λ' = -101°56'06"0 , β' = +33°10'47"0

     166.08022 ENTER^
    -33.41009 ENTER^
   -101.56060    ENTER^
     33.10470 XEQ "TGC"   >>>>   D = 12138.70252 km                                                      ---Execution time = 11.8s---

-Exact D = 12138.702497....

Notes:

-If | λ - λ' | = 160° , the maximum relative error is about 2.27 10^-8 and the root-mean-square relative error is ~ 0.63 10^-8
-If | λ - λ' | = 120° , the root-mean-square relative error is ~ 0.46 10^-8

-Here are a few maximum errors ( in centimeters ):

| λ - λ' | = 0°

( 0° , 51° ) ( 0° , -51° ) m = -11.5 cm
( 90° , 90° ) ( 90° , 43° ) M = +4.3 cm

| λ - λ' | = 120°

( -45°, +19° ) ( +75° , -30° ) m = -11.9 cm
( +15°, +36° ) ( +135° , -25° ) M = +15.3 cm

| λ - λ' | = 160°

( -65°, +62° ) ( +95° , -56° ) M = +42.2 cm
( +10°, +63° ) ( +170° , -55° ) m = -41.4 cm

| λ - λ' | = 168°

( +6°, +55° ) ( +174°, -49° ) m = -82.2 cm
( -69°, +58° ) ( +99° , -54° ) M = +76.9 cm

d) Program#4

Formula:

   Dist/(a' µ) ~ 1 + (f/2) ( -A + B (Sin µ)/µ )
                + (f²/4) [ A ( -4 + µ Cot µ ) + B ( 2 (Sin µ)/µ + µ / Sin µ ) + (A²/4) ( 17 + (Sin µ) (Cos µ)/µ - 4 µ Cot µ ) - A.B ( 2 (Sin µ)/µ + µ / Sin µ ) - ( B²/2 ) (Sin µ) (Cos µ)/µ ]
                  + (f³/8) { ( 1 - A ) [ ( -A ( 1 - A ) (4/3) µ² + A ( 21.A - 13 ) + ( B + A Cos µ ) [ ( µ / Sin µ ) ( 6 - 13.A) + 5.B + ( µ² / Sin² µ ) ( B + ( 3 A - 2 ) Cos µ ) ) ]
      + ( 3.A² - 9.B² ) ( Sin µ ) ( Cos µ ) / µ ] + B³ [ ( Sin µ ) / µ + 2 ( Sin³ µ ) / µ ] }

+   (f'/2) [ -A"+ B" (Sin µ)/µ ]
+ ( f f'/4 ) (1-A) [ µ² ( -2 AB" - 4 AA" ) + 5 BB" ]   f/2 ~ 1/596.51444 & f'/2 ~ 1/364465 & a' = 6378.137 km

     Sin² µ = Sin² ( λ - λ' ) Cos² β' + [ Cos β Sin β' - Cos β' Sin β Cos ( λ - λ' ) ]²
      Cos µ = Sin β Sin β' + Cos β Cos β' Cos ( λ - λ' )

A = [ Sin² β + Sin² β' - 2 Sin β Sin β' Cos µ ] / Sin² µ
B = [ ( Sin² β + Sin² β' ) Cos µ - 2 Sin β Sin β' ] / Sin² µ

   -A" = 2 ( 1 - A' ) - A
   -B" = -2 B' - B with:

   A' = ( cos² β Sin² λ + cos² β' Sin² λ' - 2 cos β Sin λ cos β' Sin λ' Cos µ ) / Sin² µ
   B' = [ ( cos² β Sin² λ + cos² β' Sin² λ' ) Cos µ - 2 cos β Sin λ cos β' Sin λ' ] / Sin² µ

Data Registers: R00 thru R12: temp R00 = µ R01 = -A + B ( Sin µ ) / µ R02 = - B ( Sin µ ) / µ R08 = -A" R09 = -B" R12 = 5.B.B" ( 1 - A )
Flags: /
Subroutines: /

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

63 ST+ 10
64 SQRT
65 RCL 08
66 RCL 00
67 RCL 09
68 +
69 ST* 08
70 STO 09
71 RCL 03
72 ST- 08
73 *
74 -
75 ST+ X
76 X<> 09
77 STO 03
78 R-P
79 RCL 08
80 +
81 ST+ X
82 STO 08
83 X<>Y
84 D-R
85 STO 00
86 R^
87 /
88 STO 06
89 RCL 01
90 RCL 03
91 ST* 01
92 RCL 05
93 ST- 01

94 *
95 -
96 STO 02
97 ST+ 09
98 RCL 01
99 ST+ 04
100 ST+ 08
101 RCL 03
102 *
103 +
104 STO 05
105 LASTX
106 RCL 03
107 +
108 ST+ X
109 +
110 *
111 RCL 01
112 13
113 *
114 STO T
115 -
116 6
117 -
118 *
119 RCL 02
120 5
121 *
122 STO 12
123 +
124 RCL 05

125 *
126 +
127 RCL 01
128 X^2
129 STO 07
130 21
131 *
132 +
133 RCL 07
134 RCL 02
135 X^2
136 ST* 10
137 ST+ 10
138 STO 11
139 ST+ X
140 -
141 RCL 03
142 RCL 06
143 ST/ 02
144 ST* 05
145 /
146 ST* 11
147 *
148 ST+ 07
149 RCL 11
150 -
151 RCL 01
152 RCL 04
153 ST* 12
154 *
155 4

156 ST/ 07
157 *
158 ST+ 07
159 RCL 00
160 X^2
161 *
162 STO 11
163 3
164 ST* Z
165 /
166 596.51444
167 1/X
168 ST* 10
169 RDN
170 +
171 +
172 *
173 RCL 05
174 RCL 02
175 ST- 01
176 ST* 10
177 ST+ X
178 +
179 -
180 RCL 04
181 *
182 RCL 10
183 -
184 RCL 07
185 +
186 *

187 RCL 01
188 +
189 *
190 RCL 09
191 ST* 12
192 2
193 /
194 RCL 08
195 +
196 RCL 11
197 *
198 RCL 12
199 -
200 R^
201 *
202 RCL 08
203 -
204 RCL 09
205 RCL 06
206 /
207 +
208 364465
209 /
210 -
211 RCL 00
212 ST* Y
213 +
214 6378.137
215 *
216 END

( 290 bytes / SIZE 013 )

STACK	INPUTS	OUTPUTS
T	λ ( ° ' " )	/
Z	β ( ° ' " )	/
Y	λ' ( ° ' " )	/
X	β' ( ° ' " )	D ( km )

Where D is the geodesic distance on the triaxial ellipsoid

Example1: λ = +166°08'02"2 , β = -33°41'00"9
λ' = -101°56'06"0 , β' = +33°10'47"0

     166.08022 ENTER^
    -33.41009 ENTER^
   -101.56060    ENTER^
     33.10470 XEQ "TGC"   >>>>   D = 12138.70249 km                                                      ---Execution time = 12.4s---

-Exact D = 12138.702497....

Example2: λ = 0° , β = -16°
λ' = +177° , β' = +10°

0 ENTER^
     -16 ENTER^
   177 ENTER^
     10 R/S >>>>   D = 19266.90923 km

-Exact D = 19266.909329....
-In this example, the great elliptic arc length = 19267.285727...

Remark:

-If we replace lines 211 to 216 with:

RCL 00 ST* Y ST* Z ST+ Z + 6378.137 ST* Z * END

Y-output = Geodesic distance on the ellipsoid of revolution ( WGS84 )

-For instance:

     166.08022 ENTER^
    -33.41009 ENTER^
   -101.56060    ENTER^
     33.10470 XEQ "TGC"   >>>>   D = 12138.70249 km
X<>Y d = 12138.68505 km ( WGS84 )

Notes:

-If | λ - λ' | = 160° , the maximum relative error is about 2.24 10^-8

-If | λ - λ' | = 80° , the root-mean-square relative error is ~ 0.32 10^-8
-If | λ - λ' | = 120° , the root-mean-square relative error is ~ 0.37 10^-8
-If | λ - λ' | = 140° , the root-mean-square relative error is ~ 0.47 10^-8
-If | λ - λ' | = 150° , the root-mean-square relative error is ~ 0.52 10^-8
-If | λ - λ' | = 160° , the root-mean-square relative error is ~ 0.60 10^-8
-If | λ - λ' | = 168° , the root-mean-square relative error is ~ 0.82 10^-8
-If | λ - λ' | = 170° , the root-mean-square relative error is ~ 0.96 10^-8

-Here are a few maximum errors ( in centimeters ):

| λ - λ' | = 0°

( 0° , 50° ) ( 0° , -90° ) m = -7.5 cm
( 90° , 79° ) ( 90° , 42° ) M = +3.7 cm

| λ - λ' | = 120°

( -45°, +20° ) ( +75° , -31° ) m = -13.4 cm
( +30°, +21° ) ( +150° , +20° ) M = +11.5 cm

| λ - λ' | = 140°

( -55°, +25° ) ( +85° , -32° ) m = -21.3 cm
( +20°, +17° ) ( +160° , +17° ) M = +19.6 cm

| λ - λ' | = 150°

( -60°, +28° ) ( +90° , -32° ) m = -26.3 cm
( 0°, +32° ) ( +150° , -27° ) M = +24.8 cm

| λ - λ' | = 160°

( -80°, +62° ) ( +80° , -54° ) M = +40.5 cm
( 10°, +55° ) ( +170° , -63° ) m = -41.5 cm

| λ - λ' | = 168°

( -69°, +58° ) ( +99° , -54° ) M = +77.5 cm
( 6° , +59° ) ( +174° , -54° ) m = -77.9 cm

| λ - λ' | = 170°

( -70° , +54° ) ( +100° , -50° ) M = +96.5 cm
( -10° , +50° ) ( +160° , -54° ) m = -98.6 cm

e) Program#5 ( With Great Elliptic Arc Length )

-In most cases, the following program is even more accurate !
-This method is not completely a generalization of Andoyer's formula:

-It first uses the formula listed in "Great Elliptic Arc Length for the HP41" paragraph 2°)b):

   m = [ y² (a²/b² - 1)/2 + z² (a²/c² - 1)/2 ] / ( Sin² µ )
With n = [ y'² (a²/b² - 1)/2 + z'² (a²/c² - 1)/2 ] / ( Sin² µ )
p = 2 [ y.y' (a²/b² - 1)/2 + z.z' (a²/c² - 1)/2 ] / ( Sin² µ )

R/a = 1 / SQRT [ 1 + m Sin² (µ-t) + n Sin² t + p Sin t Sin (µ-t) ] = 1 / SQRT ( 1 + q ) where q is a small number

-After some calculations, it yields:

   1 + q = ( 1 + m Sin² µ ) [ 1 + e Sin² t + g Sin t Cos t ] = ( 1 + m Sin² µ ) ( 1 + e/2 ) [ 1 + E Cos 2.t + F Sin 2.t ] = ( 1 + m Sin² µ ) ( 1 + e/2 ) [ 1 + G Sin ( 2.t + D ) ]

= ( 1 + m Sin² µ ) ( 1 + e/2 ) ( 1 + h )

-The formula of order 3 to compute the arc length = §₀^µ [ R² + ( dR/dt )² ]^1/2 dt is

  [ R² + ( dR/dt )² ]^1/2 ~ a ( 1 + m Sin² µ )^-1/2 (1 + e/2 )^-1/2 [ 1 - h/2 + (3/8) h² - (5/16) h³ + (1/8) h'² - (5/16) h h'² ] where h' = dh/dt

-And the great elliptic arc length is:

  Δ / a   = ( 1 + m Sin² µ )^-1/2 ( 1 + e/2 )^-1/2 { µ + ( G / 4 )    [ - Cos D + Cos ( 2.µ + D ) ] + ( G² / 32 ) [ 14 µ - ( Sin D ) ( Cos D ) + Sin ( 2.µ + D ) Cos ( 2.µ + D ) ]
+ ( G³ / 32 ) [ -5 Cos D + 5 Cos ( 2.µ + D ) - 5 Cos³ D + 5 Cos³ ( 2.µ + D ) ] }

-Then, the difference between geodesic distance & great elliptic arc lentgh is given in "Terrestrial Geodesic Distance(3)"

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

with a' = (a+b)/2 + (a-b)/2 Cos (λ+λ')

-The angular separation µ is computed with Sin µ & Cos µ

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

01 LBL "TGC"
02 DEG
03 HR
04 STO 02
05 RDN
06 HR
07 STO 06
08 STO 15
09 X<> Z
10 HR
11 ST+ 15
12 X<>Y
13 HR
14 1
15 STO 07
16 P-R
17 X<>Y
18 STO 14
19 RDN
20 P-R
21 STO 11
22 X<>Y
23 STO 13
24 RCL 06
25 RCL 02
26 RCL 07
27 P-R
28 X<>Y
29 STO 12
30 RDN
31 P-R
32 STO 09
33 RCL 11
34 *
35 X<>Y
36 STO 10
37 RCL 13
38 *
39 STO 05
40 +
41 RCL 12
42 RCL 14
43 *
44 STO 00
45 +
46 STO 04

47 RCL 10
48 X^2
49 2195 E-8
50 STO 02
51 ST* 05
52 *
53 6750546 E-9
54 STO 03
55 ST* 00
56 RCL 12
57 X^2
58 *
59 +
60 STO 01
61 RCL 13
62 X^2
63 RCL 02
64 *
65 RCL 14
66 X^2
67 RCL 03
68 *
69 +
70 STO 06
71 RCL 10
72 RCL 11
73 *
74 RCL 09
75 RCL 13
76 *
77 -
78 X^2
79 RCL 12
80 RCL 13
81 *
82 RCL 10
83 RCL 14
84 *
85 -
86 X^2
87 +
88 RCL 11
89 RCL 12
90 *
91 RCL 09
92 RCL 14

93 *
94 -
95 X^2
96 +
97 STO 09
98 SQRT
99 STO 11
100 ST/ 12
101 ST/ 14
102 RCL 04
103 R-P
104 X<>Y
105 STO 08
106 RCL 00
107 RCL 05
108 +
109 ST+ X
110 RCL 09
111 ST/ 01
112 ST/ 06
113 /
114 STO 02
115 RCL 04
116 *
117 RCL 01
118 -
119 RCL 08
120 ST+ X
121 STO 05
122 COS
123 RCL 06
124 *
125 -
126 RCL 02
127 RCL 04
128 ST+ X
129 RCL 06
130 *
131 -
132 RCL 06
133 RCL 11
134 ST* Z
135 X^2
136 *
137 RCL 07
138 +

139 STO 03
140 ST+ X
141 ST/ Z
142 /
143 RCL 07
144 RCL Z
145 -
146 ST* 03
147 ST/ Z
148 /
149 R-P
150 STO 01
151 X<>Y
152 STO 02
153 RCL 07
154 P-R
155 STO 00
156 X<>Y
157 STO 06
158 RCL 02
159 RCL 05
160 +
161 RCL 07
162 P-R
163 STO 10
164 ST* Y
165 X^2
166 LASTX
167 ST* Y
168 +
169 RCL 00
170 ST* 06
171 ST- 10
172 X^2
173 LASTX
174 ST* Y
175 +
176 -
177 RCL 01
178 STO T
179 *
180 5
181 *
182 +
183 RCL 06
184 -

185 RCL 08
186 D-R
187 STO 00
188 14
189 *
190 +
191 *
192 8
193 /
194 RCL 10
195 +
196 *
197 4
198 /
199 RCL 00
200 +
201 RCL 03
202 SQRT
203 /
204 6378.172
205 STO 16
206 *
207 STO 10
208 RCL 12
209 X^2
210 RCL 14
211 ST* 12
212 X^2
213 +
214 STO 05
215 RCL 04
216 ST* 05
217 STO 08
218 STO 13
219 RCL 12
220 ST+ X
221 ST- 05
222 *
223 -
224 ST- 07
225 ST* 13
226 STO 14
227 RCL 00
228 RCL 11
229 /

230 STO 06
231 ST/ 08
232 STO 12
233 RCL 05
234 RCL 13
235 +
236 ST* 12
237 RCL 04
238 RCL 07
239 ST+ X
240 *
241 -
242 *
243 RCL 14
244 13
245 *
246 -
247 6
248 +
249 RCL 12
250 *
251 5
252 ST* 13
253 RCL 08
254 6
255 *
256 -
257 RCL 05
258 *
259 RCL 07
260 4
261 *
262 RCL 06
263 /
264 -
265 RCL 13
266 +
267 RCL 05
268 *
269 +
270 RCL 00
271 ST+ X
272 X^2
273 RCL 07
274 ST* 00

275 *
276 3
277 ST* 08
278 /
279 RCL 08
280 7
281 +
282 RCL 14
283 *
284 -
285 4
286 +
287 RCL 14
288 ST- 12
289 *
290 -
291 RCL 15
292 COS
293 .035
294 ST- 16
295 *
296 RCL 16
297 +
298 STO 01
299 6356.752314
300 -
301 RCL 01
302 ST* 00
303 ST+ X
304 /
305 STO Z
306 *
307 RCL 05
308 RCL 06
309 /
310 -
311 RCL 12
312 +
313 *
314 *
315 RCL 00
316 *
317 RCL 10
318 +
319 END

( 407 bytes / SIZE 017 )

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

Where D is the geodesic distance on the triaxial ellipsoid & d =R10 = great elliptic arc length

Example: λ = +166°08'02"2 , β = -33°41'00"9
λ' = -101°56'06"0 , β' = +33°10'47"0

     166.08022 ENTER^
    -33.41009 ENTER^
   -101.56060    ENTER^
     33.10470 XEQ "TGC"   >>>>   D = 12138.70250 km                                                      ---Execution time = 18.7s---
      LASTX = R10 = d = 12138.70371 km

-Exact D = 12138.702497....

  Notes:

-With free42 , replace line 49 by 2195022 E-11 and line 53 by 675054581 E-11 ( exact values are: a²/b²-1 & a²/c²-1 )

-With the great elliptic arc length, maximum error is about 1 millimeter
-But with the geodesic distance - except if λ - λ' = 0 - errors may be much larger:

-If | λ - λ' | = 160° , the maximum relative error is about 2.31 10^-8 and the root-mean-square relative error is ~ 0.53 10^-8
-If | λ - λ' | = 120° , the root-mean-square relative error is ~ 0.22 10^-8

-Here are a few maximum errors ( in centimeters ):

| λ - λ' | = 0°

maximum error ~ 1 mm

| λ - λ' | = 120°

( 0°, +43° ) ( +120° , -60° ) m = -10.4 cm
( -30°, +68° ) ( +90° , -51° ) M = +10.4 cm

| λ - λ' | = 160°

( -65°, +58° ) ( +95° , -50° ) M = +42.4 cm
( -5°, +51° ) ( +155° , -58° ) m = -42.8 cm

| λ - λ' | = 168°

( -69°, +54° ) ( +99° , -49° ) M = +75.1 cm
( -9° , +54° ) ( +159° , -58° ) m = -76.7 cm

f) My favorite Programs

f1) Program#1

Formula:

   Dist/(a' µ) ~ 1 + (f/2) ( -A + B (Sin µ)/µ )
                + (f²/4) [ A ( -4 + µ Cot µ ) + B ( 2 (Sin µ)/µ + µ / Sin µ ) + (A²/4) ( 17 + (Sin µ) (Cos µ)/µ - 4 µ Cot µ ) - A.B ( 2 (Sin µ)/µ + µ / Sin µ ) - ( B²/2 ) (Sin µ) (Cos µ)/µ ]
                  + (f³/8) { ( 1 - A ) [ ( -A ( 1 - A ) (4/3) µ² + A ( 21.A - 13 ) + ( B + A Cos µ ) [ ( µ / Sin µ ) ( 6 - 13.A) + 5.B + ( µ² / Sin² µ ) ( B + ( 3 A - 2 ) Cos µ ) ) ]
      + ( 3.A² - 9.B² ) ( Sin µ ) ( Cos µ ) / µ ] + B³ [ ( Sin µ ) / µ + 2 ( Sin³ µ ) / µ ] }
   + (f⁴/16) ( 1 - A ) ( µ³ / Sin³ µ ) ( B + A Cos µ ) [ 2.B + ( 4.A - 2 ) Cos µ ]

+   (f'/2) [ -A"+ B" (Sin µ)/µ ]
+ ( f f'/4 ) (1-A) [ µ² ( -2 AB" - 4 AA" ) + 5 BB" ]
+ ( f ² f'/8 ) (1-A) ( µ² / Sin² µ ) ( A² + B² ) ( -28.A" - 24.B" ) f/2 ~ 1/596.51444 & f'/2 ~ 1/364465 & a' = 6378.137 km

     Sin² µ = Sin² ( λ - λ' ) Cos² β' + [ Cos β Sin β' - Cos β' Sin β Cos ( λ - λ' ) ]²
      Cos µ = Sin β Sin β' + Cos β Cos β' Cos ( λ - λ' )

A = [ Sin² β + Sin² β' - 2 Sin β Sin β' Cos µ ] / Sin² µ
B = [ ( Sin² β + Sin² β' ) Cos µ - 2 Sin β Sin β' ] / Sin² µ

   -A" = 2 ( 1 - A' ) - A
   -B" = -2 B' - B with:

   A' = ( cos² β Sin² λ + cos² β' Sin² λ' - 2 cos β Sin λ cos β' Sin λ' Cos µ ) / Sin² µ
   B' = [ ( cos² β Sin² λ + cos² β' Sin² λ' ) Cos µ - 2 cos β Sin λ cos β' Sin λ' ] / Sin² µ

Data Registers: R00 thru R14: temp R00 = µ R01 = -A + B ( Sin µ ) / µ R02 = - B ( Sin µ ) / µ R08 = -A" R09 = -B" ( Sin µ ) / µ R12 = 5.B.B" ( 1 - A )
Flags: /
Subroutines: /

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

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

69 RCL 03
70 ST- 08
71 *
72 -
73 ST+ X
74 X<> 09
75 STO 03
76 R-P
77 RCL 08
78 +
79 ST+ X
80 STO 08
81 X<>Y
82 D-R
83 STO 00
84 R^
85 /
86 STO 06
87 RCL 01
88 RCL 03
89 ST* 01
90 RCL 05
91 ST- 01
92 *
93 -
94 STO 02
95 ST+ 09
96 RCL 01
97 ST+ 04
98 ST+ 08
99 RCL 03
100 *
101 +
102 STO 05

103 LASTX
104 RCL 03
105 +
106 ST+ X
107 +
108 STO 07
109 RCL 05
110 +
111 *
112 596.51444
113 STO 10
114 /
115 RCL 07
116 +
117 *
118 RCL 01
119 13
120 *
121 STO T
122 -
123 6
124 -
125 *
126 RCL 02
127 5
128 *
129 STO 12
130 +
131 RCL 05
132 *
133 +
134 RCL 01
135 X^2
136 STO 07

137 STO 14
138 21
139 *
140 +
141 RCL 07
142 RCL 02
143 X^2
144 STO 11
145 ST* 13
146 ST+ 13
147 ST+ 14
148 ST+ X
149 -
150 RCL 03
151 RCL 06
152 ST/ 02
153 ST* 05
154 /
155 ST* 11
156 *
157 ST+ 07
158 RCL 11
159 -
160 RCL 04
161 4
162 ST/ 07
163 *
164 ST* 14
165 RCL 01
166 *
167 ST+ 07
168 RCL 00
169 X^2
170 *

171 STO 11
172 3
173 ST* Z
174 /
175 +
176 +
177 RCL 10
178 ST/ 13
179 /
180 RCL 05
181 RCL 02
182 ST- 01
183 ST* 13
184 ST+ X
185 +
186 -
187 RCL 04
188 ST* 12
189 *
190 RCL 07
191 +
192 RCL 13
193 -
194 RCL 10
195 ST/ 14
196 /
197 RCL 01
198 +
199 RCL 08
200 RCL 09
201 ST* 12
202 2
203 /
204 +

205 ST* 11
206 LASTX
207 +
208 6
209 *
210 RCL 08
211 +
212 RCL 14
213 *
214 RCL 06
215 ST/ 09
216 X^2
217 *
218 RCL 12
219 +
220 RCL 11
221 -
222 RCL 10
223 ST/ Z
224 /
225 RCL 08
226 RCL 09
227 -
228 +
229 364465
230 /
231 +
232 RCL 00
233 ST* Y
234 +
235 6378.137
236 *
237 END

( 319 bytes / SIZE 015 )

STACK	INPUTS	OUTPUTS
T	λ ( ° ' " )	/
Z	β ( ° ' " )	/
Y	λ' ( ° ' " )	/
X	β' ( ° ' " )	D ( km )

Where D is the geodesic distance on the triaxial ellipsoid

Example1: λ = +166°08'02"2 , β = -33°41'00"9
λ' = -101°56'06"0 , β' = +33°10'47"0

     166.08022 ENTER^
    -33.41009 ENTER^
   -101.56060    ENTER^
     33.10470 XEQ "TGC"   >>>>   D = 12138.70249 km                                                      ---Execution time = 13.3s---

-Exact D = 12138.702497....

Example2: λ = -16° , β = -16°
λ' = +157° , β' = +13°

     -16 ENTER^ ENTER^
   157 ENTER^
     13   R/S >>>>   D = 19205.13652 km

-Exact D = 19205.13653....
-In this example, the great elliptic arc length = 19205.463988...

Remarks:

-If we replace lines 232 to 237 with:

RCL 00 ST* Y ST* Z ST+ Z + 6378.137 ST* Z * END

Y-output = Geodesic distance on the ellipsoid of revolution ( WGS84 )

-For instance:

     166.08022 ENTER^
    -33.41009 ENTER^
   -101.56060    ENTER^
     33.10470 XEQ "TGC"   >>>>   D = 12138.70249 km
X<>Y d = 12138.68505 km ( WGS84 )

-If you prefer:

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

replace lines 112-229-235 with 596.51284 364101 6378.13661

-The result in example 1 becomes: 12138.70176 km

Notes:

-If | λ - λ' | = 160° , the maximum relative error is about 2 10^-8

-If | λ - λ' | = 120° , the root-mean-square relative error is ~ 0.37 10^-8
-If | λ - λ' | = 160° , the root-mean-square relative error is ~ 0.59 10^-8
-If | λ - λ' | = 168° , the root-mean-square relative error is ~ 0.75 10^-8
-If | λ - λ' | = 170° , the root-mean-square relative error is ~ 0.84 10^-8

-Here are a few maximum errors ( in centimeters ):

| λ - λ' | = 0°

( 0° , 50° ) ( 0° , -90° ) m = -7.5 cm
( 90° , 79° ) ( 90° , 42° ) M = +3.7 cm

| λ - λ' | = 90°

( +30°, +43° ) ( +120°, -78° ) m = -8.6 cm
( -30°, +64° ) ( +60° , -16° ) M = +7.2 cm

| λ - λ' | = 120°

( -45°, +21° ) ( +75° , -32° ) m = -13.7 cm
( +30°, +21° ) ( +150° , +20° ) M = +11.9 cm

| λ - λ' | = 160°

( -80°, +34° ) ( +80° , -34° ) m= -37.1 cm
( 10°, +35° ) ( +170° , -35° ) M = +36.2 cm

| λ - λ' | = 168°

( -84°, +59° ) ( +84° , -53° ) M = +62.2 cm
( 6° , +58° ) ( +174° , -52° ) m = -61.9 cm

| λ - λ' | = 170°

( +5° , +52° ) ( +175° , -57° ) m = -73.9 cm
( -85° , +53° ) ( +85° , -58° ) M = +74.7 cm

-If | λ - λ' | = 172° , the error can reach 104 cm ( perhaps a little more... )

-But the error may be relatively small.
-For instance, if λ = β = 0 , λ' = 172° here are a few errors ( in millimeters )

β'	0°	1°	2°	3°	4°	5°	6°	7°	8°	9°	10°	11°	12°	13°	14°	15°	20°	30°	40°	50°	60°	70°	80°	90°
error	~0	46	37	-182	-484	-646	-591	-390	-147	66	221	316	363	376	368	347	203	47	2	-12	-21	-39	-37	-35

-If we compare with the program listed in §4°)d), the precision is ~ the same if | λ - λ' | <= 120°
-But if | λ - λ' | = 160° , 170° ... we have a better accuracy... even better than the program in § 4°)e) !

f2) Program(s)#2

-Unfortunately with the program above, if µ is very small, there are sometimes very large roundoff-errors.

-For instance,

if   λ = 56° , β = 26°    λ' = 56°00'00"36 ,  β' = 26°00'00"36 "TGC" above gives D = 0.014967941 km but exact distance = 0.014957847... km
if   λ = 56° , β = 26°    λ' = 56°00'00"0036 ,  β' = 26°00'00"0036 it yields: D = 909806.6455 km (!!!) completely wrong.

-If we want to calculate very small distances, we must compute µ A A' B B' with other formulae:

   H = (λ - λ')/2 I = (λ + λ')/2 F = (β + β ')/2 G = (β - β ')/2

S = Sin² (µ/2) = Sin² G Cos² H + Cos² F Sin² H
C = Cos² (µ/2) = Cos² G Cos² H + Sin² F Sin² H

A = ( Sin² G Cos² F ) / S + ( Cos² G Sin² F ) / C
B = ( Sin² G Cos² F ) / S - ( Cos² G Sin² F ) / C

-These formulae to compute S C A B may be found here: Jean Meeus - "Astronomical Algorithms" - Willmann-Bell - ISBN 0-943396-35-2

-After some calculations, I've found the following formulae to compute A' & B'

A' = ( Cos I Sin H Cos F Cos G - Sin I Cos H Sin F Sin G )² / S + ( Sin I Cos H Cos F Cos G - Cos I Sin H Sin F Sin G )² / C
B' = ( Cos I Sin H Cos F Cos G - Sin I Cos H Sin F Sin G )² / S - ( Sin I Cos H Cos F Cos G - Cos I Sin H Sin F Sin G )² / C

   A" = 2 ( A' - 1 ) + A and B" = 2 B' + B

-Of course, the program becomes longer:

Data Registers: R00 thru R15: temp R00 = µ R01 = A - B ( Sin µ ) / µ R02 = B ( Sin µ ) / µ R08 = A" R09 = B" ( Sin µ ) / µ R12 = 5.B.B" ( 1 - A )
Flags: /
Subroutines: /

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

38 STO 07
39 RCL 11
40 RCL 04
41 P-R
42 ST* 00
43 ST* 07
44 ST* 10
45 RDN
46 STO 09
47 *
48 ST* 12
49 X^2
50 +
51 X<> 10
52 X^2
53 RCL 06
54 X^2
55 RCL 07
56 X^2
57 +
58 STO 08
59 /
60 STO 02
61 RCL 03
62 ST* 07
63 RCL 09
64 ST* 06
65 *
66 X^2
67 RCL 10
68 /
69 ST- 02
70 +
71 STO 01
72 RCL 14
73 RCL 04
74 P-R

75 ST* 12
76 RCL 07
77 *
78 X<>Y
79 ST* 00
80 RCL 06
81 *
82 -
83 X^2
84 RCL 08
85 /
86 STO 09
87 RCL 00
88 RCL 12
89 -
90 X^2
91 RCL 10
92 /
93 ST- 09
94 +
95 ST+ X
96 RCL 15
97 ST* 09
98 -
99 +
100 X<> 08
101 SQRT
102 RCL 10
103 SQRT
104 R-P
105 RDN
106 ST+ X
107 D-R
108 STO 00
109 LASTX
110 RCL 04
111 P-R

112 STO 03
113 ST* T
114 RDN
115 ST/ Y
116 X^2
117 ST+ X
118 STO 13
119 X<>Y
120 STO 06
121 RCL 02
122 ST+ 09
123 R^
124 +
125 STO 05
126 LASTX
127 RCL 03
128 -
129 ST+ X
130 +
131 STO 07
132 RCL 05
133 +
134 *
135 596.51444
136 STO 10
137 /
138 RCL 07
139 +
140 *
141 RCL 01
142 ST- 04
143 13
144 *
145 STO 07
146 -
147 6
148 +

149 *
150 RCL 02
151 5
152 *
153 STO 12
154 +
155 RCL 05
156 *
157 RCL 07
158 -
159 RCL 01
160 X^2
161 STO 07
162 STO 14
163 21
164 *
165 +
166 RCL 07
167 RCL 02
168 X^2
169 STO 11
170 ST* 13
171 ST+ 13
172 ST+ 14
173 ST+ X
174 -
175 RCL 03
176 RCL 06
177 ST/ 02
178 ST* 05
179 /
180 ST* 11
181 *
182 ST+ 07
183 RCL 11
184 -
185 RCL 04
186 4

187 ST/ 07
188 *
189 ST* 14
190 RCL 01
191 *
192 ST- 07
193 RCL 00
194 X^2
195 *
196 STO 11
197 3
198 ST* Z
199 /
200 -
201 +
202 RCL 10
203 ST/ 13
204 /
205 RCL 05
206 RCL 02
207 ST- 01
208 ST* 13
209 ST+ X
210 +
211 +
212 RCL 04
213 ST* 12
214 *
215 RCL 07
216 +
217 RCL 13
218 +
219 RCL 10
220 ST/ 14
221 /
222 RCL 01
223 -
224 RCL 08

225 RCL 09
226 ST* 12
227 RCL 15
228 /
229 +
230 ST* 11
231 LASTX
232 +
233 6
234 *
235 RCL 08
236 +
237 RCL 14
238 *
239 RCL 06
240 ST/ 09
241 X^2
242 *
243 RCL 12
244 -
245 RCL 11
246 +
247 RCL 10
248 ST/ Z
249 /
250 RCL 08
251 RCL 09
252 -
253 +
254 364465
255 /
256 -
257 RCL 00
258 ST* Y
259 +
260 6378.137
261 *
262 END

( 340 bytes / SIZE 016 )

STACK	INPUTS	OUTPUTS
T	λ ( ° ' " )	/
Z	β ( ° ' " )	/
Y	λ' ( ° ' " )	/
X	β' ( ° ' " )	D ( km )

Where D is the geodesic distance on the triaxial ellipsoid

Example1: λ = +166°08'02"2 , β = -33°41'00"9
λ' = -101°56'06"0 , β' = +33°10'47"0

     166.08022 ENTER^
    -33.41009 ENTER^
   -101.56060    ENTER^
     33.10470 XEQ "TGC"   >>>>   D = 12138.70249 km                                                      ---Execution time = 14.2s---

-Exact D = 12138.702497....

Example2: λ = -16° , β = -16°
λ' = +157° , β' = +13°

     -16 ENTER^ ENTER^
   157 ENTER^
     13   R/S >>>>   D = 19205.13652 km

-Exact D = 19205.13653....
-In this example, the great elliptic arc length = 19205.463988...

Example3:      λ = 56° , β = 26°    λ' = 56°00'00"36 ,  β' = 26°00'00"36

   56    ENTER^
26    ENTER^
56.000036 ENTER^
26.000036 R/S >>>> D = 0.014957847 km

-Exact distance = 0.014957847... km

Example4:      λ = 56° , β = 26°    λ' = 56°00'00"0036 ,  β' = 26°00'00"0036

      56       ENTER^
26 ENTER^
56.00000036 ENTER^
26.00000036 R/S >>>> D = 0.000149578497 km

-Exact distance = 0.000149578497... km

-If the coordinates are in decimal degrees, we can save a few bytes.
-Here is a version with

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

Data Registers: R00 thru R14: temp R00 = µ R01 = A - B ( Sin µ ) / µ R02 = B ( Sin µ ) / µ R08 = A" R09 = B" ( Sin µ ) / µ R12 = 5.B.B" ( 1 - A )
Flags: /
Subroutines: /

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

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

73 RCL 06
74 *
75 +
76 X^2
77 RCL 08
78 /
79 STO 09
80 RCL 00
81 RCL 12
82 +
83 X^2
84 RCL 10
85 /
86 ST- 09
87 +
88 ST+ X
89 RCL 14
90 ST* 09
91 -
92 +
93 X<> 08
94 SQRT
95 RCL 10
96 SQRT
97 R-P
98 RDN
99 ST+ X
100 D-R
101 STO 00
102 LASTX
103 RCL 04
104 P-R
105 STO 03
106 R^
107 ST- 04
108 *

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

145 +
146 RCL 05
147 *
148 RCL 07
149 -
150 RCL 01
151 X^2
152 STO 07
153 STO 11
154 21
155 *
156 +
157 RCL 07
158 RCL 02
159 X^2
160 STO 06
161 ST* 13
162 ST+ 11
163 ST+ 13
164 ST+ X
165 -
166 RCL 03
167 R^
168 STO 03
169 ST/ 02
170 ST* 05
171 /
172 ST* 06
173 *
174 ST+ 07
175 RCL 06
176 -
177 RCL 04
178 4
179 ST/ 07
180 *

181 ST* 11
182 RCL 01
183 *
184 ST- 07
185 RCL 00
186 X^2
187 *
188 STO 06
189 3
190 ST* Z
191 /
192 -
193 +
194 RCL 10
195 ST/ 13
196 /
197 RCL 05
198 RCL 02
199 ST- 01
200 ST* 13
201 ST+ X
202 +
203 +
204 RCL 04
205 ST* 12
206 *
207 RCL 07
208 +
209 RCL 13
210 +
211 RCL 10
212 ST/ 11
213 /
214 RCL 01
215 -
216 RCL 08
217 RCL 09

218 ST* 12
219 RCL 14
220 /
221 +
222 ST* 06
223 LASTX
224 +
225 6
226 *
227 RCL 08
228 +
229 RCL 11
230 *
231 RCL 03
232 ST/ 09
233 X^2
234 *
235 RCL 12
236 -
237 RCL 06
238 +
239 RCL 10
240 ST/ Z
241 /
242 RCL 08
243 RCL 09
244 -
245 +
246 364101
247 /
248 -
249 RCL 00
250 ST* Y
251 +
252 6378.13661
253 *
254 END

( 336 bytes / SIZE 015 )

STACK	INPUTS	OUTPUTS
T	λ ( ° )	/
Z	β ( ° )	/
Y	λ' ( ° )	/
X	β' ( ° )	D ( km )

Where D is the geodesic distance on the triaxial ellipsoid

Example1: λ = +166°1339444 , β = -33°68358333
λ' = -101°935 , β' = +33°17972222

     166.1339444 ENTER^
    -33.68358333 ENTER^
   -101.9350000    ENTER^
     33.17972222 XEQ "TGC"   >>>>   D = 12138.70176 km                                                      ---Execution time = 14.1s---

Example2: λ = -16° , β = -16°
λ' = +157° , β' = +13°

     -16 ENTER^ ENTER^
   157 ENTER^
     13   R/S >>>>   D = 19205.13532 km

Remark:

-If we replace lines 249 to 254 with:

RCL 00 ST* Y ST* Z ST+ Z + 6378.13661 ST* Z * END

Y-output = Geodesic distance on the ellipsoid of revolution

-For instance:

     166.1339444 ENTER^
    -33.68358333 ENTER^
   -101.9350000    ENTER^
     33.17972222 XEQ "TGC"   >>>>   D = 12138.70176 km
X<>Y d = 12138.68430 km

-With the golden search to minimize the sum of 2 geodesic distances, the maximum error becomes much smaller:

01 LBL "DDD"
02 STO 19
03 RDN
04 STO 18
05 RDN
06 STO 17
07 X<>Y
08 STO 16
09 R^
10 X<>Y
11 -
12 R^
13 COS
14 P-R
15 RCL Z
16 COS
17 +
18 R-P
19 CLX

20 RCL 16
21 +
22 STO 20
23 GTO 00
24 LBL 12
25 STO 28
26 RCL 20
27 X<>Y
28 RCL 18
29 RCL 19
30 XEQ "TGC"
31 STO 24
32 RCL 16
33 RCL 17
34 RCL 20
35 RCL 28
36 XEQ "TGC"
37 RCL 24
38 +

39 RTN
40 LBL 00
41 90
42 ENTER
43 STO 22
44 CHS
45 STO 21
46 -
47 5
48 SQRT
49 1
50 +
51 2
52 /
53 STO 27
54 /
55 RCL 21
56 +
57 STO 25

58 XEQ 12
59 STO 26
60 LBL 10
61 RCL 25
62 RCL 21
63 -
64 RCL 27
65 /
66 RCL 21
67 +
68 STO 23
69 XEQ 12
70 STO 24
71 VIEW 24
72 RCL 26
73 X<Y?
74 GTO 02
75 RCL 23
76 X<> 25
77 STO 22

78 RCL 24
79 STO 26
80 GTO 04
81 LBL 02
82 RCL 22
83 STO 21
84 RCL 23
85 STO 22
86 LBL 04
87 E-4
88 RCL 21
89 RCL 23
90 -
91 ABS
92 X>Y?
93 GTO 10
94 RCL 23
95 RCL 24
96 CLD
97 END

( 168 bytes / SIZE 029 )

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

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

Example: With (L,b) = (0°,0°) & (L',b') = (179.85°,0°)

        0         ENTER^ ENTER^
    179.85    ENTER^
        0         XEQ "DDD" >>>>   D = 20001.90494 km
                                     X<>Y b" ~ -75.52°

-Exact distance = 20001.904992...

Notes:

-Very slow with HP41 but only a few seconds with V41.
-Here is another method to minimize the sum of 2 geodesic distances:

01 LBL "DDD"
02 STO 19
03 RDN
04 STO 18
05 RDN
06 STO 17
07 RDN
08 STO 16
09 RDN
10 X<> T
11 -
12 R^
13 COS
14 P-R
15 X<>Y
16 X<> Z
17 COS
18 ST+ Y
19 X<> L
20 SIN
21 R^
22 SIN
23 +

24 X<> Z
25 X<>Y
26 R-P
27 RCL Z
28 X<>Y
29 R-P
30 X<> Z
31 RCL 16
32 +
33 X<>Y
34 STO 21
35 X<>Y
36 STO 20
37 10
38 STO 28
39 RCL 21
40 STO 26
41 XEQ 10
42 STO 27
43 LBL 01
44 RCL 26
45 RCL 28
46 +

47 XEQ 10
48 STO 29
49 LBL 02
50 RCL 26
51 RCL 28
52 -
53 XEQ 10
54 STO 30
55 LBL 03
56 VIEW 27
57 RCL 29
58 X>Y?
59 X<>Y
60 RCL 27
61 X>Y?
62 X<>Y
63 RCL 27
64 X=Y?
65 GTO 05
66 CLX
67 RCL 29
68 X=Y?
69 GTO 04

70 RCL 30
71 X<> 27
72 STO 29
73 RCL 28
74 ST- 26
75 GTO 02
76 LBL 04
77 X<> 27
78 STO 30
79 RCL 28
80 ST+ 26
81 RCL 26
82 +
83 XEQ 10
84 STO 29
85 RCL 30
86 GTO 03
87 LBL 05
88 3
89 ST/ 28
90 RCL 28
91 ABS
92 E-4

93 X<Y?
94 GTO 01
95 GTO 11
96 LBL 10
97 STO 24
98 RCL 16
99 RCL 17
100 RCL 20
101 RCL 24
102 XEQ "TGC"
103 STO 25
104 RCL 20
105 RCL 24
106 RCL 18
107 RCL 19
108 XEQ "TGC"
109 RCL 25
110 +
111 RTN
112 LBL 11
113 RCL 27
114 CLD
115 END

( 203 bytes / SIZE 031 )

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

Where D is the geodesic distance on the triaxial ellipsoid

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

        0         ENTER^ ENTER^
    179.85    ENTER^
        0         XEQ "DDD" >>>>   D = 20001.90494 km

-Exact distance = 20001.904992...
-The error becomes ~ 5 cm instead of 19 km

Example2: λ = 5° , β = 52°
λ' = 175° , β' = -57°

     5 ENTER^
52 ENTER^
   175 ENTER^
    -57 R/S >>>>   D = 19160.59422 km ---Execution time = 11m21s---

-Exact distance = 19160.5941938...
-The error becomes ~ 3 cm instead of 74 cm

5°) Transformation of Coordinates

a) Geocentric Coordinates <> Geodetic Coordinates

-The formulae between the geodetic longitude L & geodetic latitude B and the cartesian coordinates x , y , z of a point on a triaxial ellipsoid are:

    x = k a² Cos B Cos L
    y = k b² Cos B Sin L
    z = k c² Sin B

where k = ( a² Cos² B Cos² L + b² Cos² B Sin² L + c² Sin² B )^-1/2

Data Registers: /
Flags: /
Subroutines: /

01 LBL "LCLD"
02 DEG
03 HR
04 1
05 P-R
06 X<>Y

07 6356752.314
08 X^2
09 /
10 X<> Z
11 HR
12 X<>Y

13 P-R
14 6378102
15 X^2
16 ST/ Z
17 CLX
18 6378172

19 X^2
20 /
21 R-P
22 X<>Y
23 RDN
24 R-P

25 X<>Y
26 HMS
27 R^
28 HMS
29 END

( 62 bytes / SIZE 000 )

STACK	INPUTS	OUTPUTS
Y	λ ( ° ' " )	β' ( ° ' " )
X	β ( ° ' " )	λ' ( ° ' " )

Example: Geocentric longitude & latitude: λ = 124°41' β = 57°16'

124.41 ENTER^
57.16 XEQ "LCLD" >>>> λ' = 124°40'57"881 = geodetic longitude ---Execution time = 4s---
X<>Y β' = 57°26'28"9454 = geodetic latitude

01 LBL "LDLC"
02 DEG
03 HR
04 1
05 P-R
06 X<>Y

07 6356752.314
08 X^2
09 *
10 X<> Z
11 HR
12 X<>Y

13 P-R
14 6378172
15 X^2
16 *
17 X<>Y
18 6378102

19 X^2
20 *
21 X<>Y
22 R-P
23 X<>Y
24 HMS

25 RDN
26 R-P
27 X<>Y
28 HMS
29 R^
30 END

( 62 bytes / SIZE 000 )

STACK	INPUTS	OUTPUTS
Y	λ ( ° ' " )	β' ( ° ' " )
X	β ( ° ' " )	λ' ( ° ' " )

Example: Geodetic longitude & latitude: λ = 124°40'57"881 β = 57°26'28"9454

124.4057881 ENTER^
57.26289454 XEQ "LDLC" >>>> λ' = 124°40'59"999 = geocentric longitude ---Execution time = 4s---
X<>Y β' = 57°16 = geocentric latitude

b) WGS84 Coordinates -> Triaxial Ellipsoid Geocentric Coordinates

-We simply add 14°92911 to the WGS84 longitude λ to obtain λ' = triaxial ellipsoid geocentric latitude.
-With the WGS84 latitude β , we have Tan β' = (6356.752314/6378.137)² Tan β ( β' = Triaxial ellipsoid geocentric latitude )

Data Registers: /
Flags: /
Subroutines: /

01 LBL "GD2GC3"
02 DEG
03 HR
04 1
05 P-R
06 .99330562
07 /
08 R-P
09 RDN
10 HMS
11 X<>Y
12 HR
13 14.92911
14 +
15 HMS
16 END

( 42 bytes / SIZE 000 )

STACK	INPUTS	OUTPUTS
Y	λ ( ° ' " )	β' ( ° ' " )
X	β ( ° ' " )	λ' ( ° ' " )

Example: λ = 124°41 ' β = 57°16' ( WGS84 geocentric longitude & geodetic latitude )

124.41 ENTER^
57.16 XEW "GD2GC3" >>>> λ' = 139°36'44"796 ---Execution time = 2s---
X<>Y β' = 57°05'28"9352

Remark:

-The values of a & b that are used in all these programs are those given by WGS84.

-According to G. Panou, R. Korakitis, and G. Pantazis ( cf "Fitting a triaxial ellipsoid to a geoid model" ( 2020 ) )

a = 6378.17192 km
b = 6378.10206 km
c = 6356.75217 km

and the equatorial longitude of the major semi-axis = –14.9367 degrees

-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

hp41programs

Terrestrial Geodesic Distance(4) for the HP-41