gravity3axial

# Gravity on a Triaxial Ellipsoid for the HP-41

Overview

1°)  Calculating  ga , gb , gc
2°)  Generalized Somigliana's Formula

a)  Geodetic Longitudes & Latitudes
b)  Geocentric Longitudes & Latitudes

-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

-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

a)  Geodetic Longitudes & Latitudes

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

-Here is another variant with

a = 6378172 m            ga = 9.780378635 m/s2
b = 6378102 m            gb = 9.780272308 m/s2
c = 6356752.314 m     gc = 9.832184675 m/s2

Data Registers:            R00 thru  R05: temp
Flags: /
Subroutines:  /

 01 LBL "GRV"  02 DEG  03 ST+ X  04 STO 00            05 X<>Y  06 HR  07 1  08 P-R  09 R^  10 HR  11 14.92911  12 +  13 X<>Y 14 P-R  15 X^2  16 6378172  17 STO 01  18 STO 04            19 *  20 ST* 04  21 9.780378635  22 *  23 X<>Y  24 X^2  25 6378102  26 ST+ 01 27 STO 05  28 *  29 ST* 05  30 9.780272308  31 *  32 +  33 X<>Y  34 X^2  35 STO 02            36 6356752.314  37 STO 03  38 *  39 ST* 03 40 9.832184675  41 *  42 +  43 RCL 03  44 RCL 04            45 RCL 05  46 +  47 +  48 SQRT  49 /  50 RCL 00  51 RCL 01  52 /  53 STO 00 54 ABS  55 3  56 *  57 RCL 02            58 74.56430504  59 /  60 +  61 2.013605194  62 -  63 RCL 00  64 *  65 *  66 +  67 END

( 159 bytes / SIZE 006 )

 STACK INPUTS OUTPUTS Z L ( ° ' " ) g(0) Y B ( ° ' " ) g(0) X h ( m ) g(h)

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.800514846  m/s2                                      ---Execution time = 6s---
RDN    g(0)  = 9.800721604  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.790658424  m/s2
RDN       g(0)    = 9.795921698  m/s2

-Of course, if the geodetic longitude is measured from the equatorial major axis ( longitude = 14.92911 W ), delete lines 11-12

b)  Geocentric Longitudes & Latitudes

-In the following programs, the longitudes are measured from the equatorial major axis ( longitude = 14.92911 W ) and:

a = 6378172 m            ga = 9.780378635 m/s2
b = 6378102 m            gb = 9.780272308 m/s2
c = 6356752.314 m     gc = 9.832184675 m/s2

Data Registers:            R00 thru  R05: temp
Flags: /
Subroutines:  /

 01 LBL "GRV"  02 DEG  03 ST+ X  04 STO 00            05 X<>Y  06 HR  07 1  08 P-R  09 R^  10 HR  11 X<>Y  12 P-R  13 6356752.314  14 STO 03  15 X^2  16 ST/ T  17 CLX 18 6378102  19 STO 05            20 X^2  21 ST/ Z  22 CLX  23 6378172  24 STO 01  25 X^2  26 /  27 R-P  28 R^  29 X<>Y  30 R-P  31 SIGN  32 P-R  33 RCL Z  34 X<>Y 35 P-R  36 X^2  37 RCL 01  38 STO 04            39 *  40 ST* 04  41 9.780378635  42 *  43 X<>Y  44 X^2  45 RCL 05  46 ST+ 01  47 *  48 ST* 05  49 9.780272308  50 *  51 + 52 X<>Y  53 X^2  54 STO 02  55 RCL 03            56 *  57 ST* 03  58 9.832184675  59 *  60 +  61 RCL 03  62 RCL 04  63 RCL 05  64 +  65 +  66 SQRT  67 /  68 RCL 00 69 RCL 01  70 /  71 STO 00            72 ABS  73 3  74 *  75 RCL 02  76 74.56430504  77 /  78 +  79 2.013605194  80 -  81 RCL 00  82 *  83 *  84 +  85 END

( 173 bytes / SIZE 006 )

 STACK INPUTS OUTPUTS Z L ( ° ' " ) g(0) Y B ( ° ' " ) g(0) X h ( m ) g(h)

Example:         L = -101°56'06"  ,    B = +33°10'47" N    h = 1706 m

-101.5606    ENTER^
33.1047      ENTER^
1706          R/S     >>>>   g(1709) = 9.790658215  m/s2                                      ---Execution time = 8.7s---
RDN       g(0)    = 9.795921489  m/s2

Notes:

-If M is the point with altitude = h, the longitude L & latitude B are the coordinates of the point P with altitude = 0  ( MP normal to the ellipsoid )

-We can simplify the program if longitude & latitude are expressed in decimal degrees and if h is always 0:

Data Registers:            R00 thru  R02: temp
Flags: /
Subroutines:  /

 01 LBL "GRV"  02 DEG  03 1  04 P-R  05 X<>Y  06 6356752.314  07 STO 02            08 X^2  09 /  10 R^  11 RCL Z  12 P-R 13 6378102  14 STO 01            15 X^2  16 ST/ Z  17 CLX  18 6378172  19 STO 00  20 X^2  21 /  22 R-P  23 R^  24 X<>Y 25 R-P  26 SIGN  27 P-R  28 RCL Z  29 X<>Y  30 P-R  31 X^2  32 RCL 00            33 *  34 ST* 00  35 9.780378635  36 * 37 X<>Y  38 X^2  39 RCL 01            40 *  41 ST* 01  42 9.780272308  43 *  44 +  45 X<>Y  46 X^2  47 RCL 02  48 * 49 ST* 02  50 9.832184675  51 *  52 +  53 RCL 00            54 RCL 01  55 RCL 02  56 +  57 +  58 SQRT  59 /  60 END

( 126 bytes / SIZE 003 )

 STACK INPUTS OUTPUTS Y L ( deg ) / X B ( deg ) g(0)

Example:         L = -101°935  ,    B = +33°17972222 N

-101.935       ENTER^
33.17972222   XEQ "GRV"     >>>>   g(0)    = 9.795921489  m/s2                 ---Execution time = 7.2s---

Reference:

[1]   Michèle Caputo - "Some Space Gravity Formulas and the Dimensions and the Mass of the Earth"