Terrestrial Geodesic Distance ( 3-axial Ellipsoid
)
"TGD" computes the geodesic distance between 2 points ( L , b ) (
L' , b' )
where L , L' are the geodetic longitudes & b , b'
are the geodetic latitudes.
This program employs the Andoyer's formula of order 3 to evaluate
the geodesic distance on the ellipsoid of revolution WGS84
and the difference between the great elliptic arc distances on WGS84
& the triaxial ellipsoid to estimate the geodesic distance on the triaxial
ellipsoid.
The semi-axes are:
a = 6378.172 km
b = 6378.102 km
and the longitude of the major axis = L0 = -14°92911
c = 6356.75231425 km
-Of course, you can change these values in variables 'A' 'B' 'C'
'L0'
-> Longitudes are positive East.
Andoyer's Formula of Order 3:
With L , L' = longitudes &
B , B' = parametric latitudes , f = (a-b) / a
Tan B = (1-f) Tan b where b = geodetic latitude
H = ( L' - L )/2 ; F = ( B + B' )/2 ;
G = ( B' - B )/2
S = sin2 G cos2 H + cos2
F sin2 H , µ = 2 Arc Sin S1/2
we define A = ( sin2 G cos2 F )/S
+ ( sin2 F cos2 G )/( 1-S )
and
B = ( sin2 G cos2 F )/S - ( sin2
F cos2 G )/( 1-S )
-After "some" calculations, it yields:
DistWGS84 = a { µ - (f/2) ( A.µ +
B Sin µ )
+ (f2/4) [ A µ2 Cot µ + B µ2
/ Sin µ + A2 ( µ + Sin µ Cos µ
- 4 µ2 Cot µ ) / 4 - A.B µ2
/ Sin µ - ( B2/2 ) Sin µ Cos µ ]
+ (f3/8) [ A ( 2 µ2 Cot µ -
2 µ3 Cot2 µ - 8 µ3
/ 6 ) + A2 ( -5 µ2 Cot µ + 5 µ3
Cot2 µ + 8 µ3 / 3 + (1/2) Sin µ
Cos µ + µ /2 )
+ A3 ( 3 µ2 Cot µ - 3 µ3
Cot2 µ - 8 µ3 / 6 + 2 Sin µ
Cos3 µ - (9/4) Sin µ Cos µ - µ
/ 4 + 2 Sin3 µ Cos3 µ + 2 Sin5
µ Cos µ )
+ B ( 2 µ2 / Sin µ - 2 µ3
Cos µ / Sin2 µ ) + A.B ( -5 µ2
/ Sin µ + 6 µ3 Cos µ / Sin2
µ + µ Cos µ )
+ A2 B ( 3 µ2 / Sin µ - 4 µ3
Cos µ / Sin2 µ - 2 Sin3
µ Cos2 µ - µ Cos µ - 2 Sin5
µ + 2 Sin µ - (3/2) Sin µ Cos2
µ )
+ B2 ( - Sin µ Cos µ + µ3
/ Sin2 µ + µ ) + A B2 (
(1/2) Sin µ Cos µ - µ - µ3 /
Sin2 µ )
+ B3 ( (2/3) Sin3 µ - (1/2) Sin µ
) ] }
-Cf http://hp41programs.yolasite.com/geodesic2.php
for more details about this method.
| STACK |
INPUTS |
OUTPUT |
| Level 4 |
L ( ° ' " ) |
/ |
| Level 3 |
b ( ° ' "
) |
/ |
| Level 2 |
L' ( °
' " ) |
/ |
| Level 1 |
b' ( ° ' "
) |
Distance ( km ) |
Example: Calculate the geodesic distance
between the Sydney Observatory ( Long = 151°12'17"8 E , Lat = 33°51'41"1
S )
and the Palomar Observatory ( Long = 116°51'50"4 W ,
Lat = 33°21'22"4 N )
151.12178 ENTER
-33.51411 ENTER
-116.51504 ENTER
33.21224
TGD >>>> D = 12138.6575674
km
---Execution time = 2.8s
with an HP48G/GX or 4.3s with an HP48S/SX---
Notes:
-If you key in LASTARG
just after executing "TGD" , Level 2 = great elliptic arc distance
= 12138.6587539 in this example
-It is calculated by Andoyer's formula of order 2
-The root-mean-square error is about 2 cm if 80° < | L'
- L | < 120° and the maximum error is about 13 cm
-The root-mean-square error is about 4 cm if 110°
< | L' - L | < 150° and the maximum error is about
34 cm
-When the difference in the longitudes is 165°,
the error can reach about 75 cm.
-With L = 0° , b = 2° & L' = 175° , b'
= -4° error is about 2 meters.
-This program does not give accurate results for nearly antipodal points.
-The method gives a good precision because the equator flattening is very
small.
-Otherwise, it would not work.
-"TGD" calls R->P & P->R to perform
rectangular <> polar conversions
-With the HP48 version, I've used externals that check if levels 1&2
are 2 real numbers before executing the conversions.
-I don't know the similar externals for the HP50, so the HP50 version simply
employs RPL programs.
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
-Many thanks to Professor Richard
H. Rapp who sent me reference [3] to find the typo in Sodano's
formula !