Gravity on a Triaxial Ellipsoid for the HP-41
Overview
1°) Calculating ga
, gb , gc
2°) Generalized Somigliana's Formula
-Somigliana's formula works on an ellipsoid of revolution.
-In reference [1], M. Caputo gives a generalization of this formula that
also works for triaxial ellipsoids.
x2/a2 + y2/b2 + z2/c2 = 1
-The unique assumption is that the gravity vector must always be normal
to the ellipsoid:
-In other words, the ellipsoid is an equipotential surface.
-But to apply this formula, we first have to calculate the gravity at the
points of intersection of the axes and the ellipsoid itself
1°) Calculating ga , gb
, gc
-The formulas are given in reference [1] ( warning: M. Caputo denotes a1 = b , a2 = a , a3 = c ). With small differences in the notations:
Let E2 = (b2 - c2) / c2 and n = (a2 - b2) / b2
Ai,j = A'i,j + n A"i,j
A'1,1 = A'1,2 = A'2,1 =
A'2,2 = [ 0.75 / (b2 - c2)5/2
] [ Arc tan E - (E/3) (5.E2 + 3) / (1+E2)2
]
A"1,1 = [ 5 / 16 / (b2 - c2)7/2
] [ -Arc tan E + ( E/15/(1+E2) ) (20 - (5-13.E4)/(1+E2)2
) ] A"2,1 = A"1,2
= 3 A"1,1 ; A"2,2 = 5 A"1,1
A'1,3 = A'2,3 = [ 3 / (b2
- c2)5/2 ] [ -Arc tan E + ( E/3/(1+E2)
) (2E2+3) ]
A"1,3 = [ 15 / 8 / (b2 - c2)7/2
] [ Arc tan E - (E/30) ( 25 + (5-9.E4) / (1+E2)2
) ]
A"2,3 = 3 A"1,1
-Then, we must calculate k1 & k2 with
D.K1 / w2
= - A'1,1 - n [ A'1,1 + 6 A"1,1 + (c2/b2)
A'2,3 ]
w2 = angular velocity of the Earth
= 7.292115 E-5 rd/s
D.K2 / w2
= - A'1,1 - n [ A'1,1 - (c2/b2)
A'2,3 ]
where D = 4 A'1,1 [ 2 A'1,1 - (c2/b2) A'2,3 ] - 2 n (c2/b2) A'2,3 ( A'1,1 + 6 A"1,1 ) + 4 n A'1,1 [ 2 A'1,1 + 12 A"1,1 - (c2/b2) A"2,3 ]
-Finally, with GM = Earth gravitational constant = 3.9860044188 E14 m3/s2
• ga/a = ( GM + 4 K2
/ a2 ) / ( abc ) - 2 ( A1,2
K1 + 3 A2,2 K2 ) - w2
• gb/b = ( GM + 4 K1
/ b2 ) / ( abc ) - 2 ( 3 A1,1
K1 + A2,1 K2 ) - w2
• gc/c =
GM / ( abc ) - 2 ( A1,3 K1 + A2,3 K2
)
-Unfortunately, the expressions involving Arc Tan E would produce
low accuracy because of cancellation of leading digits.
-In the following program, I've used truncated series to get a better precision,
namely:
Arc tan E - (E/3) (5.E2 + 3) / (1+E2)2
= (8/15) E5 - (8/7) E7 + (16/9) E9
- (80/33) E11 + (40/13) E13 - (56/15) E15
+ ...
-Arc tan E + ( E/15/(1+E2) ) (20 - (5-13.E4)/(1+E2)2
) = -(16/35) E7 + (64/45) E9 - (32/11) E11
+ (64/13) E13 - (112/15) E15 + (896/85) E17
- ....
-Arc tan E + ( E/3/(1+E2) ) (2E2+3) = (2/15)
E5 - (4/21) E7 + (2/9) E9 - (8/33) E11
+ (10/39) E13 - (4/15) E15 + ...
Arc tan E - (E/30) ( 25 + (5-9.E4) / (1+E2)2
) = -(8/105) E7 + (8/45) E9 - (16/55) E11
+ (16/39) E13 - (8/15) E15 + (56/85) E17
- ....
-For the Earth, E2 = (b2 - c2) /
c2 ~ 0.0067 so these terms are sufficient. More
terms should be used for larger E-values.
Data Registers: R00 = n = (a2 - b2) / b2
R01 = a R04 =
ga R07 thru R14:
temp
R02 = b R05 =
gb
R03 = c R06 =
gc
Flags: /
Subroutines: /
01 LBL "GABC" 02 STO 03 03 RDN 04 STO 02 05 X<>Y 06 ENTER^ 07 STO 01 08 RCL 02 09 ST+ Z 10 - 11 * 12 RCL 02 13 X^2 14 / 15 STO 00 16 RCL 02 17 RCL 02 18 RCL 03 19 ST+ Z 20 - 21 * 22 STO 07 23 RCL 03 24 X^2 25 / 26 ENTER^ 27 ENTER^ 28 STO 04 29 14 30 * 31 CHS 32 15 33 / 34 1.3 35 1/X 36 + 37 * 38 1.65 39 1/X 40 - 41 * 42 4 43 9 44 / 45 + 46 * 47 3.5 48 1/X 49 - 50 * 51 7.5 52 1/X 53 + 54 * 55 * 56 X<>Y 57 SQRT 58 STO 05 59 * 60 3 61 * 62 RCL 07 63 2.5 64 Y^X 65 STO 11 66 / 67 STO 09 68 CHS |
69 STO 13 70 STO 14 71 CLX 72 56 73 * 74 85 75 / 76 7 77 15 78 / 79 - 80 * 81 3.25 82 1/X 83 + 84 * 85 5.5 86 1/X 87 - 88 * 89 4 90 45 91 / 92 + 93 * 94 35 95 1/X 96 - 97 * 98 * 99 * 100 RCL 05 101 * 102 5 103 * 104 RCL 02 105 X^2 106 * 107 RCL 07 108 RCL 11 109 * 110 STO 12 111 / 112 STO 10 113 CLX 114 3.75 115 / 116 CHS 117 3.9 118 1/X 119 + 120 * 121 8 122 33 123 / 124 - 125 * 126 4.5 127 1/X 128 + 129 * 130 5.25 131 1/X 132 - 133 * 134 7.5 135 1/X 136 + |
137 * 138 * 139 RCL 05 140 * 141 3 142 * 143 RCL 11 144 / 145 STO 11 146 CLX 147 7 148 * 149 85 150 / 151 15 152 1/X 153 - 154 * 155 19.5 156 1/X 157 + 158 * 159 27.5 160 1/X 161 - 162 * 163 45 164 1/X 165 + 166 * 167 105 168 1/X 169 - 170 * 171 * 172 * 173 RCL 05 174 * 175 15 176 * 177 RCL 02 178 X^2 179 * 180 RCL 12 181 / 182 STO 12 183 RCL 10 184 6 185 * 186 RCL 09 187 + 188 STO 08 189 RCL 11 190 RCL 03 191 RCL 02 192 / 193 X^2 194 STO 07 195 * 196 2 197 / 198 + 199 RCL 00 200 * 201 ST- 13 202 RCL 09 203 RCL 11 204 RCL 07 |
205 * 206 2 207 / 208 - 209 RCL 00 210 * 211 ST- 14 212 RCL 08 213 ST+ X 214 RCL 12 215 3 216 * 217 RCL 07 218 * 219 - 220 RCL 09 221 * 222 ST+ X 223 RCL 08 224 RCL 11 225 * 226 RCL 07 227 * 228 - 229 RCL 00 230 ST* 10 231 ST* 12 232 * 233 RCL 09 234 ST+ X 235 RCL 11 236 RCL 07 237 * 238 - 239 RCL 09 240 * 241 ST+ X 242 + 243 ST+ X 244 7.292115 E-5 245 X^2 246 STO 08 247 / 248 ST/ 13 249 ST/ 14 250 3.986004419 E14 251 STO 06 252 RCL 13 253 4 254 * 255 RCL 02 256 X^2 257 / 258 + 259 STO 05 260 RCL 06 261 RCL 14 262 4 263 * 264 RCL 01 265 X^2 266 / 267 + 268 RCL 01 269 RCL 02 270 RCL 03 271 * 272 * |
273 ST/ 06 274 ST/ 05 275 / 276 STO 04 277 RCL 09 278 RCL 10 279 + 280 3 281 * 282 RCL 13 283 * 284 RCL 10 285 3 286 * 287 RCL 09 288 + 289 STO 07 290 RCL 14 291 * 292 + 293 ST+ X 294 RCL 08 295 + 296 ST- 05 297 RCL 07 298 RCL 13 299 * 300 RCL 10 301 5 302 * 303 RCL 09 304 + 305 3 306 * 307 RCL 14 308 * 309 + 310 ST+ X 311 RCL 08 312 + 313 ST- 04 314 RCL 11 315 RCL 12 316 + 317 RCL 13 318 * 319 RCL 12 320 3 321 * 322 RCL 11 323 + 324 RCL 14 325 * 326 + 327 ST+ X 328 ST- 06 329 RCL 01 330 ST* 04 331 RCL 02 332 ST* 05 333 RCL 03 334 ST* 06 335 RCL 06 336 RCL 05 337 RCL 04 338 END |
( 442 bytes / SIZE 015 )
STACK | INPUTS | OUTPUTS |
Z | a | gc |
Y | b | gb |
X | c | ga |
where a >= b >= c are expressed in meters & g in m/s2
Example:
a = 6378171.645 m
b = 6378101.575 m
c = 6356751.868 m
6378171.645 ENTER^
6378101.575 ENTER^
6356751.868 XEQ "GABC" >>>>
ga = 9.780379978 m/s2
---Execution time = 18s---
RDN gb = 9.780273552 m/s2
RDN gc = 9.832185873 m/s2
Notes:
-"GABC" does not check that a >= b >= c
-You could add X>Y? SF 99 after lines 03
and 01.
-This method supposes that the equatorial flattening is enough small to
use a linear approximation in k = (a2 - b2) /
b2
-For the Earth, k ~ 0.000022
-For larger k-values, the results could be disappointing.
-The mean equatorial radius of the Earth is 6378136.61 m
and the polar flattening is f = 1/298.256421
which leads to a polar radius = c = 6356751.868 m
-Other choices are clearly possible. For instance, if you prefer the WGS 84 constants:
polar radius = c = 6356752.314 m
mean equatorial radius = 6378137 m
2°) Generalized Somigliana's Formula
-At the surface of the ellipsoid ( altitude h = 0 ), the gravity at a point whose longitude = L and latitude = B is given by
g(0) = ( a ga Cos2 L' Cos2 B + b gb Sin2 L' Cos2 B + c gc Sin2 B ) / d
with d2 = a2 Cos2 L' Cos2 B + b2 Sin2 L' Cos2 B + c2 Sin2 B and L' = L + 14°92911
( The longitude of the major equatorial axis is 14°92911
W )
-If h # 0, an approximate formula is used:
g(h) = g(0) [ 1 - 2 (h/a') ( 1 + f + m - 2 f Sin2
B ) + 3 sign(h) (h2/a'2) ]
with a' = (a+b)/2 and m = a.b.c w2 / GM
w = 7.292115
E-5 rd/s = angular velocity of the Earth
GM = 3.986004419 E14 m3/s2 = Earth gravitational
constant
Data Registers: R00 = 2.h/(a+b)
R01 = a R04 thru
R07: temp
R02 = b
R03 = c
Flags: /
Subroutine: "S-R" Spherical-Rectangular
conversion ( cf "Spherical Coordinates for the HP-41" paragraph 1°-a)
)
01 LBL "GRV" 02 DEG 03 ST+ X 04 STO 00 05 RDN 06 HR 07 X<>Y 08 HR 09 14.92911 10 + 11 1 12 XEQ "S-R" 13 X^2 14 6378171.645 |
15 STO 01 16 * 17 STO 04 18 9.780379982 19 * 20 X<>Y 21 X^2 22 6378101.575 23 STO 02 24 * 25 STO 05 26 9.780273549 27 * 28 + |
29 X<>Y 30 X^2 31 STO 07 32 6356751.868 33 STO 03 34 * 35 STO 06 36 9.832185871 37 * 38 + 39 RCL 02 40 RCL 01 41 + 42 ST/ 00 |
43 X<> L 44 RCL 04 45 * 46 RCL 02 47 RCL 05 48 * 49 + 50 RCL 03 51 RCL 06 52 * 53 + 54 SQRT 55 / 56 RCL 00 |
57 ABS 58 3 59 * 60 RCL 07 61 74.56410525 62 / 63 + 64 2.013605211 65 - 66 RCL 00 67 * 68 * 69 + 70 END |
( 172 bytes / SIZE 008 )
STACK | INPUTS | OUTPUTS |
Z | L ( ° ' " ) | g(0) |
Y | B ( ° ' " ) | g(0) |
X | h ( m ) | g(h) |
---Execution time = 7s---
Example1: L = 77°03'56" W , B = 38°55'17"2 N , h = 67m ( U.S. Naval Observatory , Washington D.C. )
-77.0356 ENTER^
38.55172 ENTER^
67
XEQ "GRV" >>>> g(67) = 9.800516081
m/s2
RDN g(0) = 9.800722840 m/s2
Example2: L = 116°51'50"4 W , B = 33°21'22"4 N h = 1706 m ( Mount Palomar Observatory )
-116.51504 ENTER^
33.21224 ENTER^
1706
R/S >>>> g(1709) = 9.790659652
m/s2
RDN g(0) = 9.795922927
m/s2
Notes:
-The longitudes are positive East.
-Latitudes & longitudes are both supposed geodetic, but using a geocentric
longitude will not change the result significantly.
-The values of ga , gb , gc
( lines 18-26-36 ) are obtained by GABC, rewritten for the HP48:
ga = 9.780379982
gb = 9.780273549
Replace lines 18-26-36 by other numbers if you prefer
gc = 9.832185871
-Unfortunately, the Earth gravitational field is much more complex than
that of a triaxial ellipsoid,
so this program only gives slightly better results than the classical
Somigliana's formula.
-The 3 parameters ga , gb , gc may also be chosen differently, but they should verify the Pizzetti theorem:
ga/a + gb/b + gc/c = 3.GM / (abc) - 2 w2
-The following variant uses 2 M-Code routines DOT ( dot product ) and
S-R ( Spherical-Rectangular conversion )
01 LBL "GRV" 02 DEG 03 ST+ X 04 STO 00 05 RDN 06 HR 07 X<>Y 08 HR 09 14.92911 10 + 11 1 12 S-R 13 X^2 |
14 6378171.645 15 STO 01 16 * 17 STO O 18 X<>Y 19 X^2 20 6378101.575 21 STO 02 22 * 23 STO N 24 R^ 25 X^2 26 STO 04 |
27 6356751.868 28 STO 03 29 * 30 STO M 31 9.780379982 32 9.780273549 33 9.832185871 34 DOT 35 RCL 01 36 RCL 02 37 RCL 03 38 DOT 39 SQRT |
40 ST/ T 41 R^ 42 RCL 00 43 RCL 01 44 RCL 02 45 + 46 / 47 STO 00 48 ABS 49 3 50 * 51 RCL 04 52 74.56410525 |
53 / 54 + 55 2.013605211 56 - 57 RCL 00 58 * 59 * 60 + 61 CLA 62 END |
( 167 bytes / SIZE 005 )
STACK | INPUTS | OUTPUTS |
Z | L ( ° ' " ) | g(0) |
Y | B ( ° ' " ) | g(0) |
X | h ( m ) | g(h) |
---Execution time = 6s---
Example1: L = 77°03'56" W , B = 38°55'17"2 N , h = 67m ( U.S. Naval Observatory , Washington D.C. )
-77.0356 ENTER^
38.55172 ENTER^
67
XEQ "GRV" >>>> g(67) = 9.800516081
m/s2
RDN g(0) = 9.800722840 m/s2
Example2: L = 116°51'50"4 W , B = 33°21'22"4 N h = 1706 m ( Mount Palomar Observatory )
-116.51504 ENTER^
33.21224 ENTER^
1706
R/S >>>> g(1709) = 9.790659652
m/s2
RDN g(0) = 9.795922927
m/s2
Notes:
-The results are identical !
-Lines 11-12 may also be replaced by X<>Y
1 P-R RCL Z X<>Y P-R
( 66 lines / 171 bytes )
-The alpha "register" is cleared.
Reference:
[1] Michèle Caputo - "Some Space Gravity Formulas and
the Dimensions and the Mass of the Earth"