Terrestrial Geodesic Distance(3) for the HP-41
Overview
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
2°) Triaxial Ellipsoid
a) Andoyer2 + Andoyer2
b) Andoyer2 + Andoyer3
c) Geocentric Longitudes & Geodetic Latitudes
d) Geocentric Longitudes & Latitudes
e) Andoyer3 + Andoyer3
Latest Update: §1°) b4)
-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 )
-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 )
-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 )
+ (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 )
-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
-We can add other terms.
-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 )
-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 )
-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 )
-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
-Yet other added terms:
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 02 66 ST+ X 67 RCL 01 68 4 |
69 * 70 1 71 + 72 RCL 03 73 * 74 + 75 RCL 06 76 * 77 596.51444 78 STO 10 79 / 80 RCL 02 81 - 82 RCL 01 83 1 84 - 85 RCL 03 86 * 87 - 88 RCL 01 89 RCL 03 90 * 91 STO 12 92 RCL 02 93 + 94 STO 05 95 * 96 RCL 06 97 ST* Y 98 X^2 99 * 100 RCL 10 101 / 102 STO 11 |
103 RCL 07 104 X^2 105 24 106 * 107 18 108 - 109 RCL 02 110 STO 09 111 ST* 09 112 * 113 RCL 01 114 33 115 * 116 25 117 - 118 RCL 03 119 * 120 + 121 RCL 02 122 * 123 RCL 12 124 X^2 125 19 126 * 127 + 128 2 129 ST* 09 130 X^2 131 - 132 RCL 02 133 * 134 RCL 01 135 15 136 * |
137 11 138 - 139 RCL 01 140 X^2 141 ST- 09 142 * 143 RCL 03 144 * 145 - 146 RCL 06 147 / 148 STO 07 149 RCL 08 150 ST* 09 151 RCL 06 152 RCL 05 153 RCL 03 154 RCL 04 155 ST+ X 156 * 157 - 158 * 159 RCL 01 160 7 161 * 162 2 163 - 164 STO 08 165 + 166 * 167 RCL 02 168 3 169 * 170 STO 03 |
171 - 172 RCL 05 173 * 174 RCL 08 175 2.5 176 * 177 RCL 00 178 ST+ X 179 X^2 180 RCL 04 181 * 182 3 183 / 184 - 185 RCL 02 186 18 187 * 188 RCL 06 189 / 190 + 191 RCL 01 192 * 193 + 194 RCL 11 195 - 196 RCL 04 197 * 198 RCL 07 199 + 200 RCL 10 201 / 202 RCL 05 203 RCL 03 204 RCL 06 |
205 ST* Z 206 / 207 STO 03 208 ST+ X 209 + 210 RCL 04 211 * 212 + 213 RCL 09 214 RCL 01 215 X^2 216 - 217 15 218 * 219 RCL 01 220 ST+ 03 221 4 222 ST/ Z 223 * 224 + 225 + 226 RCL 10 227 / 228 RCL 03 229 - 230 RCL 10 231 / 232 RCL 00 233 ST* Y 234 + 235 6378.137 236 * 237 END |
( 300 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 )
-33.51411 ENTER^
-116.51504 ENTER^
33.21224 XEQ "DGT" >>>> D = 12138.68432 km ---Execution time = 12.5s---
-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
-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
-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
-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
-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
-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
-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
-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
-33.86141667 ENTER^
-116.8640000 ENTER^
33.35622222 XEQ "DGC" >>>> D = 12157.20890 km ---Execution time = 10s---
Advantages:
239 bytes instead of 246 & 10 seconds instead of 11 seconds
Disadvantage:
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 )
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---
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---
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---
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] !