Template

# Terrestrial Geodesic Distance(3) for the HP-41

Overview

1°)  Ellipsoid of Revolution ( WGS84 )

a) Andoyer3 ( 2 programs )
b) Andoyer2 + A few terms of order 3 ( without calculating parametric latitudes )

b1) Program#1
b2) Program#2
b3) Program(s)#3
b4) Program#4

c) Vincenty's Formula
d) Geocentric Latitudes

d1) Program(s)#1
d2) Program#2
d3) Program#3
d4) Program(s)#4
d5) Program(s)#5 - My Favorite Programs
d6) Program#6 ( latitudes geocentric or geodetic )

2°)  Triaxial Ellipsoid

a) Andoyer2 + Andoyer2
b) Andoyer2 + Andoyer3
c) Geocentric Longitudes & Geodetic Latitudes
d) Geocentric Longitudes & Latitudes
e) Andoyer3 + Andoyer3

Latest Updates:   §1°) b4) & 1°) d6)

-With the triaxial ellipsoid, latitudes and longitudes are geodetic... except in paragraph 2°)c) & 2°)d)

-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

1°)  Ellipsoid of Revolution  ( WGS84 )

a)  Andoyer3

a = 6378.137 km
b = 6356.752314 km

-Like the program listed in "Terrestrial Geodesic Distance(2)",  this program uses Andoyer's formula of order 3.
-But it is re-written:

Dist = a µ { 1 - (f/2) ( A + B ( Sin µ ) / µ )                                ( µ , A & B employ parametric latitudes )

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

+ (f3/8) [ (1-A) { A (2-3A) ( µ Cot µ - µ2 Cot2 µ ) - A (1-A) (8/6) µ2 + B [ (2-3A) µ / Sin µ - (2-4A) µ2 ( Cos µ ) / Sin2 µ  + B ( µ2 / Sin2 µ + 1 ) ] }
+    B3 ( (2/3) ( Sin3 µ ) / µ - (1/2) ( Sin µ ) / µ ) + B2 (A/2-1) ( Sin µ Cos µ ) / µ
+ A2 B ( Sin µ Cos2 µ ) / ( 2µ ) + A B (1-A) Cos µ + (A2/2) (1-A/2) ( 1 + ( Sin µ Cos µ ) / µ ) ]

-Thus, we save 16 bytes !

-And if it's re-written:

D / ( a µ ) = 1 - (f/2) ( A + B ( Sin µ ) / µ )
+ (f2/4) [ ( 1 - A ) ( µ / Sin µ ) ( B + A Cos µ ) + (1/2) [ ( A2/2 ) ( 1 + (Sin µ  Cos µ) / µ ) - B2 (Sin µ  Cos µ) / µ ) ) ]
+ (f3/8) { ( 1 - A ) { B2 + A [ B Cos µ - ( 1 - A ) (8/6) µ2 ] + ( B + A Cos µ ) ( µ / Sin µ ) [ - ( 3 A - 2 ) + ( µ / Sin µ ) ( ( 3 A - 2 ) Cos µ + B ) ] }
+ [ B3 ( (2/3) Sin2 µ - 1/2 )  + ( A/2 - 1 ) B2 Cos µ ] (Sin µ) / µ ) - (A2/2) [ ( A/2 - 1 ) ( 1 + (Sin µ  Cos µ) / µ ) - B (Sin µ  Cos2 µ) / µ ] }

we save 13 extra bytes !!

Data Registers:   R00 = µ   R01 = A   R02 = B   R03 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 05  17 P-R  18 LASTX  19 298.25722  20 STO 04  21 ST+ 04  22 1/X  23 -  24 STO 06            25 /  26 R-P  27 X<> 01  28 RCL 05 29 P-R  30 RCL 06  31 /  32 R-P  33 RDN  34 +  35 2  36 /  37 ST- Y  38 RCL 05  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 05            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 05  63 R^  64 STO 01  65 -  66 /  67 ST- 02  68 +  69 ST- 05  70 X<> 01  71 SQRT  72 ASIN  73 ST+ X  74 D-R  75 STO 00            76 LASTX  77 1  78 P-R  79 STO 03  80 RDN  81 STO 09  82 /  83 STO 06  84 RCL 01 85 R^  86 *  87 RCL 02  88 +  89 *  90 STO 08  91 RCL 01  92 3  93 *  94 2  95 -  96 STO 07            97 RCL 03  98 *  99 RCL 02 100 + 101 RCL 06 102 * 103 RCL 07 104 - 105 * 106 RCL 02 107 RCL 03 108 * 109 RCL 00 110 ST+ X 111 X^2 112 RCL 05 113 * 114 3 115 / 116 - 117 RCL 01 118 * 119 + 120 RCL 02 121 X^2 122 STO 07 123 + 124 RCL 05 125 ST* 08 126 * 127 RCL 09         128 X^2 129 1.5 130 ST/ Y 131 FRC 132 - 133 RCL 02 134 * 135 RCL 01 136 2 137 ST- Y 138 / 139 STO 09 140 RCL 03 141 * 142 + 143 RCL 07 144 * 145 RCL 06 146 / 147 + 148 RCL 03 149 RCL 06 150 / 151 ST* 03 152 ST* 07 153 1 154 + 155 STO 05 156 RCL 09         157 * 158 RCL 02 159 RCL 03 160 * 161 - 162 RCL 01 163 X^2 164 2 165 / 166 ST* 05 167 * 168 - 169 RCL 04 170 / 171 RCL 05 172 RCL 07 173 - 174 2 175 / 176 + 177 RCL 08 178 + 179 RCL 04 180 / 181 RCL 02         182 RCL 06 183 / 184 RCL 01 185 + 186 - 187 RCL 04 188 / 189 RCL 00 190 ST* Y 191 + 192 6378.137 193 * 194 END

( 239 bytes / SIZE 010 )

 STACK INPUTS OUTPUTS T λ ( ° ' " ) / Z β ( ° ' " ) / Y λ' ( ° ' " ) / X β' ( ° ' " ) 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 = 12.9s---

-Here are a few root-mean-square & maximal errors in centimeters, found with free42.

 | λ'-λ | 90° 120° 150° 160° 165° 170° 171° 172° 173° 174° 175° 176° rms ~ 0 0.01 0.1 0.4 1 3 4 5.5 7.8 11.7 19 34 max ~ 0 0.06 0.8 3 7.7 27 39 52 82 136 205 439

-I've computed the errors with β & β' integers.
-So the maximal errors may be a little larger especially if | λ'-λ | = 173° , 174° ....

-We can obtain more accurate results with very large µ
if we compute µ with the same formulae used in Vinenty's method  ( Sin µ  &  Cos µ  give accurate µ )
-This program = 246 bytes instead of 239 but it's faster !

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

 01 LBL "AND"  02 DEG  03 HR  04 STO 06  05 X<> T  06 HMS-  07 HR  08 STO 01  09 X<>Y  10 HR  11 1  12 STO 05  13 P-R  14 LASTX  15 298.25722  16 STO 04            17 ST+ 04  18 1/X  19 -  20 STO 07  21 /  22 R-P  23 RCL 01  24 RCL Z  25 X<> 06  26 RCL 05  27 P-R  28 RCL 07 29 /  30 R-P  31 CLX  32 SIGN  33 P-R  34 X<>Y  35 STO 00  36 STO 02  37 ST+ 02  38 STO 03  39 X^2  40 STO 01  41 RDN  42 P-R  43 STO 07            44 RCL 06  45 RCL 05  46 P-R  47 ST* 03  48 ST* 07  49 RDN  50 ST* 00  51 ST* 02  52 ST* Y  53 X^2  54 ST+ 01  55 X<> 03  56 - 57 X^2  58 X<>Y  59 X^2  60 +  61 ST/ 01  62 ST/ 02  63 SQRT  64 STO 09  65 RCL 00  66 RCL 07  67 +  68 STO 03  69 R-P  70 X<>Y  71 D-R  72 STO 00            73 RCL 09  74 /  75 STO 06  76 RCL 01  77 RCL 02  78 RCL 03  79 ST* 01  80 *  81 -  82 X<> 01  83 RCL 02  84 - 85 STO 02  86 RCL 01  87 ST- 05  88 RCL 03  89 *  90 +  91 *  92 STO 08  93 RCL 01  94 3  95 *  96 2  97 -  98 STO 07            99 RCL 03 100 * 101 RCL 02 102 + 103 R^ 104 * 105 RCL 07 106 - 107 * 108 RCL 02 109 RCL 03 110 * 111 RCL 00 112 ST+ X 113 X^2 114 RCL 05 115 * 116 3 117 / 118 - 119 RCL 01 120 * 121 + 122 RCL 02 123 X^2 124 STO 07 125 + 126 RCL 05         127 ST* 08 128 * 129 RCL 09 130 X^2 131 1.5 132 ST/ Y 133 FRC 134 - 135 RCL 02 136 * 137 RCL 01 138 2 139 ST- Y 140 / 141 STO 09 142 RCL 03 143 * 144 + 145 RCL 07 146 * 147 RCL 06 148 / 149 + 150 RCL 03         151 RCL 06 152 / 153 ST* 03 154 ST* 07 155 1 156 + 157 STO 05 158 RCL 09 159 * 160 RCL 02 161 RCL 03 162 * 163 - 164 RCL 01 165 X^2 166 2 167 / 168 ST* 05 169 * 170 - 171 RCL 04 172 / 173 RCL 05 174 RCL 07 175 - 176 2 177 / 178 + 179 RCL 08         180 + 181 RCL 04 182 / 183 RCL 02 184 RCL 06 185 / 186 RCL 01 187 + 188 - 189 RCL 04 190 / 191 RCL 00 192 ST* Y 193 + 194 6378.137 195 * 196 END

( 246 bytes / SIZE 010 )

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

Example1:    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 = 12.4s---

Example2:    ( 0° , 0° )  ( 180° , 0°00'01" )

0       ENTER^  ENTER^
180     ENTER^
0.0001  XEQ "DGT"   >>>>   D= 20003.90075 km

( exact result = 20003.90074.... )

b)  Andoyer2 + A few Terms of order 3 ( Without Calculating Parametric Latitudes )

b1)  Program#1

-If less accurate results are acceptable, we can improve a little Andoyer formula of order 2 after adding some terms of order 3 like  µ^3 ....
-These routines don't calculate parametric latitudes

Formulae:

-With   H = ( λ' - λ )/2  ;   F = ( β + β' )/2  ;  G = ( β' - β )/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 )                    calculated with geodetic latitudes.
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 ( 1 - A ) (8/6) µ3 + ( B + A Cos µ ) ( µ3 / Sin2 µ ) ( B + ( 3 A - 2 ) Cos µ ) ] ( 1 - A ) }

Data Registers:   R00 = µ  R01 = A  R02 = B  R03 to R07:  temp
Flags: /
Subroutines: /

 01 LBL "DGT"  02 DEG  03 HR  04 X<> T  05 HMS-  06 HR  07 2  08 /  09 SIN  10 X^2  11 RDN  12 HR  13 +  14 2  15 /  16 ST- Y  17 1  18 STO 04            19 P-R 20 X^2  21 STO 02  22 X<>Y  23 X^2  24 STO 01  25 RDN  26 X<>Y  27 RCL 04            28 P-R  29 X^2  30 ST* 01  31 RDN  32 X^2  33 ST* 02  34 STO T  35 -  36 *  37 +  38 ST/ 02 39 RCL 02  40 RCL 01  41 RCL 04  42 R^  43 STO 01  44 -  45 /  46 ST- 02  47 +  48 ST- 04  49 X<> 01  50 SQRT  51 ASIN  52 ST+ X  53 D-R  54 STO 00            55 LASTX  56 1  57 P-R 58 STO 03  59 RDN  60 /  61 STO 06  62 X^2  63 RCL 01  64 R^  65 *  66 RCL 02  67 +  68 STO 05  69 ST* Y  70 RCL 03  71 RCL 04            72 ST+ X  73 *  74 -  75 *  76 RCL 00 77 ST+ X  78 X^2  79 RCL 01  80 *  81 RCL 04  82 *  83 3  84 /  85 -  86 596.51444  87 STO 07  88 /  89 RCL 05            90 RCL 02  91 3  92 *  93 RCL 06  94 ST/ 03  95 ST* Z 96 /  97 STO 05  98 ST+ X  99 + 100 + 101 RCL 04 102 * 103 RCL 03 104 RCL 02 105 X^2 106 ST+ X 107 RCL 01 108 X^2 109 STO 03         110 - 111 * 112 RCL 03 113 - 114 15 115 * 116 RCL 01 117 4 118 ST/ Z 119 * 120 + 121 + 122 RCL 07 123 / 124 RCL 01 125 RCL 05 126 + 127 - 128 RCL 07         129 / 130 RCL 00 131 ST* Y 132 + 133 6378.137 134 * 135 END

( 178 bytes / SIZE 008 )

 STACK INPUTS OUTPUTS T λ ( ° ' " ) / Z β ( ° ' " ) / Y λ' ( ° ' " ) / X β' ( ° ' " ) 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 "DGT"   >>>>   D =  12138.68424 km                                          ---Execution time = 8.6s---

Notes:

-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 λ = 0°  β = 5°  λ' = 175°  β' = 3°  this program gives  18968.13253 km  whereas Vincenty's formula gives 18968.13070 km
-So, an error of 1.83 m: slightly less than the similar program listed in "Terrestrial Geodesic Distance(2) for the HP41" !

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

-The second row = maximum error of Andoyer's formula of order 2 without any term of order 3
-Third row = maximum error with "DGT" 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.

 | λ'-λ | 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.17 0.24 0.34 0.52 0.87 0.98 1.11 1.28 1.51 1.82 3.56 And3 ~ 0 0.008 0.01 0.03 0.08 0.27 0.39 0.52 0.82 1.36 2.05 4.39

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

-Surprisingly, if λ' - λ > 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 !

-Here are also a few root-mean-square errors in centimeters.

 | λ'- λ | 20° 45° 90° 120° 150° 160° 170° 171° 175° And2afew3 4 5 5.7 5.4 6.7 10.3 19.3 20.9 31.3

b2)  Program#2

-The following program employs the formula:

Dist ~   a µ { 1 - (f/2) ( A + 3.B (Sin µ)/µ )
+ (f2/4) [ A ( 4 + µ Cot µ ) + B ( 6(Sin µ)/µ + µ / Sin µ ) - A2 ( 15 + 15 (Sin µ) (Cos µ)/µ + 4 µ Cot µ ) / 4 - A.B ( 6(Sin µ)/µ + µ / Sin µ ) + ( B2/2 ) 15 (Sin µ) (Cos µ)/µ ]
+ (f3/8) [ ( -A ( 1 - A ) (8/6) µ2 + ( B + A Cos µ ) ( µ2 / Sin2 µ ) ( B + ( 3 A - 2 ) Cos µ ) ) ( 1 - A )
+ B2 (3A-3) + B ( - 3A(1-A) Cos µ + (1-A) (7A-2) µ/Sin µ ) + A(1-A) (7A-2) µ Cot µ ] }

Data Registers:   R00 = µ  R01 = A  R02 = B  R03 to R08:  temp
Flags: /
Subroutines: /

 01 LBL "DGT"  02 DEG  03 HR  04 X<> T  05 HMS-  06 HR  07 2  08 /  09 SIN  10 X^2  11 RDN  12 HR  13 +  14 2  15 /  16 ST- Y  17 1  18 STO 04            19 P-R  20 X^2  21 STO 02 22 X<>Y  23 X^2  24 STO 01  25 RDN  26 X<>Y  27 RCL 04  28 P-R  29 X^2  30 ST* 01  31 RDN  32 X^2  33 ST* 02  34 STO T  35 -  36 *  37 +  38 ST/ 02  39 RCL 02            40 RCL 01  41 RCL 04  42 R^ 43 STO 01  44 -  45 /  46 ST- 02  47 +  48 ST- 04  49 X<> 01  50 SQRT  51 ASIN  52 ST+ X  53 D-R  54 STO 00  55 LASTX  56 1  57 P-R  58 STO 03  59 STO 08            60 RDN  61 /  62 STO 06  63 ST/ 08 64 RCL 01  65 R^  66 *  67 RCL 02  68 +  69 STO 05  70 RCL 03  71 RCL 04  72 ST+ X  73 *  74 -  75 *  76 RCL 01  77 7  78 *  79 2  80 -  81 +  82 RCL 06            83 *  84 RCL 02 85 3  86 *  87 STO 03  88 -  89 RCL 05  90 *  91 RCL 00  92 ST+ X  93 X^2  94 RCL 01  95 *  96 RCL 04  97 *  98 3  99 / 100 - 101 596.51444 102 STO 07       103 / 104 RCL 05 105 RCL 03 106 RCL 06 107 ST* Z 108 / 109 STO 05 110 ST+ X 111 + 112 + 113 RCL 04 114 * 115 RCL 08 116 RCL 02 117 X^2 118 ST+ X 119 RCL 01 120 X^2 121 STO 03       122 - 123 * 124 RCL 03 125 - 126 15 127 * 128 RCL 01 129 4 130 ST/ Z 131 * 132 + 133 + 134 RCL 07 135 / 136 RCL 01 137 RCL 05 138 + 139 - 140 RCL 07 141 / 142 RCL 00       143 ST* Y 144 + 145 6378.137 146 * 147 END

( 189 bytes / SIZE 009 )

 STACK INPUTS OUTPUTS T λ ( ° ' " ) / Z β ( ° ' " ) / Y λ' ( ° ' " ) / X β' ( ° ' " ) 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 "DGT"   >>>>   D =  12138.68424 km                                          ---Execution time = 9s---

-Here are a few root-mean-square & maximal errors in centimeters.

 | λ'-λ | 20° 45° 90° 120° 150° 160° 170° 171° 175° rms 3.7 4.2 5 5.8 6.3 6.3 6.8 7.2 19.0 max 9.5 12.4 12.9 13.1 14.1 14.1 26.8 37.6 201

b3)  Program#3

-Yet other added terms ( computed with least-squares approximations )

Dist ~   a µ { 1 - (f/2) ( A + 3.B (Sin µ)/µ )
+ (f2/4) [ A ( 4 + µ Cot µ ) + B ( 6(Sin µ)/µ + µ / Sin µ ) - A2 ( 15 + 15 (Sin µ) (Cos µ)/µ + 4 µ Cot µ ) / 4 - A.B ( 6(Sin µ)/µ + µ / Sin µ ) + ( B2/2 ) 15 (Sin µ) (Cos µ)/µ ]
+ (f3/8) [ ( -A ( 1 - A ) (8/6) µ2 + ( B + A Cos µ ) ( µ2 / Sin2 µ ) ( B + ( 3 A - 2 ) Cos µ ) ) ( 1 - A )
+ B2 (3A-3) + B ( - 3A(1-A) Cos µ + (1-A) (7A-2) µ/Sin µ ) + A(1-A) (7A-2) µ Cot µ
+ ( 1-A ) ( A ( 7A-2 )sqrt(7) + ( 11B2 - 8A2 ) (Sin µ) (Cos µ)/µ  ) ] }

Data Registers:   R00 = µ  R01 = A  R02 to R08:  temp
Flags: /
Subroutines: /

 01 LBL "DGT"  02 DEG  03 HR  04 X<> T  05 HMS-  06 HR  07 2  08 /  09 SIN  10 X^2  11 RDN  12 HR  13 +  14 2  15 /  16 ST- Y  17 1  18 STO 04            19 P-R  20 X^2  21 STO 02  22 X<>Y 23 X^2  24 STO 01  25 RDN  26 X<>Y  27 RCL 04  28 P-R  29 X^2  30 ST* 01  31 RDN  32 X^2  33 ST* 02  34 STO T  35 -  36 *  37 +  38 ST/ 02  39 RCL 02            40 RCL 01  41 RCL 04  42 R^  43 STO 01  44 - 45 /  46 ST- 02  47 +  48 X<> 01  49 SQRT  50 ASIN  51 ST+ X  52 D-R  53 STO 00  54 LASTX  55 RCL 04  56 P-R  57 STO 03  58 RDN  59 /  60 STO 06            61 RCL 01  62 ST- 04  63 R^  64 *  65 RCL 02  66 + 67 STO 05  68 RCL 03  69 RCL 04  70 ST+ X  71 *  72 -  73 *  74 RCL 01  75 7  76 *  77 2  78 -  79 STO 07  80 +  81 *  82 RCL 02            83 3  84 *  85 STO 08  86 -  87 RCL 05  88 * 89 RCL 07  90 7  91 SQRT  92 *  93 RCL 00  94 ST+ X  95 X^2  96 RCL 04  97 *  98 3  99 / 100 - 101 RCL 01         102 * 103 + 104 RCL 08 105 RCL 02 106 ST* Y 107 X^2 108 STO 07 109 RCL 01 110 X^2 111 STO 02 112 - 113 ST+ 07 114 8 115 * 116 + 117 RCL 03 118 RCL 06         119 ST* 05 120 ST/ 08 121 / 122 ST* 07 123 * 124 + 125 596.51444 126 STO 03 127 / 128 RCL 05 129 + 130 RCL 08 131 ST+ X 132 + 133 RCL 04 134 * 135 RCL 01 136 ST+ 08 137 RCL 07         138 RCL 02 139 - 140 ST+ Y 141 4 142 ST* Z 143 / 144 - 145 + 146 RCL 03 147 / 148 RCL 08 149 - 150 RCL 03 151 / 152 RCL 00 153 ST* Y 154 + 155 6378.137 156 * 157 END

( 202 bytes / SIZE 009 )

 STACK INPUTS OUTPUTS T λ ( ° ' " ) / Z β ( ° ' " ) / Y λ' ( ° ' " ) / X β' ( ° ' " ) 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 "DGT"   >>>>   D =  12138.68429 km                                          ---Execution time = 9.4s---

Notes:

-If you want to save B in register R02, replace lines 111 & 138 by STO 09 & RCL 09
-Here are a few root-mean-square & maximal errors in centimeters.

 | λ'-λ | 20° 45° 90° 120° 150° 160° 170° 171° 175° rms 3.1 2.5 2.2 2.8 3.6 3.9 4.8 5.3 17.9 max 7.4 6.3 5.9 7.6 8.9 9.1 24.4 35 198

-Like with "AND", we can obtain more accurate results with very large µ
if we compute µ with the same formulae used in Vinenty's method  ( Sin µ  &  Cos µ  give accurate µ )
-This program = 209 bytes instead of 202 but it's faster !

Data Registers:   R00 = µ  R01 = A  R02 to R08:  temp
Flags: /
Subroutines: /

 01 LBL "DGT"  02 DEG  03 HR  04 STO 06  05 X<> T  06 HMS-  07 HR  08 X<>Y  09 HR  10 1  11 STO 04  12 P-R  13 X<>Y  14 STO 00            15 STO 02  16 ST+ 02  17 STO 03  18 X^2  19 STO 01  20 RDN  21 P-R  22 STO 05  23 RCL 06 24 RCL 04  25 P-R  26 ST* 03  27 ST* 05  28 RDN  29 ST* 00  30 ST* 02  31 ST* Y  32 X^2  33 ST+ 01  34 X<> 03  35 -  36 X^2  37 X<>Y  38 X^2  39 +  40 ST/ 01  41 ST/ 02  42 SQRT  43 STO 06            44 RCL 00  45 RCL 05  46 + 47 STO 03  48 R-P  49 X<>Y  50 D-R  51 STO 00  52 RCL 06  53 /  54 STO 06  55 RCL 01  56 RCL 02  57 RCL 03  58 ST* 01  59 *  60 -  61 X<> 01  62 RCL 02  63 -  64 STO 02            65 RCL 01  66 ST- 04  67 RCL 03  68 *  69 + 70 STO 05  71 LASTX  72 RCL 03  73 -  74 ST+ X  75 +  76 *  77 RCL 01  78 7  79 *  80 2  81 -  82 STO 07  83 +  84 *  85 RCL 02  86 3  87 *  88 STO 08            89 -  90 RCL 05  91 *  92 RCL 07 93 7  94 SQRT  95 *  96 RCL 00  97 ST+ X  98 X^2  99 RCL 04 100 * 101 3 102 / 103 - 104 RCL 01         105 * 106 + 107 RCL 08 108 RCL 02 109 ST* Y 110 X^2 111 STO 07 112 RCL 01 113 X^2 114 STO 02 115 - 116 ST+ 07 117 8 118 * 119 + 120 RCL 03 121 RCL 06 122 ST* 05 123 ST/ 08 124 / 125 ST* 07 126 * 127 + 128 596.51444 129 STO 03         130 / 131 RCL 05 132 + 133 RCL 08 134 ST+ X 135 + 136 RCL 04 137 * 138 RCL 01 139 ST+ 08 140 RCL 07 141 RCL 02 142 - 143 ST+ Y 144 4 145 ST* Z 146 / 147 - 148 + 149 RCL 03         150 / 151 RCL 08 152 - 153 RCL 03 154 / 155 RCL 00 156 ST* Y 157 + 158 6378.137 159 * 160 END

( 209 bytes / SIZE 009 )

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

Where  D is the geodesic distance on the WGS84 ellipsoid

Example1:    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 "DGT"   >>>>   D =  12138.68429 km                                          ---Execution time = 9s---

Example2:    ( 0° , 0° )  ( 180° , 0°00'01" )

0       ENTER^  ENTER^
180     ENTER^
0.0001  XEQ "DGT"   >>>>   D= 20003.90072 km

( exact result = 20003.90074.... )

b4)  Program#4

Dist ~   a µ { 1 - (f/2) ( A + 3.B (Sin µ)/µ )
+ (f2/4) [ A ( 4 + µ Cot µ ) + B ( 6(Sin µ)/µ + µ / Sin µ ) - A2 ( 15 + 15 (Sin µ) (Cos µ)/µ + 4 µ Cot µ ) / 4 - A.B ( 6(Sin µ)/µ + µ / Sin µ ) + ( B2/2 ) 15 (Sin µ) (Cos µ)/µ ]
+ (f3/8) [ ( 1 - A ) { -A ( 1 - A ) (4/3) µ2 + 18 A B (Sin µ)/µ + A (17.5A-5) +
( B + A Cos µ ) [ -3 B + ( µ / Sin µ ) ( 7 A - 2 ) + ( µ2 / Sin2 µ ) ( B + ( 3 A - 2 ) Cos µ ) ] }
+ ( 11 - 15 A ) A2 (Sin µ)(Cos µ)/µ - 4 B (Sin µ)/µ + 19 A2 B (Sin µ)(Cos2 µ)/µ
+ ( 33A-25 ) B2 (Sin µ)(Cos µ)/µ + B3 (Sin µ) ( 24 Sin2 µ - 18 ) / µ  ]
+ (f4/16) ( 1 - A ) ( µ3 / Sin3 µ ) ( B + A Cos µ ) [ B + ( A - 1 ) Cos µ ]
+ (f5/32) ( 1 - A ) ( µ4 / Sin4 µ ) ( B + A Cos µ ) [ -2.B + ( -4.A - 1 ) Cos µ ]

Data Registers:   R00 = µ  R01 = A  R02 = B  R03 to R12:  temp
Flags: /
Subroutines: /

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

( 274 bytes / SIZE 013 )

 STACK INPUTS OUTPUTS T λ ( ° ' " ) / Z β ( ° ' " ) / Y λ' ( ° ' " ) / X β' ( ° ' " ) 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 "DGT"   >>>>   D =  12138.68432 km                                          ---Execution time = 11.4s---

-Here are a few root-mean-square & maximal errors in centimeters.

 | λ'-λ | 20° 50° 90° 120° 150° 160° 164° 168° 170° rms 0.8 0.6 0.8 0.9 0.9 1.0 1.1 1.9 3.0 max 2.2 2.2 2.2 2.2 2.3 2.4 3.1 8.5 16.2

c)  Vincenty's Formula

-A similar program is listed in "Terrestrial Geodesic Distance(1)" but the following one employs the complete formula.
-Not very useful with an HP41 which works with 10 digit numbers, but interesting with free42.

-Cf reference [5] :  "VCT" calculates A & B with shorter formulae...

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

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

( 297 bytes / SIZE 015 )

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

where   Az1  =  forward azimuth  ( from P1 to P2 )
and     Az2  =  back azimuth  ( from P2 to P1 )

Example:

151.12178    ENTER^
-33.51411    ENTER^
-116.51504   ENTER^
33.21224    XEQ "VCT"  >>>>      D   =  12138.68432 km                                          ---Execution time = 42s---
RDN    Az1  =  62°19'38"5790
RDN    Az2  = -118°18'05"388

Notes:

-The accuracy is excellent, provided the 2 points are not nearly antipodal.

-Line 147 X#Y?  may produce an infinte loop
-So, it's better to replace it by the M-Code routine  X#Y??
-Or replace lines 147-148 by   -  ABS  E-7  X<Y?  GTO 01

-In fact, WGS84 flattening is  1 / 298.257223563
-So, we can replace line 23 with  298.257223563
-It doesn't really change the results with an HP41, but with free42 or DM42, it yields:

12138.6843157568...     62°19'38"5789628...    -118°18'05"3880993...

-If you don't want to calculate the Azimuths, delete lines  224 to 255.

d)  Geocentric Latitudes

d1)  Program(s)#1

-If the latitudes are geocentric, we can use the programs listed in paragraphs 1°)a) & 1°)c):

-With "AND", simply replace lines 25 & 31 by  *  instead of  /  ( lines 21 & 29 in the second "AND" program )
-Likewise in "VCT"  replace lines 28 & 39 by  *  instead of  /

-But if we prefer a formula without computing parametric latitudes, we can use the following one.

Dist ~   a µ { 1 + (f/2) ( -A + B (Sin µ)/µ )
+ (f2/4) [ A ( -4 + µ Cot µ ) + B ( 2 (Sin µ)/µ + µ / Sin µ ) + (A2/4) ( 17 + (Sin µ) (Cos µ)/µ - 4 µ Cot µ ) - A.B ( 2 (Sin µ)/µ + µ / Sin µ ) - ( B2/2 ) (Sin µ) (Cos µ)/µ ]
+ (f3/8)  ( 1 - A ) [ ( -A ( 1 - A ) (4/3) µ2 + ( B + A Cos µ ) ( µ2 / Sin2 µ ) ( B + ( 3 A - 2 ) Cos µ ) )
+  ( 7.6 B2 + 7.6 A.B Cos µ ) + (µ/Sin µ) ( B (5.6 - 12.2 A) + A (Cos µ) (7.6-15.2 A) ) ] }

Data Registers:   R00 = µ  R01 = A  R02 = B  R03 to R06:  temp
Flags: /
Subroutines: /

 01 LBL "DGC"  02 DEG  03 HR  04 X<> T  05 HMS-  06 HR  07 2  08 /  09 SIN  10 X^2  11 RDN  12 HR  13 +  14 2  15 /  16 ST- Y  17 1  18 STO 04            19 P-R  20 X^2  21 STO 02 22 X<>Y  23 X^2  24 STO 01  25 RDN  26 X<>Y  27 RCL 04  28 P-R  29 X^2  30 ST* 01  31 RDN  32 X^2  33 ST* 02  34 STO T  35 -  36 *  37 +  38 ST/ 02  39 RCL 02            40 RCL 01  41 RCL 04  42 R^ 43 STO 01  44 -  45 /  46 ST- 02  47 +  48 X<> 01  49 SQRT  50 ASIN  51 ST+ X  52 D-R  53 STO 00  54 LASTX  55 RCL 04  56 P-R  57 STO 03  58 RDN  59 /  60 RCL 01            61 ST- 04  62 R^  63 * 64 RCL 02  65 +  66 STO 05  67 RCL 03  68 RCL 04  69 ST+ X  70 *  71 -  72 *  73 RCL 05  74 ST* Y  75 7.6  76 *  77 STO 06  78 RCL 02            79 ST+ X  80 -  81 ST+ Y  82 ST+ X  83 RCL 02  84 ST* 06 85 +  86 RCL 01  87 *  88 -  89 *  90 RCL 06  91 +  92 RCL 00  93 ST+ X  94 X^2  95 RCL 01  96 *  97 RCL 04  98 *  99 3 100 / 101 - 102 596.51444 103 STO 06         104 / 105 RCL 05 106 RCL 02 107 R^ 108 ST/ 03 109 ST* Z 110 / 111 STO 05 112 ST+ X 113 + 114 + 115 RCL 04 116 * 117 17 118 RCL 01 119 X^2 120 ST* Y 121 RCL 02         122 X^2 123 ST+ X 124 - 125 RCL 03 126 * 127 + 128 RCL 01 129 ST- 05 130 4 131 ST/ Z 132 * 133 - 134 + 135 RCL 06 136 / 137 RCL 05 138 + 139 RCL 06 140 / 141 RCL 00         142 ST* Y 143 + 144 6378.137 145 * 146 END

( 197 bytes / SIZE 007 )

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

Where  D is the geodesic distance on the WGS84 ellipsoid

Example:     Long = 151°12'17"8 E , Lat = 33°51'41"1 S
Long = 116°51'50"4 W ,  Lat = 33°21'22"4 N

151.12178  ENTER^
-33.51411   ENTER^
-116.51504  ENTER^
33.21224   XEQ "DGC"   >>>>   D =  12157.20894 km                                          ---Execution time = 9.1s---

-Exact result = 12157.208910....

Notes:

-Here are a few root-mean-square & maximal errors in centimeters.

 | λ'-λ | 20° 45° 90° 120° 150° 160° 164° 168° 171° 175° 176° rms 3.5 3.2 3.6 3.7 3.0 2.9 3.2 4.2 6.3 20.8 36.3 max 8.4 7.1 8.4 9.9 8.4 10.3 14.5 23.0 36.7 206 439

-The following program computes µ with Sin µ & Cos µ to get a better precision with nearly antipodal points:

Data Registers:   R00 = µ  R01 = A  R02 = B  R03 to R06:  temp
Flags: /
Subroutines: /

 01 LBL "DGC"  02 DEG  03 HR  04 STO 06  05 X<> T  06 HMS-  07 HR  08 X<>Y  09 HR  10 1  11 STO 04            12 P-R  13 X<>Y  14 STO 00  15 STO 02  16 ST+ 02  17 STO 03  18 X^2  19 STO 01  20 RDN  21 P-R 22 STO 05  23 RCL 06  24 RCL 04            25 P-R  26 ST* 03  27 ST* 05  28 RDN  29 ST* 00  30 ST* 02  31 ST* Y  32 X^2  33 ST+ 01  34 X<> 03  35 -  36 X^2  37 X<>Y  38 X^2  39 +  40 ST/ 01  41 ST/ 02  42 SQRT 43 STO 06  44 RCL 00            45 RCL 05  46 +  47 STO 03  48 R-P  49 X<>Y  50 D-R  51 STO 00  52 RCL 06  53 /  54 RCL 01  55 RCL 02  56 RCL 03  57 ST* 01  58 *  59 -  60 X<> 01  61 RCL 02  62 -  63 STO 02 64 RCL 01  65 ST- 04  66 RCL 03            67 *  68 +  69 STO 05  70 LASTX  71 RCL 03  72 -  73 ST+ X  74 +  75 *  76 RCL 05  77 ST* Y  78 7.6  79 *  80 STO 06  81 RCL 02  82 ST+ X  83 -  84 ST+ Y 85 ST+ X  86 RCL 02            87 ST* 06  88 +  89 RCL 01  90 *  91 -  92 *  93 RCL 06  94 +  95 RCL 00  96 ST+ X  97 X^2  98 RCL 01  99 * 100 RCL 04 101 * 102 3 103 / 104 - 105 596.51444 106 STO 06 107 / 108 RCL 05 109 RCL 02         110 R^ 111 ST/ 03 112 ST* Z 113 / 114 STO 05 115 ST+ X 116 + 117 + 118 RCL 04 119 * 120 17 121 RCL 01 122 X^2 123 ST* Y 124 RCL 02 125 X^2 126 ST+ X 127 - 128 RCL 03 129 * 130 + 131 RCL 01         132 ST- 05 133 4 134 ST/ Z 135 * 136 - 137 + 138 RCL 06 139 / 140 RCL 05 141 + 142 RCL 06 143 / 144 RCL 00 145 ST* Y 146 + 147 6378.137 148 * 149 END

( 204 bytes / SIZE 007 )

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

Where  D is the geodesic distance on the WGS84 ellipsoid

Example1:     Long = 151°12'17"8 E , Lat = 33°51'41"1 S
Long = 116°51'50"4 W ,  Lat = 33°21'22"4 N

151.12178  ENTER^
-33.51411   ENTER^
-116.51504  ENTER^
33.21224   XEQ "DGC"   >>>>   D =  12157.20894 km                                          ---Execution time = 8.6s---

Example2:    ( 0° , 0° )  ( 180° , 0°00'01" )

0       ENTER^  ENTER^
180     ENTER^
0.0001  XEQ "DGT"   >>>>   D= 20003.90051 km

( exact result = 20003.900536... )

d2)  Program#2

Formula:

D/(a.µ) ~   1 + (f/2) ( -A + B (Sin µ)/µ )
+ (f2/4) [ A ( -4 + µ Cot µ ) + B ( 2 (Sin µ)/µ + µ / Sin µ ) + (A2/4) ( 17 + (Sin µ) (Cos µ)/µ - 4 µ Cot µ ) - A.B ( 2 (Sin µ)/µ + µ / Sin µ ) - ( B2/2 ) (Sin µ) (Cos µ)/µ ]
+ (f3/8) ( 1 - A ) { -A ( 1 - A ) (4/3) µ2 +  A ( 22.A - 11 )  + ( B + A Cos µ ) [ 7 B + ( µ / Sin µ ) ( 6 - 12.A) + ( µ2 / Sin2 µ ) ( B + ( 3 A - 2 ) Cos µ ) ) ] }

Data Registers:   R00 = µ  R01 = A  R02 = B  R03 to R07:  temp
Flags: /
Subroutines: /

 01 LBL "DGC"  02 DEG  03 HR  04 STO 06  05 X<> T  06 HMS-  07 HR  08 X<>Y  09 HR  10 11  11 STO 07  12 SIGN  13 STO 04            14 P-R  15 X<>Y  16 STO 00  17 STO 02  18 ST+ 02  19 STO 03  20 X^2  21 STO 01 22 RDN  23 P-R  24 STO 05  25 RCL 06  26 RCL 04            27 P-R  28 ST* 03  29 ST* 05  30 RDN  31 ST* 00  32 ST* 02  33 ST* Y  34 X^2  35 ST+ 01  36 X<> 03  37 -  38 X^2  39 X<>Y  40 X^2  41 +  42 ST/ 01 43 ST/ 02  44 SQRT  45 STO 06  46 RCL 00  47 RCL 05            48 +  49 STO 03  50 R-P  51 X<>Y  52 D-R  53 STO 00  54 RCL 06  55 /  56 RCL 01  57 RCL 02  58 RCL 03  59 ST* 01  60 *  61 -  62 X<> 01  63 RCL 02 64 -  65 STO 02  66 RCL 01  67 ST- 04  68 RCL 03            69 *  70 +  71 STO 05  72 LASTX  73 RCL 03  74 -  75 ST+ X  76 +  77 *  78 RCL 01  79 ST* 07  80 RCL 04  81 -  82 ST* 07  83 6  84 * 85 -  86 *  87 RCL 02  88 7  89 *  90 +  91 RCL 05  92 *  93 RCL 07  94 +  95 RCL 00  96 ST+ X  97 X^2  98 RCL 01  99 * 100 RCL 04         101 * 102 3 103 / 104 - 105 596.51444 106 STO 06 107 / 108 RCL 05 109 RCL 02 110 R^ 111 ST/ 03 112 ST* Z 113 / 114 STO 05 115 ST+ X 116 + 117 + 118 RCL 04         119 * 120 17 121 RCL 01 122 X^2 123 ST* Y 124 RCL 02 125 X^2 126 ST+ X 127 - 128 RCL 03 129 * 130 + 131 RCL 01 132 ST- 05 133 4 134 ST/ Z 135 * 136 - 137 + 138 RCL 06 139 / 140 RCL 05         141 + 142 RCL 06 143 / 144 RCL 00 145 ST* Y 146 + 147 6378.137 148 * 149 END

( 200 bytes / SIZE 008 )

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

Where  D is the geodesic distance on the WGS84 ellipsoid

Example:     Long = 151°12'17"8 E , Lat = 33°51'41"1 S
Long = 116°51'50"4 W ,  Lat = 33°21'22"4 N

151.12178  ENTER^
-33.51411   ENTER^
-116.51504  ENTER^
33.21224   XEQ "DGC"   >>>>   D =  12157.20894 km                                          ---Execution time = 8.6 s---

-Here are a few root-mean-square & maximal errors in centimeters.

 | λ'-λ | 20° 45° 90° 120° 160° 164° 168° 171° 175° 176° rms 3.3 2.4 2.7 3.0 2.9 3.5 5.2 8.4 27.3 45.7 max 7.2 5.5 4.8 5.9 6.4 10.8 24.0 50.0 220 456

-If we replace line 10 with 12 ( instead of 11 ) and line 88 with 8 ( instead of 7 ) the errors become:

 | λ'-λ | 20° 45° 90° 120° 160° 164° 168° 171° 175° 176° rms 3.2 2.5 2.9 3.2 3.0 3.4 4.9 8.2 27.1 45.5 max 7.0 5.5 5.5 6.8 6.6 9.3 22.3 48.3 219 455

-Perhaps a little better precision.

d3)  Program#3

Formula:

D/(a.µ) ~   1 + (f/2) ( -A + B (Sin µ)/µ )
+ (f2/4) [ A ( -4 + µ Cot µ ) + B ( 2 (Sin µ)/µ + µ / Sin µ ) + (A2/4) ( 17 + (Sin µ) (Cos µ)/µ - 4 µ Cot µ ) - A.B ( 2 (Sin µ)/µ + µ / Sin µ ) - ( B2/2 ) (Sin µ) (Cos µ)/µ ]
+ (f3/8) { ( 1 - A ) [ ( -A ( 1 - A ) (4/3) µ2 +  A ( 21.A - 13 )  + ( B + A Cos µ ) [ ( µ / Sin µ ) ( 6 - 13.A) + 5.B + ( µ2 / Sin2 µ ) ( B + ( 3 A - 2 ) Cos µ ) ) ]
+  ( 3.A2 - 9.B2 ) ( Sin µ ) ( Cos µ ) / µ ] + B3 [ ( Sin µ ) / µ + 2 ( Sin3 µ ) / µ ] }

Data Registers:   R00 = µ  R01 = A  R02 = -A + B ( Sin µ ) / µ   R03 to R09:  temp
Flags: /
Subroutines: /

 01 LBL "DGC"  02 DEG  03 HR  04 STO 06  05 X<> T  06 HMS-  07 HR  08 X<>Y  09 HR  10 1  11 STO 04  12 P-R  13 X<>Y  14 STO 00  15 STO 02  16 ST+ 02  17 STO 03  18 X^2  19 STO 01            20 RDN  21 P-R  22 STO 05  23 RCL 06  24 RCL 04 25 P-R  26 ST* 03  27 ST* 05  28 RDN  29 ST* 00  30 ST* 02  31 ST* Y  32 X^2  33 ST+ 01  34 X<> 03  35 -  36 X^2  37 X<>Y  38 X^2  39 +  40 ST/ 01  41 ST/ 02  42 STO 08            43 ST+ 08  44 SQRT  45 STO 06  46 RCL 00  47 RCL 05  48 + 49 STO 03  50 R-P  51 X<>Y  52 D-R  53 STO 00  54 RCL 06  55 /  56 STO 06  57 RCL 01  58 RCL 02  59 RCL 03  60 ST* 01  61 *  62 -  63 X<> 01  64 RCL 02  65 -  66 STO 02            67 RCL 01  68 ST- 04  69 RCL 03  70 *  71 +  72 STO 05 73 LASTX  74 RCL 03  75 -  76 ST+ X  77 +  78 *  79 RCL 01  80 13  81 *  82 STO 07  83 -  84 6  85 +  86 X<>Y  87 ST/ 03  88 *  89 RCL 02            90 5  91 *  92 +  93 RCL 05  94 *  95 RCL 07  96 - 97 RCL 01  98 X^2  99 STO 07 100 21 101 * 102 + 103 RCL 07 104 RCL 02 105 ST* 08 106 ST+ 08 107 X^2 108 ST* 08 109 ST- Y 110 - 111 STO 09 112 LASTX 113 - 114 3 115 * 116 RCL 03         117 ST* 09 118 * 119 + 120 RCL 00 121 ST+ X 122 X^2 123 RCL 01 124 * 125 RCL 04 126 * 127 3 128 / 129 - 130 RCL 04 131 * 132 RCL 08 133 RCL 06 134 ST/ 02 135 ST* 05 136 / 137 + 138 596.51444 139 STO 06         140 / 141 RCL 05 142 RCL 02 143 ST+ X 144 + 145 RCL 04 146 * 147 + 148 RCL 07 149 RCL 09 150 + 151 RCL 01 152 ST- 02 153 RCL 04 154 * 155 4 156 ST/ Z 157 * 158 - 159 + 160 RCL 06 161 / 162 RCL 02 163 + 164 RCL 06         165 / 166 RCL 00 167 ST* Y 168 + 169 6378.137 170 * 171 END

( 225 bytes / SIZE 010 )

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

Where  D is the geodesic distance on the WGS84 ellipsoid

Example:     Long = 151°12'17"8 E , Lat = 33°51'41"1 S
Long = -116°51'50"4 W ,  Lat = 33°21'22"4 N

151.12178  ENTER^
-33.51411   ENTER^
-116.51504  ENTER^
33.21224   XEQ "DGC"   >>>>   D =  12157.20890 km                                          ---Execution time = 9.4s---

-Here are a few root-mean-square & maximal errors in centimeters.

 | λ'-λ | 20° 45° 90° 120° 160° 164° 168° 171° 175° 176° rms 1.3 1.1 1.2 1.3 1.9 2.1 3.0 5.4 22.9 41.0 max 3.8 3.8 4.0 4.2 4.7 7.0 17.0 41.2 208 457

d4)  Program(s)#4

Formula:

D/(a.µ) ~   1 + (f/2) ( -A + B (Sin µ)/µ )
+ (f2/4) [ A ( -4 + µ Cot µ ) + B ( 2 (Sin µ)/µ + µ / Sin µ ) + (A2/4) ( 17 + (Sin µ) (Cos µ)/µ - 4 µ Cot µ ) - A.B ( 2 (Sin µ)/µ + µ / Sin µ ) - ( B2/2 ) (Sin µ) (Cos µ)/µ ]
+ (f3/8) { ( 1 - A ) [ ( -A ( 1 - A ) (4/3) µ2 +  A ( 21.A - 13 )  + ( B + A Cos µ ) [ ( µ / Sin µ ) ( 6 - 13.A) + 5.B + ( µ2 / Sin2 µ ) ( B + ( 3 A - 2 ) Cos µ ) ) ]
+  ( 3.A2 - 9.B2 ) ( Sin µ ) ( Cos µ ) / µ ] + B3 [ ( Sin µ ) / µ + 2 ( Sin3 µ ) / µ ] }
+ (f4/16) ( 1 - A ) ( µ3 / Sin3 µ ) ( B + A Cos µ ) [ 2.B + ( 4.A - 2 ) Cos µ ]

-We may have a better precision when µ is very small if we use the following formulae to compute µ   A  &  B:

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

S  = Sin2 (µ/2)  = Sin2 G Cos2 H  + Cos2 F Sin2 H
C = Cos2 (µ/2) = Cos2 G Cos2 H + Sin2 F Sin2 H

A = ( Sin2 G Cos2 F ) / S + ( Cos2 G Sin2 F ) / C
B = ( Sin2 G Cos2 F ) / S - ( Cos2 G Sin2 F ) / C

-The following program takes the coordinates in decimal degrees:

Data Registers:   R00 = µ  R01 = A - B ( Sin µ ) / µ   R02 = B ( Sin µ ) / µ   R03 to R10:  temp
Flags: /
Subroutines: /

 01 LBL "DGC"  02 DEG  03 X<> T  04 -  05 X<> Z  06 -  07 2  08 ST/ Z  09 /  10 ST+ Z  11 1  12 STO 04  13 P-R  14 STO 03  15 RDN  16 STO 06            17 STO 10  18 X<> L  19 P-R  20 ST* 06  21 R^  22 *  23 X^2  24 X<>Y  25 STO 07 26 R^  27 RCL 04  28 P-R  29 ST* 07  30 ST* 10  31 RDN  32 ST* 03  33 *  34 X^2  35 +  36 X<> 10  37 X^2  38 RCL 06  39 X^2  40 RCL 07  41 X^2  42 +  43 STO 01            44 /  45 STO 02  46 RCL 03  47 X^2  48 RCL 10  49 /  50 ST- 02 51 +  52 X<> 01  53 SQRT  54 RCL 10  55 SQRT  56 R-P  57 X<>Y  58 ST+ X  59 D-R  60 STO 00  61 LASTX  62 RCL 04  63 P-R  64 STO 03  65 RDN  66 ST/ Y  67 X^2  68 ST+ X  69 STO 09            70 X<>Y  71 STO 06  72 RCL 02  73 RCL 01  74 ST- 04  75 RCL 03 76 *  77 +  78 STO 05  79 LASTX  80 RCL 03  81 -  82 ST+ X  83 +  84 STO 07  85 RCL 05  86 +  87 *  88 596.51444  89 STO 10  90 /  91 RCL 07  92 +  93 *  94 RCL 01            95 13  96 *  97 STO 07  98 -  99 6 100 + 101 * 102 RCL 02 103 5 104 * 105 + 106 RCL 05 107 * 108 RCL 07 109 - 110 RCL 01 111 X^2 112 STO 07 113 21 114 * 115 + 116 RCL 07         117 RCL 02 118 X^2 119 STO 08 120 ST* 09 121 ST+ 09 122 ST+ X 123 - 124 RCL 03 125 RCL 06 126 ST/ 02 127 ST* 05 128 / 129 ST* 08 130 * 131 ST+ 07 132 RCL 08 133 - 134 RCL 04 135 4 136 ST/ 07 137 * 138 RCL 01 139 * 140 ST- 07   141 RCL 00         142 X^2 143 * 144 3 145 ST* Z 146 / 147 - 148 + 149 RCL 10 150 ST/ 09 151 / 152 RCL 05 153 RCL 02 154 ST- 01 155 ST* 09 156 ST+ X 157 + 158 + 159 RCL 04 160 * 161 RCL 07 162 + 163 RCL 09 164 + 165 RCL 10 166 / 167 RCL 01         168 - 169 RCL 10 170 / 171 RCL 00 172 ST* Y 173 + 174 6378.137 175 * 176 END

( 233 bytes / SIZE 011 )

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

Where  D is the geodesic distance on the WGS84 ellipsoid

Example:     Long = 151°2049444 E , Lat = 33°86141667 S
Long = 116.8640000 W ,  Lat = 33°35622222 N

151.2049444   ENTER^
-33.86141667   ENTER^
-116.8640000    ENTER^
33.35622222   XEQ "DGC"   >>>>   D =  12157.20890 km                                          ---Execution time = 10.1s---

-Here are a few root-mean-square & maximal errors in centimeters.

 | λ'-λ | 20° 45° 90° 120° 160° 164° 168° 171° 175° 176° rms 1.3 1.1 1.2 1.3 1.9 2.1 2.6 4.1 16.7 30.3 max 3.8 3.8 4.0 4.2 4.7 4.7 10.1 23.6 133 280

-We have a better precision if  B3 [ ( Sin µ ) / µ + 2 ( Sin3 µ ) / µ ]  is replaced with   B3 [ -7 ( Sin µ ) / µ + 10 ( Sin3 µ ) / µ ] + A2 B [ 7 ( Sin µ ) / µ - 7 ( Sin3 µ ) / µ ]

Formula:

D/(a.µ) ~   1 + (f/2) ( -A + B (Sin µ)/µ )
+ (f2/4) [ A ( -4 + µ Cot µ ) + B ( 2 (Sin µ)/µ + µ / Sin µ ) + (A2/4) ( 17 + (Sin µ) (Cos µ)/µ - 4 µ Cot µ ) - A.B ( 2 (Sin µ)/µ + µ / Sin µ ) - ( B2/2 ) (Sin µ) (Cos µ)/µ ]
+ (f3/8) { ( 1 - A ) [ ( -A ( 1 - A ) (4/3) µ2 +  A ( 21.A - 13 )  + ( B + A Cos µ ) [ ( µ / Sin µ ) ( 6 - 13.A) + 5.B + ( µ2 / Sin2 µ ) ( B + ( 3 A - 2 ) Cos µ ) ) ]
+  ( 3.A2 - 9.B2 ) ( Sin µ ) ( Cos µ ) / µ ] + B3 [ -7 ( Sin µ ) / µ + 10 ( Sin3 µ ) / µ ] + A2 B [ 7 ( Sin µ ) / µ - 7 ( Sin3 µ ) / µ ] }
+ (f4/16) ( 1 - A ) ( µ3 / Sin3 µ ) ( B + A Cos µ ) [ 2.B + ( 4.A - 2 ) Cos µ ]

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

S  = Sin2 (µ/2)  = Sin2 G Cos2 H  + Cos2 F Sin2 H
C = Cos2 (µ/2) = Cos2 G Cos2 H + Sin2 F Sin2 H

A = ( Sin2 G Cos2 F ) / S + ( Cos2 G Sin2 F ) / C
B = ( Sin2 G Cos2 F ) / S - ( Cos2 G Sin2 F ) / C

-This program takes the coordinates in decimal degrees:

Data Registers:   R00 = µ  R01 = A - B ( Sin µ ) / µ   R02 = B ( Sin µ ) / µ   R03 to R10:  temp
Flags: /
Subroutines: /

 01 LBL "DGC"  02 DEG  03 X<> T  04 -  05 X<> Z  06 -  07 2  08 ST/ Z  09 /  10 ST+ Z  11 1  12 STO 04  13 P-R  14 STO 03  15 RDN  16 STO 05  17 STO 06            18 X<> L  19 P-R  20 ST* 06  21 R^  22 *  23 X^2  24 X<>Y  25 STO 07  26 R^ 27 RCL 04  28 P-R  29 ST* 05  30 ST* 07  31 RDN  32 ST* 03  33 *  34 X^2  35 +  36 RCL 05  37 X^2  38 RCL 06  39 X^2  40 RCL 07  41 X^2  42 +  43 STO 01  44 /  45 STO 02            46 RCL 03  47 X^2  48 R^  49 /  50 ST- 02  51 +  52 X<> 01 53 SQRT  54 X<>Y  55 SQRT  56 R-P  57 X<>Y  58 ST+ X  59 D-R  60 STO 00  61 STO 09  62 LASTX  63 RCL 04  64 P-R  65 STO 03  66 RCL 01  67 *  68 STO 07            69 X^2  70 X<>Y  71 ST/ 09  72 X^2  73 10  74 *  75 7  76 ST* Z  77 -  78 RCL 02 79 X^2  80 STO 08  81 *  82 +  83 X<> 09  84 RCL 02  85 RCL 07  86 +  87 STO 05  88 LASTX  89 RCL 03  90 -  91 ST+ X  92 +  93 STO 07  94 RCL 05            95 +  96 *  97 596.51444  98 STO 10  99 / 100 RCL 07 101 + 102 * 103 RCL 01 104 ST- 04 105 13 106 * 107 STO 06 108 - 109 6 110 + 111 * 112 RCL 02 113 5 114 * 115 + 116 RCL 05 117 * 118 RCL 06 119 - 120 RCL 01         121 X^2 122 STO 07 123 21 124 * 125 + 126 RCL 07 127 RCL 08 128 ST+ X 129 - 130 RCL 03 131 R^ 132 ST/ 02 133 ST* 05 134 / 135 ST* 08 136 * 137 ST+ 07 138 RCL 08 139 - 140 RCL 04 141 4 142 ST/ 07 143 * 144 RCL 01 145 * 146 ST- 07 147 RCL 00         148 X^2 149 * 150 3 151 ST* Z 152 / 153 - 154 + 155 RCL 10 156 ST/ 09 157 / 158 RCL 05 159 RCL 02 160 ST- 01 161 ST* 09 162 ST+ X 163 + 164 + 165 RCL 04 166 * 167 RCL 07 168 + 169 RCL 09 170 + 171 RCL 10 172 / 173 RCL 01 174 - 175 RCL 10         176 / 177 RCL 00 178 ST* Y 179 + 180 6378.137 181 * 182 END

( 238 bytes / SIZE 011 )

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

Where  D is the geodesic distance on the WGS84 ellipsoid

Example:     Long = 151°2049444 E , Lat = 33°86141667 S
Long = 116.8640000 W ,  Lat = 33°35622222 N

151.2049444   ENTER^
-33.86141667   ENTER^
-116.8640000    ENTER^
33.35622222   XEQ "DGC"   >>>>   D =  12157.20890 km                                          ---Execution time = 10.5s---

-Here are a few root-mean-square & maximal errors in centimeters.

 | λ'-λ | 1° 20° 45° 90° 120° 160° 164° 168° 170° 172° rms 0.9 0.9 0.9 0.8 0.9 1.2 1.4 2.1 3.0 5.0 max 2.4 2.4 2.4 2.4 2.3 2.4 4.1 10.0 16.8 34.0

d5)  Program(s)#5 - My Favorite Programs

Formula:

D/(a.µ) ~   1 + (f/2) ( -A + B (Sin µ)/µ )
+ (f2/4) [ A ( -4 + µ Cot µ ) + B ( 2 (Sin µ)/µ + µ / Sin µ ) + (A2/4) ( 17 + (Sin µ) (Cos µ)/µ - 4 µ Cot µ ) - A.B ( 2 (Sin µ)/µ + µ / Sin µ ) - ( B2/2 ) (Sin µ) (Cos µ)/µ ]
+ (f3/8) { ( 1 - A ) [ ( -A ( 1 - A ) (4/3) µ2 +  A ( 21.A - 13 )  + ( B + A Cos µ ) [ ( µ / Sin µ ) ( 6 - 13.A) + 5.B + ( µ2 / Sin2 µ ) ( B + ( 3 A - 2 ) Cos µ ) ) ]
+  ( 3.A2 - 9.B2 ) ( Sin µ ) ( Cos µ ) / µ ] + B3 [ -7 ( Sin µ ) / µ + Pi2 ( Sin3 µ ) / µ ] + A2 B [ 7 ( Sin µ ) / µ - 7 ( Sin3 µ ) / µ ] }
+ (f4/16) ( 1 - A ) ( µ3 / Sin3 µ ) ( B + A Cos µ ) [ 2.B + ( 4.A - 2 ) Cos µ ]
+ (f5/32) ( 1 - A ) ( µ4 / Sin4 µ ) ( B + A Cos µ ) [ 8.B + ( 16.A - 3 ) Cos µ ]

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

S  = Sin2 (µ/2)  = Sin2 G Cos2 H  + Cos2 F Sin2 H
C = Cos2 (µ/2) = Cos2 G Cos2 H + Sin2 F Sin2 H

A = ( Sin2 G Cos2 F ) / S + ( Cos2 G Sin2 F ) / C
B = ( Sin2 G Cos2 F ) / S - ( Cos2 G Sin2 F ) / C

-The following program takes the coordinates in decimal degrees:

Data Registers:   R00 = µ  R01 = A - B ( Sin µ ) / µ   R02 = B ( Sin µ ) / µ   R03 to R09:  temp
Flags: /
Subroutines: /

 01 LBL "DGC"  02 DEG  03 ST- Z  04 X<> T  05 -  06 X<>Y  07 2  08 ST/ Z  09 /  10 ST+ Z  11 1  12 STO 04  13 P-R  14 STO 03  15 RDN  16 STO 05  17 STO 06  18 X<> L  19 P-R  20 ST* 06  21 R^  22 *  23 X^2  24 X<>Y  25 STO 07            26 R^  27 RCL 04 28 P-R  29 ST* 05  30 ST* 07  31 RDN  32 ST* 03  33 *  34 X^2  35 +  36 RCL 05  37 X^2  38 RCL 06  39 X^2  40 RCL 07  41 X^2  42 +  43 STO 01  44 /  45 STO 02  46 RCL 03            47 X^2  48 R^  49 /  50 ST- 02  51 +  52 X<> 01  53 SQRT  54 X<>Y 55 SQRT  56 R-P  57 X<>Y  58 ST+ X  59 D-R  60 STO 00  61 LASTX  62 RCL 04  63 P-R  64 STO 03  65 RCL 01  66 ST- 04  67 *  68 STO 06  69 X^2  70 X<>Y  71 ST/ Z  72 PI  73 *  74 X^2  75 7  76 ST* Z  77 -  78 RCL 02            79 ST+ 06  80 X^2  81 STO 08 82 *  83 +  84 X<> 06  85 STO 05  86 RCL 03  87 RCL 04  88 *  89 -  90 ST+ X  91 STO 07  92 RCL 03  93 +  94 4  95 *  96 RCL 03            97 +  98 *  99 596.51444 100 STO 09 101 ST/ Z 102 / 103 RCL 07 104 + 105 * 106 RCL 05 107 - 108 RCL 07 109 + 110 * 111 RCL 01 112 13 113 * 114 STO 07 115 - 116 6 117 + 118 * 119 RCL 02 120 5 121 * 122 + 123 RCL 05 124 * 125 RCL 07 126 - 127 RCL 01         128 X^2 129 STO 07 130 21 131 * 132 + 133 RCL 07 134 RCL 08 135 ST+ X 136 - 137 RCL 03 138 R^ 139 ST/ 02 140 ST* 05 141 / 142 ST* 08 143 * 144 ST+ 07 145 RCL 08 146 - 147 RCL 04 148 4 149 ST/ 07 150 * 151 RCL 01         152 * 153 ST- 07 154 RCL 00 155 X^2 156 * 157 3 158 ST* Z 159 / 160 - 161 + 162 RCL 09 163 ST/ 06 164 / 165 RCL 05 166 + 167 RCL 02 168 ST- 01 169 ST* 06 170 ST+ X 171 + 172 RCL 04 173 * 174 RCL 07 175 + 176 RCL 06         177 + 178 RCL 09 179 / 180 RCL 01 181 - 182 RCL 09 183 / 184 RCL 00 185 ST* Y 186 + 187 6378.137 188 * 189 END

( 246 bytes / SIZE 010 )

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

Where  D is the geodesic distance on the WGS84 ellipsoid

Example:     Long = 151°2049444 E , Lat = 33°86141667 S
Long = 116.8640000 W ,  Lat = 33°35622222 N

151.2049444   ENTER^
-33.86141667   ENTER^
-116.8640000    ENTER^
33.35622222   XEQ "DGC"   >>>>   D =  12157.20890 km                                          ---Execution time = 11s---

-Here are a few root-mean-square & maximal errors in centimeters.

 | λ'-λ | 1° 20° 45° 90° 120° 160° 164° 168° 170° 172° rms 0.9 0.9 0.9 0.8 0.9 1.2 1.4 2.0 2.9 4.8 max 2.4 2.4 2.4 2.4 2.3 2.4 3.9 9.1 15.1 28.8

-We can save a few bytes if we compute  µ , A & B with the formulae:

Sin2 µ = Sin2 ( λ - λ' ) Cos2 β + [ Cos β' Sin β - Cos β Sin β' Cos ( λ - λ' ) ]2
Cos µ = Sin β Sin β' + Cos β Cos β' Cos ( λ - λ' )

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

Data Registers:   R00 = µ  R01 = A - B ( Sin µ ) / µ   R02 = B ( Sin µ ) / µ   R03 to R09:  temp
Flags: /
Subroutines: /

 01 LBL "DGC"  02 DEG  03 X<> T  04 -  05 X<>Y  06 1  07 STO 04  08 P-R  09 X<>Y  10 STO 00  11 STO 01  12 ST+ 01  13 STO 03  14 X^2  15 STO 02  16 RDN  17 P-R  18 STO 05  19 RCL Z  20 RCL 04            21 P-R  22 ST* 03  23 ST* 05  24 RDN  25 ST* 00  26 ST* 01 27 ST* Y  28 X^2  29 ST+ 02  30 X<> 03  31 -  32 X^2  33 X<>Y  34 X^2  35 +  36 ST/ 01  37 ST/ 02  38 SQRT  39 STO 07  40 RCL 00  41 RCL 05  42 +  43 STO 03            44 R-P  45 X<>Y  46 D-R  47 STO 00  48 RCL 07  49 /  50 RCL 02  51 RCL 03  52 ST* 02 53 RCL 01  54 ST- 02  55 *  56 -  57 STO 01  58 ST- 04  59 RCL 03  60 *  61 STO 06  62 X^2  63 RCL 07  64 PI  65 *  66 X^2  67 7  68 ST* Z  69 -  70 RCL 02            71 ST+ 06  72 X^2  73 STO 08  74 *  75 +  76 X<> 06  77 STO 05  78 RCL 03 79 RCL 04  80 *  81 -  82 ST+ X  83 STO 07  84 RCL 03  85 +  86 4  87 *  88 RCL 03  89 +  90 *  91 596.51444  92 STO 09            93 ST/ Z  94 /  95 RCL 07  96 +  97 *  98 RCL 07  99 + 100 RCL 05 101 - 102 * 103 RCL 01 104 13 105 * 106 STO 07 107 - 108 6 109 + 110 * 111 RCL 02 112 5 113 * 114 + 115 RCL 05 116 * 117 RCL 07 118 - 119 RCL 01 120 X^2 121 STO 07         122 21 123 * 124 + 125 RCL 07 126 RCL 08 127 ST+ X 128 - 129 RCL 03 130 R^ 131 ST/ 02 132 ST* 05 133 / 134 ST* 08 135 * 136 ST+ 07 137 RCL 08 138 - 139 RCL 04 140 4 141 ST/ 07 142 * 143 RCL 01 144 * 145 ST- 07 146 RCL 00         147 X^2 148 * 149 3 150 ST* Z 151 / 152 - 153 + 154 RCL 09 155 ST/ 06 156 / 157 RCL 05 158 + 159 RCL 02 160 ST- 01 161 ST* 06 162 ST+ X 163 + 164 RCL 04 165 * 166 RCL 07 167 + 168 RCL 06 169 + 170 RCL 09         171 / 172 RCL 01 173 - 174 RCL 09 175 / 176 RCL 00 177 ST* Y 178 + 179 6378.137 180 * 181 END

( 239 bytes / SIZE 010 )

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

Where  D is the geodesic distance on the WGS84 ellipsoid

Example:     Long = 151°2049444 E , Lat = 33°86141667 S
Long = 116.8640000 W ,  Lat = 33°35622222 N

151.2049444   ENTER^
-33.86141667   ENTER^
-116.8640000    ENTER^
33.35622222   XEQ "DGC"   >>>>   D =  12157.20890 km                                          ---Execution time = 10s---

239 bytes instead of 246 & 10 seconds instead of 11 seconds

If µ is very small, the result may be less accurate... and if µ is very very small, the result may be completely wrong:

Examples:

( λ , β ) = ( 0 , 26° )   ( λ' , β' ) = ( 0°0001 , 26°0001 )  ->   0.014967864 km    ( exact distance = 0.014957872... )

( λ , β ) = ( 0 , 26° )   ( λ' , β' ) = ( 0°000001 , 26°000001 )  ->   100992.4327 km  ( !!! )  ( exact distance = 0.00014957874... )

Conclusion:

-If µ is not very small, this program is preferable
-Otherwise, I choose the 1st program in this paragraph ( 246 bytes )

d6)  Program#6  ( latitudes geocentric or geodetic )

-Same formula.

Data Registers:   R00 = µ  R01 = A - B ( Sin µ ) / µ   R02 = B ( Sin µ ) / µ   R03 to R09:  temp

Flag:  CF 01  with geocentric latitudes
SF 01  with geodetic latitudes

Subroutines: /

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

( 270 bytes / SIZE 010 )

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

Where  D is the geodesic distance on the WGS84 ellipsoid

Example:     Long = 151°12'17"8 E , Lat = 33°51'41"1 S
Long = -116°51'50"4 W ,  Lat = 33°21'22"4 N

a)  Geocentric latitudes:    CF 01

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

b)  Geodetic latitudes:    SF 01

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

Note:

-If the latitudes are geodetic, they are changed into geocentric latitudes.

2°)  Triaxial Ellipsoid

a)  Andoyer2 + Andoyer2

-Like in "Terrestrial Geodesic Distance(2) for the HP41" , the distances of the 2 points M & N  on the great ellipses passing by M & N  are:

Δ  on the triaxial model
Δ' on the WGS84 ellipsoid

-Then, d = the geodesic distance 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 well approximated by the simple relation:

D ~  Δ + ( d - Δ' )

-The following programs employ similar methods used in "Terrestrial Geodesic Distance(2) for the HP41" to calculate great elliptic arc lengths on a triaxial ellipsoid.

-But to compute the difference between geodesic distance & great elliptic arc distance on an ellipsoid of revolution, we use formulae with geocentric latitudes.
-It is less complicated and thus, we save 1 data register and about 29 bytes & 4 seconds in execution time !

-Latitudes & longitudes are geodetic.

-The semi-axes of the triaxial ellipsoid 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.

-With a formula of order 2, the difference between geodesic distance & great elliptic arc length on an ellipsoid of revolution is

(d - Δ') / a'  =  µ ( 1 - A ) { (f2/4) [ A ( µ Cot µ - 1 ) + B ( µ / Sin µ - (Sin µ) / µ ) ]      with  a' = (a+b)/2

-It's the same formula if A , B & µ  are computed with geodetic, parametric or geocentric latitudes.
-But since µ , Sin µ , Cos µ are calculated in the 1st part of these programs, we can compute A , B more easily:

A = [ Sin2 g + Sin2 g' - 2 (Sin g)(Sin g') Cos µ ] / Sin2 µ                  where   g & g'  are geocentric latitudes and  µ = angle ( OM , ON )
B = [ ( Sin2 g + Sin2 g' ) Cos µ - 2 (Sin g)(Sin g') ] / Sin2 µ

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

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

( 323 bytes / SIZE 015 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

Example1:    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 = 19s---

LASTX = R10 =  12138.65874 km = great elliptic arc distance

Notes:

-When line 94 is executed,  the semi axes of the great ellipse =  R11/X & R11/Y.
-With the example above,  6378.151658 km & 6368.318461 km
-When the program stops, R02 = µ = 109°0855113

-We can perhaps decrease roundoff-errors if we add    ENTER^   SIGN   ST* Z   *   after line 81  ( just before  R-P )

Example2:       L = -95°  b = -25°      L' = 65°  b' = 39°

-95  ENTER^
-25  ENTER^
65  ENTER^
39   XEQ"TGD"   >>>>     D = 17582.16899  km
LASTX   d = 17582.26541  km

-The exact D-value is   17582.1683844
-So, error ~ 61 cm.

-This is about the maximum error if  | λ' - λ | <= 160°  ( except for nearly antipodal points )

b)  Andoyer2 + Andoyer3

-If  A , B & µ  are computed with geocentric latitudes, the difference between geodesic distance & great elliptic arc length is ( formula of order 3 ):

(d - Δ') / a'  =  µ ( 1 - A ) { (f2/4) [ A ( µ Cot µ - 1 ) + B ( µ / Sin µ - (Sin µ) / µ ) ]                              with  a' = (a+b)/2
+ (f3/8) [ B2 ( µ2 / Sin2 µ - 6 ( Sin µ  Cos µ ) / µ + 5 )
+ B (  ( -4 + 4 A ) ( Sin µ ) / µ + 5 A Cos µ + ( 4 A - 2 ) µ2 Cos µ / Sin2 µ + ( -13 A + 6 ) µ / Sin µ )
+ A ( 3A ( Sin µ  Cos µ ) / µ + ( -13 A + 6 ) µ Cot µ - ( 1 - A ) ( 8 µ2 / 6 ) + ( 3 A - 2 )  µ2 Cot2 µ - 4 + 7 A ) ] }

-We can save 19 bytes if this formula is re-written:

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

A = [ Sin2 g + Sin2 g' - 2 (Sin g)(Sin g') Cos µ ] / Sin2 µ             where   g & g'  are geocentric latitudes and  µ = angle ( OM , ON )
B = [ ( Sin2 g + Sin2 g' ) Cos µ - 2 (Sin g)(Sin g') ] / Sin2 µ

-In the program below, instead of  a' = (a+b)/2  ,  it uses   a' = (a+b)/2 + (a-b)/2  Cos (L+L')
-The precision is better.

-And the semi-axes of the great ellipse is also computed more quickly ( line 121 to 140 )

-We find an angle D such that    ( T + S + 2U ) Sin2 (W/2) + E ( Sin µ ) ( Sin ( µ - D ) )    has an extremum
-And the semi-axes of the great ellipse are

a1 = ( Sin W ) / ( ( T + S + 2U ) Sin2 (W/2) - E Cos2 (D/2) ) -1/2
a2 = ( Sin W ) / ( ( T + S + 2U ) Sin2 (W/2) + E Sin2 (D/2) ) -1/2

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

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

( 397 bytes / SIZE 015 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

Example1:    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 = 20s---

LASTX = R10 =  12138.65874 km = great elliptic arc distance

Notes:

-When line 140 is executed,  the semi axes of the great ellipse =  R11/X & R11/Y.
-With the example above,  6378.151658 km & 6368.318461 km
-When the program stops, R02 = µ = 109°0855113

Example2:       λ = -95°  b = -25°      λ' = 65°  b' = 39°

-95  ENTER^
-25  ENTER^
65  ENTER^
39   XEQ"TGD"   >>>>     D = 17582.16859  km
LASTX   d = 17582.26544  km

-The exact D-value is   17582.1683844...
-Error ~ 21 cm.

Notes:

-Here are a few comparisons between the errors with & without lines  62 to 66
-These results are obtained with free42

 ( λ , β )  ( λ' , β' ) error without errors with (-95°,-40°) (65°,50°) 52 cm 32 cm (-5°, -40° ) (155°,50) -54 cm -34 cm (-95°,-20°) (65°,10°) 43 cm 5 cm (-5°, -20°)  (155°,10°) -41 cm -4 cm (-96°,49°) (66°,-41°) 68 cm 47 cm ( 2° , 49°) (164°,-41°) -64 cm -43 cm

Note:

-If you want to use Andoyer's formula of order 3 for the great elliptic arc distance,
replace lines 145 to 172 by lines 145 to 200 of the version listed in paragraph 2°)e).

c)  Geocentric Longitudes & Geodetic Latitudes

-In fact, the longitudes & latitudes we read here and there are computed with WGS84 ellipsoid of revolution.
-So, the longitudes are geocentric and the latitudes are geodetic.

-The following program computes the geodesic distance on WGS84 if SF 00 & triaxial ellipsoid if CF 00

Data Registers:   R00 thru R15: temp

Flag:                     CF 00 = triaxial ellipsoid
SF 00 = biaxial ellipsoid
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 ENTER  19 STO 01  20 .07  21 -  22 STO 02  23 +  24 2  25 /  26 STO 15  27 FS? 00  28 STO 01  29 FS? 00  30 STO 02  31 21.384686  32 -  33 STO 03  34 RCL 15  35 /  36 X^2  37 STO 04  38 RCL 11            39 RCL 12  40 XEQ 01  41 STO 07  42 STO 12  43 RCL 09  44 STO 06  45 RCL 08 46 STO 05  47 RCL 13  48 ST+ 11  49 RCL 14  50 XEQ 01  51 STO 14  52 RCL 07  53 -  54 X^2  55 RCL 09  56 RCL 06  57 -  58 X^2  59 RCL 08  60 RCL 05  61 -  62 X^2  63 +  64 +  65 SQRT  66 STO 00  67 RCL 01  68 ST/ 05  69 ST/ 08  70 RCL 02  71 ST/ 06  72 ST/ 09  73 -  74 RCL 11  75 COS  76 *  77 2  78 ST/ 00  79 /  80 ST+ 15  81 RCL 03            82 ST/ 07  83 ST/ 10  84 RCL 08  85 X^2  86 RCL 09  87 X^2  88 RCL 10  89 X^2  90 + 91 +  92 STO 13  93 RCL 05  94 ST* 08  95 X^2  96 RCL 06  97 ST* 09  98 X^2  99 RCL 07 100 ST* 10 101 X^2 102 + 103 + 104 ST+ 13 105 - 106 RCL 00 107 ASIN 108 ST+ X 109 STO 02 110 1 111 STO 07 112 P-R 113 STO 04 114 RDN 115 STO 11 116 * 117 RCL 08 118 RCL 09 119 RCL 10 120 + 121 + 122 ST+ X 123 STO 09 124 RCL 13 125 ST+ 09 126 RCL 04         127 * 128 - 129 R-P 130 X<>Y 131 STO 08 132 2 133 / 134 X<>Y 135 SQRT 136 P-R 137 X^2 138 CHS 139 X<>Y 140 X^2 141 RCL 00 142 X^2 143 RCL 09 144 * 145 ST+ Z 146 + 147 STO 05 148 SQRT 149 STO 06 150 X<>Y 151 X<=0? 152 X<> 05 153 SQRT 154 ST- Y 155 ST+ X 156 / 157 STO 05 158 RCL 02 159 D-R 160 STO 00 161 RCL 08 162 ST+ X 163 COS 164 RCL 04 165 RCL 11 166 STO 10 167 * 168 * 169 - 170 4 171 / 172 RCL 05         173 * 174 RCL 08 175 COS 176 RCL 10 177 * 178 - 179 + 180 RCL 05 181 * 182 + 183 RCL 06 184 / 185 ST* 10 186 GTO 02 187 LBL 01 188 1 189 P-R 190 RCL 04 191 / 192 R-P 193 SIGN 194 P-R 195 X<>Y 196 STO 10 197 RDN 198 P-R 199 STO 08 200 X<>Y 201 STO 09 202 R^ 203 RTN 204 LBL 02 205 RCL 12 206 X^2 207 RCL 14 208 ST* 12 209 X^2 210 + 211 STO 05 212 RCL 04 213 STO 09         214 STO 13 215 ST* 05 216 RCL 12 217 ST+ X 218 ST- 05 219 * 220 - 221 RCL 11 222 X^2 223 ST/ 05 224 / 225 ST- 07 226 ST* 13 227 STO 14 228 RCL 00 229 RCL 11 230 / 231 STO 06 232 ST/ 09 233 STO 12 234 RCL 05 235 RCL 13 236 + 237 ST* 12 238 RCL 04 239 RCL 07 240 ST+ X 241 * 242 - 243 * 244 RCL 14 245 13 246 * 247 - 248 6 249 + 250 RCL 12 251 * 252 5 253 ST* 13 254 RCL 09 255 6 256 * 257 - 258 RCL 05         259 * 260 RCL 07 261 4 262 * 263 RCL 06 264 / 265 - 266 RCL 13 267 + 268 RCL 05 269 * 270 + 271 RCL 00 272 ST+ X 273 X^2 274 RCL 07 275 ST* 00 276 * 277 3 278 ST* 09 279 / 280 RCL 09 281 7 282 + 283 RCL 14 284 * 285 - 286 4 287 + 288 RCL 14 289 ST- 12 290 * 291 - 292 RCL 01 293 RCL 03 294 - 295 RCL 01 296 ST* 00 297 ST+ X 298 / 299 STO Z 300 * 301 RCL 05         302 RCL 06 303 / 304 - 305 RCL 12 306 + 307 * 308 * 309 RCL 00 310 * 311 RCL 10 312 + 313 END

( 394 bytes / SIZE 016 )

 STACK INPUTS OUTPUTS T λ ( ° ' " ) / Z β ( ° ' " ) / Y λ' ( ° ' " ) / X β' ( ° ' " ) 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 )

CF 00  = Triaxial Ellipsoid

151.12178  ENTER^
-33.51411   ENTER^
-116.51504  ENTER^
33.21224   XEQ "TGD"   >>>>   D = 12138.70175 km                                                      ---Execution time = 21s---
LASTX = R10 =  12138.70296 km = great elliptic arc distance

SF 00  = Biaxial Ellipsoid

151.12178  ENTER^
-33.51411   ENTER^
-116.51504  ENTER^
33.21224   XEQ "TGD"   >>>>   D = 12138.68432 km
LASTX = R10 =  12138.68552 km = great elliptic arc distance

Note:

-If you want to use Andoyer's formula of order 3 for the great elliptic arc distance,
replace lines 158 to 185 by lines 145 to 200 of the version listed in paragraph 2°)e).

d)  Geocentric Longitudes & Latitudes

-We can simplify the program if all the coordinates are geocentric.

Data Registers:   R00 thru R14: temp

Flag:                     CF 00 = triaxial ellipsoid
SF 00 = biaxial ellipsoid
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 FC? 00  26 GTO 00  27 RCL 01  28 RCL 02  29 +  30 2  31 /  32 STO 01  33 STO 02  34 LBL 00  35 RCL 11  36 RCL 12  37 1  38 P-R  39 X<>Y  40 STO 07            41 STO 12  42 RDN  43 P-R 44 STO 05  45 X<>Y  46 STO 06  47 RCL 13  48 ST+ 11  49 RCL 14  50 1  51 P-R  52 X<>Y  53 STO 10  54 STO 14  55 RDN  56 P-R  57 STO 08  58 RCL 05  59 -  60 X^2  61 X<>Y  62 STO 09  63 RCL 06  64 -  65 X^2  66 +  67 RCL 14  68 RCL 07  69 -  70 X^2  71 +  72 SQRT  73 2  74 /  75 STO 00  76 RCL 01  77 ST/ 05  78 ST/ 08  79 RCL 02            80 ST+ 01  81 ST/ 06  82 ST/ 09  83 -  84 RCL 11  85 COS  86 * 87 ST+ 01  88 RCL 03  89 ST/ 07  90 ST/ 10  91 RCL 08  92 X^2  93 RCL 09  94 X^2  95 RCL 10  96 X^2  97 +  98 +  99 STO 13 100 RCL 05 101 ST* 08 102 X^2 103 RCL 06 104 ST* 09 105 X^2 106 RCL 07 107 ST* 10 108 X^2 109 + 110 + 111 ST+ 13 112 - 113 RCL 00 114 ASIN 115 ST+ X 116 STO 02 117 1 118 STO 07 119 P-R 120 STO 04 121 RDN 122 STO 11         123 * 124 RCL 08 125 RCL 09 126 RCL 10 127 + 128 + 129 ST+ X 130 STO 09 131 RCL 13 132 ST+ 09 133 RCL 04 134 * 135 - 136 R-P 137 X<>Y 138 STO 08 139 2 140 ST/ 01 141 / 142 X<>Y 143 SQRT 144 P-R 145 X^2 146 CHS 147 X<>Y 148 X^2 149 RCL 00 150 X^2 151 RCL 09 152 * 153 ST+ Z 154 + 155 STO 05 156 SQRT 157 STO 06 158 X<>Y 159 X<=0? 160 X<> 05 161 SQRT 162 ST- Y 163 ST+ X 164 / 165 STO 05         166 RCL 02 167 D-R 168 STO 00 169 RCL 08 170 ST+ X 171 COS 172 RCL 04 173 RCL 11 174 STO 10 175 * 176 * 177 - 178 4 179 / 180 RCL 05 181 * 182 RCL 08 183 COS 184 RCL 10 185 * 186 - 187 + 188 RCL 05 189 * 190 + 191 RCL 06 192 / 193 ST* 10 194 RCL 12 195 X^2 196 RCL 14 197 ST* 12 198 X^2 199 + 200 STO 05 201 RCL 04 202 STO 09 203 STO 13 204 ST* 05 205 RCL 12 206 ST+ X 207 ST- 05 208 * 209 - 210 RCL 11         211 X^2 212 ST/ 05 213 / 214 ST- 07 215 ST* 13 216 STO 14 217 RCL 00 218 RCL 11 219 / 220 STO 06 221 ST/ 09 222 STO 12 223 RCL 05 224 RCL 13 225 + 226 ST* 12 227 RCL 04 228 RCL 07 229 ST+ X 230 * 231 - 232 * 233 RCL 14 234 13 235 * 236 - 237 6 238 + 239 RCL 12 240 * 241 5 242 ST* 13 243 RCL 09 244 6 245 * 246 - 247 RCL 05 248 * 249 RCL 07 250 4 251 * 252 RCL 06         253 / 254 - 255 RCL 13 256 + 257 RCL 05 258 * 259 + 260 RCL 00 261 ST+ X 262 X^2 263 RCL 07 264 ST* 00 265 * 266 3 267 ST* 09 268 / 269 RCL 09 270 7 271 + 272 RCL 14 273 * 274 - 275 4 276 + 277 RCL 14 278 ST- 12 279 * 280 - 281 RCL 01 282 RCL 03 283 - 284 RCL 01 285 ST* 00 286 ST+ X 287 / 288 STO Z 289 * 290 RCL 05 291 RCL 06 292 / 293 - 294 RCL 12         295 + 296 * 297 * 298 RCL 00 299 * 300 RCL 10 301 + 302 END

( 379 bytes / SIZE 015 )

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

Where  D = the geodesic distance & d = great elliptic arc length.

Example:    Calculate the geodesic distance between the Sydney Observatory ( Long = 151°12'17"8 E , Lat = 33°41'00"8954 S )    ( geocentric longitudes & latitudes )
and the Palomar Observatory ( Long = 116°51'50"4 W ,  Lat = 33°10'46"9543 N )

CF 00  = Triaxial Ellipsoid

151.12178        ENTER^
-33.41008954   ENTER^
-116.51504        ENTER^
33.10469543   XEQ "TGD"   >>>>   D = 12138.70176 km                                                                  ---Execution time = 18s---
LASTX = R10 =  12138.70297 km = great elliptic arc distance

SF 00  = Biaxial Ellipsoid

151.12178        ENTER^
-33.41008954   ENTER^
-116.51504        ENTER^
33.10469543   XEQ "TGD"   >>>>   D = 12138.68431 km
LASTX = R10 =  12138.68551 km = great elliptic arc distance

Notes:

-If  β = geocentric latitude and g = geodetic latitude

Tan β = ( 6356.7523142 / 6378.1372 )  Tan g

-Of course, if the angles are expressed in decimal degrres, delete lines 15-12-06-03  ( HR )

-If you want to use Andoyer's formula of order 3 for the great elliptic arc distance,
replace lines 166 to 193 by lines 145 to 200 of the version listed below.

e)  Andoyer3+ Andoyer3

-Though it seems impossible to find a formula of order 3 without calculating parametric latitudes in general case, I've found such a formula if the longitudes are the same.
-So, to compute the great elliptic arc distance with the geocentric coordinates, here is the formula employed by the following program ( lines 145 to 200 ):

D/a = W + (f/2) ( -W + B.SinW )  + (f 2 /4) [ W/4 + ( 1-2B2 ) ( SinW ) ( CosW )/4 ] + (f 3/8) [ W/4 + ( 1 - 2B2 + 30.B.CosW ) ( SinW ) ( CosW )/4 - (5/2) B3 Sin (3W) ]

a = major axis of the great ellipse , f = eccentricity
2F = b1+b2 where b1 & b2 = geocentric "latitudes"   B = Cos(2F)
W = | b1 - b2 |

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

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

( 433 bytes / SIZE 015 )

 STACK INPUTS OUTPUTS T λ ( ° ' " ) / Z β ( ° ' " ) / Y λ' ( ° ' " ) / X β' ( ° ' " ) 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 = 20.8s---

LASTX = R10 =  12138.65875 km = great elliptic arc distance

Notes:

-The great elliptic arc distance is almost perfect ( without roudoff-errors ):
-The precision is of the order of 1 millimeter !

-Unfortunately, the precision in the geodesic distance is not so good:

If | λ - λ' | = 120° ,  the maximum error is about 12 cm with ( 0° , 46° ) ( 120° ,-65° )  and root-mean-square  ~ 4 cm
If | λ - λ' | = 160° ,  the maximum error is about 40 cm with ( 0° , 60° ) ( 160° ,-52° )  and root-mean-square  ~ 9 cm

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" -
[5]  https://en.wikipedia.org/w/index.php?title=Vincenty%27s_formulae&oldid=1081502579

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