Geodesic2

# Great Elliptic Arc Length for the HP-41

Overview

1°)  Ellipsoid of Revolution
2°)  Triaxial Ellipsoid

a)  With Andoyer's Formula
b)  Numerical Integration ( 2 programs )

1°)  Ellipsoid of Revolution

-This program employs WGS84 ellipsoid:  a = 6378.137 km  &  flattening  f = 1/298.25722  ( lines 166 & 18 )

-The latitudes are geodetic.

"GEL" employs a formula of order 3 given in reference 1 ( page 69 )

Formula:

With      L , L' = longitudes  &  p , p' = parametric latitudes  ,  f = earth flattening  ,  e2  = f ( 2 - f )

Tan p = (1-f) Tan l  where  l = geodetic latitude

H = ( L' - L )/2  ;   F = ( p + p' )/2  ;  G = ( p' - p )/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 )

Dist / a  =  µ - (e2/4) ( A µ + B Sin µ ) - (e4/128) [ A2 ( 6 µ - Sin 2µ ) + 8 A B Sin µ + 2 B2 Sin 2µ ]
- (e6/1536) [ 3 A3 ( 10 µ - 3 Sin 2µ ) + 3 A2 B ( 15 Sin µ - Sin 3µ ) + 18 A B2 Sin 2µ  + 4 B3 Sin 3µ ]

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

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

( 206 bytes / SIZE 008 )

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

Example:
Calculate the geodesic distance between the Sydney Observatory ( Long = 151°12'17"8 E , Lat = 33°51'41"1 S )
and the Palomar Observatory ( Long = 116°51'50"4 W ,  Lat = 33°21'22"4 N )

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

2°)  Triaxial Ellipsoid

a)  With Andoyer's Formula

-The latitudes L and the longitudes B are geodetic.

-The semi-axes are:

a = 6378.172 km
b = 6378.102 km                 and the longitude of the equatorial major-axis = 14°929 West
c = 6356.752314 km

-So, the mean radius of the equator and the flattening correspond to WGS84.

-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.

-A point P on the arc MN may be determined by the vectorial equality:

With  OM / OM = ( x , y , z )   & ON / ON = ( x' , y' ,z' )      ( the unit vectors defined by OM & ON )     we have

OP = ( Sin W ) / ( A2 + B2 + C2 )1/2    for a point  P  on the ellipsoid  and   with  µ = ( OM , OP )

A = [ x Sin ( W - µ ) + x' Sin µ ] / a
B = [ y Sin ( W - µ ) + y' Sin µ ] / b           W = the angle ( OM , ON )
C = [ z Sin ( W - µ ) + z' Sin µ ] / c

N  P
/   /
/  /
/ / µ
O//-----------------------------M

>>> The value of µ for the axes of the great ellipse (E) is found by equating to zero the derivative of OP with respect to µ

-After some calculations, it yields:    ( Sin 2µ ) / ( Cos 2µ ) = P / Q  with

P  =  -2 [ ( x Sin W ) ( x Cos W - x' ) / a2 + ( y Sin W ) ( y Cos W - y' ) / b2 + ( z Sin W ) ( z Cos W - z' ) / c2 ]
Q =  [ x2 Sin2 W - ( x Cos W - x' )2 ] / a2 + [ y2 Sin2 W - ( y Cos W - y' )2 ] / b2 + [ z2 Sin2 W - ( z Cos W - z' )2 ] / c2

>>> However, it's simpler & faster to compute

A2 + B2 + C2 = T Sin2 ( W - µ ) + S Sin2 µ + 2 U Sin µ Sin ( W - µ )

where  S = x'2 / a2 + y'2 / b2 + z'2 / c2  ,   T =  x2 / a2 + y2 / b2 + z2 / c2  ,   U =  xx' / a2 + yy' / b2 + zz' / c2

and   P = 2 U - 2 T Cos W  ,  Q = T Sin W + ( 2 U Cos W - S - T Cos2 W ) / Sin W

-Then, the distance on the great ellipse of the triaxial ellipsoid - between M & N is computed by Andoyer's formula of order 2, re-written with geocentric "latitudes":

D/a' = W + (f/2) (-W+(Cos2F)(SinW)) +(f2/4) [W/4+(SinW)(CosW)/4-(Cos2(2F))(SinW)(CosW)/2]

a' = major axis of the great ellipse , f = eccentricity
W = angular separation
2F = b1+b2 where b = geocentric "latitudes"

Data Registers:     R01 = a  R02 = b  R03 = c                                                                                     R00 to R14:  temp

R06 & R07 ( or R07 & R06 ) = major & minor semi-axes of the great ellipse.
Flags: /
Subroutines: /

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

( 289 bytes / SIZE 015 )

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

Example:
Calculate the geodesic distance between the Sydney Observatory ( Long = 151°12'17"8 E , Lat = 33°51'41"1 S )
and the Palomar Observatory ( Long = 116°51'50"4 W ,  Lat = 33°21'22"4 N )

151.12178  ENTER^
-33.51411   ENTER^
-116.51504  ENTER^
33.21224   XEQ "GEL"   >>>>   D = 12138.68576 km                                                      ---Execution time = 17s---

b)  Numerical Integration

-The arc length:    §0W [ R2 + ( dR/dµ )2 ] 1/2 dµ    is calculated by Newton-Cotes formula of order 9

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

-Line 48 is a three-byte GTO 03

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

( 349 bytes / SIZE 023 )

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

Example:
Calculate the geodesic distance between the Sydney Observatory ( Long = 151°12'17"8 E , Lat = 33°51'41"1 S )
and the Palomar Observatory ( Long = 116°51'50"4 W ,  Lat = 33°21'22"4 N )

151.12178  ENTER^
-33.51411   ENTER^
-116.51504  ENTER^
33.21224   XEQ "GEL"   >>>>   D = 12138.65876 km                                                      ---Execution time = 56s---

Notes:

W is computed from  Cos W = x x' + y y' + z z'  so the precision may be small if W is very small

-With the formula:   W = 2 Arc Sin ( || u - v || / 2 )  the accuracy is better with small W-values.
-But the precision will be even better if we use

|| u x v || = Sin W  &  u.v  = Cos W

-Gauss-Legendre 10-point formula is employed.
-The interval  [0,W]  is divided in n subintervals

Data Registers:         •  R00 = n                                          ( Register R00 is to be initialized before executing "GEL" )

R05 = a  R06 = b  R07 = c              R01 to R34:  temp
Flags: /
Subroutines: /

-Line 76 is a three-byte GTO 03

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

( 436 bytes / SIZE 035 )

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

•   With  n = 1  STO 00

151.12178  ENTER^
-33.51411   ENTER^
-116.51504  ENTER^
33.21224   XEQ "GEL"   >>>>   D = 12138.65876 km                                                      ---Execution time = 63s---

-Same result with  n = 2

Notes:

-For the Earth, n = 1 is enough for a good precision.
-If the ellipsoid has larger flattenings, test this program with different n-values to obtain the best accuracy.
( modify lines 07 & 17 to 24 to store the different axis-values )

Reference:

[1]  Paul D. Thomas - "MATHEMATICAL MODELS FOR NAVIGATION SYSTEMS" - US Naval Oceanographic Office.