Template

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

Overview

1°)  Geodetic Longitudes & Latiudes

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

2°)  Geocentric Longitudes & Geodetic Latitudes

a)  Program#1
b)  Program#2

3°)  Geocentric Longitudes & Parametric Latitudes

4°)  Geocentric Longitudes & Latitudes

a)  Program(s)#1
b)  Program(s)#2
c)  Program#3

d)  Program#4
e)  Program#5 ( with great elliptic arc length )
f)  My favorite Programs

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

5°)  Transformation of Coordinates

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

Latest Update:
paragraph 4°)a)

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

-The semi-axes are:

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

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

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

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

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

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

1°)  Geodedic Longitudes & Latitudes

a)  Program#1

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

Dist ~   a { µ - (f/2) ( A.µ + 3.B Sin µ )
+ (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 )

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

-The angle µ  is computed with:

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

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

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

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

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

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

( 248 bytes / SIZE 010 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

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

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

-Exact result = 12138.657552.....

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

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

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

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

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

-Exact D = 18968.180062....

Notes:

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

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

| λ - λ' | = 120°

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

| λ - λ' | = 160°

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

| λ - λ' | = 164°

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

| λ - λ' | = 168°

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

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

b)  Program#2

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

Dist ~   a { µ - (f/2) ( A.µ + 3.B Sin µ )
+ (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 ) µ3 + ( B + A Cos µ ) ( µ3 / Sin2 µ ) ( B + ( 3 A - 2 ) Cos µ ) ] ( 1 - A )

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

-The angle µ  is computed with:

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

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

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

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

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

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

( 243 bytes / SIZE 011 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

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

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

-Exact result = 12138.657552.....

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

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

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

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

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

-Exact D = 18968.180062....

Notes:

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

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

| λ - λ' | = 120°

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

| λ - λ' | = 160°

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

| λ - λ' | = 168°

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

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

c)  Program#3

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

-Formula:

Dist ~   a' µ { 1 + (f/2) ( -A - 3.B (Sin µ)/µ )
+ (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) µ + ( 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 µ)/µ  ) ]

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

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

A" = [ cos2 β ( Sin2 λ - Cos2 λ )  + cos2 β' ( Sin2 λ' - Cos2 λ' ) + 2  cos β cos β' Cos ( λ + λ' )  Cos µ ] / Sin2 µ
B" = [ ( cos2 β ( Sin2 λ - Cos2 λ )  + cos2 β' ( Sin2 λ' - Cos2 λ' ) ) Cos  µ  + 2 cos β cos β' Cos ( λ + λ' ) ] / Sin2 µ

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

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

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

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

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

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

( 283 bytes / SIZE 013 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

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

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

-Exact  D = 12138.657552....

Notes:

-The precision becomes:

| λ - λ' | = 120°

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

| λ - λ' | = 160°

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

| λ - λ' | = 168°

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

| λ - λ' | = 170°

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

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

d)  Program#4

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

-Formula:

Dist ~   a' µ { 1 + (f/2) ( -A - 3.B (Sin µ)/µ )
+ (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) µ + ( 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 µ)/µ  ) ]

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

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 µ

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

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

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

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

( 288 bytes / SIZE 012 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

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

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

-Exact  D = 12138.657552....

Notes:

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

| λ - λ' | = 160°

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

| λ - λ' | = 162°

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

2°)  Geocentric Longitudes & Geodetic Latitudes

a)  Program#1

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

-Formula:

Dist ~   a' µ { 1 + (f/2) ( -A - 3.B (Sin µ)/µ )
+ (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) µ + ( 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 µ)/µ  ) ]

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

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

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

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

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

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

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

( 286 bytes / SIZE 012 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

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

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

-Exact  D = 12138.701763....

Notes:

-The precision is:

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

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

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

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

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

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

b)  Program#2

-Formula:

Dist ~   a' µ { 1 + (f/2) ( -A - 3.B (Sin µ)/µ )
+ (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) µ + ( 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 µ)/µ  ) ]

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

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

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

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

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

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

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

( 286 bytes / SIZE 012 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

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

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

-Exact  D = 12138.701763....

Notes:

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

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

| λ - λ' | = 120°

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

| λ - λ' | = 160°

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

| λ - λ' | = 164°

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

| λ - λ' | = 168°

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

3°)  Geocentric Longitudes & Parametric Latitudes

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

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

Formula:

Dist / ( a' µ )  ~ 1 - (f/2) ( A + B ( Sin µ ) / µ )
+ (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 µ) / µ ] }

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

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 µ

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

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

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

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

( 329 bytes / SIZE 013 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

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

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

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

-After replacing lines 23 & 31 with  /

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

Notes:

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

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

| λ - λ' | = 0°

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

| λ - λ' | = 120°

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

| λ - λ' | = 160°

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

| λ - λ' | = 168°

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

4°)  Geocentric Longitudes & Latitudes

a)  Program(s)#1

Formula:

Dist/(a µ) ~  1 + (f/2) ( -A + B (Sin µ)/µ )
+ (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 µ ) / µ ] }

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

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 µ

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

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

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

( 282 bytes / SIZE 013 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

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

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

-Exact D = 16362.072084...

Notes:

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

| λ - λ' | = 120°

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

| λ - λ' | = 140°

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

| λ - λ' | = 160°

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

| λ - λ' | = 168°

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

-Here is another program slightly more accurate:

Formula:

Dist/(a µ) ~  1 + (f/2) ( -A + B (Sin µ)/µ )
+ (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 µ ) / µ ] }

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

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

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 µ

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

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

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

( 284 bytes / SIZE 012 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

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

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

-Exact D = 16362.072084...

Notes:

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

| λ - λ' | = 120°

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

| λ - λ' | = 140°

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

| λ - λ' | = 160°

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

| λ - λ' | = 168°

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

-The maximum errors are smaller in the following programs.

b)  Program(s)#2

Formula:

Dist ~   a' µ { 1 + (f/2) ( -A + B (Sin µ)/µ )
+ (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) ) ]

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

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 µ

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

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

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

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

( 277 bytes / SIZE 010 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

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

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

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

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

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

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

Notes:

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

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

| λ - λ' | = 0°

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

| λ - λ' | = 120°

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

| λ - λ' | = 160°

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

| λ - λ' | = 168°

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

Remarks:

-If we replace lines 201 to 206 with:

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

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

-For instance:

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

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

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

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

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

( 267 bytes / SIZE 010 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

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

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

-Exact D = 12138.702497....

c)  Program#3

Formula:

Dist ~   a' µ { 1 + (f/2) ( -A + B (Sin µ)/µ )
+ (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 ( 24 A - 12 ) + ( B + A Cos µ ) [ 8 B + ( µ / Sin µ ) ( 6 - 12 A ) + ( µ2 / Sin2 µ ) ( B + ( 3 A - 2 ) Cos µ ) ) ] ]

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

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 µ

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

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

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

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

( 272 bytes / SIZE 010 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

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

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

-Exact D = 12138.702497....

Notes:

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

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

| λ - λ' | = 0°

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

| λ - λ' | = 120°

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

| λ - λ' | = 160°

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

| λ - λ' | = 168°

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

d)  Program#4

Formula:

Dist/(a' µ) ~  1 + (f/2) ( -A + B (Sin µ)/µ )
+ (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 µ ) / µ ] }

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

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 µ

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

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

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

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

( 290 bytes / SIZE 013 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

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

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

-Exact D = 12138.702497....

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

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

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

Remark:

-If we replace lines 211 to 216 with:

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

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

-For instance:

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

Notes:

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

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

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

| λ - λ' | = 0°

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

| λ - λ' | = 120°

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

| λ - λ' | = 140°

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

| λ - λ' | = 150°

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

| λ - λ' | = 160°

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

| λ - λ' | = 168°

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

| λ - λ' | = 170°

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

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

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

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

m = [ y2 (a2/b2 - 1)/2 + z2 (a2/c2 - 1)/2 ] / ( Sin2 µ )
With    n  = [ y'2 (a2/b2 - 1)/2 + z'2 (a2/c2 - 1)/2 ] / ( Sin2 µ )
p  = 2 [ y.y' (a2/b2 - 1)/2 + z.z' (a2/c2 - 1)/2 ] / ( Sin2 µ )

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

-After some calculations, it yields:

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

= ( 1 + m Sin2 µ ) ( 1 + e/2 ) ( 1 + h )

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

[ R2 + ( dR/dt )2 ] 1/2  ~   a  ( 1 + m Sin2 µ ) -1/2 (1 + e/2 )-1/2 [ 1 - h/2 + (3/8) h2 - (5/16) h3 + (1/8) h'2 - (5/16) h h'2 ]             where  h' = dh/dt

-And the great elliptic arc length is:

Δ / a   =   ( 1 + m Sin2 µ )-1/2 ( 1 + e/2 )-1/2 {  µ + ( G / 4 )    [ - Cos D + Cos ( 2.µ + D ) ]  + ( G2 / 32 ) [ 14 µ - ( Sin D ) ( Cos D ) + Sin ( 2.µ + D ) Cos ( 2.µ + D ) ]
+ ( G3 / 32 ) [ -5 Cos D + 5 Cos ( 2.µ + D ) - 5 Cos3 D + 5 Cos3 ( 2.µ + D ) ] }

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

(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 ] }

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

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

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

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

( 407 bytes / SIZE 017 )

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

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

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

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

-Exact D = 12138.702497....

Notes:

-With free42 , replace line 49 by  2195022 E-11  and line 53 by  675054581 E-11  ( exact values are:  a2/b2-1  &  a2/c2-1 )

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

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

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

| λ - λ' | = 0°

maximum error ~ 1 mm

| λ - λ' | = 120°

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

| λ - λ' | = 160°

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

| λ - λ' | = 168°

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

f)  My favorite Programs

f1)  Program#1

Formula:

Dist/(a' µ) ~  1 + (f/2) ( -A + B (Sin µ)/µ )
+ (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 µ ]

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

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 µ

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

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

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

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

( 319 bytes / SIZE 015 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

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

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

-Exact D = 12138.702497....

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

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

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

Remarks:

-If we replace lines 232 to 237 with:

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

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

-For instance:

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

-If you prefer:

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

replace lines 112-229-235 with   596.51284   364101   6378.13661

-The result in example 1 becomes:  12138.70176 km

Notes:

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

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

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

| λ - λ' | = 0°

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

| λ - λ' | = 90°

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

| λ - λ' | = 120°

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

| λ - λ' | = 160°

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

| λ - λ' | = 168°

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

| λ - λ' | = 170°

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

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

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

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

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

f2)  Program(s)#2

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

-For instance,

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

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

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

S  = 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

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

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

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

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

-Of course, the program becomes longer:

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

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

( 340 bytes / SIZE 016 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

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

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

-Exact D = 12138.702497....

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

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

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

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

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

-Exact distance = 0.014957847... km

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

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

-Exact distance = 0.000149578497... km

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

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

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

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

( 336 bytes / SIZE 015 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

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

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

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

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

Remark:

-If we replace lines 249 to 254 with:

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

Y-output = Geodesic distance on the ellipsoid of revolution

-For instance:

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

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

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

( 168 bytes / SIZE 029 )

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

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

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

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

-Exact distance = 20001.904992...

Notes:

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

 01 LBL "DDD"  02 STO 19  03 RDN  04 STO 18  05 RDN  06 STO 17  07 RDN  08 STO 16            09 RDN  10 X<> T  11 -  12 R^  13 COS  14 P-R  15 X<>Y  16 X<> Z  17 COS  18 ST+ Y  19 X<> L  20 SIN  21 R^  22 SIN  23 + 24 X<> Z  25 X<>Y  26 R-P  27 RCL Z  28 X<>Y  29 R-P  30 X<> Z  31 RCL 16            32 +  33 X<>Y  34 STO 21  35 X<>Y  36 STO 20  37 10  38 STO 28  39 RCL 21  40 STO 26  41 XEQ 10  42 STO 27  43 LBL 01  44 RCL 26  45 RCL 28  46 + 47 XEQ 10  48 STO 29  49 LBL 02  50 RCL 26  51 RCL 28  52 -  53 XEQ 10  54 STO 30            55 LBL 03  56 VIEW 27  57 RCL 29  58 X>Y?  59 X<>Y  60 RCL 27  61 X>Y?  62 X<>Y  63 RCL 27  64 X=Y?  65 GTO 05  66 CLX  67 RCL 29  68 X=Y?  69 GTO 04 70 RCL 30  71 X<> 27  72 STO 29  73 RCL 28  74 ST- 26  75 GTO 02  76 LBL 04  77 X<> 27  78 STO 30  79 RCL 28            80 ST+ 26  81 RCL 26  82 +  83 XEQ 10  84 STO 29  85 RCL 30  86 GTO 03  87 LBL 05  88 3  89 ST/ 28  90 RCL 28  91 ABS  92  E-4 93 X

( 203 bytes / SIZE 031 )

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

Where  D is the geodesic distance on the triaxial ellipsoid

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

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

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

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

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

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

5°)  Transformation of Coordinates

a)  Geocentric Coordinates <> Geodetic Coordinates

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

x = k a2 Cos B Cos L
y = k b2 Cos B Sin L
z = k c2 Sin B

where    k = ( a2 Cos2 B Cos2 L + b2 Cos2 B Sin2 L + c2 Sin2 B ) -1/2

Data Registers: /
Flags: /
Subroutines: /

 01 LBL "LCLD"  02 DEG  03 HR  04 1  05 P-R  06 X<>Y 07 6356752.314  08 X^2  09 /  10 X<> Z  11 HR  12 X<>Y 13 P-R  14 6378102  15 X^2  16 ST/ Z            17 CLX  18 6378172 19 X^2  20 /  21 R-P  22 X<>Y            23 RDN  24 R-P 25 X<>Y            26 HMS  27 R^  28 HMS  29 END

( 62 bytes / SIZE 000 )

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

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

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

 01 LBL "LDLC"  02 DEG  03 HR  04 1  05 P-R  06 X<>Y 07 6356752.314  08 X^2  09 *  10 X<> Z            11 HR  12 X<>Y 13 P-R  14 6378172  15 X^2  16 *  17 X<>Y            18 6378102 19 X^2  20 *  21 X<>Y            22 R-P  23 X<>Y  24 HMS 25 RDN  26 R-P  27 X<>Y            28 HMS  29 R^  30 END

( 62 bytes / SIZE 000 )

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

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

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

b)  WGS84 Coordinates -> Triaxial Ellipsoid Geocentric Coordinates

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

Data Registers: /
Flags: /
Subroutines: /

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

( 42 bytes / SIZE 000 )

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

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

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

Remark:

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

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

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

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

-IERS2003 suggests other numbers:  mean equatorial radius = 6378.13661 km & flattening = 1/298.256421 so, polar radius = 6356.751868 km

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

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