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 Distance ( Ellipsoid of Revolution )

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

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

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

0°)  Method

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

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

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

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

D ~  Δ + ( d - Δ' )

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

1°)  With Andoyer's Formula of order 2

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

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

OP = ( Sin W ) / ( 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 X#0?  85 /  86 ATAN  87 STO 08  88 2  89 ST/ 01  90 ST/ 11  91 /  92 STO 05  93 XEQ 02  94 STO 06  95 RCL 05  96 90  97 +  98 XEQ 02  99 STO 07 100 RCL 02 101 D-R 102 RCL 08 103 COS 104 STO 05 105 X^2 106 ST+ X 107 RCL 04 108 RCL 15 109 ST* 05 110 * 111 ST* Y 112 - 113 - 114 4 115 / 116 RCL 06           117 ST/ 15 118 RCL 07 119 ST- Y 120 ST+ X 121 / 122 STO 07 123 * 124 RCL 05 125 - 126 + 127 RCL 07 128 * 129 + 130 ST* 15 131 RCL 14 132 ENTER 133 GTO 03 134 LBL 00 135 RCL 05 136 * 137 RCL 09 138 RCL 06 139 * 140 + 141 RCL 10 142 RCL 07 143 * 144 + 145 RTN 146 LBL 01 147 1 148 P-R 149 X<>Y 150 RCL 03 151 X^2 152 * 153 STO 10           154 RDN 155 P-R 156 RCL 01 157 X^2 158 * 159 STO Z 160 X^2 161 X<>Y 162 RCL 02 163 X^2 164 * 165 STO 09 166 X^2 167 RCL 10 168 X^2 169 + 170 + 171 SQRT 172 ST/ 09 173 ST/ 10 174 / 175 STO 08 176 RTN 177 LBL 02 178 RCL 02 179 2 180 / 181 + 182 ENTER 183 SIN 184 ENTER 185 X^2 186 RCL 13 187 * 188 RCL 02 189 R^ 190 - 191 SIN 192 ST* Z 193 X^2 194 RCL 00           195 * 196 + 197 X<>Y 198 RCL 09 199 * 200 + 201 SQRT 202 RTN 203 LBL 03 204 RCL 12 205 + 206 2 207 / 208 ST- Y 209 1 210 P-R 211 X^2 212 STO T 213 RDN 214 X^2 215 STO 08 216 SIGN 217 P-R 218 X^2 219 ST* 08 220 CLX 221 + 222 X^2 223 ST* T 224 - 225 RCL 11 226 COS 227 X^2 228 * 229 - 230 / 231 STO 02 232 RCL 08           233 1 234 LASTX 235 STO 08 236 - 237 / 238 ST- 02 239 + 240 X<> 08 241 SQRT 242 ASIN 243 ST+ X 244 D-R 245 STO 00 246 LASTX 247 1 248 P-R 249 RDN 250 / 251 STO 07 252 R^ 253 * 254 RCL 08 255 ST* Y 256 - 257 RCL 07 258 ENTER 259 1/X 260 - 261 RCL 02 262 * 263 + 264 1 265 RCL 08 266 - 267 * 268 RCL 03           269 RCL 01 270 ST- Y 271 ST+ X 272 / 273 X^2 274 * 275 RCL 00 276 * 277 RCL 01 278 * 279 RCL 15 280 + 281 END

( 358 bytes / SIZE 016 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

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

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

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

Notes:

-Longitudes are geodetic and positive East.

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

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

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

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

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

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

2°)  With Andoyer's Formula of order 3

a) Ellipsoid of Revolution

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

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

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

+  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 X#0?  85 /  86 ATAN  87 STO 08  88 2  89 ST/ 01  90 ST/ 11  91 /  92 STO 05  93 XEQ 02  94 STO 06  95 RCL 05  96 90  97 +  98 XEQ 02  99 STO 07 100 RCL 02           101 D-R 102 RCL 08 103 COS 104 STO 05 105 X^2 106 ST+ X 107 RCL 04 108 RCL 15 109 ST* 05 110 * 111 ST* Y 112 - 113 - 114 4 115 / 116 RCL 06 117 ST/ 15 118 RCL 07 119 ST- Y 120 ST+ X 121 / 122 STO 07 123 * 124 RCL 05 125 - 126 + 127 RCL 07 128 * 129 + 130 ST* 15 131 GTO 03 132 LBL 00 133 RCL 05 134 * 135 RCL 09 136 RCL 06 137 * 138 + 139 RCL 10 140 RCL 07 141 * 142 + 143 RTN 144 LBL 01 145 1 146 P-R 147 X<>Y 148 RCL 03 149 X^2 150 * 151 STO 10           152 RDN 153 P-R 154 RCL 01 155 X^2 156 * 157 STO Z 158 X^2 159 X<>Y 160 RCL 02 161 X^2 162 * 163 STO 09 164 X^2 165 RCL 10 166 X^2 167 + 168 + 169 SQRT 170 ST/ 09 171 ST/ 10 172 / 173 STO 08 174 RTN 175 LBL 02 176 RCL 02 177 2 178 / 179 + 180 ENTER 181 SIN 182 ENTER 183 X^2 184 RCL 13 185 * 186 RCL 02 187 R^ 188 - 189 SIN 190 ST* Z 191 X^2 192 RCL 00 193 * 194 + 195 X<>Y 196 RCL 09 197 * 198 + 199 SQRT 200 RTN 201 LBL 03 202 RCL 14 203 1 204 STO 08           205 STO 09 206 P-R 207 RCL 01 208 RCL 03 209 - 210 RCL 01 211 / 212 STO 04 213 SIGN 214 LASTX 215 - 216 STO 05 217 / 218 R-P 219 X<> 12 220 RCL 09 221 P-R 222 RCL 05 223 / 224 R-P 225 RDN 226 + 227 2 228 ST/ 04 229 / 230 ST- Y 231 RCL 09 232 P-R 233 X^2 234 STO T 235 RDN 236 X^2 237 STO 14 238 SIGN 239 P-R 240 X^2 241 ST* 14 242 CLX 243 + 244 X^2 245 ST* T 246 - 247 RCL 11 248 COS 249 X^2 250 * 251 - 252 / 253 STO 02 254 RCL 14 255 RCL 09 256 LASTX 257 STO 14           258 - 259 / 260 ST- 02 261 + 262 ST- 09 263 X<> 14 264 SQRT 265 ASIN 266 ST+ X 267 D-R 268 STO 00 269 LASTX 270 RCL 08 271 P-R 272 STO 03 273 STO 06 274 STO 13 275 RDN 276 STO 05 277 / 278 STO 07 279 ST/ Z 280 X^2 281 ST* 03 282 RCL 08 283 + 284 R^ 285 ST- 08 286 ST+ X 287 - 288 RCL 02 289 * 290 2 291 RCL 07 292 ST- 03 293 ST/ Y 294 - 295 RCL 06 296 - 297 RCL 03 298 ST+ 03 299 ST* 13 300 4 301 * 302 - 303 RCL 09 304 ST* Z 305 * 306 RCL 03 307 - 308 RCL 14           309 * 310 - 311 RCL 03 312 - 313 RCL 02 314 * 315 RCL 13 316 RCL 00 317 ST* 01 318 ST+ X 319 X^2 320 3 321 ST* Z 322 / 323 STO 11 324 + 325 STO 12 326 RCL 08 327 - 328 RCL 14 329 * 330 RCL 12 331 ST+ X 332 - 333 RCL 13 334 + 335 RCL 08 336 + 337 RCL 14 338 * 339 RCL 13 340 ST+ X 341 + 342 RCL 11 343 + 344 RCL 14 345 * 346 - 347 RCL 04 348 * 349 RCL 06 350 RCL 07 351 * 352 RCL 14 353 ST* Y 354 - 355 RCL 07 356 ENTER 357 1/X 358 - 359 RCL 02 360 * 361 + 362 RCL 09           363 * 364 + 365 RCL 04 366 X^2 367 * 368 RCL 01 369 * 370 RCL 15 371 + 372 END

( 466 bytes / SIZE 016 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

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

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

LASTX = R10 =  12138.65875 km = great elliptic arc distance

Notes:

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

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

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

c) Triaxial Ellipsoid  ( Without Parametric Coordinates )

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

(d - Δ') / a  =  (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 µ )

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

 01 LBL "TGD"  02 DEG  03 HR  04 STO 14  05 RDN  06 HR  07 14.92911  08 STO 11  09 +  10 STO 13  11 RDN  12 HR  13 STO 12  14 X<>Y  15 HR  16 ST+ 11  17 6378.172  18 STO 01  19 .07  20 -  21 STO 02  22 21.349686  23 -  24 STO 03  25 RCL 11  26 RCL 12  27 XEQ 01  28 STO 05  29 RCL 09  30 STO 06  31 RCL 10  32 STO 07  33 RCL 13  34 ST- 11  35 RCL 14  36 XEQ 01  37 XEQ 00  38 STO 04  39 RCL 01  40 ST/ 05  41 ST/ 08  42 RCL 02  43 ST+ 01  44 ST/ 06  45 ST/ 09  46 RCL 03            47 ST/ 07  48 ST/ 10  49 RCL 08  50 X^2  51 RCL 09  52 X^2 53 RCL 10  54 X^2  55 +  56 +  57 STO 13  58 RCL 05  59 X^2  60 RCL 06  61 X^2  62 RCL 07  63 X^2  64 +  65 +  66 STO 00  67 -  68 RCL 04  69 ACOS  70 STO 02  71 SIN  72 STO 15  73 *  74 RCL 08  75 XEQ 00  76 ST+ X  77 STO 09  78 RCL 00  79 RCL 13  80 +  81 RCL 04  82 *  83 -  84 X#0?  85 /  86 ATAN  87 STO 08  88 2  89 ST/ 01  90 ST/ 11  91 /  92 STO 05  93 XEQ 02  94 STO 06  95 RCL 05  96 90  97 +  98 XEQ 02  99 STO 07 100 RCL 02           101 D-R 102 RCL 08 103 COS 104 STO 05 105 X^2 106 ST+ X 107 RCL 04 108 RCL 15 109 ST* 05 110 * 111 ST* Y 112 - 113 - 114 4 115 / 116 RCL 06 117 ST/ 15 118 RCL 07 119 ST- Y 120 ST+ X 121 / 122 STO 07 123 * 124 RCL 05 125 - 126 + 127 RCL 07 128 * 129 + 130 ST* 15 131 RCL 14 132 ENTER 133 GTO 03 134 LBL 00 135 RCL 05 136 * 137 RCL 09 138 RCL 06 139 * 140 + 141 RCL 10 142 RCL 07 143 * 144 + 145 RTN 146 LBL 01 147 1 148 P-R 149 X<>Y 150 RCL 03           151 X^2 152 * 153 STO 10 154 RDN 155 P-R 156 RCL 01 157 X^2 158 * 159 STO Z 160 X^2 161 X<>Y 162 RCL 02 163 X^2 164 * 165 STO 09 166 X^2 167 RCL 10 168 X^2 169 + 170 + 171 SQRT 172 ST/ 09 173 ST/ 10 174 / 175 STO 08 176 RTN 177 LBL 02 178 RCL 02 179 2 180 / 181 + 182 ENTER 183 SIN 184 ENTER 185 X^2 186 RCL 13 187 * 188 RCL 02 189 R^ 190 - 191 SIN 192 ST* Z 193 X^2 194 RCL 00 195 * 196 + 197 X<>Y 198 RCL 09 199 * 200 + 201 SQRT 202 RTN 203 LBL 03 204 RCL 12           205 + 206 2 207 / 208 ST- Y 209 1 210 STO 07 211 P-R 212 X^2 213 STO T 214 RDN 215 X^2 216 STO 08 217 SIGN 218 P-R 219 X^2 220 ST* 08 221 CLX 222 + 223 X^2 224 ST* T 225 - 226 RCL 11 227 COS 228 X^2 229 * 230 - 231 / 232 STO 02 233 RCL 08 234 RCL 07 235 LASTX 236 STO 08 237 - 238 / 239 ST- 02 240 + 241 ST- 07 242 X<> 08 243 SQRT 244 ASIN 245 ST+ X 246 D-R 247 STO 00 248 LASTX 249 1 250 P-R 251 STO 09 252 STO 12 253 STO 13           254 RDN 255 / 256 STO 06 257 ST/ 12 258 X^2 259 ST* 09 260 3 261 - 262 RCL 12 263 ST+ X 264 + 265 RCL 02 266 * 267 3 268 RCL 06 269 ST* Y 270 RCL 09 271 + 272 ST+ X 273 STO 10 274 STO 11 275 ST+ X 276 + 277 ST- 09 278 RCL 13 279 3 280 * 281 - 282 8 283 RCL 06 284 / 285 - 286 RCL 07 287 ST* Z 288 * 289 RCL 10 290 + 291 4 292 STO 14 293 RCL 06 294 / 295 ST- 11 296 - 297 RCL 08 298 * 299 + 300 RCL 11 301 - 302 RCL 02 303 * 304 RCL 09 305 RCL 13           306 ST* 10 307 * 308 RCL 12 309 + 310 9 311 + 312 RCL 00 313 ST+ X 314 X^2 315 3 316 / 317 ST- 14 318 - 319 RCL 07 320 * 321 RCL 10 322 - 323 RCL 14 324 ST- 10 325 + 326 RCL 08 327 * 328 RCL 10 329 + 330 RCL 08 331 * 332 - 333 RCL 01 334 RCL 03 335 - 336 RCL 01 337 ST+ X 338 / 339 STO 04 340 * 341 RCL 06 342 RCL 13 343 * 344 RCL 08 345 ST* Y 346 - 347 RCL 06 348 ENTER 349 1/X 350 - 351 RCL 02 352 * 353 + 354 RCL 07 355 * 356 + 357 RCL 04           358 X^2 359 * 360 RCL 00 361 * 362 RCL 01 363 * 364 RCL 15 365 + 366 END

( 458 bytes / SIZE 016 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

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

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

LASTX = R15 =  12138.65875 km = great elliptic arc distance

Note:

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

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

a)  Ellipsoid of Revolution

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

Formulae:

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

S = 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 X#0?  85 /  86 ATAN  87 STO 08  88 2  89 ST/ 01  90 ST/ 11  91 /  92 STO 05            93 XEQ 02  94 STO 06  95 RCL 05  96 90 97 +  98 XEQ 02  99 STO 07 100 RCL 02 101 D-R 102 RCL 08 103 COS 104 STO 05 105 X^2 106 ST+ X 107 RCL 04 108 RCL 15 109 ST* 05 110 * 111 ST* Y 112 - 113 - 114 4 115 / 116 RCL 06 117 ST/ 15 118 RCL 07 119 ST- Y 120 ST+ X 121 / 122 STO 07 123 * 124 RCL 05 125 - 126 + 127 RCL 07 128 * 129 + 130 ST* 15 131 RCL 14 132 ENTER 133 GTO 03 134 LBL 00 135 RCL 05 136 * 137 RCL 09 138 RCL 06 139 * 140 + 141 RCL 10           142 RCL 07 143 * 144 + 145 RTN 146 LBL 01 147 1 148 P-R 149 X<>Y 150 RCL 03 151 X^2 152 * 153 STO 10 154 RDN 155 P-R 156 RCL 01 157 X^2 158 * 159 STO Z 160 X^2 161 X<>Y 162 RCL 02 163 X^2 164 * 165 STO 09 166 X^2 167 RCL 10 168 X^2 169 + 170 + 171 SQRT 172 ST/ 09 173 ST/ 10 174 / 175 STO 08 176 RTN 177 LBL 02 178 RCL 02 179 2 180 / 181 + 182 ENTER 183 SIN 184 ENTER 185 X^2 186 RCL 13 187 * 188 RCL 02           189 R^ 190 - 191 SIN 192 ST* Z 193 X^2 194 RCL 00 195 * 196 + 197 X<>Y 198 RCL 09 199 * 200 + 201 SQRT 202 RTN 203 LBL 03 204 RCL 12 205 + 206 2 207 / 208 ST- Y 209 1 210 STO 07 211 P-R 212 X^2 213 STO T 214 RDN 215 X^2 216 STO 08 217 SIGN 218 P-R 219 X^2 220 ST* 08 221 CLX 222 + 223 X^2 224 ST* T 225 - 226 RCL 11 227 COS 228 X^2 229 * 230 - 231 / 232 STO 02 233 RCL 08 234 RCL 07           235 LASTX 236 STO 08 237 - 238 / 239 ST- 02 240 + 241 ST- 07 242 X<> 08 243 SQRT 244 ASIN 245 ST+ X 246 D-R 247 STO 00 248 LASTX 249 1 250 P-R 251 STO 04 252 STO 05 253 STO 09 254 RDN 255 / 256 ST* 05 257 STO 06 258 X^2 259 ST* 09 260 RCL 02 261 * 262 RCL 06 263 ST- 09 264 RCL 09 265 ST* 04 266 ST+ 09 267 4 268 * 269 + 270 RCL 07 271 ST* Z 272 * 273 RCL 09 274 + 275 RCL 08 276 * 277 + 278 RCL 09 279 - 280 RCL 02 281 * 282 RCL 04           283 RCL 00 284 ST+ X 285 X^2 286 3 287 ST* Z 288 / 289 STO Z 290 + 291 STO 09 292 RCL 08 293 * 294 RCL 09 295 ST+ X 296 - 297 RCL 04 298 ST+ Z 299 ST+ Z 300 + 301 RCL 08 302 ST* 05 303 * 304 + 305 RCL 08 306 ST- 05 307 * 308 - 309 RCL 01 310 RCL 03 311 - 312 RCL 01 313 ST+ X 314 / 315 STO 04 316 * 317 RCL 05 318 RCL 06 319 ENTER 320 1/X 321 - 322 RCL 02 323 * 324 + 325 RCL 07 326 * 327 + 328 RCL 04 329 X^2 330 * 331 RCL 00           332 * 333 RCL 01 334 * 335 RCL 15 336 + 337 END

( 427 bytes / SIZE 016 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

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

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

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

Notes:

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

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

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

c)  Triaxial Ellipsoid  ( With Parametric Coordinates )

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

d - Δ' = a (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 X#0?  85 /  86 ATAN  87 STO 08  88 2  89 ST/ 01  90 ST/ 11  91 /  92 STO 05  93 XEQ 02  94 STO 06  95 RCL 05            96 90  97 +  98 XEQ 02  99 STO 07 100 RCL 02 101 D-R 102 RCL 08 103 COS 104 STO 05 105 X^2 106 ST+ X 107 RCL 04 108 RCL 15 109 ST* 05 110 * 111 ST* Y 112 - 113 - 114 4 115 / 116 RCL 06 117 ST/ 15 118 RCL 07 119 ST- Y 120 ST+ X 121 / 122 STO 07 123 * 124 RCL 05 125 - 126 + 127 RCL 07 128 * 129 + 130 ST* 15 131 GTO 03 132 LBL 00 133 RCL 05 134 * 135 RCL 09 136 RCL 06 137 * 138 + 139 RCL 10           140 RCL 07 141 * 142 + 143 RTN 144 LBL 01 145 1 146 P-R 147 X<>Y 148 RCL 03 149 X^2 150 * 151 STO 10 152 RDN 153 P-R 154 RCL 01 155 X^2 156 * 157 STO Z 158 X^2 159 X<>Y 160 RCL 02 161 X^2 162 * 163 STO 09 164 X^2 165 RCL 10 166 X^2 167 + 168 + 169 SQRT 170 ST/ 09 171 ST/ 10 172 / 173 STO 08 174 RTN 175 LBL 02 176 RCL 02 177 2 178 / 179 + 180 ENTER 181 SIN 182 ENTER 183 X^2 184 RCL 13 185 * 186 RCL 02 187 R^ 188 - 189 SIN 190 ST* Z 191 X^2 192 RCL 00           193 * 194 + 195 X<>Y 196 RCL 09 197 * 198 + 199 SQRT 200 RTN 201 LBL 03 202 RCL 14 203 1 204 STO 07 205 P-R 206 RCL 01 207 RCL 03 208 - 209 RCL 01 210 / 211 STO 04 212 SIGN 213 LASTX 214 - 215 STO 05 216 / 217 R-P 218 X<> 12 219 RCL 07 220 P-R 221 RCL 05 222 / 223 R-P 224 RDN 225 + 226 2 227 ST/ 04 228 / 229 ST- Y 230 RCL 07 231 P-R 232 X^2 233 STO T 234 RDN 235 X^2 236 STO 08 237 SIGN 238 P-R 239 X^2 240 ST* 08 241 CLX 242 + 243 X^2 244 ST* T 245 - 246 RCL 11           247 COS 248 X^2 249 * 250 - 251 / 252 STO 02 253 RCL 08 254 RCL 07 255 LASTX 256 STO 08 257 - 258 / 259 ST- 02 260 + 261 ST- 07 262 X<> 08 263 SQRT 264 ASIN 265 ST+ X 266 D-R 267 STO 00 268 LASTX 269 1 270 P-R 271 STO 05 272 STO 09 273 STO 13 274 RDN 275 / 276 ST* 05 277 STO 06 278 X^2 279 ST* 09 280 RCL 02 281 * 282 RCL 06 283 ST- 09 284 RCL 09 285 ST+ 09 286 ST* 13 287 4 288 * 289 + 290 RCL 07 291 ST* Z 292 * 293 RCL 09 294 + 295 RCL 08           296 * 297 + 298 RCL 09 299 - 300 RCL 02 301 * 302 RCL 13 303 RCL 00 304 ST* 01 305 ST+ X 306 X^2 307 3 308 ST* Z 309 / 310 STO Z 311 + 312 STO 09 313 RCL 08 314 * 315 RCL 09 316 ST+ X 317 - 318 RCL 13 319 ST+ Z 320 ST+ Z 321 + 322 RCL 08 323 ST* 05 324 * 325 + 326 RCL 08 327 ST- 05 328 * 329 - 330 RCL 04 331 * 332 RCL 05 333 RCL 06 334 ENTER 335 1/X 336 - 337 RCL 02 338 * 339 + 340 RCL 07 341 * 342 + 343 RCL 04           344 X^2 345 * 346 RCL 01 347 * 348 RCL 15 349 + 350 END

( 444 bytes / SIZE 016 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

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

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

Notes:

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

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

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

4°)  With Vincenty's Formula

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

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

-Line 47 is a three-byte GTO 14

 01 LBL "TGDV"  02 DEG  03 HR  04 STO 21          05 RDN  06 HR  07 14.92911  08 STO 00  09 +  10 STO 20  11 RDN  12 HR  13 STO 19  14 X<>Y  15 HR  16 RCL 00  17 +  18 STO 18            19 6378.172  20 STO 01  21 .07  22 -  23 STO 02  24 6356.752314  25 STO 03  26 XEQ 01  27 STO 17  28 RCL 01  29 RCL 02  30 +  31 2  32 /  33 STO 01  34 STO 02  35 XEQ 01  36 ST- 17  37 RCL 01  38 STO 06  39 RCL 03  40 STO 07  41 RCL 21  42 RCL 19  43 RCL 20  44 RCL 18  45 XEQ 13  46 ENTER  47 GTO 14  48 LBL 01  49 RCL 18  50 RCL 19  51 XEQ 02  52 RCL 08  53 STO 05  54 RCL 09  55 STO 06  56 RCL 10  57 STO 07  58 RCL 20  59 RCL 21  60 XEQ 02  61 RCL 05  62 RCL 08  63 *  64 RCL 06 65 RCL 09  66 *  67 RCL 07   68 RCL 10          69 *  70 +  71 +  72 STO 00  73 ACOS  74 STO 04  75 RCL 01  76 ST/ 05  77 ST/ 08  78 RCL 02  79 ST/ 06  80 ST/ 09  81 RCL 03            82 ST/ 07  83 ST/ 10  84 RCL 05  85 X^2  86 RCL 06  87 X^2  88 RCL 07  89 X^2  90 +  91 +  92 STO 11  93 RCL 08  94 X^2  95 RCL 09  96 X^2  97 RCL 10  98 X^2  99 + 100 + 101 X<> 10 102 RCL 07 103 * 104 RCL 06 105 RCL 09 106 * 107 RCL 05 108 RCL 08 109 * 110 + 111 + 112 STO 09 113 RCL 00 114 * 115 ST+ X 116 RCL 10 117 - 118 RCL 11 119 RCL 00 120 X^2 121 * 122 - 123 RCL 11 124 RCL 04 125 SIN 126 STO 08 127 ST/ Z 128 * 129 + 130 RCL 09 131 RCL 11 132 RCL 00         133 * 134 - 135 ST+ X 136 X<>Y 137 X#0? 138 / 139 ATAN 140 2 141 / 142 90 143 STO 00 144 X<>Y 145 X<0? 146 + 147 STO 05 148 XEQ 03 149 STO 06           150 RCL 05 151 RCL 00 152 + 153 XEQ 03 154 STO 07 155 RCL 06 156 X>Y? 157 GTO 00 158 X<> 07 159 STO 06 160 RCL 00 161 ST+ 05 162 LBL 00 163 PI 164 R-D 165 RCL 05 166 X=Y? 167 CLX 168 STO 05 169 CHS 170 1 171 P-R 172 RCL 07 173 RCL 06 174 / 175 X^2 176 STO 09 177 * 178 R-P 179 CLX 180 RCL 04 181 RCL 05 182 - 183 1 184 P-R 185 RCL 09 186 * 187 R-P 188 CLX 189 ENTER 190 LBL 13 191 - 192 360 193 MOD 194 PI 195 R-D 196 STO 16         197 X>Y? 198 CLX 199 ST+ X 200 - 201 STO 11 202 STO 04 203 CLX 204 RCL 07 205 CHS 206 RCL 06 207 ST+ Y 208 / 209 STO 08 210 SIGN 211 LASTX 212 - 213 STO 15           214 P-R 215 LASTX 216 / 217 R-P 218 RDN 219 STO 09 220 CLX 221 RCL 15 222 P-R 223 LASTX 224 / 225 R-P 226 X<>Y 227 STO 10 228 LBL 12 229 RCL 16 230 RCL 04 231 X>Y? 232 ACOS 233 RCL 10 234 COS 235 STO 05 236 STO 12 237 RCL 04 238 SIN 239 * 240 STO 14 241 RCL 09 242 COS 243 ST* 05 244 ST* 14 245 RCL 10 246 SIN 247 STO 13 248 * 249 RCL 09 250 SIN 251 ST* 13 252 RCL 12 253 * 254 RCL 04 255 COS 256 ST* 05 257 * 258 - 259 R-P 260 RCL 05         261 RCL 13  262 + 263 R-P 264 X<> 14 265 X<>Y 266 STO 05 267 SIN 268 X#0? 269 / 270 ASIN 271 STO 12 272 COS 273 X^2 274 RCL 05 275 COS 276 RCL 13           277 ENTER 278 + 279 R^ 280 X=0? 281 SIGN 282 / 283 - 284 STO 13 285 CLX 286 3 287 * 288 4 289 X<>Y 290 - 291 RCL 08 292 * 293 4 294 + 295 * 296 RCL 08 297 * 298 16 299 / 300 RCL 05 301 SIN 302 RCL 13 303 RCL Z 304 * 305 * 306 R-D 307 RCL 05 308 + 309 RCL 12 310 SIN 311 * 312 RCL 08 313 * 314 1 315 R^ 316 - 317 * 318 RCL 11 319 + 320 ENTER 321 X<> 04 322 - 323 ABS 324  E-7 325 XY 393 RCL 03         394 X^2 395 * 396 STO 10 397 X^2 398 STO 04 399 RDN 400 P-R 401 RCL 01 402 X^2 403 * 404 STO 08 405 X^2 406 ST+ 04 407 X<>Y 408 RCL 02 409 X^2 410 * 411 STO 09           412 X^2 413 RCL 04 414 + 415 SQRT 416 ST/ 08 417 ST/ 09 418 ST/ 10 419 RTN 420 LBL 03 421 ENTER 422 SIN 423 STO 07 424 X^2 425 RCL 10 426 * 427 X<>Y 428 RCL 04 429 - 430 SIN 431 ST* 07 432 X^2 433 RCL 11 434 * 435 + 436 RCL 07 437 RCL 09 438 ST+ X 439 * 440 - 441 SQRT 442 RCL 08 443 X<>Y 444 / 445 RTN 446 LBL 14 447 RCL 17 448 + 449 END

( 564 bytes / SIZE 022 )

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

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

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

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

Notes:

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

5°)  Minimizing the Sum of 2 Distances

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

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

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

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

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

 01 LBL "D+D"  02 STO 26          03 RDN  04 STO 25   05 RDN  06 STO 24   07 X<>Y  08 STO 23  09 R^  10 X<>Y  11 HMS-  12 HR  13 RCL 26  14 HR  15 COS  16 P-R  17 RCL 24  18 HR  19 COS  20 +  21 R-P  22 CLX 23 RCL 23   24 HR  25 +  26 HMS  27 STO 27          28 GTO 00  29 LBL 12  30 STO 35  31 RCL 27  32 X<>Y  33 RCL 25  34 RCL 26  35 XEQ "TGD"  36 STO 31  37 RCL 23  38 RCL 24  39 RCL 27  40 RCL 35  41 XEQ "TGD"  42 RCL 31  43 +  44 RTN 45 LBL 00  46 90  47 ENTER  48 STO 29          49 CHS  50 STO 28   51 -  52 5  53 SQRT  54 1  55 +  56 2  57 /  58 STO 34  59 /  60 RCL 28  61 HR  62 +  63 HMS  64 STO 32  65 XEQ 12  66 STO 33 67 LBL 10  68 RCL 32   69 RCL 28          70 HMS-  71 HR  72 RCL 34   73 /  74 RCL 28  75 HR  76 +  77 HMS  78 STO 30  79 XEQ 12  80 STO 31  81 VIEW 31  82 RCL 33  83 X 32  87 STO 29  88 RCL 31 89 STO 33   90 GTO 04  91 LBL 02  92 RCL 29          93 STO 28   94 RCL 30  95 STO 29  96 LBL 04  97  E-4  99 RCL 28 100 RCL 30 101 HMS- 102 HR 103 ABS 104 X>Y? 105 GTO 10 106 RCL 30 107 RCL 31 108 CLD 109 END

( 180 bytes / SIZE 036 )

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

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

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

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

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

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

Notes:

-The HP41 displays the successive approximations.

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

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

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

6°)  Other Formulae for the Great Elliptic Distance

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

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

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

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

Δ  ~  a [ W - {  ( f 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 Distance

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

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

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

 01 LBL "TGD+"  02 DEG  03 HR  04 STO 01  05 X<> T  06 HMS-  07 HR  08 2  09 /  10 SIN  11 X^2  12 STO 00  13 X<>Y  14 HR  15 1  16 STO 11  17 P-R  18 LASTX  19 298.25722  20 STO 04  21 ST+ 04  22 1/X  23 -  24 STO 05  25 /  26 R-P  27 X<> 01  28 RCL 11  29 P-R  30 RCL 05            31 /  32 R-P  33 RDN  34 +  35 2 36 /  37 ST- Y  38 RCL 11  39 P-R  40 X^2  41 STO 02  42 X<>Y  43 X^2  44 STO 01  45 RDN  46 X<>Y  47 RCL 11  48 P-R  49 X^2  50 ST* 01  51 RDN  52 X^2  53 ST* 02  54 STO Z  55 -  56 RCL 00  57 *  58 +  59 ST/ 02  60 RCL 02  61 RCL 01  62 RCL 11  63 R^  64 STO 01            65 -  66 /  67 ST- 02  68 +  69 X<> 01  70 SQRT 71 ASIN  72 ST+ X  73 D-R  74 STO 00  75 LASTX  76 RCL 11  77 P-R  78 STO 03  79 STO 06  80 STO 08  81 RDN  82 STO 05  83 /  84 STO 07  85 ST/ 08  86 X^2  87 ST* 03  88 RCL 11  89 +  90 STO Y  91 RCL 08  92 ST+ 11  93 ST- Z  94 2  95 /  96 STO 09  97 -  98 RCL 01  99 * 100 - 101 RCL 05           102 X^2 103 1.5 104 STO 10 105 ST/ Y 106 FRC 107 - 108 RCL 07 109 / 110 RCL 02 111 * 112 + 113 RCL 02 114 * 115 RCL 06 116 X^2 117 RCL 07 118 ST- 03 119 ST+ X 120 / 121 RCL 03 122 4 123 ST/ 11 124 * 125 RCL 06 126 + 127 RCL 07 128 + 129 STO 05 130 - 131 RCL 01 132 * 133 RCL 05 134 + 135 RCL 03 136 ST+ X 137 STO 05           138 + 139 RCL 01 140 * 141 + 142 RCL 05 143 - 144 RCL 02 145 * 146 RCL 03 147 RCL 06 148 * 149 STO 03 150 3 151 * 152 RCL 00 153 X^2 154 ST+ X 155 RCL 10 156 / 157 STO 10 158 + 159 RCL 11 160 + 161 STO 05 162 RCL 01 163 * 164 RCL 05 165 ST+ X 166 - 167 RCL 03 168 ST+ 10 169 ST+ 10 170 + 171 RCL 01 172 * 173 RCL 10           174 + 175 RCL 01 176 * 177 - 178 RCL 04 179 / 180 RCL 11 181 RCL 06 182 RCL 07 183 * 184 STO Z 185 - 186 RCL 01 187 * 188 + 189 RCL 02 190 ST* 09 191 RCL 07 192 ST- 09 193 * 194 - 195 RCL 01 196 * 197 + 198 RCL 02 199 RCL 09 200 * 201 - 202 RCL 04 203 / 204 RCL 02 205 ST* 06 206 RCL 07 207 / 208 STO 09           209 RCL 01 210 + 211 STO 10 212 - 213 RCL 08 214 3 215 - 216 4 217 / 218 RCL 01 219 * 220 RCL 09 221 - 222 RCL 01 223 ST* Y 224 + 225 RCL 06 226 2 227 / 228 RCL 09           229 ST* Y 230 - 231 - 232 RCL 04 233 / 234 RCL 10 235 - 236 X<>Y 237 RCL 04 238 ST/ Z 239 / 240 RCL 00 241 ST* Z 242 ST* Y 243 ST+ Z 244 + 245 6378.137 246 ST* Z 247 * 248 END

( 310 bytes / SIZE 012 )

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

Where  D = geodesic distance & d = great elliptic distance

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

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

Notes:

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

Conclusion:

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

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

Remark:

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

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

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

References:

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

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