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)  F
ormula 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 ) / ( A2 + B2 + C2 )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' ) / a2 + ( y Sin W ) ( y Cos W - y' ) / b2 + ( z Sin W ) ( z Cos W - z' ) / c2 ]
Q =  [ x2 Sin2 W - ( x Cos W - x' )2 ] / a2 + [ y2 Sin2 W - ( y Cos W - y' )2 ] / b2 + [ z2 Sin2 W - ( z Cos W - z' )2 ] / c2

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

A2 + B2 + C2 = T Sin2 ( W - µ ) + S Sin2 µ + 2 U Sin µ Sin ( W - µ )
where  S = x'2 / a2 + y'2 / b2 + z'2 / c2  ,   T =  x2 / a2 + y2 / b2 + z2 / c2  ,   U =  xx' / a2 + yy' / b2 + zz' / c2
and   P = 2 U - 2 T Cos W  ,  Q = T Sin W + ( 2 U Cos W - S - T Cos2 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:       A2 + B2 + C2 = T Sin2 ( W/2 - µ ) + S Sin2 ( 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))  +(f2 /4) [W/4+(SinW)(CosW)/4-(Cos2 (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 = sin2 G  cos2 H  +  cos2 F  sin2 H  ;   µ = Arcsin S1/2  ;  R = ( S.(1-S) )1/2

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

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

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

and   f = flattening = 1/298.25722

And we also have:

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

( here,   µ = 2 Arcsin S1/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 µ )

+  f2 [ p ((m-1)/2) µ2 Csc µ + m ((1-m)/2) µ2 Cot µ + ((m2-8.p2)/16) Sin µ Cos µ + (m2/16) µ - (m2/8) Sin µ Cos3 µ + (mp/2) Sin µ Cos2 µ ]

+  f3 [ µ2 Csc µ (p/4) (-2+5m-3m2) + µ3 Csc µ Cot µ (p/2) (1-4m+3m2) - µ2 Cot µ (m/4) (-2+5m-3m2) + µ3 (m/6) (-1+2m-m2)
+ µ3 Cot2 µ (m/2) (-1+3m-2m2) + Sinµ Cos µ (-8 p2 + m2 - 8 p2m - m3/2 )/16 + µ  (8 p2 + m2 - 8 p2m - m3/2 )/16
+ µ3 Csc2 µ (p2/2) (1-m) + Sin µ Cos3 µ (m2/16) (-2+m) + µ Cos2 µ (m2/4) (1-m) + Sin µ Cos2 µ (pm/2) (1+m) + p2m Sin3 µ Cos µ
+ µ Cos µ (3pm/4) (-1+m) + (pm2/2) Sin5 µ + Sin µ (p/2) (p2-m2) + (m3/12) Sin3 µ Cos3 µ - (2p3/3) Sin3 µ ]  }

where    µ = 2 Arc Sin [ sin2 ((B'-B)/2)  cos2 ((L'-L)/2)  +  cos2 ((B+B')/2)  sin2 ((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) ]2 / Sin2 µ   ,   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 = sin2 G  cos2 H  +  cos2 F  sin2 H ,  µ = 2 Arc Sin S1/2

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

-After "some" calculations, it yields:

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

+ (f2/4) [ A µ2 Cot µ + B µ2 / Sin µ + A2 ( µ + Sin µ Cos µ - 4 µ2 Cot µ ) / 4 - A.B µ2 / Sin µ - ( B2/2 ) Sin µ Cos µ ]

+ (f3/8) [ A ( 2 µ2 Cot µ - 2 µ3 Cot2 µ - 8 µ3 / 6 ) + A2 ( -5 µ2 Cot µ + 5 µ3 Cot2 µ + 8 µ3 / 3 + (1/2) Sin µ Cos µ + µ /2 )
+ A3 ( 3 µ2 Cot µ - 3 µ3 Cot2 µ - 8 µ3 / 6 - (1/4) Sin µ Cos µ - µ / 4 )
+ B ( 2 µ2 / Sin µ - 2 µ3 Cos µ / Sin2 µ ) + A.B ( -5 µ2 / Sin µ + 6 µ3 Cos µ / Sin2 µ + µ Cos µ )
+ A2 B ( 3 µ2 / Sin µ - 4 µ3 Cos µ / Sin2 µ - µ Cos µ + (1/2) Sin µ Cos2 µ )
+ B2 ( - Sin µ Cos µ + µ3 / Sin2 µ  + µ ) + A B2 ( (1/2) Sin µ Cos µ - µ3 / Sin2 µ - µ )
+ B3 ( (2/3) Sin3 µ - (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 bparam Y b' bparam 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  =  (f2/4) µ [ (A2-A) ( 1 - µ Cot µ ) - B.(A-1) ( µ / Sin µ - (Sin µ) / µ ) ]                 where  µ , A , B are computed with parametric latitudes

+ (f3/8) [ A ( 2 µ2 Cot µ - 2 µ3 Cot2 µ - 8 µ3 / 6 ) + A2 ( -5 µ2 Cot µ + 5 µ3 Cot2 µ + 8 µ3 / 3 + Sin µ Cos µ - µ  )
+ A3 ( 3 µ2 Cot µ - 3 µ3 Cot2 µ - 8 µ3 / 6 - Sin µ Cos µ + µ )
+ B ( 2 µ2 / Sin µ - 2 µ3 Cos µ / Sin2 µ ) + A.B ( -5 µ2 / Sin µ + 6 µ3 Cos µ / Sin2 µ + µ Cos µ - 2 Sin µ )
+ A2 B ( 3 µ2 / Sin µ - 4 µ3 Cos µ / Sin2 µ - µ Cos µ + 2 Sin µ )
+ B2 ( - 2 Sin µ Cos µ + µ3 / Sin2 µ  + µ ) + A B2 ( 2 Sin µ Cos µ - µ3 / Sin2 µ - µ )

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  =  (f2/4) µ [ (A2-A) ( 1 - µ Cot µ ) - B.(A-1) ( µ / Sin µ - (Sin µ) / µ ) ]                 where  µ , A , B are computed with geodetic latitudes

+ (f3/8) [ A ( 4 µ - 2 µ2 Cot µ - 2 µ3 Cot2 µ - 8 µ3 / 6 )
+ A2 ( 9 µ2 Cot µ + 5 µ3 Cot2 µ + 8 µ3 / 3 - Sin µ Cos µ - 13 µ  )
+ A3 ( -7 µ2 Cot µ - 3 µ3 Cot2 µ - 8 µ3 / 6 + Sin µ Cos µ + 9 µ )
+ B ( -2 µ2 / Sin µ - 2 µ3 Cos µ / Sin2 µ + 4 Sin µ )
+ A.B ( 9 µ2 / Sin µ + 6 µ3 Cos µ / Sin2 µ - 3 µ Cos µ - 12 Sin µ )
+ A2 B ( -7 µ2 / Sin µ - 4 µ3 Cos µ / Sin2 µ + 3 µ Cos µ + 8 Sin µ )
+ B2 ( 2 Sin µ Cos µ + µ3 / Sin2 µ  - 3 µ )
+ A B2 ( -2 Sin µ Cos µ - µ3 / Sin2 µ + 3 µ )

-We may simplify this formula like this:

(d - Δ') / a  =  µ ( 1 - A ) { (f2/4) [ A ( µ Cot µ - 1 ) + B ( µ / Sin µ - (Sin µ) / µ ) ]
+ (f3/8) [ B2 ( µ2 / Sin2 µ + 2 ( Sin µ  Cos µ ) / µ - 3 )
+ B (  ( 4 - 8 A ) ( Sin µ ) / µ - 3 A Cos µ + ( 4 A - 2 ) µ2 Cos µ / Sin2 µ + ( 7 A - 2 ) µ / Sin µ )
+ A ( -A ( Sin µ  Cos µ ) / µ + ( 7 A - 2 ) µ Cot µ - 8 ( 1 - A ) ( µ2 / 6 ) + ( 3 A - 2 )  µ2 Cot2 µ + 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 = sin2 G  cos2 H  +  cos2 F  sin2 H  ;   µ = 2 Arcsin S1/2

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

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

+ (f2/4) [ A ( 4.µ + µ2 Cot µ ) + B ( 6.Sin µ + µ2 / Sin µ ) - A2 ( 15.µ + 15 Sin µ Cos µ + 4 µ2 Cot µ ) / 4 - A.B ( 6.Sin µ + µ2 / Sin µ ) + ( B2/2 ) 15 Sin µ Cos µ ]

+ (f3/8) [ A ( 2 µ2 Cot µ - 2 µ3 Cot2 µ - 8 µ3 / 6 ) + A2 ( -5 µ2 Cot µ + 5 µ3 Cot2 µ + 8 µ3 / 3 )
+ A3 ( 3 µ2 Cot µ - 3 µ3 Cot2 µ - 8 µ3 / 6 )
+ B ( 2 µ2 / Sin µ - 2 µ3 Cos µ / Sin2 µ ) + A.B ( -5 µ2 / Sin µ + 6 µ3 Cos µ / Sin2 µ )
+ A2 B ( 3 µ2 / Sin µ - 4 µ3 Cos µ / Sin2 µ  )
+ B2 ( µ3 / Sin2 µ  ) + A B2 ( - µ3 / Sin2 µ )

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 = sin2 G  cos2 H  +  cos2 F  sin2 H  ;   µ = 2 Arcsin S1/2

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

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

+ a (f3/8) [ A ( 2 µ2 Cot µ - 2 µ3 Cot2 µ - 8 µ3 / 6 ) + A2 ( -5 µ2 Cot µ + 5 µ3 Cot2 µ + 8 µ3 / 3 )
+ A3 ( 3 µ2 Cot µ - 3 µ3 Cot2 µ - 8 µ3 / 6 )
+ B ( 2 µ2 / Sin µ - 2 µ3 Cos µ / Sin2 µ ) + A.B ( -5 µ2 / Sin µ + 6 µ3 Cos µ / Sin2 µ )
+ A2 B ( 3 µ2 / Sin µ - 4 µ3 Cos µ / Sin2 µ  )
+ B2 ( µ3 / Sin2 µ  ) + A B2 ( - µ3 / Sin2 µ )

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 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 (f2/4) µ [ (A2-A) ( 1 - µ Cot µ ) - B.(A-1) ( µ / Sin µ - (Sin µ) / µ ) ]

+ a (f3/8) [ A ( 2 µ2 Cot µ - 2 µ3 Cot2 µ - 8 µ3 / 6 ) + A2 ( -5 µ2 Cot µ + 5 µ3 Cot2 µ + 8 µ3 / 3 )
+ A3 ( 3 µ2 Cot µ - 3 µ3 Cot2 µ - 8 µ3 / 6 )
+ B ( 2 µ2 / Sin µ - 2 µ3 Cos µ / Sin2 µ ) + A.B ( -5 µ2 / Sin µ + 6 µ3 Cos µ / Sin2 µ )
+ A2 B ( 3 µ2 / Sin µ - 4 µ3 Cos µ / Sin2 µ  )
+ B2 ( µ3 / Sin2 µ  ) + A B2 ( - µ3 / Sin2 µ )

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 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 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 XY 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 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 y2 + f y'2 + f ' z2 + f ' z'2 )  ( W/2 - (1/4) Sin 2.W ) + ( f y.y' + f ' z.z' ) ( Sin W - W Cos W ) } / Sin2 W ]

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

with       R / Sin W = 1 / ( A2 + B2 + C2 )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 + Δ1 - Δ2
-It might also be calculated by  D ~ d . Δ12   ( 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 = [ y2 ( 2.f + 3.f2 ) + z2 ( 2.f ' + 3.f '2 ) ] / ( 2 Sin2 W )
With    n  = [ y'2 ( 2.f + 3.f2 ) + z'2 ( 2.f ' + 3.f '2 ) ] / ( 2 Sin2 W )         [ In fact,  2.f + 3.f2 & 2.f ' + 3.f '2  may be replaced by  (a2/b2 - 1)/2  and  (a2/c2 - 1)/2 ]
p  = [ y.y' ( 2.f + 3.f2 ) + z.z' ( 2.f ' + 3.f '2 ) ] / ( Sin2 W )

Δ / (a.Sin W)  ~   W / Sin W + ( p2/2 - m - n ) ( W/Sin W - Cos W ) / 2  + (p/2) ( W Cot W - 1 )
+  ( m2 + n2 ) [ (13/16) ( W/Sin W - Cos W ) + (1/8) Sin2 W Cos W ) ]
+ 3 ( p2/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) Sin4 W + ( Cos W ) { (-13/8) ( W/Sin W - Cos W ) + (3/4) Sin2 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' ( ° ' " ) D2 (km ) X b' ( ° ' " ) D3 ( km ) L / d3 (km)

Where  D3 = geodesic distance ( triaxial ellipsoid )             & d3 = great elliptic distance ( triaxial ellipsoid ) = R20
D2 = geodesic distance ( biaxial ellipsoid ) = R00  & d2 =  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"   >>>>     D3 = 12138.65754 km                                                       ---Execution time = 72s---
X<>Y    D2  = 12138.68432 km = R00

and    LASTX = d3 = 12138.65876 km = R20  &   RCL 21 = d2 = 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] !