hp41programs

gravity3axial

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.274 m
    b = 6378101.946 m
    c = 6356751.868 m

    6378171.274   ENTER^
    6378101.946   ENTER^
    6356751.868   XEQ "GABC"   >>>>   ga = 9.780379416  m/s2                ---Execution time = 18s---
                                                      RDN   gb = 9.780274115  m/s2
                                                      RDN   gc = 9.832185879  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
-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:

      ga = 9.780378119  m/s2
      gb = 9.780272819  m/s2
      gc = 9.832184675  m/s2
 

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.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/s2
                                          RDN    g(0)  = 9.800723037  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.790660016  m/s2
                                     RDN       g(0)    = 9.795923291  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"
-But because of roundoff-errors, the last decimals are doubtful:
-An HP-48 - which works with 12 digits - returns

    ga = 9.780379417
    gb = 9.780274111            Replace lines 18-26-36 by these numbers if you prefer
    gc = 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 ga , gb , gc may also be chosen differently to get a better fit, 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.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/s2
                                          RDN    g(0)  = 9.800723037  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.790660016  m/s2
                                     RDN       g(0)    = 9.795923291  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"