Geodesics on a Triaxial Ellipsoid for the HP-41
Overview
1°) Jacobi's Method
a) Program#1
b) Program#2
c) Program(s)#3
2°) Longitude & Latitude >>> Cartesian Coordinates
a) Geodetic Coordinates
b) Astronomical Coordinates
c) Geocentric Coordinates
3°) Cartesian Coordinates <> Jacobi
Parameters ( u , v )
4°) Approximate Methods
a) Program#1
b) Program#2
c) Program#3
d) Program#4
e) For the Earth only ( 2 programs )
Latest Update: A second program in §4e)
-The routines listed in paragraph 1 employ the rigorous method found by
Jacobi in 1838.
-In paragraph 4, the results are only approximate but still very good
for the Earth:
-With the last 3 programs, the accuracy is of the order of 20 cm, except
for nearly antipodal points.
1°) Jacobi's Method
a) Program#1
-The geodesic distance between 2 points on a scalene - or triaxial - ellipsoid
is much more difficult to calculate
than the geodesic distance on an ellpsoid of revolution.
-The program hereunder uses the formulae given by Jacobi in reference
[1]
-Given an ellipsoid x2 / a + y2 / b + z2 / c = 1 with a > b > c ( here, a , b , c are the squares of the semi-axes ),
and 2 points M ( u1 , v1 ) & N ( u2 , v2 ) where x y z and u v are related by:
x = [ a / ( a - c ) ]1/2 Cos u
( a - b Sin2 v - c Cos2 v )1/2
y = b1/2 Sin u
Cos v
z = [ c / ( a - c ) ]1/2 Sin v
( a Sin2 u + b Cos2 u - c )1/2
-First, we have to determine the parameter Beta such that ( § = Integral )
§u1u2 ( a Sin2 u + b Cos2 u )1/2 ( a Sin2 u + b Cos2 u - c ) -1/2 ( (a-b) Sin2 u + Beta ) -1/2 du
= +/- §v1v2 ( b Sin2 v + c Cos2 v )1/2 ( a - b Sin2 v - c Cos2 v ) -1/2 ( (b-c) Cos2 v - Beta ) -1/2 dv
-This is done by the routine "BETA"
-Then, the geodesic distance s between the 2 points may be computed by 2 formulae:
s = §u1u2 [ P + R ( dv/du )2 ]1/2 du (1)
s = §v1v2 [ P ( du/dv )2 +R ]1/2 dv (2)
where the geodesic parameters are given by
P = ( ¶x / ¶u
)2 + ( ¶y / ¶u )2 + ( ¶z
/ ¶u )2
R = ( ¶x / ¶v )2 + ( ¶y
/ ¶v )2 + ( ¶z / ¶v )2
( the other geodesic parameter Q = ( ¶x / ¶u )( ¶x / ¶v ) + ( ¶y / ¶u )( ¶y / ¶v ) + ( ¶z / ¶u )( ¶z / ¶v ) = 0 here )
-It yields:
P = ( a / (a-c) ) Sin2 u ( a - b Sin2 v - c Cos2 v ) + b Cos2 u Cos2 v + ( c (a-b)2 / (a-c) ) Sin2 u Cos2 u Sin2 v ( a Sin2 u + b Cos2 u - c ) -1
R = ( c / (a-c) ) Cos2 v ( a Sin2 u + b Cos2 u - c ) + b Sin2 u Sin2 v + ( a (b-c)2 / (a-c) ) Sin2 v Cos2 v Cos2 u ( a - b Sin2 v - c Cos2 v ) -1
-"DIST" calculate s with Gauss-Legendre 10-point formula.
-The intervall of integration is divided in n subintervals ( n must
be stored in R00 )
-For each u-value, we have to calculate the corresponding v-value:
-So we must solve a non-linear equation of the form §v1v
f(v) dv - §u1u g(u) du = 0
-This is done by applying again Gauss-Lendre integration + Newton's method
!
-So you can guess that a good emulator in turbo mode is not superfluous...
Data Registers: • R00 = n ( Registers R00 thru R07 and R40 thru R49 are to be initialized before executing "BETA" & "DIST" )
• R01 = a
• R04 = u1
• R06 = u2
R09 to R39: temp
• R02 = b
• R05 = v1
• R07 = v2
• R03 = c
>>> R08 = Beta
• R40 = 0.2955242247
• R45 = 0.6794095683
• R41 = 0.1488743390
• R46 = 0.1494513492
• R42 = 0.2692667193
• R47 = 0.8650633667
which are the coefficients for Gauss-Legendre 10-point formula
• R43 = 0.4333953941
• R48 = 0.06667134431
• R44 = 0.2190863625
• R49 = 0.9739065285
Flags: F01 F02 F06 F07 F09 F10
CF 01 if ( u1 - u2
).( v1 - v2 ) > 0
SF 01 if ( u1 - u2
).( v1 - v2 ) < 0
CF 02 if you want to compute the geodesic distance with formula
(1)
SF 02 if you want to compute the geodesic distance with
formula (2)
Subroutines: /
| 01 LBL "BETA" 02 DEG 03 CF 10 04 CF 07 05 49.039 06 STO 10 07 RDN 08 STO 21 09 X<>Y 10 STO 09 11 XEQ 00 12 STO 20 13 LBL 01 14 RCL 21 15 VIEW 21 16 XEQ 00 17 ENTER 18 ENTER 19 X<> 20 20 - 21 X#0? 22 / 23 RCL 09 24 RCL 21 25 STO 09 26 STO T 27 - 28 * 29 + 30 STO 21 31 X#Y? 32 GTO 01 33 STO 08 34 RTN 35 LBL "DIST" 36 LBL 12 37 CF 10 38 RCL 04 39 RCL 06 40 + 41 RCL 05 42 RCL 07 43 + 44 2 45 ST/ Z 46 / 47 FS? 02 48 X<>Y 49 STO 12 50 49.039 51 STO 10 52 RCL 04 53 FS? 02 54 RCL 05 55 STO 11 56 RCL 06 57 FS? 02 58 RCL 07 59 E-8 60 STO 22 61 SF 07 62 XEQ 14 63 ABS 64 D-R 65 CLD 66 RTN 67 GTO 12 68 LBL 00 69 STO 08 70 RCL 06 71 RCL 04 72 STO 11 |
73 CF 06 74 XEQ 14 75 STO 19 76 RCL 07 77 RCL 05 78 STO 11 79 SF 06 80 XEQ 14 81 FS? 01 82 CHS 83 ST- 19 84 RCL 19 85 RTN 86 GTO 00 87 LBL 13 88 CLX 89 STO 28 90 CLX 91 RCL 21 92 STO 23 93 - 94 RCL 00 95 STO 24 96 ST+ X 97 / 98 STO 25 99 LBL 06 100 ST+ 23 101 RCL 10 102 STO 26 103 LBL 07 104 RCL 23 105 STO 27 106 RCL 25 107 RCL IND 26 108 * 109 ST+ 27 110 - 111 XEQ 02 112 X<> 27 113 XEQ 02 114 RCL 27 115 + 116 DSE 26 117 RCL IND 26 118 * 119 ST+ 28 120 DSE 26 121 GTO 07 122 RCL 25 123 ST+ 23 124 DSE 24 125 GTO 06 126 ST* 28 127 RCL 28 128 RTN 129 LBL 14 130 CLX 131 STO 18 132 CLX 133 RCL 11 134 STO 13 135 - 136 RCL 00 137 STO 14 138 ST+ X 139 / 140 STO 15 141 LBL 08 142 ST+ 13 143 RCL 10 144 STO 16 |
145 LBL 09 146 RCL 13 147 STO 17 148 RCL 15 149 RCL IND 16 150 * 151 ST+ 17 152 - 153 FS? 07 154 XEQ 05 155 FC? 07 156 XEQ 02 157 X<> 17 158 FS? 07 159 XEQ 05 160 FC? 07 161 XEQ 02 162 RCL 17 163 + 164 DSE 16 165 RCL IND 16 166 * 167 ST+ 18 168 DSE 16 169 GTO 09 170 RCL 15 171 ST+ 13 172 DSE 14 173 GTO 08 174 ST* 18 175 RCL 18 176 RTN 177 LBL 02 178 FS? 06 179 GTO 03 180 SIN 181 X^2 182 RCL 01 183 RCL 02 184 - 185 * 186 STO Y 187 RCL 02 188 + 189 RCL X 190 RCL 03 191 - 192 / 193 X<>Y 194 RCL 08 195 + 196 / 197 X<0? 198 SF 09 199 FS? 10 200 ABS 201 SQRT 202 RTN 203 LBL 03 204 COS 205 X^2 206 RCL 02 207 RCL 03 208 - 209 * 210 ENTER 211 CHS 212 RCL 02 213 + 214 RCL 01 215 RCL Y 216 - |
217 / 218 X<>Y 219 RCL 08 220 - 221 / 222 X<0? 223 SF 09 224 FS? 10 225 ABS 226 SQRT 227 RTN 228 LBL 05 229 STO 31 230 FC? 02 231 RCL 04 232 FS? 02 233 RCL 05 234 STO 21 235 CF 06 236 FS? 02 237 SF 06 238 XEQ 13 239 FS? 01 240 CHS 241 STO 20 242 RCL 05 243 FS? 02 244 RCL 04 245 STO 21 246 SF 10 247 RCL 12 248 STO 30 249 SF 06 250 FS? 02 251 CF 06 252 LBL 10 253 CF 09 254 RCL 30 255 FC? 02 256 XEQ 03 257 FS? 02 258 XEQ 02 259 STO 29 260 RCL 30 261 ENTER 262 XEQ 13 263 RCL 20 264 - 265 RCL 29 266 / 267 ST- 30 268 VIEW 30 269 RCL 30 270 X=0? 271 SIGN 272 / 273 ABS 274 RCL 22 275 X<Y? 276 GTO 10 277 FS? 09 278 SF 41 279 RCL 31 280 FS? 02 281 X<> 30 282 STO 31 283 STO 12 284 1 285 P-R 286 X^2 287 STO 32 288 RCL 02 |
289 * 290 X<>Y 291 X^2 292 STO 33 293 RCL 01 294 * 295 + 296 STO 36 297 RCL 03 298 - 299 RCL 30 300 FC? 02 301 STO 12 302 COS 303 X^2 304 STO 34 305 * 306 RCL 03 307 * 308 RCL 01 309 RCL 03 310 - 311 / 312 RCL 30 313 SIN 314 X^2 315 STO 35 316 RCL 33 317 * 318 RCL 02 319 * 320 + 321 RCL 34 322 RCL 03 323 * 324 RCL 35 325 RCL 02 326 * 327 + 328 STO 37 329 RCL 01 330 X<>Y 331 - 332 RCL 34 333 RCL 35 334 * 335 X<>Y 336 / 337 RCL 32 338 * 339 RCL 01 340 * 341 RCL 02 342 RCL 03 343 - 344 X^2 345 * 346 RCL 01 347 RCL 03 348 - 349 / 350 + 351 STO 38 352 RCL 32 353 RCL 33 354 * 355 RCL 35 356 * 357 RCL 01 358 RCL 02 359 - 360 X^2 |
361 * 362 RCL 03 363 * 364 RCL 01 365 RCL 03 366 - 367 / 368 RCL 36 369 RCL 03 370 - 371 / 372 RCL 32 373 RCL 34 374 * 375 RCL 02 376 * 377 + 378 RCL 01 379 RCL 37 380 - 381 RCL 33 382 * 383 RCL 01 384 ST* Y 385 RCL 03 386 - 387 / 388 + 389 STO 39 390 RCL 36 391 RCL 01 392 RCL 37 393 - 394 * 395 RCL 02 396 RCL 03 397 - 398 RCL 34 399 * 400 RCL 08 401 - 402 * 403 RCL 37 404 RCL 36 405 RCL 03 406 - 407 * 408 RCL 01 409 RCL 02 410 - 411 RCL 33 412 * 413 RCL 08 414 + 415 * 416 FS? 02 417 X<>Y 418 / 419 RCL 39 420 RCL 38 421 FS? 02 422 X<>Y 423 RCL Z 424 * 425 + 426 SQRT 427 RTN 428 END |
( 642 bytes / SIZE 050 )
| STACK | BETA-INPUT | BETA-OUTPUT | DIST-OUTPUT |
| Y | beta1 | / | / |
| X | beta2 | beta | s |
Where beta1 & beta2 are 2 initial guesses, beta the solution and s = geodesic distance.
Example1: Ellipsoid x2 / 41 + y2 / 37 + z2 / 35 = 1
u1 = +10°
u2 = +75°
( all in decimal degrees )
v1 = -15°
v2 = +61°
41 STO 01 10 STO
04 75 STO 06
37 STO 02 -15 STO 05
61 STO 07
35 STO 03
-If you choose n = 4, 4 STO 00
CF 01 since ( u1 - u2 ) & ( v1 - v2 ) have the same sign
and if your initial guesses for beta are 0.1 and 0.2
0.1 ENTER^
0.2 XEQ "BETA" >>>>
beta = 0.1986540209 = R08
-With n = 8 STO 00, the same guesses give beta = 0.1986540199
which is stored in R08
-Let's keep this value for beta in R08.
-Then, simply press R/S or XEQ "DIST" to get the geodesic distance
-With n = 2 STO 00 and CF 02, it yields
s = 8.594822585
-With n = 4 STO 00
R/S >>>>
s = 8.594822582
-You could want to check if using the 2nd formula confirms this value:
2 STO 00 SF 02 R/S >>>>
s = 8.594822577
4 STO 00 R/S >>>>
s = 8.594822577
-So the exact result is probably about 8.594822580
Example2: On the Earth: a triaxial ellipsoid model with semi-axes 6378.172 km , 6378.102 km , 6356.752 km ,
and the equatorial major axis in the direction 14°93W , 165°07E
• Store the squares of these semi-axes into R01-R02-R03
-Calculate the geodesic distance between the U.S. Naval Observatory at
Washington (D.C.) L1 = 77°03'56"0 W ;
b1 = 38°55'17"2 N
and Cape Town Observatory ( South Africa ) L2
= 18°28'35"7 E ; b2 = 33°56'02"5 S
-With the programs listed in paragraph 2°)b) and 3°), we find the following Jacobi parameters:
u1 = -62°16090881
u2 = 33°42630683
store these 4 numbers R04 R06
respectively
v1 = +38°84397819
v2 = -33°88879132
into
R05 R07
-Set flag F01 SF 01 since ( u1 - u2 ).( v1 - v2 ) < 0
-You will find, with n = 4 STO 00 , that beta = 136050.1048 and s = 12709.56546 km
Notes:
-For n = 1, the execution time t for "DIST" is about 3 seconds with V41
in turbo mode
-For n = 2, t = 8s
-For n = 4, t = 26s
-Multiply these values by 600 to get the execution time with a real HP-41...
-In fact, getting Beta is seldom easy !
-Generally, you will have to calculate f(beta) for many beta-values
before finding a small interval (b,b') for which f(b)
< 0 & f(b') > 0
-For that purpose, XEQ 00 with CF 07 CF 10 and R10 = 49.039 with your guess in X-register.
-Flag F10 is used to avoid a DATA ERROR line 198 or 223:
-When the Newton's method is applied - lines 252 to 276 - it may happen
that the iterations lead to a negative radicand.
-However, the last iteration must not be the result of a wrong calculation:
-In this case, line 278 would produce a NONEXISTENT message.
-The method also works if a = b > c or if a > b = c , but it's not the fastest way to solve the problem in this case...
Problem:
-This program cannot find the geodesic distance between Washington and Paris because the geodesic looks like this:
*
*
*
*
*
*
Paris
*
Washington
-In other words, the derivative dv / du changes its sign
-The program hereunder overcomes this limitation.
b) Program#2
-So, we have to compute §v1v0 f(v) dv - §v0v2 f(v) dv
with f(v) = ( b Sin2 v + c Cos2 v )1/2 ( a - b Sin2 v - c Cos2 v ) -1/2 ( (b-c) Cos2 v - Beta ) -1/2
where the denominator ( (b-c) Cos2 v - Beta ) vanishes between v1 & v2 i-e ( (b-c) Cos2 v0 - Beta ) = 0
-In order to remove this singularity without producing another singularity for v = 0, I've made the following change of variable:
Cos w = [ ( Cos2 v - K ) / ( 1 - K
) ] 1/2
with K = beta / ( b - c )
Sin w = Sin v / ( 1 - K ) 1/2
-The integral becomes: ( b-c ) -1/2 §w1w2 g(w) dw
with g(w) = [ b - ( b-c ).{ K + ( 1-K) Cos2
w } ]1/2 [ ( a-b ) - ( b-c ).{ K + ( 1-K) Cos2
w } ] -1/2 [ K + ( 1-K) Cos2 w } ] -1/2
-Flags F01-F02-F04 determine the form of the geodesic between the 2 points:
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
* *
*
CF 01 CF 04
SF 01 CF 04
CF 01-CF 02-SF04
SF 01-CF 02-SF 04
CF 02 or SF 02 <-- It's
your choice --> CF 02 or SF 02
-Thus, F02 must be cleared if F04 is set ( add FS? 04 CF 02 inside the listing below if you prefer )
-Of course, it's not always obvious to know the good case in advance.
-So, you'll probably have to make several tests before finding the right
configuration.
Data Registers: • R00 = n ( Registers R00 thru R07 and R42 thru R51 are to be initialized before executing "BETA" & "DIST" )
• R01 = a
• R04 = u1
• R06 = u2
R09 to R38: temp
• R02 = b
• R05 = v1
• R07 = v2
• R03 = c
>>> R08 = Beta
R39 = K , R40 = w1 , R41 = w2
• R42 = 0.2955242247
• R47 = 0.6794095683
• R43 = 0.1488743390
• R48 = 0.1494513492
• R44 = 0.2692667193
• R49 = 0.8650633667
which are the coefficients for Gauss-Legendre 10-point formula
• R45 = 0.4333953941
• R50 = 0.06667134431
• R46 = 0.2190863625
• R51 = 0.9739065285
Flags: F01 F02 F04 F06 F07 F09 F10
F01 & F04 indicate the kind of geodesic ( see above )
CF 02 if you want to compute the geodesic distance with formula
(1)
SF 02 if you want to compute the geodesic distance with
formula (2)
Subroutines: /
| 01 LBL "BETA" 02 DEG 03 CF 10 04 CF 07 05 51.041 06 STO 10 07 RDN 08 STO 21 09 X<>Y 10 STO 09 11 XEQ 00 12 STO 20 13 LBL 01 14 RCL 21 15 VIEW 21 16 XEQ 00 17 ENTER 18 ENTER 19 X<> 20 20 - 21 X#0? 22 / 23 RCL 09 24 RCL 21 25 STO 09 26 STO T 27 - 28 * 29 + 30 STO 21 31 X#Y? 32 GTO 01 33 STO 08 34 RTN 35 LBL "DIST" 36 LBL 12 37 CF 10 38 RCL 05 39 RCL 07 40 + 41 RCL 04 42 RCL 06 43 + 44 RCL 40 45 RCL 41 46 + 47 2 48 ST/ T 49 ST/ Z 50 / 51 FS? 02 52 X<> Z 53 STO 12 54 51.041 55 STO 10 56 RCL 04 57 FS? 02 58 RCL 05 59 STO 11 60 RCL 06 61 FS? 02 62 RCL 07 63 E-8 64 STO 22 65 SF 07 66 XEQ 14 67 D-R 68 ABS 69 CLD 70 RTN 71 GTO 12 72 LBL 00 73 STO 08 74 RCL 02 75 RCL 03 76 - 77 / 78 STO 39 79 RCL 06 80 RCL 04 81 STO 11 82 CF 06 83 XEQ 14 84 STO 19 85 RCL 07 |
86 RCL 05
87 FS? 02 88 GTO 11 89 X<>Y 90 1 91 P-R 92 X^2 93 RCL 39 94 - 95 SQRT 96 FS? 04 97 CHS 98 R-P 99 X<>Y 100 STO 41 101 RCL 05 102 1 103 P-R 104 X^2 105 RCL 39 106 - 107 SQRT 108 R-P 109 RDN 110 STO 40 111 LBL 11 112 STO 11 113 SF 06 114 XEQ 14 115 FS? 01 116 CHS 117 ST- 19 118 RCL 19 119 RTN 120 GTO 00 121 LBL 13 122 CLX 123 STO 28 124 CLX 125 RCL 21 126 STO 23 127 - 128 RCL 00 129 STO 24 130 ST+ X 131 / 132 STO 25 133 LBL 06 134 ST+ 23 135 RCL 10 136 STO 26 137 LBL 07 138 RCL 23 139 STO 27 140 RCL 25 141 RCL IND 26 142 * 143 ST+ 27 144 - 145 XEQ 02 146 X<> 27 147 XEQ 02 148 RCL 27 149 + 150 DSE 26 151 RCL IND 26 152 * 153 ST+ 28 154 DSE 26 155 GTO 07 156 RCL 25 157 ST+ 23 158 DSE 24 159 GTO 06 160 ST* 28 161 RCL 28 162 RTN 163 LBL 14 164 CLX 165 STO 18 166 CLX 167 RCL 11 168 STO 13 169 - 170 RCL 00 |
171 STO 14
172 ST+ X 173 / 174 STO 15 175 LBL 08 176 ST+ 13 177 RCL 10 178 STO 16 179 LBL 09 180 RCL 13 181 STO 17 182 RCL 15 183 RCL IND 16 184 * 185 ST+ 17 186 - 187 FS? 07 188 XEQ 05 189 FC? 07 190 XEQ 02 191 X<> 17 192 FS? 07 193 XEQ 05 194 FC? 07 195 XEQ 02 196 RCL 17 197 + 198 DSE 16 199 RCL IND 16 200 * 201 ST+ 18 202 DSE 16 203 GTO 09 204 RCL 15 205 ST+ 13 206 DSE 14 207 GTO 08 208 ST* 18 209 RCL 18 210 RTN 211 LBL 02 212 FS? 06 213 GTO 03 214 SIN 215 X^2 216 RCL 01 217 RCL 02 218 - 219 * 220 STO Y 221 RCL 02 222 + 223 RCL X 224 RCL 03 225 - 226 / 227 X<>Y 228 RCL 08 229 + 230 / 231 X<0? 232 SF 09 233 FS? 10 234 ABS 235 SQRT 236 RTN 237 LBL 03 238 FS? 02 239 GTO 04 240 COS 241 X^2 242 1 243 RCL 39 244 - 245 * 246 RCL 39 247 + 248 STO Y 249 RCL 03 250 RCL 02 251 - 252 * 253 RCL 02 254 + 255 ENTER |
256 CHS 257 RCL 01 258 + 259 / 260 X<>Y 261 / 262 RCL 02 263 RCL 03 264 - 265 / 266 X<0? 267 SF 09 268 FS? 10 269 ABS 270 SQRT 271 RTN 272 LBL 04 273 COS 274 X^2 275 RCL 02 276 RCL 03 277 - 278 * 279 ENTER 280 CHS 281 RCL 02 282 + 283 RCL 01 284 RCL Y 285 - 286 / 287 X<>Y 288 RCL 08 289 - 290 / 291 X<0? 292 SF 09 293 FS? 10 294 ABS 295 SQRT 296 RTN 297 LBL 05 298 STO 31 299 FC? 02 300 RCL 04 301 FS? 02 302 RCL 05 303 STO 21 304 CF 06 305 FS? 02 306 SF 06 307 XEQ 13 308 FS? 01 309 CHS 310 STO 20 311 RCL 40 312 FS? 02 313 RCL 04 314 STO 21 315 SF 10 316 RCL 12 317 STO 30 318 SF 06 319 FS? 02 320 CF 06 321 LBL 10 322 CF 09 323 RCL 30 324 FC? 02 325 XEQ 03 326 FS? 02 327 XEQ 02 328 STO 29 329 RCL 30 330 ENTER 331 XEQ 13 332 RCL 20 333 - 334 RCL 29 335 / 336 ST- 30 337 VIEW 30 338 RCL 30 339 X=0? 340 SIGN |
341 / 342 ABS 343 RCL 22 344 X<Y? 345 GTO 10 346 FS? 09 347 SF 41 348 RCL 30 349 STO 12 350 FS? 02 351 X<> 31 352 STO 30 353 FS? 02 354 GTO 05 355 SIN 356 X^2 357 1 358 RCL 39 359 - 360 * 361 SQRT 362 ASIN 363 STO 30 364 LBL 05 365 RCL 31 366 1 367 P-R 368 X^2 369 STO 32 370 RCL 02 371 * 372 X<>Y 373 X^2 374 STO 33 375 RCL 01 376 * 377 + 378 STO 36 379 RCL 03 380 - 381 RCL 30 382 COS 383 X^2 384 STO 34 385 * 386 RCL 03 387 * 388 RCL 01 389 RCL 03 390 - 391 / 392 RCL 30 393 SIN 394 X^2 395 STO 35 396 RCL 33 397 * 398 RCL 02 399 * 400 + 401 RCL 34 402 RCL 03 403 * 404 RCL 35 405 RCL 02 406 * 407 + 408 STO 37 409 RCL 01 410 X<>Y 411 - 412 RCL 34 413 RCL 35 414 * 415 X<>Y 416 / 417 RCL 32 418 * 419 RCL 01 420 * 421 RCL 02 422 RCL 03 423 - 424 X^2 425 * |
426 RCL 01 427 RCL 03 428 - 429 / 430 + 431 STO 38 432 RCL 32 433 RCL 33 434 * 435 RCL 35 436 * 437 RCL 01 438 RCL 02 439 - 440 X^2 441 * 442 RCL 03 443 * 444 RCL 01 445 RCL 03 446 - 447 / 448 RCL 36 449 RCL 03 450 - 451 / 452 RCL 32 453 RCL 34 454 * 455 RCL 02 456 * 457 + 458 RCL 01 459 RCL 37 460 - 461 RCL 33 462 * 463 RCL 01 464 ST* Y 465 RCL 03 466 - 467 / 468 + 469 STO 35 470 RCL 36 471 RCL 01 472 RCL 37 473 - 474 * 475 RCL 02 476 RCL 03 477 - 478 RCL 34 479 * 480 RCL 08 481 - 482 * 483 RCL 37 484 RCL 36 485 RCL 03 486 - 487 * 488 RCL 01 489 RCL 02 490 - 491 RCL 33 492 * 493 RCL 08 494 + 495 * 496 FS? 02 497 X<>Y 498 / 499 RCL 35 500 RCL 38 501 FS? 02 502 X<>Y 503 RCL Z 504 * 505 + 506 SQRT 507 RTN 508 END |
( 746 bytes / SIZE 052 )
| STACK | BETA-INPUT | BETA-OUTPUT | DIST-OUTPUT |
| Y | beta1 | / | / |
| X | beta2 | beta | s |
Where beta1 & beta2 are 2 initial guesses, beta the solution and s = geodesic distance.
Example1&2: The same as in paragraph 1°) a)
-With the first example CF 01 CF 02 ( or SF 02 ) and CF 04
, you'll get the same results for beta & s
-With the second example SF 01 CF 02 ( or SF 02 ) and CF
04, same results too.
Example3: The same ellipsoidal model of the Earth: semi-axes 6378.172 km , 6378.102 km , 6356.752 km ,
with the equatorial major axis in the direction 14°93W , 165°07E
-Calculate the geodesic distance between the U.S. Naval Observatory at
Washington (D.C.) L1 = 77°03'56"0 W ,
b1 = 38°55'17"2 N
and the "Observatoire de Paris" L2 = 2°20'13"8
E , b2 = 48°50'11"2 N
-The routines in paragraph 2°) b) and 3°) return the Jacobi parameters
u1 = -62°16090881
u2 = +17°30207991
store these 4 numbers R04 R06
respectively
v1 = +38°84397819
v2 = +48°83869959
into
R05 R07
-With n = 4 STO 00 , CF 01 CF 02 SF 04 , it yields beta = 101348.0008 and the geodesic distance s = 6181.626001 km
Notes:
-With an oblate ellipsoid of revolution with a = b = 6378.137 km
& c = 6356.752 km ,
it yields 12709.52608 km in example 2 and 6181.621921
km in example 3
which may be found more easily by Vincenty's formula.
-This routine is not perfect !
-The small number in R22 = E-8 may sometimes produce an infinite
loop inside LBL 10
-Replace it by a less small number ( line 63 ) or
-Stop the routine
-Set flag 21 , R/S
-Place 0 in X-register when the program stops ( line 337
)
-Clear flag 21 and R/S again.
-Moreover, the change of variable w cannot be used if the 2 points are on the equator... except for nearly antipodal points !
-For example, if L1 = b1 = 0 and L1 = 179°51' E , b2 = 0 you first find
u1 = 14°93015654
u2 = 194°7801551
v1 = 0
v2 = 0
-With n = 8 STO 00 it yields beta = 16815.92370
& the geodesic distance s = 20001.90110 km
c) Program(s)#3
-Here, the w-values are computed by the classical Runge-Kutta method of
order 4 ( LBL 03 )
and the geodesic distance itself by the Newton-Côtes 7-point
formula ( LBL 01 )
-We assume that the parameter "beta" has been already calculated by
the first part of the program above,
but "FBET" may also be used ( with a root-finding program ) to
solve f(beta) = 0
-Moreover, another version of "BETA" is also listed at the end of this paragraph ( see below )
-"DGD" uses the same change of argument of paragraph 1°)b) and flags
F01 & F04 must be set or cleared as before:
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
* *
*
CF 01-CF 04
SF 01-CF 04
CF 01-SF04
SF 01-SF 04
-Flag F02 is unused here, so "DGD" cannot be used to apply the second
formula for the distance.
-Here, we compute:
s = §u1u2
[ P + R ( dv/du )2 ]1/2 du
Data Registers: • R00 = n ( Registers R00 thru R08 are to be initialized before executing "DGD" - R00 to R07 for "FBET" )
• R01 = a
• R04 = u1
• R06 = u2
R09 = s
R11 = u R13 to R21: temp
• R02 = b
• R05 = v1
• R07 = v2
R10 = h
R12 = w
• R03 = c
>>> • R08 = Beta
R22 = K , R23 = w2
Flags: F00 F01 F04 F10
- F01 & F04 indicate the kind of geodesic ( see
above )
Subroutines: /
| 01 LBL "DGD" 02 CF 00 03 LBL 00 04 DEG 05 RCL 06 06 RCL 04 07 STO 11 08 - 09 RCL 00 10 STO 21 11 6 12 * 13 / 14 STO 10 15 RCL 08 16 RCL 02 17 RCL 03 18 - 19 / 20 STO 22 21 RCL 05 22 1 23 P-R 24 X^2 25 RCL Z 26 - 27 SQRT 28 R-P 29 X<>Y 30 STO 12 31 RCL 07 32 1 33 P-R 34 X^2 35 RCL 22 36 - 37 SQRT 38 FS? 04 39 CHS 40 R-P 41 X<>Y 42 STO 23 43 FS? 00 44 RTN 45 CLX 46 STO 09 47 SF 10 48 LBL 01 49 XEQ 02 50 41 51 FC?C 10 52 ST+ X 53 * 54 ST+ 09 55 XEQ 03 56 XEQ 02 57 216 58 * |
59 ST+ 09 60 XEQ 03 61 XEQ 02 62 27 63 * 64 ST+ 09 65 XEQ 03 66 XEQ 02 67 272 68 * 69 ST+ 09 70 XEQ 03 71 XEQ 02 72 27 73 * 74 ST+ 09 75 XEQ 03 76 XEQ 02 77 216 78 * 79 ST+ 09 80 XEQ 03 81 DSE 21 82 GTO 01 83 XEQ 02 84 41 85 * 86 RCL 09 87 + 88 RCL 10 89 * 90 140 91 / 92 D-R 93 ABS 94 STO 09 95 RTN 96 STO 00 97 GTO 00 98 LBL 02 99 RCL 11 100 1 101 P-R 102 X^2 103 STO 13 104 RCL 02 105 * 106 X<>Y 107 X^2 108 STO 14 109 RCL 01 110 * 111 + 112 STO 17 113 RCL 03 114 - 115 STO 18 116 RCL 12 |
117 SIN 118 1 119 RCL 22 120 - 121 SQRT 122 * 123 STO 16 124 ST* 16 125 ASIN 126 COS 127 X^2 128 STO 15 129 RCL 03 130 * 131 RCL 16 132 RCL 02 133 * 134 + 135 STO 19 136 RCL 01 137 X<>Y 138 - 139 STO 20 140 RCL 01 141 RCL 01 142 RCL 03 143 - 144 / 145 RCL 14 146 * 147 * 148 RCL 02 149 RCL 13 150 * 151 RCL 15 152 * 153 + 154 RCL 01 155 RCL 02 156 - 157 X^2 158 RCL 01 159 RCL 03 160 ST* Z 161 - 162 / 163 RCL 14 164 * 165 RCL 13 166 * 167 RCL 16 168 * 169 RCL 18 170 / 171 + 172 RCL 02 173 RCL 03 174 - |
175 X^2 176 RCL 01 177 ST* Y 178 RCL 03 179 - 180 / 181 RCL 15 182 * 183 RCL 16 184 * 185 RCL 13 186 * 187 RCL 20 188 ST* 17 189 / 190 RCL 02 191 RCL 14 192 * 193 RCL 16 194 * 195 + 196 RCL 01 197 RCL 03 198 ST- Y 199 X<>Y 200 / 201 RCL 15 202 * 203 RCL 18 204 * 205 + 206 RCL 01 207 RCL 02 208 - 209 RCL 14 210 * 211 RCL 08 212 + 213 RCL 19 214 * 215 ST* 18 216 CLX 217 RCL 02 218 RCL 03 219 - 220 RCL 15 221 * 222 RCL 08 223 - 224 RCL 17 225 * 226 RCL 18 227 / 228 * 229 + 230 SQRT 231 RTN 232 LBL 03 |
233 RCL 12 234 RCL 11 235 XEQ 04 236 STO 17 237 2 238 / 239 RCL 12 240 + 241 RCL 11 242 RCL 10 243 2 244 / 245 + 246 STO 18 247 XEQ 04 248 ST+ 17 249 ST+ 17 250 2 251 / 252 RCL 12 253 + 254 RCL 18 255 XEQ 04 256 ST+ 17 257 ST+ 17 258 RCL 12 259 + 260 RCL 10 261 RCL 11 262 + 263 STO 11 264 XEQ 04 265 RCL 17 266 + 267 6 268 / 269 ST+ 12 270 RTN 271 LBL 04 272 1 273 P-R 274 X^2 275 STO 13 276 RCL 02 277 * 278 X<>Y 279 X^2 280 STO 14 281 RCL 01 282 * 283 + 284 STO 16 285 RCL 03 286 - 287 ST/ 16 288 R^ 289 COS 290 X^2 |
291 1 292 RCL 22 293 STO T 294 - 295 * 296 + 297 STO 15 298 RCL 02 299 RCL 03 300 - 301 * 302 RCL 02 303 X<>Y 304 ST- Y 305 RCL 01 306 RCL 02 307 - 308 + 309 RCL 15 310 * 311 RCL 02 312 RCL 03 313 - 314 * 315 X<>Y 316 / 317 RCL 16 318 * 319 RCL 01 320 RCL 02 321 - 322 RCL 14 323 * 324 RCL 08 325 + 326 / 327 SQRT 328 RCL 10 329 * 330 FS? 01 331 CHS 332 RTN 333 LBL "FBET" 334 STO 08 335 SF 00 336 XEQ 00 337 6 338 ST* 21 339 LBL 05 340 XEQ 03 341 DSE 21 342 GTO 05 343 RCL 12 344 RCL 23 345 - 346 END |
( 475 bytes / SIZE 024 )
| STACK | DGD-INPUT | DGD-OUTPUT | FBET-INPUT | FBET-OUTPUT |
| Y | / | / | / | / |
| X | / | s | beta | f(beta) |
Where s = geodesic distance and f(beta) = R12 - R23 = 0 if we achieve a perfect accuracy
Example: The same ellipsoidal model of the Earth: semi-axes 6378.172 km , 6378.102 km , 6356.752 km ,
with the equatorial major axis in the direction 14°93W , 165°07E
-Calculate the geodesic distance between the U.S. Naval Observatory at
Washington (D.C.) L1 = 77°03'56"0 W ,
b1 = 38°55'17"2 N
and the "Observatoire de Paris" L2 = 2°20'13"8
E , b2 = 48°50'11"2 N
-The same Jacobi parameters
u1 = -62°16090881
u2 = +17°30207991
store these 4 numbers R04 R06
respectively
v1 = +38°84397819
v2 = +48°83869959
into
R05 R07
-If we already know that beta = 101348.0008 STO 08
-With n = 1 STO 00 , CF 01 SF 04
XEQ "DGD" >>>> s = 6181.533192 = R09 & R12 = 108°0852527 instead of R23 = 108°0854380
-To try another n-value, simply place it in X-register and R/S
n = 2 R/S >>>> s = 6181.625172 = R09 & R12 = 108°0854216 ---Execution time = 4mn36s---
n = 4 R/S >>>> s = 6181.625979 = R09 & R12 = 108°0854368 ---Execution time = 8mn59s---
n = 8 R/S >>>> s = 6181.626012 = R09 & R12 = 108°0854379
-So the results obtained with n = 8 or even n = 4 are correct: the errors are smaller than a few centimeters !
Notes:
-The geodesic distance is obtained much faster than with the first 2 programs
!
-So, this program is better for a real HP-41... and V41 too !
-However, calculating beta will be easier with "BETA" listed in paragraph
1°)-b)
-Nevertheless, you can also use "FBET" to solve f(beta) = 0
-In this case, it's unuseful to compute the distance s ( flag F00 is
set line 335 )
-For example 101348.0008 XEQ "FBET" gives
-0.000001200 with n = 4 in R00
and -0.000000100 with n = 8
in R00
-So, a root-finding program could call "FBET" as a subroutine to get
the correct beta-value.
Variant:
-Instead of using the classical RK4, we can also use a Runge-Kutta method
of order 6.
-For instance, replace lines 232 to 270 ( LBL 03 ) by
| 232 LBL 03 233 RCL 12 234 RCL 11 235 XEQ 04 236 STO 17 237 3 238 / 239 RCL 12 240 + 241 RCL 10 242 3 243 / 244 RCL 11 245 + 246 STO 24 247 XEQ 04 248 STO 18 249 1.5 250 / 251 RCL 12 252 + 253 RCL 10 254 1.5 255 / 256 RCL 11 257 + 258 XEQ 04 259 STO 19 260 CHS 261 RCL 18 |
262 4 263 * 264 + 265 RCL 17 266 + 267 12 268 / 269 RCL 12 270 + 271 RCL 24 272 XEQ 04 273 STO 20 274 90 275 * 276 RCL 19 277 35 278 * 279 + 280 RCL 18 281 110 282 * 283 - 284 RCL 17 285 25 286 * 287 + 288 48 289 / 290 RCL 12 291 + |
292 RCL 10 293 1.2 294 / 295 RCL 11 296 + 297 XEQ 04 298 STO 24 299 10 300 / 301 RCL 20 302 2 303 / 304 + 305 RCL 19 306 8 307 / 308 - 309 RCL 18 310 55 311 * 312 RCL 17 313 18 314 * 315 - 316 120 317 / 318 - 319 RCL 12 320 + 321 RCL 10 |
322 6 323 / 324 RCL 11 325 + 326 XEQ 04 327 STO 25 328 80 329 * 330 RCL 20 331 118 332 * 333 - 334 RCL 18 335 99 336 * 337 + 338 39 339 / 340 RCL 24 341 128 342 * 343 RCL 19 344 215 345 * 346 + 347 RCL 17 348 783 349 * 350 - 351 780 |
352 / 353 + 354 RCL 12 355 + 356 RCL 10 357 RCL 11 358 + 359 STO 11 360 XEQ 04 361 RCL 17 362 + 363 13 364 * 365 RCL 24 366 RCL 25 367 + 368 32 369 * 370 + 371 RCL 19 372 RCL 20 373 + 374 55 375 * 376 + 377 .5 378 % 379 ST+ 12 380 RTN |
( 645 bytes / SIZE 026 ) for
the whole program
-With the same example and n = 1 STO 00 , CF 01 SF 04
XEQ "DGD" >>>> s = 6181.537185 = R09 & R12 = 108°0854502 instead of R23 = 108°0854380
n = 2 R/S >>>> s = 6181.625364 = R09 & R12 = 108°0854381
n = 4 R/S >>>> s = 6181.625979 = R09 & R12 = 108°0854380
n = 8 R/S >>>> s = 6181.626012 = R09 & R12 = 108°0854380
-The convergence is obviously better for w but the precision is not really increased for the geodesic distance, except for small n-values.
-Here is also another variant of "BETA" that employs Ridder's method to
solve f(beta) = 0
-You have to place in registers X & Y 2 initial guesses
b & b' such that f(b).f(b') < 0
Data Registers: • R00 = n ( Registers R00 thru R07 and R26 thru R35 are to be initialized before executing "BETA" )
• R01 = a
• R04 = u1
• R06 = u2
R10 to R25: temp
• R02 = b
• R05 = v1
• R07 = v2
• R03 = c
>>> R08 = Beta
R09 = f(Beta)
• R26 = 0.2955242247
• R31 = 0.6794095683
• R27 = 0.1488743390
• R32 = 0.1494513492
• R28 = 0.2692667193
• R33 = 0.8650633667
• R29 = 0.4333953941
• R34 = 0.06667134431
• R30 = 0.2190863625
• R35 = 0.9739065285
Flags: F01 F04 F10
F01 & F04 must be used as above.
Subroutines: /
| 01 LBL "BETA" 02 DEG 03 35.025 04 STO 25 05 RDN 06 STO 11 07 X<>Y 08 STO 12 09 XEQ 00 10 STO 20 11 E-8 12 STO 24 13 RCL 11 14 XEQ 00 15 STO 09 16 LBL 01 17 VIEW 11 18 RCL 11 19 RCL 12 20 + 21 2 22 / 23 STO 21 24 XEQ 00 25 ABS 26 RCL 24 27 X>Y? 28 GTO 05 29 LASTX 30 STO 22 31 RCL 20 32 RCL 09 33 - 34 SIGN 35 * 36 RCL 22 37 X^2 38 RCL 09 39 RCL 20 40 * |
41 - 42 SQRT 43 / 44 RCL 21 45 RCL 12 46 - 47 * 48 RCL 21 49 + 50 STO 23 51 XEQ 00 52 ABS 53 RCL 24 54 X>Y? 55 GTO 06 56 RCL 22 57 LASTX 58 * 59 X<0? 60 GTO 02 61 RCL 20 62 LASTX 63 * 64 X<0? 65 GTO 03 66 LASTX 67 X<> 09 68 STO 20 69 RCL 23 70 X<> 11 71 STO 12 72 GTO 04 73 LBL 02 74 RCL 22 75 STO 20 76 RCL 21 77 STO 12 78 LBL 03 79 LASTX 80 STO 09 |
81 RCL 23 82 STO 11 83 LBL 04 84 RCL 12 85 RCL 11 86 X#Y? 87 GTO 01 88 GTO 13 89 LBL 05 90 RCL 21 91 GTO 07 92 LBL 06 93 RCL 23 94 LBL 07 95 STO 11 96 LASTX 97 STO 09 98 LBL 13 99 RCL 09 100 RCL 11 101 STO 08 102 CLD 103 RTN 104 LBL 00 105 STO 08 106 RCL 02 107 RCL 03 108 - 109 / 110 STO 10 111 RCL 06 112 RCL 04 113 ENTER 114 CF 10 115 XEQ 14 116 STO 19 117 RCL 07 118 1 119 P-R 120 X^2 |
121 RCL 10
122 - 123 SQRT 124 FS? 04 125 CHS 126 R-P 127 X<>Y 128 RCL 05 129 1 130 P-R 131 X^2 132 RCL 10 133 - 134 SQRT 135 R-P 136 SF 10 137 XEQ 14 138 FS? 01 139 CHS 140 ST- 19 141 RCL 19 142 RTN 143 GTO 00 144 LBL 14 145 CLX 146 STO 18 147 RDN 148 STO 13 149 - 150 RCL 00 151 STO 14 152 ST+ X 153 / 154 STO 15 155 LBL 08 156 ST+ 13 157 RCL 25 158 STO 16 159 LBL 09 160 RCL 13 |
161 STO 17 162 RCL 15 163 RCL IND 16 164 * 165 ST+ 17 166 - 167 XEQ 02 168 X<> 17 169 XEQ 02 170 RCL 17 171 + 172 DSE 16 173 RCL IND 16 174 * 175 ST+ 18 176 DSE 16 177 GTO 09 178 RCL 15 179 ST+ 13 180 DSE 14 181 GTO 08 182 RCL 18 183 * 184 RTN 185 LBL 02 186 FS? 10 187 GTO 03 188 SIN 189 X^2 190 RCL 01 191 RCL 02 192 - 193 * 194 STO Y 195 RCL 02 196 + 197 RCL X 198 RCL 03 199 - 200 / |
201 X<>Y 202 RCL 08 203 + 204 / 205 SQRT 206 RTN 207 LBL 03 208 COS 209 X^2 210 1 211 RCL 10 212 STO T 213 - 214 * 215 + 216 STO Y 217 RCL 03 218 RCL 02 219 - 220 * 221 RCL 02 222 + 223 ENTER 224 CHS 225 RCL 01 226 + 227 / 228 X<>Y 229 / 230 RCL 02 231 RCL 03 232 - 233 / 234 SQRT 235 RTN 236 END |
( 335 bytes / SIZE 036 )
| STACK | INPUTS | OUTPUTS |
| Y | b | f(beta) ~ 0 |
| X | b' | beta |
Example: Again the same one
-With n = 1 STO 00 , CF 01 SF 04 we find that 100000 < beta < 110000
100000 ENTER^
110000 XEQ "BETA" >>>>
beta = 101348.0011 X<>Y
f(beta) = -4 E-9
---Execution time = 6m42s---
-With n = 2 STO 00
100000 ENTER^
110000 XEQ "BETA" >>>>
beta = 101348.0010 X<>Y
f(beta) = -1 E-9
-With n = 4 STO 00
100000 ENTER^
110000 XEQ "BETA" >>>>
beta = 101348.0008 X<>Y
f(beta) = -4 E-9
Notes:
-"BETA" uses Ridder's method as root-finding routine ( lines 01 to 103
)
-Though not very fast, Gaussian integration is however less slow than
using "FBET"
-When the program stops, register R08 = beta, ready for using "DGD"
above.
-In order to find 2 values for which the results have opposite signs,
you can XEQ 00 provided R25 = 35.025
= control number of the Gauss-Legendre coefficients.
2°) Longitude & Latitude >>>>
Cartesian Coordinates
a) Geodetic Coordinates
-The formulae between the geodetic longitude l
& latitude b and the cartesian coordinates
x , y , z
of a point on a triaxial ellipsoid are given in references [2]
& [3]
x = w Cos b Cos l
y = w ( 1 - e2 ) Cos b Sin l
In fact l + 14°93 instead
of l
z = w ( 1 - e' 2 ) Sin b
where e2 = 1 - b2/a2 , e' 2 = 1 - c2/a2 , w = a / ( 1 - e' 2 Sin2b - e2 Cos2 b Sin2l ) 1/2
-This program is written for the Earth, assuming that the semi-axes are:
a = 6378.172 km with
the equatorial major axis in the direction 14°93W
, 165°07E
b = 6378.102 km equatorial
minor semi-axis
c = 6356.752 km polar
minor semi-axis
-Change these values if you know better approximations ( lines 08-14-17-19
).
-The longitudes are measured positively Eastward from the meridian of
Greenwich ( otherwise, replace line 09 + by a - ).
Data Registers: /
Flags: /
Subroutine: "S-R" Spherical-Rectangular
conversion ( cf "Spherical Coordinates for the HP-41" paragraph 1°-a)
)
| 01 LBL "LB-XYZ" 02 DEG 03 HR 04 STO N 05 COS 06 X<>Y 07 HR 08 14.93 09 + 10 STO M 11 SIN |
12 * 13 X^2 14 6356.752 15 X^2 16 X<>Y 17 6378.102 18 X^2 19 6378.172 20 STO P 21 X^2 22 ST- Y |
23 ST- T 24 ST/ T 25 / 26 STO O 27 * 28 X<>Y 29 RCL N 30 SIN 31 X^2 32 * 33 + |
34 PI 35 SIGN 36 ST+ O 37 ST+ Z 38 + 39 SQRT 40 RCL P 41 X<>Y 42 / 43 RCL M 44 RCL N |
45 X<> Z 46 XEQ "S-R" 47 R^ 48 ST* T 49 X<> O 50 ST* Z 51 RDN 52 CLA 53 END |
( 110 bytes / SIZE 000 )
| STACK | INPUTS | OUTPUTS |
| Z | / | z |
| Y | l ( ° ' " ) | y |
| X | b ( ° ' " ) | x |
---Execution time = 6s---
Example: The Observatoire de Paris l = 2°20'13"8 E , b = 48°50'11"2 N
2.20138 ENTER^
48.50112 XEQ "LB-XYZ" >>>>
x = 4016.614852 km
RDN y = 1248.484248 km
RDN z = 4778.596570 km
Note:
-The longitudes given in almanachs are more likely geocentric longitudes.
-The program below employs this definition.
b) Astronomical Coordinates
-If the longitude l is geocentric i-e equals
the angle between the meridian plane passing through the observer and
the Prime Meridian plane,
and if the latitude b is the angle
between the local vertical in the local meridian plane and the plane
of the Equator,
the rectangular coordinates x , y , z of a point P on the triaxial
ellipsoid may be found as follows:
r2 = b2 + ( a2-b2 ) / [ 1 + (a2/b2) Tan2 L ] with r = OQ where Q = point of intersection of the equator and the local meridian
R2 = c2 + ( r2-c2 ) / [ 1 + (c2/r2) Tan2 b ] and L = l + 14°93
x = R Cos µ Cos L
y = R Cos µ Sin L
where R = OP and Tan µ = (c2/r2)
Tan b
z = R Sin µ
Data Registers: /
Flags: /
Subroutine: "S-R" Spherical-Rectangular
conversion ( cf "Spherical Coordinates for the HP-41" paragraph 1°-a)
)
| 01 LBL "LB-XYZ" 02 DEG 03 HR 04 STO N 05 X<>Y 06 HR 07 14.93 08 + 09 STO M 10 1 11 P-R 12 X^2 13 X<>Y |
14 6378.172 15 STO O 16 6378.102 17 STO P 18 / 19 * 20 X^2 21 + 22 / 23 RCL O 24 X^2 25 RCL P 26 X^2 |
27 STO T 28 - 29 * 30 + 31 RCL N 32 1 33 P-R 34 X^2 35 X<>Y 36 X^2 37 6356.752 38 X^2 39 STO P |
40 R^ 41 STO O 42 / 43 * 44 + 45 / 46 RCL O 47 RCL P 48 - 49 * 50 RCL P 51 ST/ O 52 + |
53 SQRT 54 RCL N 55 1 56 P-R 57 RCL O 58 * 59 R-P 60 X<> M 61 R^ 62 XEQ "S-R" 63 CLA 64 END |
( 121 bytes / SIZE 000 )
| STACK | INPUTS | OUTPUTS |
| Z | / | z |
| Y | l ( ° ' " ) | y |
| X | b ( ° ' " ) | x |
---Execution time = 7s---
Example: Again the Observatoire de Paris l = 2°20'13"8 E , b = 48°50'11"2 N
2.20138 ENTER^
48.50112 XEQ "LB-XYZ" >>>>
x = 4016.607084 km
RDN y = 1248.509238 km
RDN z = 4778.596570 km
-We also get R = OP = 6366.073591 in register T
Notes:
-Strictly speaking, the perpendicular to the local meridian plane is not
perpendicular to the triaxial ellipsoid,
but the difference is negligible: both versions return the same
z-values
-This would not be the case if the flattening were large.
-One could use TAN and ATAN instead of the polar
<> rectangular conversion,
but there would be an OUT OF RANGE if the latitude = 90° and
F24 is clear.
-Moreover, the HP-41 wrongly returns ATAN ( -9.999999999 E99 )
= 90° instead of -90°
-So, using P-R & R-P is preferable
-I've used this program to get the coordinates in the examples of paragraph 1
-If the latitude b is strictly geodetic, i-e if it is the angle between
the normal to the ellipsoid and the equatorial plane,
the geocentric longitude b' may be computed by the formula:
Tan b' = (c2/a2) (Tan b ) [ Cos2 L + (a4/b4) Sin2 L ] 1/2 where L is the geocentric longitude.
-This formula is used in paragraph 4°)d)
c) Geocentric
Coordinates
Formula:
x = r Cos L Cos B
y = r Sin L Cos B r > 0 ; r2 [ ( Cos L Cos B )2 / a2 + ( Sin L Cos B )2 / b2 + ( Sin B )2 / c2 ] = 1
z = r Sin B
Data Registers: R01 = a2 = 6378.172^2 R02 = b2 = 6378.102^2 R03 = c2 = 6356.752^2 R09 & R10: temp
Flags: /
Subroutines: /
| 01 LBL "LB-XYZ" 02 DEG 03 1 04 P-R 05 RCL Z 06 X<>Y 07 P-R 08 STO 09 |
09 RDN 10 STO 10 11 X<>Y 12 ENTER 13 X^2 14 6356.752 15 X^2 16 STO 03 |
17 / 18 RCL 10 19 X^2 20 6378.102 21 X^2 22 STO 02 23 / 24 + |
25 RCL 09 26 X^2 27 6378.172 28 X^2 29 STO 01 30 / 31 + |
32 SQRT 33 ST/ 09 34 ST/ 10 35 / 36 RCL 10 37 RCL 09 38 END |
( 73 bytes / SIZE 011 )
| STACK | INPUTS | OUTPUTS |
| Z | / | z |
| Y | L ( deg ) | y |
| X | B ( deg ) | x |
Example: L = 124.77° B = 37.64°
124.77 ENTER^
37.64 XEQ "LB-XYZ" >>>> x = -2876.665405 km ---Execution time = 3.3s---
RDN y = 4143.606935 km
RDN z = 3890.225649 km
Note:
-This program stores a2 , b2 & c2 in registers R01 R02 & R03
3°) Cartesian Coordinates <> Jacobi Parameters
( u , v )
-The coordinates of a point on the ellipsoid x2 / a + y2 / b + z2 / c = 1 ( a > b > c ) are seldom given by ( u , v ) where
x = [ a / ( a - c ) ]1/2 Cos u
( a - b Sin2 v - c Cos2 v )1/2
y = b1/2 Sin u
Cos v
z = [ c / ( a - c ) ]1/2 Sin v
( a Sin2 u + b Cos2 u - c )1/2
-The 2 routines below perform the conversions
-In the conversion ( x , y , z ) >>> ( u , v ) , we have to solve the equation:
( b2 - b.c ) Sin4 v + [ b.y2 + b.c - b2 - a.y2 - b.(a-c)/c z2 ] Sin2 v + b.(a-c)/c z2 = 0
-So, the 2nd part of the routine calls "P2" ( quadratic equation )
-The sign of v is determined by the sign of z ( lines 100-101-102 )
Data Registers: R00 & R04 thru R08: unused ( Registers R01 thru R03 are to be initialized before executing "UV-XYZ" & "XYZ-UV" )
• R01 = a
R09 to R11: temp
• R02 = b
• R03 = c
x2 / a + y2 / b + z2 / c = 1
( a > b > c )
Flag: F00
Subroutine: "P2" ( quadratic equations
) - cf "Polynomials for the HP-41" paragraph 1°)-a)
| 01 LBL "UV-XYZ" 02 1 03 P-R 04 STO 09 05 X^2 06 RCL 02 07 * 08 X<>Y 09 STO 10 10 X^2 11 RCL 01 12 * 13 + 14 RCL 03 15 ST- Y 16 * 17 SQRT 18 STO 11 19 SIGN 20 P-R 21 ST* 10 22 X^2 23 RCL 03 24 * 25 X<>Y 26 ST* 11 27 X^2 |
28 RCL 02 29 * 30 + 31 CHS 32 RCL 01 33 ST+ Y 34 * 35 SQRT 36 RCL 01 37 RCL 03 38 - 39 SQRT 40 ST/ 11 41 / 42 ST* 09 43 RCL 11 44 RCL 10 45 RCL 02 46 SQRT 47 * 48 RCL 09 49 RTN 50 LBL "XYZ-UV" 51 STO 09 52 RDN 53 STO 10 54 X<>Y |
55 STO 11 56 RCL 02 57 RCL 03 58 - 59 RCL 02 60 * 61 ENTER^ 62 CHS 63 RCL 02 64 RCL 01 65 - 66 RCL 10 67 X^2 68 * 69 + 70 RCL 01 71 RCL 03 72 ST- Y 73 / 74 RCL 02 75 * 76 RCL 11 77 X^2 78 * 79 ST- Y 80 R^ 81 X=0? |
82 GTO 01 83 RDN 84 XEQ "P2" 85 FS? 00 86 SF 41 87 X>Y? 88 X<>Y 89 X<0? 90 X<>Y 91 GTO 02 92 LBL 01 93 X<> Z 94 / 95 CHS 96 LBL 02 97 ENTER^ 98 SQRT 99 ASIN 100 RCL 11 101 SIGN 102 * 103 X<>Y 104 RCL 03 105 RCL 02 106 - 107 * 108 RCL 03 |
109 - 110 RCL 01 111 ST+ Y 112 * 113 SQRT 114 RCL 01 115 RCL 03 116 - 117 SQRT 118 RCL 09 119 * 120 X<>Y 121 / 122 RCL 10 123 RCL 02 124 SQRT 125 / 126 R^ 127 COS 128 X=0? 129 SIGN 130 / 131 X<>Y 132 R-P 133 RDN 134 END |
( 171 bytes / SIZE 012 )
| STACK | INPUTS1 | OUTPUTS1 | INPUTS2 | OUTPUTS2 |
| Z | / | z | z | v |
| Y | v | y | y | v |
| X | u | x | x | u |
Example: The ellipsoid x2
/ 41 + y2 / 37 + z2 / 35 = 1
u = +10° , v = -15°
41 STO 01
37 STO 02
35 STO 03
-Set the HP-41 in DEG mode: XEQ "DEG"
15 CHS ENTER^
10 XEQ "UV-XYZ" >>>>
x = 6.235047001
RDN y = 1.020269420
RDN z = -0.910302041
-With these values in registers X , Y , Z ( RDN RDN )
XEQ "XYZ-UV" ( or simply R/S ) gives again u = 10 RDN v = -15 , with small roudoff-errors in the last decimals.
Notes:
-These programs work in all angular modes.
-"XYZ-UV" does not check that ( x , y , z ) is on the ellipsoid.
-Lines 80 to 83 and 91 to 96 are only useful if b = c
-If the quadratic equation has 2 complex conjugate roots - which should
not happen - the HP-41 displays NONEXISTENT ( line 86 )
-Registers R00 & R04 to R08 are unused because they contain important
data for the main program of paragraph 1°)
-Otherwise, R09-R10-R11 could be replaced by R00-R04-R05.
4°) Approximate Methods
a) Program#1
-The center O of the ellipsoid and 2 points M & N define a plane (
at least if they are not aligned ) that intersects the ellipsoid.
and we can compute the distance of the arc of ellipse MN thus
defined.
-This will not be - in general - the geodesic distance, but it will
remain slightly larger if the flattenings are small.
-The following "TGD" uses this method for the Earth, assuming that the semi-axes are:
a = 6378.172 km with
the equatorial major axis in the direction 14°93W
, 165°07E
b = 6378.102 km equatorial
minor semi-axis
c = 6356.752 km polar
minor semi-axis
( same values as in paragraph 2 )
-A point P on the arc MN may be determined by the vectorial equality:
OP = x u + y v where u = OM / || OM || and v = ON / || ON ||
-We have x = R Sin ( W - µ ) / Sin W and y = R Sin µ / Sin W
with W = the angle ( OM , ON ) µ = the angle ( OM , OP ) and R = OP
-Writing that P(X,Y,Z) is on the ellipsoid X2/a2 + Y2/b2 + Z2/c2 = 1, we find that
R / Sin W = 1 / ( A2 + B2 + C2 )1/2 with
A = [ x1 Sin ( W - µ )
+ x2 Sin µ ] / a
B = [ y1 Sin ( W - µ
) + y2 Sin µ ] / b
where u ( x1 , y1 , z1 )
& v ( x2 , y2 , z2 )
C = [ z1 Sin ( W - µ
) + z2 Sin µ ] / c
-Then, we compute dR / dµ and the arc length is
§0W [ R2 + ( dR/dµ )2
] 1/2 dµ
Data Registers: R00 to R09 & R20 to R37: temp ( Registers R10 thru R19 are to be initialized before executing "TGD" )
R03 = n = the number of subintervalls in which the Gaussian formula is applied. n = 1 seems enough ( line 56 )
R20 = a , R21 = b , R22 = c , R23 to R28 contain the coordinates of the unit vectors u & v , R29 = W
• R10 = 0.2955242247
• R11 = 0.1488743390
• R12 = 0.2692667193
• R13 = 0.4333953941
• R14 = 0.2190863625
• R15 = 0.6794095683
• R16 = 0.1494513492
• R17 = 0.8650633667
• R18 = 0.06667134431
• R19 = 0.9739065285
Flags: /
Subroutines: "LB-XYZ" cf paragraph 2°)b)
above
"GL10" Gauss-Legendre 10-point formula cf "Gaussian Integration
for the HP-41" paragraph 1°)-b)
Note:
-If you have the M-Code routines GX & GW no register must be
initialized and you can replace
registers R20 , ......... , R38 by R10 , ..............
, R28 in the listing below
| 01 LBL "TGD" 02 DEG 03 STO 08 04 RDN 05 STO 07 06 RDN 07 XEQ "LB-XYZ" 08 STO 23 09 X<>Y 10 STO 24 11 R-P 12 RCL Z 13 STO 25 14 R-P 15 ST/ 23 16 ST/ 24 17 ST/ 25 18 RCL 07 19 RCL 08 20 XEQ "LB-XYZ" 21 STO 26 22 X<>Y 23 STO 27 24 R-P 25 RCL Z 26 STO 28 27 R-P 28 ST/ 26 29 ST/ 27 30 ST/ 28 31 RCL 23 32 RCL 26 |
33 * 3 4 RCL 24 35 RCL 27 36 * 37 + 38 RCL 25 39 RCL 28 40 * 41 + 42 ACOS 43 STO 02 44 STO 29 45 6378.172 46 STO 20 47 .07 48 - 49 STO 21 50 6356.752 51 STO 22 52 CLX 53 STO 01 54 "S" 55 ASTO 00 56 SIGN 57 LBL 00 58 STO 03 59 XEQ "GL10" 60 RCL 29 61 SIN 62 * 63 D-R 64 RTN |
65 GTO 00 66 LBL "S" 67 STO Y 68 1 69 P-R 70 STO 31 71 RDN 72 STO 30 73 RCL 26 74 * 75 RCL 29 76 RCL Z 77 - 78 1 79 P-R 80 STO 33 81 RDN 82 STO 32 83 RCL 23 84 * 85 + 86 RCL 20 87 / 88 STO 34 89 X^2 90 RCL 32 91 RCL 24 92 * 93 RCL 27 94 RCL 30 95 * 96 + |
97 RCL 21 98 / 99 STO 35 100 X^2 101 + 102 RCL 25 103 RCL 32 104 * 105 RCL 28 106 RCL 30 107 * 108 + 109 RCL 22 110 / 111 STO 36 112 X^2 113 + 114 STO 37 115 RCL 26 116 RCL 31 117 * 118 RCL 23 119 RCL 33 120 * 121 - 122 RCL 34 123 * 124 RCL 20 125 / 126 RCL 27 127 RCL 31 128 * |
129 RCL 24 130 RCL 33 131 * 132 - 133 RCL 35 134 * 135 RCL 21 136 / 137 + 138 RCL 28 139 RCL 31 140 * 141 RCL 25 142 RCL 33 143 * 144 - 145 RCL 36 146 * 147 RCL 22 148 / 149 + 150 X^2 151 RCL 37 152 3 153 Y^X 154 / 155 RCL 37 156 1/X 157 + 158 SQRT 159 RTN 160 END |
( 279 bytes / SIZE 038 )
| STACK | INPUTS | OUTPUTS |
| T | l1 ( ° ' " ) | / |
| Z | b1 ( ° ' " ) | / |
| Y | l2 ( ° ' " ) | / |
| X | b2 ( ° ' " ) | D ( km ) |
Example: Calculate
the geodesic distance between the U.S. Naval Observatory at Washington (D.C.)
and the "Observatoire de Paris"
l1 = 77°03'56"0
W , b1 = 38°55'17"2 N
l2 = 2°20'13"8
E , b2 =
48°50'11"2 N
-77.0356 ENTER^
38.55172 ENTER^
2.20138 ENTER^
48.50112 XEQ "TGD" >>>>
D = 6181.628624 km
---Execution time = 67s---
-To perform the calculation with another n-value, simply place it in X-register
and R/S
-For instance: 2 R/S >>>> D = 6181.628624
km , same result !
Notes:
-With this example, we've found in paragraph 1°)-b) s = 6181.626001
km
-So, this simplified method gives an error < 3 meters in this case.
-For the geodesic distance between Washington & Cape Town Observatories,
it yields D = 12709.56721 km
whereas the exact value was s = 12709.56545 km,
an error < 2 meters for D
-However, with ( -40° , -40° ) ( +120° , +50° ) it yields
18095.54550 km
but the exact geodesic distance is 18095.48497 km
-So, the errors may reach at least 60 meters.
-They seem to remain relatively small ( except for nearly antipodal points
)
but Andoyer's formula of order 2 ( for an ellipsoid of revolution
) remains preferable...
-If you apply this method to the first example at the top of this page
( modify the corresponding lines in "TGD" and "LB-XYZ" ),
you get, with n = 4 , D = 8.594880194 instead
of s = 8.594822580
-The difference is more important here since the flattenings are larger
too.
-Gauss-Legendre 10-point formula may of course be replaced by another
integration formula...
b) Program#2
-To get a better accuracy, we can also divide the arc MN in several parts
-In the following program, 2 points P & Q are inserted between M
& N
and we compute the sum of 3 arc lentgths MP + PQ + QN
-P and Q are chosen to minimize this sum.
Ph
Qh
M ---------------------
P0 -------------------- Q0 ---------------------N
P -h
Q -h
-If W is the angle ( OM , ON ) , P0 & Q0 are in the plane ( OMN ) and ( OM , OP0 ) = ( OP0 , OQ0 ) = ( OQ0 , ON ) = W / 3
-The points Ph , Qh are above the plane ( OMN
) with ( OP0 , OPh ) = ( OQ0
, OQh ) = h
-The points P -h , Q -h are below the plane
( OMN ) with ( OP0 , OP -h ) = ( OQ0
, OQ -h ) = h
where h is a small angle ( I've chosen h = W / 1000 , lines 64-65 )
----------------------------------------------------------------------------------------------------------------------------------------
-Let a point P ( on the ellipsoid ) with ( OM , OP' ) = µ , ( OP' , OP ) = h , P' = orthogonal projection of P on the plane (OMN)
OP = x u + y v + z w where u = OM / || OM || , v = ON / || ON || , w = u x v ( x = cross-product )
-We have, with R = OP
x = R Cos h Sin ( W - µ ) / Sin W
y = R Cos h Sin µ / Sin W
z = R Sin h / Sin W
-Writing that P(X,Y,Z) is on the ellipsoid X2/a2 + Y2/b2 + Z2/c2 = 1, we find that
R / Sin W = 1 / ( A2 + B2 + C2 )1/2 with
A = [ x1 Cos h Sin ( W - µ
) + x2 Cos h Sin µ + ( y1z2 - y2z1
) Sin h ] / a
B = [ y1 Cos h Sin ( W - µ
) + y2 Cos h Sin µ + ( z1x2 - z2x1
) Sin h ] / b
where u ( x1 , y1 , z1 )
& v ( x2 , y2 , z2 )
C = [ z1 Cos h Sin ( W - µ
) + z2 Cos h Sin µ + ( x1y2 - x2y1
) Sin h ] / c
-----------------------------------------------------------------------------------------------------------------------------------------
-Several "distances" are calculated ( all the following "arcs" are elliptic arcs ).
f(0,0) = arc(MP0) + arc(P0Q0)
+ arc(Q0N)
f(h,0) = arc(MPh) + arc(PhQ0)
+ arc(Q0N)
f(-h,0) = arc(MP-h) + arc(P-hQ0)
+ arc(Q0N)
f(0,h) = arc(MP0) + arc(P0Qh)
+ arc(QhN)
f(0,-h) = arc(MP0) + arc(P0Q-h)
+ arc(Q-hN)
f(h,h) = arc(MPh) + arc(PhQh)
+ arc(QhN)
-f is approximated by a polynomial ( in 2 variables ) of degree 2
-Finally, the minimum is found after solving a 2x2 linear system.
Data Registers: R00 = n = number of subintervals of integration
R01 = a
R05 = x1
R08 = x2
R21 = h
R48 = f(h,0) R51 = f(0,-h)
R02 = b
R06 = y1
R09 = y2
R22 = W
R49 = f(-h,0) R52 = f(h,h)
R03 = c
R07 = z1
R10 = z2
R47 = f(0,0)
R50 = f(0,h)
>>>> When the program stops, R04 = s
Flag: F00
CF 00 for the Earth: n = 1 is sufficient and R01-R02-R03
are initialized ( lines 08 to 17 )
SF 00 for another ellipsoid: you'll have to initialize
R00-R01-R02-R03
Subroutine: "S-R" Spherical-Rectangular
conversion ( cf "Spherical Coordinates for the HP-41" paragraph 1°-a)
)
-Line 259 is a three-byte GTO 10 ( not really important )
-LBL 14 is actually the routine "LB-XYZ" listed in paragraph 2°)b)
- slightly modified.
| 01 LBL "TGD2" 02 STO 10 03 RDN 04 STO 09 05 RDN 06 FS? 00 07 GTO 00 08 6378.172 09 STO 01 10 .07 11 - 12 STO 02 13 21.35 14 - 15 STO 03 16 SIGN 17 STO 00 18 RDN 19 LBL 00 20 XEQ 14 21 STO 05 22 X<>Y 23 STO 06 24 R-P 25 RCL Z 26 STO 07 27 R-P 28 ST/ 05 29 ST/ 06 30 ST/ 07 31 RCL 09 32 RCL 10 33 XEQ 14 34 STO 08 35 X<>Y 36 STO 09 37 R-P 38 RCL Z 39 STO 10 40 R-P 41 ST/ 08 42 ST/ 09 43 ST/ 10 44 LBL 10 45 RCL 05 46 STO 41 47 RCL 08 48 STO 44 49 * 50 RCL 06 51 STO 42 52 RCL 09 53 STO 45 54 * 55 + 56 RCL 07 57 STO 43 58 RCL 10 59 STO 46 60 * 61 + 62 ACOS 63 STO 22 64 E3 65 / 66 STO 21 67 RCL 22 68 3 69 / 70 XEQ 11 71 23.032009 72 REGMOVE 73 RCL 22 74 ST+ X 75 3 76 / 77 XEQ 11 78 RCL 38 79 STO 08 80 RCL 39 81 STO 09 82 RCL 40 83 STO 10 84 XEQ 12 85 STO 48 86 STO 52 87 RCL 32 88 STO 08 89 RCL 33 90 STO 09 91 RCL 34 92 STO 10 93 XEQ 12 94 STO 49 95 RCL 35 96 STO 08 97 RCL 36 |
98 STO 09 99 RCL 37 100 STO 10 101 XEQ 12 102 STO 47 103 STO 50 104 STO 51 105 RCL 26 106 STO 05 107 RCL 27 108 STO 06 109 RCL 28 110 STO 07 111 XEQ 12 112 ST+ 47 113 RCL 29 114 STO 05 115 RCL 30 116 STO 06 117 RCL 31 118 STO 07 119 XEQ 12 120 ST+ 50 121 RCL 23 122 STO 05 123 RCL 24 124 STO 06 125 RCL 25 126 STO 07 127 XEQ 12 128 ST+ 51 129 RCL 44 130 STO 08 131 RCL 45 132 STO 09 133 RCL 46 134 STO 10 135 XEQ 12 136 ST+ 51 137 RCL 29 138 STO 05 139 RCL 30 140 STO 06 141 RCL 31 142 STO 07 143 XEQ 12 144 ST+ 50 145 ST+ 52 146 RCL 26 147 STO 05 148 RCL 27 149 STO 06 150 RCL 28 151 STO 07 152 XEQ 12 153 ST+ 47 154 ST+ 48 155 ST+ 49 156 RCL 32 157 STO 08 158 RCL 33 159 STO 09 160 RCL 34 161 STO 10 162 XEQ 12 163 ST+ 49 164 RCL 38 165 STO 08 166 RCL 39 167 STO 09 168 RCL 40 169 STO 10 170 XEQ 12 171 ST+ 48 172 RCL 29 173 STO 05 174 RCL 30 175 STO 06 176 RCL 31 177 STO 07 178 XEQ 12 179 ST+ 52 180 RCL 48 181 RCL 49 182 + 183 RCL 47 184 ST+ X 185 - 186 STO 05 187 RCL 50 188 RCL 51 189 - 190 STO 06 191 * 192 RCL 52 193 RCL 48 194 - |
195 RCL 47 196 + 197 RCL 50 198 - 199 STO 07 200 RCL 48 201 RCL 49 202 - 203 STO 09 204 * 205 - 206 STO 11 207 RCL 50 208 RCL 51 209 + 210 RCL 47 211 ST+ X 212 - 213 STO 08 214 RCL 09 215 * 216 RCL 06 217 RCL 07 218 * 219 - 220 RCL 07 221 X^2 222 RCL 05 223 RCL 08 224 * 225 - 226 ST+ X 227 ST/ 11 228 / 229 STO 10 230 RCL 05 231 * 232 RCL 09 233 + 234 RCL 10 235 * 236 RCL 08 237 RCL 11 238 * 239 RCL 06 240 + 241 RCL 11 242 * 243 + 244 2 245 / 246 RCL 10 247 RCL 11 248 * 249 RCL 07 250 * 251 + 252 RCL 47 253 + 254 STO 04 255 RTN 256 STO 00 257 41.005006 258 REGMOVE 259 GTO 10 260 LBL 11 261 STO 11 262 SIN 263 STO 13 264 RCL 08 265 * 266 RCL 22 267 RCL 11 268 - 269 SIN 270 STO 14 271 RCL 05 272 * 273 + 274 STO 26 275 RCL 21 276 COS 277 STO 15 278 * 279 RCL 06 280 RCL 10 281 * 282 RCL 07 283 RCL 09 284 * 285 - 286 RCL 21 287 SIN 288 STO 12 289 * 290 ST+ Z 291 - |
292 STO 23 293 X<>Y 294 STO 29 295 RCL 06 296 RCL 14 297 * 298 RCL 09 299 RCL 13 300 * 301 + 302 STO 27 303 RCL 15 304 * 305 RCL 07 306 RCL 08 307 * 308 RCL 05 309 RCL 10 310 * 311 - 312 RCL 12 313 * 314 ST+ Z 315 - 316 STO 24 317 X<>Y 318 STO 30 319 RCL 07 320 RCL 14 321 * 322 RCL 10 323 RCL 13 324 * 325 + 326 STO 28 327 RCL 15 328 * 329 RCL 05 330 RCL 09 331 * 332 RCL 06 333 RCL 08 334 * 335 - 336 RCL 12 337 * 338 ST+ Z 339 - 340 STO 25 341 X<>Y 342 STO 31 343 RCL 30 344 R-P 345 RCL 29 346 R-P 347 ST/ 29 348 ST/ 30 349 ST/ 31 350 RCL 28 351 RCL 27 352 R-P 353 RCL 26 354 R-P 355 ST/ 26 356 ST/ 27 357 ST/ 28 358 RCL 25 359 RCL 24 360 R-P 361 RCL 23 362 R-P 363 ST/ 23 364 ST/ 24 365 ST/ 25 366 RTN 367 LBL 12 368 CLX 369 STO 04 370 RCL 05 371 RCL 08 372 * 373 RCL 06 374 RCL 09 375 * 376 + 377 RCL 07 378 RCL 10 379 * 380 + 381 ACOS 382 STO 11 383 RCL 00 384 STO M 385 ST+ X 386 / 387 STO 12 388 .6 |
389 SQRT 390 * 391 STO 16 392 LBL 01 393 RCL 12 394 RCL 16 395 - 396 XEQ 13 397 ST+ 04 398 RCL 12 399 XEQ 13 400 1.6 401 * 402 ST+ 04 403 RCL 12 404 RCL 16 405 + 406 XEQ 13 407 ST+ 04 408 RCL 11 409 RCL 00 410 / 411 ST+ 12 412 DSE M 413 GTO 01 414 RCL 04 415 * 416 RCL 11 417 SIN 418 * 419 3.6 420 / 421 D-R 422 STO 04 423 RTN 424 LBL 13 425 STO 13 426 1 427 P-R 428 STO 15 429 RDN 430 STO 14 431 RCL 08 432 * 433 RCL 11 434 RCL 13 435 - 436 1 437 P-R 438 STO 18 439 RDN 440 STO 17 441 RCL 05 442 * 443 + 444 RCL 01 445 / 446 STO 19 447 X^2 448 RCL 06 449 RCL 17 450 * 451 RCL 09 452 RCL 14 453 * 454 + 455 RCL 02 456 / 457 STO 20 458 X^2 459 + 460 RCL 07 461 RCL 17 462 * 463 RCL 10 464 RCL 14 465 * 466 + 467 RCL 03 468 / 469 STO 13 470 X^2 471 + 472 STO 14 473 RCL 08 474 RCL 15 475 * 476 RCL 05 477 RCL 18 478 * 479 - 480 RCL 19 481 * 482 RCL 01 483 / 484 RCL 09 485 RCL 15 |
486 * 487 RCL 06 488 RCL 18 489 * 490 - 491 RCL 20 492 * 493 RCL 02 494 / 495 + 496 RCL 10 497 RCL 15 498 * 499 RCL 07 500 RCL 18 501 * 502 - 503 RCL 13 504 * 505 RCL 03 506 / 507 + 508 X^2 509 RCL 14 510 3 511 Y^X 512 / 513 RCL 14 514 1/X 515 + 516 SQRT 517 RTN 518 LBL 14 519 DEG 520 HR 521 STO N 522 X<>Y 523 HR 524 14.93 525 FS? 00 526 CLX 527 + 528 STO M 529 1 530 P-R 531 X^2 532 X<>Y 533 RCL 01 534 RCL 02 535 / 536 * 537 X^2 538 + 539 / 540 RCL 01 541 X^2 542 RCL 02 543 X^2 544 STO T 545 - 546 * 547 + 548 RCL N 549 1 550 P-R 551 X^2 552 X<>Y 553 X^2 554 RCL 03 555 X^2 556 R^ 557 STO O 558 / 559 * 560 + 561 / 562 RCL O 563 RCL 03 564 X^2 565 - 566 * 567 RCL 03 568 X^2 569 ST/ O 570 + 571 SQRT 572 RCL N 573 1 574 P-R 575 RCL O 576 * 577 R-P 578 X<> M 579 R^ 580 XEQ "S-R" 581 CLA 582 END |
( 835 bytes / SIZE 053 )
| STACK | INPUTS | OUTPUTS |
| T | l1 ( ° ' " ) | / |
| Z | b1( ° ' " ) | / |
| Y | l2 ( ° ' " ) | / |
| X | b2 ( ° ' " ) | D ( km ) |
Example1: Calculate
the geodesic distance between the U.S. Naval Observatory at Washington
(D.C.) and the "Observatoire de Paris"
l1 = 77°03'56"0
W , b1 = 38°55'17"2 N
l2 = 2°20'13"8
E , b2 =
48°50'11"2 N
CF 00
-77.0356 ENTER^
38.55172 ENTER^
2.20138 ENTER^
48.50112 XEQ "TGD2" >>>>
D = 6181.626264 km
---Execution time = 4mn04s---
-The error is about 26 centimeters.
Example2: Calculate the geodesic distance between the U.S. Naval Observatory at Washington (D.C.)
L1 = 77°03'56"0 W ; b1 = 38°55'17"2
N and Cape Town Observatory ( South Africa )
L2 = 18°28'35"7 E ; b2
= 33°56'02"5 S
-77.0356 ENTER^
38.55172 ENTER^
18.28357 ENTER^
-33.56025 XEQ "TGD2" >>>>
D = 12709.56616 km
-The exact value was s = 12709.56545 km, so, the error = 71 cm
Example3: With ( -40° , -40° ) ( +120° , +50° ) it yields 18095.48823 km
-Exact geodesic distance = 18095.48497 km -> Error = 3m26
-In this case, the error remains larger than 1 meter.
Notes:
-Gauss-Legendre 3-point formula is used to compute the integral.
-If you prefer the geodetic longitude, replace LBL 14 ( lines 518 to 581
) by the routine listed in paragraph 2°)a)
with small modifications since the semi-axes a , b , c are in
registers R01-R02-R03 here.
-If you have the cartesian coordinates x , y , z of the 2 points
M & N, store these values in R05 to R10
then divide the 3 registers R05-R06-R07 by the norm of the vector
OM and the 3 registers R08-R09-R10 by the norm of the vector ON
and press XEQ 10
-With the first example at the top of this page, you will get something like 8.594838705 ( exact result = 8.594822580 )
-The sum of the elliptic arcs lengths arc(MP0) + arc(P0Q0)
+ arc(Q0N) is in register R47 and in L-register at the end.
-For instance, in example 3, you get 18095.54550
and "TGD2" has diminished this value by about 57 meters.
-This program should not be used for nearly antipodal points... though it may work in some cases.
-The union of such 3 arcs is not really differentiable at 2 points, so
using a polynomial approximation remains doubtful.
-A better approach is given hereunder.
c) Program#3
-This version employs the same angles µ & h to determine the
position of a point P of the curve on the ellipsoid.
-But instead of adding several arc lengths in different planes, we search
an approximation of the function h(µ) under the form:
h(µ) = h1 Sin ( 180° µ / W ) + h2 Sin ( 360° µ / W ) ( in fact the first terms of a Fourier series )
-And the corresponding integral s = §0W [ R2 Cos2 h + R2 ( dh/dµ )2 + ( dR/dµ )2 ] 1/2 dµ
= ( Sin W ) §0W [ ( Cos2 h ) / T + (1/T) ( dh/dµ )2 + ( dr/dµ )2 ] 1/2 dµ is evaluated by a Gaussian quadrature.
where T = ( A2 + B2 + C2 ) with the same notations as in paragraph above:
r = 1/sqrt(T) dr/dµ = ¶r/¶µ + ( ¶r/¶h ) ( dh/dµ )
r = 1 / ( A2 + B2 + C2 )1/2 with
A = [ x1 Cos h Sin ( W - µ
) + x2 Cos h Sin µ + ( y1z2 - y2z1
) Sin h ] / a
B = [ y1 Cos h Sin ( W - µ
) + y2 Cos h Sin µ + ( z1x2 - z2x1
) Sin h ] / b
where u ( x1 , y1 , z1 )
& v ( x2 , y2 , z2 )
C = [ z1 Cos h Sin ( W - µ
) + z2 Cos h Sin µ + ( x1y2 - x2y1
) Sin h ] / c
the partial derivatives
¶r/¶µ
= ( -1/sqrt(T3) ) [ (A/a) { -x1 Cos h Cos ( W -
µ ) + x2 Cos h Cos µ }
+ (B/b) { -y1 Cos h Cos ( W - µ ) + y2 Cos
h Cos µ }
+ (C/c) { -z1 Cos h Cos ( W - µ ) + z2 Cos
h Cos µ } ]
and
¶r/¶h
= ( -1/sqrt(T3) ) [ (A/a) { -x1 Sin h Sin ( W -
µ ) - x2 Sin h Sin µ + ( y1z2
- y2z1 ) Cos h }
+ (B/b) { -y1 Sin h Sin ( W - µ ) - y2 Sin
h Sin µ + ( z1x2 - z2x1
) Cos h }
+ (C/c) { -z1 Sin h Sin ( W - µ ) - z2 Sin
h Sin µ + ( x1y2 - x2y1
) Cos h } ]
-With H = W / 1000 ( lines 89-90 ) , we calculate the following lengths
A = f(0,0)
B = f(H,0)
C = f(-H,0)
f = length s corresponding to h1 = -H , 0 , +H , h2
= -H , 0 , +H
D = f(0,H)
E = f(0,-H)
-Then, the extremum is estimated by the formula:
s ~ A - (B-C)2 / [ 8.(B+C-2.A) ] - (D-E)2
/ [ 8.(D+E-2.A) ]
Data Registers: R00 = n = number of subintervals of integration ( Registers R41 to R46 are to be initialized before executing "TGD3" )
R01 = a
R05 = x1
R08 = x2
R32 = W
R35 = f(-H,0)
R02 = b
R06 = y1
R09 = y2
R33 = f(0,0)
R36 = f(0,H)
R03 = c
R07 = z1
R10 = z2
R34 = f(H,0)
R37 = f(0,-H)
>>>> When the program stops, R04 = s
• R41 = 0.4679139346
• R44 = 0.6612093865
• R42 = 0.2386191861
• R45 = 0.1713244924
( coefficients for Gauss-Legendre 6-point formula )
• R43 = 0.3607615730
• R46 = 0.9324695142
Flag: F00
CF 00 for the Earth: n = 1 is sufficient and R01-R02-R03
are initialized ( lines 08 to 17 )
SF 00 for another ellipsoid: you'll have to initialize
R00-R01-R02-R03
Subroutine: "S-R" Spherical-Rectangular
conversion ( cf "Spherical Coordinates for the HP-41" paragraph 1°-a)
)
-Line 146 is a three-byte GTO 10 ( not really important )
-LBL 14 is actually the routine "LB-XYZ" listed in paragraph 2°)b)
- slightly modified.
-You can also add several lines to initialize R41 to R46 , for instance
after line 43.
| 01 LBL "TGD3" 02 STO 10 03 RDN 04 STO 09 05 RDN 06 FS? 00 07 GTO 00 08 6378.172 09 STO 01 10 .07 11 - 12 STO 02 13 21.35 14 - 15 STO 03 16 SIGN 17 STO 00 18 RDN 19 LBL 00 20 XEQ 14 21 STO 05 22 X<>Y 23 STO 06 24 R-P 25 RCL Z 26 STO 07 27 R-P 28 ST/ 05 29 ST/ 06 30 ST/ 07 31 RCL 09 32 RCL 10 33 XEQ 14 34 STO 08 35 X<>Y 36 STO 09 37 R-P 38 RCL Z 39 STO 10 40 R-P 41 ST/ 08 42 ST/ 09 43 ST/ 10 44 LBL 10 45 RCL 06 46 RCL 10 47 * 48 RCL 07 49 RCL 09 50 * 51 - 52 STO 21 53 RCL 07 54 RCL 08 55 * 56 RCL 05 57 RCL 10 58 * 59 - 60 STO 22 61 RCL 05 62 RCL 09 63 * 64 RCL 06 65 RCL 08 66 * 67 - 68 STO 23 69 46.04 70 STO 38 71 RCL 05 72 RCL 08 73 * |
74 RCL 06 75 RCL 09 76 * 77 + 78 RCL 07 79 RCL 10 80 * 81 + 82 ACOS 83 ENTER^ 84 STO 32 85 SIN 86 D-R 87 STO 31 88 CLX 89 E3 90 / 91 STO 11 92 CLX 93 STO 12 94 XEQ 12 95 STO 34 96 RCL 11 97 CHS 98 STO 11 99 XEQ 12 100 STO 35 101 CLX 102 X<> 11 103 STO 12 104 XEQ 12 105 STO 37 106 RCL 12 107 CHS 108 STO 12 109 XEQ 12 110 STO 36 111 CLX 112 STO 12 113 XEQ 12 114 STO 33 115 RCL 34 116 RCL 35 117 - 118 X^2 119 RCL 34 120 RCL 35 121 + 122 R^ 123 ST+ X 124 - 125 / 126 RCL 36 127 RCL 37 128 - 129 X^2 130 RCL 36 131 R^ 132 ST+ X 133 - 134 RCL 37 135 + 136 / 137 + 138 8 139 / 140 CHS 141 RCL 33 142 + 143 STO 04 144 RTN 145 STO 00 146 GTO 10 |
147 LBL 12 148 CLX 149 STO 04 150 STO 40 151 RCL 32 152 RCL 00 153 STO 30 154 ST+ X 155 / 156 STO 39 157 LBL 01 158 ST+ 40 159 RCL 38 160 STO 47 161 LBL 02 162 RCL 40 163 STO 48 164 RCL 39 165 RCL IND 47 166 * 167 ST+ 48 168 - 169 XEQ 13 170 X<> 48 171 XEQ 13 172 RCL 48 173 + 174 DSE 47 175 RCL IND 47 176 * 177 ST+ 04 178 DSE 47 179 GTO 02 180 RCL 39 181 ST+ 40 182 DSE 30 183 GTO 01 184 RCL 04 185 * 186 RCL 31 187 * 188 RTN 189 LBL 13 190 STO 13 191 1 192 P-R 193 STO 15 194 RDN 195 STO 14 196 RCL 08 197 * 198 RCL 32 199 RCL 13 200 - 201 1 202 P-R 203 STO 17 204 RDN 205 STO 16 206 RCL 05 207 * 208 + 209 STO 18 210 PI 211 R-D 212 RCL 13 213 * 214 RCL 32 215 / 216 STO Z 217 RCL 11 218 P-R 219 STO 24 |
220 X<> T 221 ST+ X 222 RCL 12 223 P-R 224 ST+ X 225 ST+ 24 226 RDN 227 + 228 1 229 P-R 230 STO 26 231 ST* Z 232 RDN 233 STO 25 234 RCL 21 235 * 236 + 237 RCL 01 238 / 239 STO 27 240 X^2 241 RCL 06 242 RCL 16 243 * 244 RCL 09 245 RCL 14 246 * 247 + 248 STO 19 249 RCL 26 250 * 251 RCL 22 252 RCL 25 253 * 254 + 255 RCL 02 256 / 257 STO 28 258 X^2 259 + 260 RCL 07 261 RCL 16 262 * 263 RCL 10 264 RCL 14 265 * 266 + 267 STO 20 268 RCL 26 269 * 270 RCL 23 271 RCL 25 272 * 273 + 274 RCL 03 275 / 276 STO 29 277 X^2 278 + 279 X<> 18 280 RCL 25 281 * 282 RCL 21 283 RCL 26 284 * 285 - 286 RCL 27 287 RCL 01 288 / 289 STO 27 290 * 291 RCL 19 292 RCL 25 |
293 * 294 RCL 22 295 RCL 26 296 * 297 - 298 RCL 28 299 RCL 02 300 / 301 STO 28 302 * 303 + 304 RCL 20 305 RCL 25 306 * 307 RCL 23 308 RCL 26 309 * 310 - 311 RCL 29 312 RCL 03 313 / 314 STO 29 315 * 316 + 317 RCL 24 318 PI 319 * 320 RCL 32 321 / 322 STO 14 323 * 324 STO 13 325 RCL 05 326 RCL 17 327 * 328 RCL 08 329 RCL 15 330 * 331 - 332 RCL 27 333 * 334 RCL 06 335 RCL 17 336 * 337 RCL 09 338 RCL 15 339 * 340 - 341 RCL 28 342 * 343 + 344 RCL 07 345 RCL 17 346 * 347 RCL 10 348 RCL 15 349 * 350 - 351 RCL 29 352 * 353 + 354 RCL 26 355 * 356 RCL 13 357 + 358 RCL 18 359 / 360 X^2 361 RCL 14 362 X^2 363 + 364 RCL 26 365 X^2 |
366 + 367 RCL 18 368 / 369 SQRT 370 RTN 371 LBL 14 372 DEG 373 HR 374 STO N 375 X<>Y 376 HR 377 14.93 378 FS? 00 379 CLX 380 + 381 STO M 382 1 383 P-R 384 X^2 385 X<>Y 386 RCL 01 387 RCL 02 388 / 389 * 390 X^2 391 + 392 / 393 RCL 01 394 X^2 395 RCL 02 396 X^2 397 STO T 398 - 399 * 400 + 401 RCL N 402 1 403 P-R 404 X^2 405 X<>Y 406 X^2 407 RCL 03 408 X^2 409 R^ 410 STO O 411 / 412 * 413 + 414 / 415 RCL O 416 RCL 03 417 X^2 418 - 419 * 420 RCL 03 421 X^2 422 ST/ O 423 + 424 SQRT 425 RCL N 426 1 427 P-R 428 RCL O 429 * 430 R-P 431 X<> M 432 R^ 433 XEQ "S-R" 434 CLA 435 END |
( 617 bytes / SIZE 049 )
| STACK | INPUTS | OUTPUTS |
| T | l1( ° ' " ) | / |
| Z | b1 ( ° ' " ) | / |
| Y | l2 ( ° ' " ) | / |
| X | b2 ( ° ' " ) | D ( km ) |
Example1: Calculate
the geodesic distance between the U.S. Naval Observatory at Washington
(D.C.) and the "Observatoire de Paris"
l1 = 77°03'56"0
W , b1 = 38°55'17"2 N
l2 = 2°20'13"8
E , b2 =
48°50'11"2 N
CF 00
-77.0356 ENTER^
38.55172 ENTER^
2.20138 ENTER^
48.50112 XEQ "TGD3" >>>>
D = 6181.626034 km
---Execution time = 5mn16s---
-The error is about 3 centimeters.
Example2: Calculate the geodesic distance between the U.S. Naval Observatory at Washington (D.C.)
L1 = 77°03'56"0 W ; b1 = 38°55'17"2
N and Cape Town Observatory ( South Africa )
L2 = 18°28'35"7 E ; b2
= 33°56'02"5 S
-77.0356 ENTER^
38.55172 ENTER^
18.28357 ENTER^
-33.56025 XEQ "TGD3" >>>>
D = 12709.56552 km
-The exact value was s = 12709.56545 km, error = 7 cm
Example3: With ( -40° , -40° ) ( +120° , +50° ) it yields 18095.48516 km
-Exact geodesic distance = 18095.48497 km -> Error = 19 cm
-Now the errors are always much smaller than 1 meter ( except for nearly
antipodal points ).
Notes:
-The elliptic arc length MN is in registers R33 and L
-Gauss-Legendre 6-point formula is used to compute the integral.
- n = 1 seems sufficient for the Earth.
-With n = 1, Gauss-Legendre 5-point formula would add an extra error
of a few decimeters for long geodesics.
-If you prefer the geodetic longitude, replace LBL 14 ( lines 371 to 434
) by the routine listed in paragraph 2°)a)
with small modifications since the semi-axes a , b , c are already
in registers R01-R02-R03.
-If you have the cartesian coordinates x , y , z of the 2 points
M & N, store these values in R05 to R10
then divide the 3 registers R05-R06-R07 by the norm of the vector
OM
and the 3 registers R08-R09-R10 by the norm of the vector ON,
then press XEQ 10
-You can also compute other arc lengths corresponding to different values
of h1 & h2
-Change lines 91 to 112 according to your choices and lines 114 to 143
to exploit these results.
Variants / Improvements:
• With the geocentric longitudes that are used here, LBL 14 may be simplified:
-Replace lines 404 to 432 by
RCL 03 X<> T
X<> M
X^2
*
1
ST* Z
R-P
and replace lines 21 to 43 with
STO 05 RCL 09
STO 09
RDN
RCL 10 X<>Y
STO 06 XEQ 14
STO 10
X<>Y
STO 08
STO 07 RDN
since the vectors are already unit vectors.
• If the inputs are geocentric longitudes and latitudes, the simplifications are even more important:
-Simply replace lines 371 to 435 by
LBL 14 14.93
XEQ "S-R"
DEG FS?
00 END
HR
CLX
X<>Y +
HR
1
d) Program#4
-In this variant, we search h(µ) under the form
h(µ) = h1 Sin ( 180° µ / W ) + h2 Sin ( 360° µ / W ) + h3 Sin ( 540° µ / W )
-Only 2 terms of the Fourier series if CF 01 but 3 terms if
SF 01.
-So, with SF 01 "TGD4" should return more accurate results.
Data Registers: R00 = n = number of subintervals of integration ( Registers R40 to R45 are to be initialized before executing "TGD4" )
R01 = a
R05 = x1
R08 = x2
R26 = W
R02 = b
R06 = y1
R09 = y2
R33 = f(0,0)
R35 to R39: temp
R03 = c
R07 = z1
R10 = z2
R34 = f(H,0)
>>>> When the program stops, R04 = s
• R40 = 0.4679139346
• R43 = 0.6612093865
• R41 = 0.2386191861
• R44 = 0.1713244924
( coefficients for Gauss-Legendre 6-point formula )
• R42 = 0.3607615730
• R45 = 0.9324695142
Flags: F00-F01-F03-F04
CF 00 for the Earth: n = 1 is sufficient if CF 01 and
R01-R02-R03 are initialized ( lines 08 to 17 )
SF 00 for another ellipsoid: you'll have to initialize
R00-R01-R02-R03
CF 01 only 2 terms in the Fourier series: same results
as "TGD3"
SF 01 3 terms are found
SF 04 geocentric longitudes and latitudes
CF 04 & CF 03 geocentric longitudes
& geodetic latitudes
CF 04 & SF 03 geodetic
longitudes and latitudes
Subroutine: "S-R" Spherical-Rectangular conversion ( cf "Spherical Coordinates for the HP-41" paragraph 1°-a) )
-Line 122 is a three-byte GTO 10 ( not really important )
| 01 LBL "TGD4" 02 STO 08 03 RDN 04 STO 07 05 RDN 06 FS? 00 07 GTO 00 08 6378.172 09 STO 01 10 .07 11 - 12 STO 02 13 21.35 14 - 15 STO 03 16 SIGN 17 STO 00 18 RDN 19 LBL 00 20 XEQ 14 21 STO 05 22 RDN 23 STO 06 24 X<>Y 25 X<> 07 26 RCL 08 27 XEQ 14 28 STO 08 29 RDN 30 STO 09 31 X<>Y 32 STO 10 33 LBL 10 34 RCL 06 35 RCL 10 36 * 37 RCL 07 38 RCL 09 39 * 40 - 41 STO 11 42 RCL 07 43 RCL 08 44 * 45 RCL 05 46 RCL 10 47 * 48 - 49 STO 12 50 RCL 05 51 RCL 09 52 * 53 RCL 06 54 RCL 08 55 * 56 - 57 STO 13 58 45.039 59 STO 36 60 RCL 05 61 RCL 08 62 * 63 RCL 06 64 RCL 09 65 * |
66 + 67 RCL 07 68 RCL 10 69 * 70 + 71 ACOS 72 STO 26 73 SIN 74 D-R 75 STO 27 76 CLX 77 STO 16 78 STO 17 79 STO 35 80 XEQ 12 81 STO 33 82 3 83 FC? 01 84 DSE X 85 STO 17 86 RCL 26 87 E3 88 / 89 STO 16 90 LBL 11 91 XEQ 12 92 STO 34 93 SIGN 94 CHS 95 ST* 16 96 XEQ 12 97 STO Y 98 RCL 34 99 ST+ Z 100 - 101 X^2 102 RCL 33 103 ENTER 104 + 105 R^ 106 - 107 X#0? 108 / 109 ST+ 35 110 DSE 17 111 GTO 11 112 RCL 35 113 8 114 ST/ Z 115 / 116 RCL 33 117 ST+ Z 118 + 119 STO 04 120 RTN 121 STO 00 122 GTO 10 123 LBL 12 124 CLX 125 STO 04 126 STO 28 127 RCL 26 128 RCL 00 129 STO 29 130 ST+ X |
131 / 132 STO 30 133 LBL 01 134 ST+ 28 135 RCL 36 136 STO 31 137 LBL 02 138 RCL 28 139 STO 32 140 RCL 30 141 RCL IND 31 142 * 143 ST+ 32 144 - 145 XEQ 13 146 X<> 32 147 XEQ 13 148 RCL 32 149 + 150 DSE 31 151 RCL IND 31 152 * 153 ST+ 04 154 DSE 31 155 GTO 02 156 RCL 30 157 ST+ 28 158 DSE 29 159 GTO 01 160 RCL 04 161 * 162 RCL 27 163 * 164 RTN 165 LBL 13 166 STO 18 167 1 168 P-R 169 STO 20 170 RDN 171 STO 19 172 RCL 08 173 * 174 RCL 18 175 RCL 26 176 ST/ 18 177 - 178 1 179 P-R 180 STO 22 181 RDN 182 STO 21 183 RCL 05 184 * 185 - 186 STO 23 187 PI 188 R-D 189 RCL 18 190 * 191 RCL 17 192 * 193 RCL 16 194 P-R 195 RCL 17 |
196 * 197 PI 198 * 199 RCL 26 200 / 201 STO 18 202 CLX 203 SIGN 204 P-R 205 STO 15 206 ST* Z 207 RDN 208 STO 14 209 ST* 23 210 RCL 11 211 * 212 + 213 RCL 01 214 / 215 STO 37 216 X^2 217 RCL 09 218 RCL 19 219 * 220 RCL 06 221 RCL 21 222 * 223 - 224 STO 24 225 RCL 15 226 * 227 RCL 12 228 RCL 14 229 ST* 24 230 * 231 + 232 RCL 02 233 / 234 STO 38 235 X^2 236 + 237 RCL 10 238 RCL 19 239 * 240 RCL 07 241 RCL 21 242 * 243 - 244 STO 25 245 RCL 15 246 * 247 RCL 13 248 RCL 14 249 ST* 25 250 * 251 + 252 RCL 03 253 / 254 STO 39 255 X^2 256 + 257 RCL 01 258 ST/ 37 259 RCL 02 260 ST/ 38 |
261 RCL 03 262 ST/ 39 263 R^ 264 X<> 23 265 RCL 11 266 RCL 15 267 * 268 - 269 RCL 37 270 * 271 RCL 24 272 RCL 12 273 RCL 15 274 * 275 - 276 RCL 38 277 * 278 + 279 RCL 25 280 RCL 13 281 RCL 15 282 * 283 - 284 RCL 39 285 * 286 + 287 RCL 18 288 * 289 STO 14 290 RCL 05 291 RCL 22 292 * 293 RCL 08 294 RCL 20 295 * 296 - 297 RCL 37 298 * 299 RCL 06 300 RCL 22 301 * 302 RCL 09 303 RCL 20 304 * 305 - 306 RCL 38 307 * 308 + 309 RCL 07 310 RCL 22 311 * 312 RCL 10 313 RCL 20 314 * 315 - 316 RCL 39 317 * 318 + 319 RCL 15 320 * 321 RCL 14 322 + 323 RCL 23 324 / 325 X^2 |
326 RCL 18 327 X^2 328 + 329 RCL 15 330 X^2 331 + 332 RCL 23 333 / 334 SQRT 335 RTN 336 LBL 14 337 DEG 338 HR 339 STO N 340 X<>Y 341 HR 342 14.93 343 FS? 00 344 CLX 345 + 346 FS? 04 347 GTO 09 348 FC? 03 349 GTO 08 350 1 351 P-R 352 RCL 01 353 RCL 02 354 / 355 X^2 356 * 357 R-P 358 X<>Y 359 LBL 08 360 STO M 361 1 362 P-R 363 X^2 364 X<>Y 365 RCL 01 366 RCL 02 367 / 368 X^2 369 * 370 X^2 371 + 372 SQRT 373 RCL 03 374 RCL 01 375 / 376 X^2 377 * 378 RCL N 379 1 380 P-R 381 RCL Z 382 / 383 R-P 384 X<> M 385 LBL 09 386 1 387 XEQ "S-R" 388 CLA 389 END |
( 554 bytes / SIZE 046 )
| STACK | INPUTS | OUTPUTS |
| T | l1 ( ° ' " ) | / |
| Z | b1 ( ° ' " ) | / |
| Y | l2 ( ° ' " ) | D1 ( km ) |
| X | b2 ( ° ' " ) | D2 or D3 ( km ) |
| L | / | D0 ( km ) |
CF 01 ---Execution
time = 4mn48s--- With n = 1
SF 01
---Execution time = 6mn41s---
---------
Examples 1-2-3:
-The results are almost identical to those given by "TGD3" with or without
F01 set
Example4: On the ellipsoid x2 / 41 + y2 / 37 + z2 / 35 = 1 the geocentric coordinates of the 2 points are:
l1 = 9°17'35"55660
l2 = 63°20'07"3683
b1 = -8°11'55"79344
b2 = 57°46'42"9935
41 SQRT STO 01
37 SQRT STO 02
35 SQRT STO 03
-If you choose n = 1 1 STO 00
SF 00 CF 01 SF 04
9.173555660 ENTER^
-8.115579344 ENTER^
63.20073683 ENTER^
57.46429935 XEQ "TGD4" >>>>
D = 8.594824283
-To calculate D with a different n-value, say 2, simply press
2 R/S >>>> 8.594824312
-If n = 4 , 4 R/S >>>>
8.594824330
-With SF 01 ( more accurate )
1 R/S >>>>
8.594823587
2 R/S >>>>
8.594823575
4 R/S >>>>
8.594823592
-In fact, this example is the 1st on this page and we have found the exact
value: about 8.594822580 whence an error of 1 E-6
Notes:
-We could replace line 82 by 4 - instead of 3 - to get 4 terms of the
Fourier series,
but the method to estimate the coefficients is not perfect and
could lead to disappointing results.
-If you replace line 82 by 2 ( instead of 3 ) and if you delete lines
83-84, "TGD4" becomes equivalent to "TGD3",
apart from the fact that "TGD4" is shorter and faster than "TGD3".
-Lines 348 to 359 are useful if the longitudes are geodetic ( set flag F03 ). We have:
Tan Lgeocentric = (b2/a2)
Tan Lgeodetic
Variant:
-Instead of using flag F01, you can initialize a register - say R46 - with
the number of terms
that you want to use in the Fourier series and replace lines 82
to 84 by RCL 46
e) For The Earth Only
-The integrand contains terms of the form: u Sin ( W-µ ) +
v Cos µ ....
-They may be re-written x Cos µ + x' Sin µ
-Moreover, the divisions by a , b , c can be performed at the beginning
of the program i-e before the loops
-So, we save bytes and time ( this can be applied to "TGD3" & "TGD4" too )
-Performing the angular functions in RAD mode also seems preferable.
-For the Earth, there is no need to subdivise the interval of integration
when Gauss-Legendre 6-point formula is used.
-Registers R00 thru R04 may also be re-used in the loops.
-All these simplifications are applied in "TGD5" hereunder.
Data Registers: R00 to R37: temp
( Registers R31 to R36 are to be initialized before executing "TGD5"
)
>>>> When the program stops, R01 = D
• R31 = 0.4679139346
• R34 = 0.6612093865
• R32 = 0.2386191861
• R35 = 0.1713244924
( coefficients for Gauss-Legendre 6-point formula )
• R33 = 0.3607615730
• R36 = 0.9324695142
Flags: F01-F03-F04
CF 01 2 terms in the Fourier series
SF 01 only 1 term is found
SF 04 geocentric longitudes and latitudes
CF 04 & CF 03 geocentric longitudes
& geodetic latitudes
CF 04 & SF 03 geodetic
longitudes and latitudes
Subroutine: "S-R" Spherical-Rectangular
conversion ( cf "Spherical Coordinates for the HP-41" paragraph 1°-a)
)
| 01 LBL "TGD5" 02 DEG 03 STO 08 04 RDN 05 STO 07 06 CLX 07 6378.172 08 STO 01 09 .07 10 - 11 STO 02 12 21.35 13 - 14 STO 03 15 CLX 16 14.93 17 STO 00 18 RDN 19 XEQ 05 20 STO 05 21 RDN 22 STO 06 23 X<>Y 24 X<> 07 25 RCL 08 26 XEQ 05 27 STO 08 28 RDN 29 STO 09 30 X<>Y 31 STO 10 32 RCL 06 33 * 34 RCL 07 35 RCL 09 36 * 37 - 38 STO 11 39 RCL 07 40 RCL 08 41 * 42 RCL 05 43 RCL 10 44 * 45 - 46 STO 12 47 RCL 05 48 RCL 09 49 * 50 RCL 06 51 RCL 08 52 * 53 - 54 STO 13 55 36.03 56 STO 37 57 RCL 05 58 RCL 08 59 * 60 RCL 06 61 RCL 09 62 * 63 + 64 RCL 07 |
65 RCL 10
66 * 67 + 68 STO 14 69 RAD 70 ACOS 71 STO 16 72 SIN 73 STO 15 74 LASTX 75 2 76 / 77 STO 17 78 RCL 05 79 RCL 14 80 * 81 ST- 08 82 RCL 01 83 ST/ 05 84 RCL 15 85 * 86 ST/ 08 87 ST/ 11 88 RCL 06 89 RCL 14 90 * 91 ST- 09 92 RCL 02 93 ST/ 06 94 RCL 15 95 * 96 ST/ 09 97 ST/ 12 98 RCL 07 99 RCL 14 100 * 101 ST- 10 102 RCL 03 103 ST/ 07 104 RCL 15 105 * 106 ST/ 10 107 ST/ 13 108 CLX 109 STO 18 110 STO 19 111 STO 20 112 XEQ 02 113 STO 21 114 2 115 FS? 01 116 SIGN 117 STO 18 118 RCL 16 119 E3 120 / 121 STO 19 122 LBL 01 123 XEQ 02 124 STO 22 125 SIGN 126 CHS 127 ST* 19 128 XEQ 02 |
129 STO Z 130 RCL 22 131 ST+ T 132 - 133 X^2 134 RCL 21 135 ST+ X 136 R^ 137 - 138 X#0? 139 / 140 ST+ 20 141 DSE 18 142 GTO 01 143 RCL 20 144 8 145 ST/ Z 146 / 147 RCL 21 148 ST+ Z 149 + 150 STO 01 151 DEG 152 RTN 153 LBL 02 154 CLX 155 STO 23 156 RCL 37 157 STO 24 158 LBL 03 159 RCL 17 160 ENTER^ 161 STO 25 162 RCL IND 24 163 * 164 ST- 25 165 + 166 XEQ 04 167 X<> 25 168 XEQ 04 169 RCL 25 170 + 171 DSE 24 172 RCL IND 24 173 * 174 ST+ 23 175 DSE 24 176 GTO 03 177 RCL 17 178 RCL 23 179 * 180 RTN 181 LBL 04 182 STO 14 183 RCL 18 184 RCL 16 185 / 186 STO 15 187 ST* 14 188 SIGN 189 P-R 190 STO 04 191 RCL 05 192 * |
193 X<>Y 194 STO 00 195 RCL 08 196 * 197 + 198 STO 01 199 PI 200 RCL 14 201 * 202 RCL 19 203 P-R 204 RCL 15 205 * 206 PI 207 * 208 STO 14 209 CLX 210 SIGN 211 P-R 212 STO 15 213 ST* Z 214 RDN 215 STO 26 216 ST* 01 217 RCL 11 218 * 219 + 220 STO 27 221 X^2 222 RCL 09 223 RCL 00 224 * 225 RCL 06 226 RCL 04 227 * 228 + 229 STO 02 230 RCL 15 231 * 232 RCL 12 233 RCL 26 234 ST* 02 235 * 236 + 237 STO 28 238 X^2 239 + 240 RCL 10 241 RCL 00 242 * 243 RCL 07 244 RCL 04 245 * 246 + 247 STO 03 248 RCL 15 249 * 250 RCL 13 251 RCL 26 252 ST* 03 253 * 254 + 255 STO 29 256 X^2 |
257 + 258 STO 30 259 RCL 11 260 RCL 15 261 * 262 RCL 01 263 - 264 RCL 27 265 * 266 RCL 12 267 RCL 15 268 * 269 RCL 02 270 - 271 RCL 28 272 * 273 + 274 RCL 13 275 RCL 15 276 * 277 RCL 03 278 - 279 RCL 29 280 * 281 + 282 STO 01 283 RCL 04 284 RCL 08 285 * 286 RCL 00 287 RCL 05 288 * 289 - 290 RCL 27 291 * 292 RCL 04 293 RCL 09 294 * 295 RCL 00 296 RCL 06 297 * 298 - 299 RCL 28 300 * 301 + 302 RCL 04 303 RCL 10 304 * 305 RCL 00 306 RCL 07 307 * 308 - 309 RCL 29 310 * 311 + 312 RCL 15 313 * 314 RCL 01 315 RCL 14 316 * 317 + 318 RCL 30 319 / 320 X^2 |
321 RCL 14
322 X^2 323 + 324 RCL 15 325 X^2 326 + 327 RCL 30 328 / 329 SQRT 330 RTN 331 LBL 05 332 HR 333 STO 10 334 X<>Y 335 HR 336 RCL 00 337 + 338 FS? 04 339 GTO 07 340 FC? 03 341 GTO 06 342 1 343 P-R 344 RCL 01 345 RCL 02 346 / 347 X^2 348 * 349 R-P 350 X<>Y 351 LBL 06 352 STO 09 353 1 354 P-R 355 X^2 356 X<>Y 357 RCL 01 358 RCL 02 359 / 360 X^2 361 * 362 X^2 363 + 364 SQRT 365 RCL 03 366 RCL 01 367 / 368 X^2 369 * 370 RCL 10 371 1 372 P-R 373 RCL Z 374 / 375 R-P 376 X<> 09 377 LBL 07 378 1 379 XEQ "S-R" 380 END |
( 512 bytes / SIZE 038 )
| STACK | INPUTS | OUTPUTS |
| T | l1 ( ° ' " ) | / |
| Z | b1 ( ° ' " ) | / |
| Y | l2 ( ° ' " ) | D1 ( km ) |
| X | b2 ( ° ' " ) | D2 or D1 ( km ) |
| L | / | D0 ( km ) |
CF 01 ---Execution
time = 3mn46s---
SF 01
---Execution time = 2mn19s---
Example1: Geodesic distance between the U.S. Naval Observatory at Washington (D.C.) and the "Observatoire de Paris"
l1 = 77°03'56"0
W , b1 = 38°55'17"2 N
l2 = 2°20'13"8
E , b2 =
48°50'11"2 N
CF 01 CF 03 CF 04
-77.0356 ENTER^
38.55172 ENTER^
2.20138 ENTER^
48.50112 XEQ "TGD5" >>>>
D2 = 6181.626035 km = R01
X<>Y D1 = 6181.626036 km
LASTX D0 = 6181.628623 km = R21
Example2: Geodesic distance between the U.S. Naval Observatory at Washington (D.C.)
l1 = 77°03'56"0
W , b1 = 38°55'17"2 N
and Cape Town Observatory ( South Africa )
l2 = 18°28'35"7
E , b2 = 33°56'02"5
S
CF 01 CF 03 CF 04
-77.0356 ENTER^
38.55172 ENTER^
18.28357 ENTER^
-33.56025 XEQ "TGD5" >>>>
D2 = 12709.56553 km = R01
X<>Y D1 = 12709.56694 km
LASTX D0 = 12709.56721 km = R21
Example3: With ( -40° , -40° ) ( +120° , +50° ) it yields:
• If CF 01 CF 03 CF 04 geocentric longitudes, geodetic latitudes
D2 = 18095.48517 km = R01
X<>Y D1 = 18095.49054
km
LASTX D0 = 18095.54547 km
= R21
• If CF 01 SF 03 CF 04 geodetic longitudes and latitudes
D2 = 18095.49441 km = R01
X<>Y D1 = 18095.49978
km
LASTX D0 = 18095.55469 km
= R21
• If CF 01 SF 04 geocentric longitudes and latitudes
D2 = 18099.69164 km = R01
X<>Y D1 = 18099.69700
km
LASTX D0 = 18099.75156 km
= R21
-If flag F01 is set, D2 is not computed ( faster but less accurate in general ).
Notes:
-The results should verify: D2 <= D1 <=
D0
-However, the formulae don't make the differences between
maximum and minimum: they just return an approximate extremum with
D1 & D2
-So, you could add X>0? CLX after line 139
-Lines 378-379 may be replaced by
X<>Y
RCL Z
1
X<>Y
P-R
P-R
thus avoiding the call to the subroutine "S-R" ... which may be replaced by an M-Code routine S-R !
-The longest geodesic distance is 20003.98600 km ( provided
a = 6378.172 km , b = 6378.102 km , c = 6356.752 km
)
-It is twice the geodesic distance between ( -14°55'48" , 0 ) and
( -14°55'48" , 90° )
-According to other sources, a = 6378.171 km
-Change lines 07 to 14 if you know better evaluations of the semi-axes...
-Despite the title of this paragraph, "TGD5" may also be used for other
celestial bodies, provided the flattenings remain small.
-Change lines 07 to 16 according to the characteristics of the triaxial
ellipsoid.
Variant:
-The program may of course be simplified if you have an M-Code routine
DOT ( dot product )
( cf for example "A few M-Code Routines for the HP-41" )
-For instance:
• Replace lines 283 to 311 by
RCL 10
*
RCL 06 RCL 04
-
RCL 04
-
RCL 00 *
DOT
*
RCL 09 *
RCL 05
RCL 07
RCL 04 -
ST* 00
RCL 00
*
RCL 08 X<> 00
• Replace lines 259 to 281 by
RCL 13
-
RCL 02 *
RCL 15
RCL 12 -
RCL 01
*
RCL 15 RCL 11
-
RCL 03
*
RCL 15 DOT
• Replace lines 220-237-255 with STO M STO N STO O respectively.
• Add CLA after line 151
-The first lines of the program may be simplified further if you have
also an M-Code routine CROSS ( cross product )
Remarks:
-The precision of these approximate methods seems quite acceptable, except for nearly antipodal points M & N.-But in this case, we can always add an intermediate point P and find the minimum
of MP + PN by a golden section search - though it will be very slow with a real HP-41...
-Of course, to get very accurate results, Jacobi's method remains the best !
-This new version employs Gauss-Legendre 10-point integration formula
-It uses the following semi-axes:
a = 6378.172 km with the equatorial major axis in the direction 14°92911W
b = 6378.102 km equatorial minor semi-axis
c = 6356.752314 km polar minor semi-axis ( like WGS84 )
-In this variant, we search h(µ) under the form
h(µ) = h1 Sin ( 180° µ / W ) + h2 Sin ( 360° µ / W ) + h3 Sin ( 540° µ / W ) + h4 Sin ( 720° µ / W )Data Registers: R00 to R49: temp
Flags: /
Subroutines: /
-With an HP41, round the coefficients in lines 59-61-....-77 to 10 decimals.
| 01 LBL "TGDX" 02 DEG 03 STO 08 04 RDN 05 STO 07 06 RDN 07 STO 06 08 X<>Y 09 STO 05 10 6378.172 11 STO 01 12 6378.102 13 STO 02 14 6356.752314 15 STO 03 16 14.92911 17 CHS 18 STO 00 19 RCL 05 20 RCL 06 21 XEQ 05 22 STO 05 23 RDN 24 STO 06 25 X<>Y 26 X<> 07 27 RCL 08 28 XEQ 05 29 STO 08 30 RDN 31 STO 09 32 X<>Y 33 STO 10 34 RCL 06 35 * 36 RCL 07 37 RCL 09 38 * 39 - 40 STO 11 41 RCL 07 42 RCL 08 43 * 44 RCL 05 45 RCL 10 46 * 47 - 48 STO 12 49 RCL 05 50 RCL 09 51 * 52 RCL 06 53 RCL 08 54 * 55 - 56 STO 13 57 40.03 58 STO 49 |
59 0.295524224715 60 STO 31 61 0.148874338982 62 STO 32 63 0.26926671931 64 STO 33 65 0.433395394129 66 STO 34 67 0.219086362516 68 STO 35 69 0.679409568299 70 STO 36 71 0.149451349151 72 STO 37 73 0.865063366689 74 STO 38 75 0.0666713443087 76 STO 39 77 0.973906528517 78 STO 40 79 RCL 05 80 RCL 08 81 * 82 RCL 06 83 RCL 09 84 * 85 + 86 RCL 07 87 RCL 10 88 * 89 + 90 STO 14 91 RAD 92 ACOS 93 STO 16 94 SIN 95 STO 15 96 LASTX 97 2 98 / 99 STO 17 100 RCL 05 101 RCL 14 102 * 103 ST- 08 104 RCL 01 105 ST/ 05 106 RCL 15 107 * 108 ST/ 08 109 ST/ 11 110 RCL 06 111 RCL 14 112 * 113 ST- 09 114 RCL 02 115 ST/ 06 116 RCL 15 |
117 * 118 ST/ 09 119 ST/ 12 120 RCL 07 121 RCL 14 122 * 123 ST- 10 124 RCL 03 125 ST/ 07 126 RCL 15 127 * 128 ST/ 10 129 ST/ 13 130 CLX 131 STO 18 132 STO 19 133 XEQ 02 134 STO 21 135 4 136 STO 18 137 44 138 STO 20 139 RCL 16 140 1 E3 141 / 142 STO 19 143 GTO 01 144 LBL 02 145 CLX 146 STO 23 147 RCL 49 148 STO 24 149 LBL 03 150 RCL 17 151 ENTER 152 STO 25 153 RCL IND 24 154 * 155 ST- 25 156 + 157 XEQ 04 158 X<> 25 159 XEQ 04 160 RCL 25 161 + 162 DSE 24 163 RCL IND 24 164 * 165 ST+ 23 166 DSE 24 167 GTO 03 168 RCL 17 169 RCL 23 170 * 171 RTN 172 LBL 04 173 STO 14 174 RCL 18 |
175 RCL 16 176 / 177 STO 15 178 ST* 14 179 SIGN 180 P-R 181 STO 04 182 RCL 05 183 * 184 X<>Y 185 STO 48 186 RCL 08 187 * 188 + 189 STO 45 190 PI 191 RCL 14 192 * 193 RCL 19 194 P-R 195 RCL 15 196 * 197 PI 198 * 199 STO 14 200 CLX 201 SIGN 202 P-R 203 STO 15 204 ST* Z 205 RDN 206 STO 26 207 ST* 45 208 RCL 11 209 * 210 + 211 STO 27 212 X^2 213 RCL 09 214 RCL 48 215 * 216 RCL 06 217 RCL 04 218 * 219 + 220 STO 46 221 RCL 15 222 * 223 RCL 12 224 RCL 26 225 ST* 46 226 * 227 + 228 STO 28 229 X^2 230 + 231 RCL 10 232 RCL 48 |
233 * 234 RCL 07 235 RCL 04 236 * 237 + 238 STO 47 239 RCL 15 240 * 241 RCL 13 242 RCL 26 243 ST* 47 244 * 245 + 246 STO 29 247 X^2 248 + 249 STO 30 250 RCL 11 251 RCL 15 252 * 253 RCL 45 254 - 255 RCL 27 256 * 257 RCL 12 258 RCL 15 259 * 260 RCL 46 261 - 262 RCL 28 263 * 264 + 265 RCL 13 266 RCL 15 267 * 268 RCL 47 269 - 270 RCL 29 271 * 272 + 273 STO 45 274 RCL 04 275 RCL 08 276 * 277 RCL 05 278 RCL 48 279 * 280 - 281 RCL 27 282 * 283 RCL 04 284 RCL 09 285 * 286 RCL 06 287 RCL 48 288 * 289 - 290 RCL 28 291 * |
292 + 293 RCL 04 294 RCL 10 295 * 296 RCL 07 297 RCL 48 298 * 299 - 300 RCL 29 301 * 302 + 303 RCL 15 304 * 305 RCL 14 306 RCL 45 307 * 308 + 309 RCL 30 310 / 311 X^2 312 RCL 14 313 X^2 314 + 315 RCL 15 316 X^2 317 + 318 RCL 30 319 / 320 SQRT 321 RTN 322 LBL 05 323 HR 324 X<>Y 325 HR 326 RCL 00 327 - 328 1 329 X<>Y 330 RDN 331 P-R 332 R^ 333 X<>Y 334 P-R 335 RCL 01 336 X^2 337 * 338 X^2 339 STO 35 340 X<> L 341 X<>Y 342 RCL 02 343 X^2 344 * 345 X^2 346 ST+ 35 347 X<> L 348 R^ 349 RCL 03 350 X^2 |
351 * 352 X^2 353 ST+ 35 354 X<> L 355 X<> Z 356 RCL 35 357 SQRT 358 ST/ T 359 ST/ Z 360 / 361 RTN 362 LBL 07 363 1 364 X<>Y 365 RDN 366 P-R 367 R^ 368 X<>Y 369 P-R 370 RTN 371 LBL 01 372 XEQ 02 373 STO 22 374 SIGN 375 CHS 376 ST* 19 377 XEQ 02 378 STO Z 379 RCL 22 380 ST+ T 381 - 382 X^2 383 RCL 21 384 ST+ X 385 R^ 386 - 387 8 388 * 389 X#0? 390 / 391 STO IND 20 392 DSE 20 393 DSE 18 394 GTO 01 395 DEG 396 RCL 21 397 ST+ 41 398 SIGN 399 RCL 41 400 ST+ 42 401 RCL 42 402 ST+ 43 403 RCL 43 404 ST+ 44 405 RCL 41 406 RCL 42 407 RCL 43 408 RCL 44 409 END |
( SIZE 050 )
| STACK | INPUTS | OUTPUT |
| T | L ( ° ' " ) | D1 |
| Z | b ( ° ' " ) | D2 |
| Y | L' ( ° ' " ) | D3 |
| X | b' ( ° ' " ) | D |
| L |
/ |
D0 |
where L , b , L' , b' are geodetic coordinates
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 "TGDX" >>>> D = 12138.657569 km
RDN 12138.657625
RDN 12138.657625
RDN 12138.658752
LASTX 12138.658756 = great elliptic arc length
-The exact result ( computed with Jacobi's method ) is D = 12138.657551942....
Notes:
-If you want to use this program with geocentric coordinates expressed in decimal degrees measured from the equatorial major axis
replace lines 21 & 28 ( XEQ 05 ) by XEQ 07 and line 16 by 0
-I've used free42 to compute these values.
References:
[1] Jacobi - "De la ligne
géodésique sur un ellipsoïde, et des différents
usages d’une transformation analytique remarquable"
[2] J. Feltens - Vector method to compute the Cartesian (X, Y,
Z) to geodetic (f, l,
h) transformation on a triaxial ellipsoid - J Geod (2009) 83:129–137
[3] Marcin Ligas - Cartesian to geodetic coordinates conversion
on a triaxial ellipsoid - J Geod (2012) 86:249–256
-With many thanks to Charles Karney for the link he sent me to get Jacobi's method:
http://lists.maptools.org/pipermail/proj/2012-October/006448.html
[4] If you want to visualize the geodesics on a triaxial ellipsoid,
go to this excellent
page