gravity3axial

Gravity on a Triaxial Ellipsoid for the HP-41

Overview

1°) Calculating g_a , g_b , g_c
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.

x²/a² + y²/b² + z²/c² = 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 g_a , g_b , g_c

-The formulas are given in reference [1] ( warning: M. Caputo denotes a₁ = b , a₂ = a , a₃ = c ). With small differences in the notations:

Let E² = (b² - c²) / c² and n = (a² - b²) / b²

A_i,j = A'_i,j + n A"_i,j

A'_1,1 = A'_1,2 = A'_2,1 = A'_2,2 = [ 0.75 / (b² - c²)^5/2 ] [ Arc tan E - (E/3) (5.E² + 3) / (1+E²)² ]
A"_1,1 = [ 5 / 16 / (b² - c²)^7/2 ] [ -Arc tan E + ( E/15/(1+E²) ) (20 - (5-13.E⁴)/(1+E²)² ) ] A"_2,1 = A"_1,2 = 3 A"_1,1 ; A"_2,2 = 5 A"_1,1

A'_1,3 = A'_2,3 = [ 3 / (b² - c²)^5/2 ] [ -Arc tan E + ( E/3/(1+E²) ) (2E²+3) ]
A"_1,3 = [ 15 / 8 / (b² - c²)^7/2 ] [ Arc tan E - (E/30) ( 25 + (5-9.E⁴) / (1+E²)² ) ] A"_2,3 = 3 A"_1,1

-Then, we must calculate k₁ & k₂ with

D.K₁ / w² = - A'_1,1 - n [ A'_1,1 + 6 A"_1,1 + (c²/b²) A'_2,3 ] w² = angular velocity of the Earth = 7.292115 E-5 rd/s
D.K₂ / w² = - A'_1,1 - n [ A'_1,1 - (c²/b²) A'_2,3 ]

where D = 4 A'_1,1 [ 2 A'_1,1 - (c²/b²) A'_2,3 ] - 2 n (c²/b²) A'_2,3 ( A'_1,1 + 6 A"_1,1 ) + 4 n A'_1,1 [ 2 A'_1,1 + 12 A"_1,1 - (c²/b²) A"_2,3 ]

-Finally, with GM = Earth gravitational constant = 3.9860044188 E14 m³/s²

    • g_a/a = ( GM + 4 K₂ / a² ) / ( abc ) - 2 ( A_1,2 K₁ + 3 A_2,2 K₂ ) - w²
    • g_b/b = ( GM + 4 K₁ / b² ) / ( abc ) - 2 ( 3 A_1,1 K₁ + A_2,1 K₂ ) - w²
    • g_c/c =      GM / ( abc ) - 2 ( A_1,3 K₁ + A_2,3 K₂ )

-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.E² + 3) / (1+E²)² = (8/15) E⁵ - (8/7) E⁷ + (16/9) E⁹ - (80/33) E¹¹ + (40/13) E¹³ - (56/15) E¹⁵ + ...
-Arc tan E + ( E/15/(1+E²) ) (20 - (5-13.E⁴)/(1+E²)² ) = -(16/35) E⁷ + (64/45) E⁹ - (32/11) E¹¹ + (64/13) E¹³ - (112/15) E¹⁵ + (896/85) E¹⁷ - ....

-Arc tan E + ( E/3/(1+E²) ) (2E²+3) = (2/15) E⁵ - (4/21) E⁷ + (2/9) E⁹ - (8/33) E¹¹ + (10/39) E¹³ - (4/15) E¹⁵ + ...
Arc tan E - (E/30) ( 25 + (5-9.E⁴) / (1+E²)² ) = -(8/105) E⁷ + (8/45) E⁹ - (16/55) E¹¹ + (16/39) E¹³ - (8/15) E¹⁵ + (56/85) E¹⁷ - ....

-For the Earth, E² = (b² - c²) / c² ~ 0.0067 so these terms are sufficient. More terms should be used for larger E-values.

Data Registers: R00 = n = (a² - b²) / b²

                                        R01 = a           R04 = g_a         R07 thru R14: temp
                                        R02 = b           R05 = g_b
                                        R03 = c           R06 = g_c
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	g_c
Y	b	g_b
X	c	g_a

where a >= b >= c are expressed in meters & g in m/s²

Example:

    a = 6378171.274 m
    b = 6378101.946 m
    c = 6356751.868 m

    6378171.274   ENTER^
    6378101.946   ENTER^
    6356751.868   XEQ "GABC"   >>>>   g_a = 9.780379416 m/s²                ---Execution time = 18s---
                                                      RDN   g_b = 9.780274115 m/s²
                                                      RDN   g_c = 9.832185879 m/s²

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 = (a² - b²) / b²
-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
-The equatorial flattening is about 1/92000
-So, I've added and subtracted 6378136.61 / 184000 to the mean equatorial radius of the Earth to get a = 6378171.274 m & b = 6378101.946 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 which leads - with the same method as above - to
a = 6378171.664 m
b = 6378102.336 m

-With these semi-axes, "GABC" returns:

      g_a = 9.780378119 m/s²
      g_b = 9.780272819 m/s²
      g_c = 9.832184675 m/s²

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 g_a Cos² L' Cos² B + b g_b Sin² L' Cos² B + c g_c Sin² B ) / d

with d² = a² Cos² L' Cos² B + b² Sin² L' Cos² B + c² Sin² 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 Sin² B ) + 3 sign(h) (h²/a'²) ] with a' = (a+b)/2 and m = a.b.c w² / GM

w = 7.292115 E-5 rd/s = angular velocity of the Earth
GM = 3.986004419 E14 m³/s² = 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.274

15 STO 01
16 *
17 STO 04
18 9.780379416
19 *
20 X<>Y
21 X^2
22 6378101.946
23 STO 02
24 *
25 STO 05
26 9.780274115
27 *
28 +

29 X<>Y
30 X^2
31 STO 07
32 6356751.868
33 STO 03
34 *
35 STO 06
36 9.832185879
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.800516278 m/s²
                                          RDN    g(0) = 9.800723037 m/s²

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.790660016 m/s²
                                     RDN       g(0)    = 9.795923291 m/s²

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 g_a , g_b , g_c ( lines 18-26-36 ) are obtained by "GABC"
-But because of roundoff-errors, the last decimals are doubtful:
-An HP-48 - which works with 12 digits - returns

    g_a = 9.780379417
    g_b = 9.780274111            Replace lines 18-26-36 by these numbers if you prefer
    g_c = 9.832185874

-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 g_a , g_b , g_c may also be chosen differently to get a better fit, but they should verify the Pizzetti theorem:

g_a/a + g_b/b + g_c/c = 3.GM / (abc) - 2 w²

-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.274
15 STO 01
16 *
17 STO O
18 X<>Y
19 X^2
20 6378101.946
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.780379416
32 9.780274115
33 9.832185879
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.800516278 m/s²
                                          RDN    g(0) = 9.800723037 m/s²

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.790660016 m/s²
                                     RDN       g(0)    = 9.795923291 m/s²

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"

hp41programs

Gravity on a Triaxial Ellipsoid for the HP-41