Geogravity

# Earth's Gravity for the HP-41

Overview

1°)  Somigliana's Formula
2°)  Somigliana's Formula + Taylor expansion
3°)  Ellipsoidal Gravity: Rigorous Formulae
4°)  Via Earth's Gravitational Potential

1°)  Somigliana's Formula

-The gravity at the surface of an ellipsoid of revolution is given by the exact formula:

g = gE ( 1 + k sin2 b ) ( 1 - e2 sin2 b ) -1/2

gE = equatorial gravity = 9.780325336 m/s2
where       k  = ( b' gP/a/gE ) - 1 = 0.00193185265241   with   gP = polar gravity = 9.8321849378 m/s2
a   =  6378137 m =  equatorial radius of the Earth  ,  b' = 6.356752.314 m = polar radius

and         e  =  eccentricity of the Earth's meridian = 0.08181919084
b  =  geographic ( geodetic ) latitude

-If an accuracy of  10-6 m/s2 is enough and in order to save bytes, Somigliana's formula may be re-written:

g = 9.780325 + 0.051631 sin2 b + 0.000228 sin4 b + . . .

-A correction of  - 3.086 10 -6 h ( m/s2 ) is added where h = altitude ( in meters ).
-This is not very accurate but it is acceptable if h is small

Data Registers: /
Flags: /
Subroutines: /

 01  LBL "GRV"  02  DEG  03  X<>Y  04  HR  05  SIN  06  X^2  07  RCL X  08  228  09  *  10  51631  11  +  12  *  13  9780325  14  +  15  X<>Y  16  3.086  17  *  18  -  19   E6  20  /  21  END

( 47 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Y b ( ° ' " ) h X h ( m ) g ( m/s2 )

Example:     b = 38°55'17"2     h = 67m

38.55172   ENTER^
67         XEQ "GRV"  >>>>   g = 9.800533 m/s2        ---Execution time = 2s---

Notes:

-The longitude is not taken into account.
-Line 15 corresponds to the "free-air" correction.

-Another correction ( the so-called "Bouguer anomaly" ) may also be added.
-It takes into accout the extra-mass for an altitude h > 0. One usually takes:  rock density = 2.67  which leads to add  1.119 10 -6 h
-In this case,  replace line 15 by 1.967
-In the example above, it yields g = 9.800608 m/s2

2°)  Somigliana's Formula + Taylor Expansion

-A more accurate result is obtained if we use a Taylor series ( truncated to the terms in h2 ) to take the altitude into account:

g(h) = g(0) [ 1 - (2 h/a ) ( 1 + f + m - 2 f sin2 b ) + 3 h2/a2 ]

where       m =  ( omega )2 a2 b' / GM = 0.00344978650684

omega = Earth's angular velogity = 7.292115 10 -5 rad/s
with        a   =  6378137 m =  equatorial radius of the Earth  ,  b' = 6.356752.314 m = polar radius
GM =  3.986004418 1014  m3/s2 = geocentric gravitational constant

-"GRV2" uses  g(h) = g(0) [ 1 - ( h/a ) ( 2.013605194 - 0.013411243 sin2 b ) + 3 h2/a2 ]

Data Registers:  R00 & R01: temp
Flags: /
Subroutines: /

 01  LBL "GRV2" 02  DEG 03  X<>Y 04  HR 05  SIN 06  X^2 07  ENTER^ 08  STO 00 09  228 10  * 11  51631 12  + 13  * 14  9780325 15  + 16   E6 17  / 18  X<>Y 19  6378137       20  / 21  STO 01 22  3 23  * 24  RCL 00 25  74.5643 26  / 27  2.013605     28  - 29  + 30  RCL 01        31  * 32  1 33  + 34  * 35  END

( 76 bytes / size 002 )

 STACK INPUTS OUTPUTS Y b ( ° ' " ) g0 ( m/s2 ) X h ( m ) gh ( m/s2 )

Example1:     b = 38°55'17"2     h = 67m

38.55172   ENTER^
67         XEQ "GRV2"  >>>>   gh = 9.800533 m/s2                   ---Execution time = 3s---
RDN    g0 = 9.800739 m/s2

Example2:     b = 38°55'17"2     h = 23456 m

38.55172   ENTER^
23456         R/S        >>>>   gh = 9.728752 m/s2
RDN    g0 = 9.800739 m/s2

Note:

-We can also use Somigliana's original formula to get even more accurate results at low altitudes:

 01  LBL "GRV2B" 02  DEG 03  6378137 04  / 05  STO 00 06  X<>Y 07  HR 08  SIN 09  X^2 10  STO 01 11  .0188941474 12  * 13  9.780325336 14  + 15  1 16  RCL 01 17  .00669438 18  * 19  - 20  SQRT 21  / 22  RCL 00          23  3 24  * 25  .0134112427 26  RCL 01 27  * 28  + 29  2.013605194 30  - 31  RCL 00          32  * 33  1 34  + 35  * 36  END

( 100 bytes / SIZE 002 )

 STACK INPUTS OUTPUTS Y b ( ° ' " ) g0 ( m/s2 ) X h ( m ) gh ( m/s2 )

Example:     b = 38°55'17"2     h = 67m

38.55172   ENTER^
67         XEQ "GRV2B"  >>>>   gh = 9.800532949 m/s2                   ---Execution time = 4s---
RDN    g0 = 9.800739708 m/s2

3°)  Ellipsoidal Gravity: Rigorous Formulae

-Let  b = geographic ( geodetic ) latitude and  h = altitude
-The gravity g above an ellipsoid of revolution may be computed by the rigorous formulae ( cf reference  ):

g  =  ( p2+ q2 )1/2

where   p =  { - [ GM + ( omega )2 a2 E ( q' / 2q0 ) ( sin2 F - 1/3 ) ] / ( u2 + E2 ) + ( omega )2 u cos2 F } / w
and     q =   [ a2 ( q / q0 ) ( u2 + E2 ) -1/2 - ( u2 + E2 )1/2 ] ( omega )2 ( sin F cos F ) / w

u  =  { [ ( x2 + y2 + z2 - E2 ) / 2 ] [ 1 + { 1 + 4 E2 z2 / ( x2 + y2 + z2 - E2 )2 }1/2 ] }1/2
with    E =   a ( 2 f - f 2 )1/2  ,   sin F = z / u   ,   cos F = ( x2 + y2 )1/2 / ( u2 + E2 )1/2
w =  [ ( u2 + E2 sin2 F )1/2 / ( u2 + E2 )1/2 ]1/2

q  = (1/2) [ ( 1 + 3 u2 / E2 ) Atan ( E / u )  - 3 E / u ]                                                            ( x2 + y2 )1/2  = ( N + h ) cos b
and    q0 = (1/2) [ ( 1 + 3 b' 2 / E2 ) Atan ( E / b' ) - 3 E / b' ]         ( b' = polar radius )                          z            = [ ( b'/a )2 N + h ] sin b
q' =   3 ( 1 + u2 / E2 ) [ 1 - ( u / E ) Atan ( E / u ) ] - 1                                                                 N          = a ( 1 - e2 sin2 b ) -1/2

-If, however, we use ATAN to calculate  q , q0 and q' , the results will not be accurate because of cancellation of leading digits.
-So, the following program employs ascending series to compute:

Atan t + 3 [ ( Atan t ) / t - 1 ] / t      =  4 t3  Sumk=0,1,2,.....  ( k + 1 ) ( - t2 )k / [ ( 2 k + 3 ) ( 2 k + 5 ) ]
3 ( 1 + 1 / t2 ) [ 1 - ( Atan t ) / t ] - 1  =  6 t2  Sumk=0,1,2,.....  ( - t2 )k / [ ( 2 k + 3 ) ( 2 k + 5 ) ]

-Since the argument t never exceeds 0.0821 for the Earth, only a few terms must be calculated.

Data Registers:   R00 = a = Equatorial radius               R03 to R13: temp
R01 = f = Earth's flattening
R02 = ( omega )^2  where  omega = angular velocity of the Earth.

Flag:  F10 is cleared at the end
Subroutines: /

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

( 242 bytes / SIZE 014 )

 STACK INPUTS OUTPUTS Z / B ( deg ) Y b ( ° ' " ) R ( m ) X h ( m ) g ( m/s2 )

where   R = distance from the center of the Earth
and      B = geocentric latitude

Example1:     b = 38°55'17"2     h = 23456 m

38.55172   ENTER^
23456      XEQ "GRV3"  >>>>   g = 9.728750374 m/s2                               ---Execution time = 22s---
RDN    R = 6393195.112 m
RDN    B = 38°73415897

Example2:     b = 38°55'17"2     h = 12345678 m

38.55172   ENTER^
12345678   XEQ "GRV3"  >>>>   g = 1.078713288 m/s2                            ---Execution time = 19s---
RDN   R = 18715394.62 m
RDN   B = 38°85746766

Notes:

-If you don't want to calculate R & B, delete lines 136-134-133-132

-If you replace line 12 by  298.2572236 ( WGS84 ), you'll get:

with b =  0°  &  h = 0 ,  g = 9.780325335     ( almost exactly  gE = 9.7803253359 )
with b = 90° &  h = 0 ,  g = 9.832184940     ( almost exactly  gP = 9.8321849379 )

-When h = 0, these formulas are equivalent to Somigliana's formula.

-More recent evaluations suggest  a = 6378136.61 m  and  1/f = 298.256421
-With  b = 38°55'17"2  &  h = 23456 m , it yields  g = 9.728751603 m/s2

-One can use an M-Code routine to compute the ascending series, for instance:

Delete lines 139 to 177
replace line 67 by CHS  and lines  68-72-78 by @GR where @GR is listed hereunder:

( the arguments are always positive.  CHS is used to set CPU-flag 9 )

092   "R"
007   "G"
000   "@"
2A0   SETDEC              this routine does not check for alpha data
0F8   READ 3(X)
128   WRIT 4(L)
244   CLRF 9
2FE   ?C#0 MS
013   JNC +02
248   SETF 9
05E   C=0 MS          C= | C |
10E   A=C ALL
135   ?NCXQ             C=
060   184D                A*C
070   N=C ALL
04E   C=0 ALL            C
35C   PT=12                =
090    LD@PT-2          2
268    WRIT 9(Q)
0B0   C=N ALL
10E   A=C ALL
068   WRIT 1(Z)
04E   C=0 ALL            C
35C   PT=12
050    LD@PT-1          =
150    LD@PT-5
226    C=C+1 S&X     15
261    ?NCXQ            C=
060    1898                A/C
0E8    WRIT 3(X)
0B0    C=N ALL         loop begins here
10E    A=C ALL
078    READ 1(Z)
135    ?NCXQ            C=
060    184D                A*C
2BE   C=-C-1 MS     C=-C
068    WRIT 1(Z)
278    READ 9(Q)
028    WRIT 0(T)
00E   A=0 ALL            A
35C   PT=12                =
162    A=A+1@PT      1
01D   ?NCXQ            C=
060    1807                A+C
268    WRIT 9(Q)
10E   A=C ALL
01D   ?NCXQ            C=
060    1807                A+C
10E   A=C ALL
135    ?NCXQ            C=
060    184D                A*C
00E   A=0 ALL            A
1BE   A=A-1 MS         =
35C   PT=12
162    A=A+1@PT      -1
01D   ?NCXQ            C=
060    1807                A+C
10E   A=C ALL
078    READ 1(Z)
0AE   A<>C ALL
261    ?NCXQ            C=
060    1898                A/C
10E   A=C ALL
046   C=0 S&X           C
270   RAMSLCT         =
038   READDATA      T
0AE   A<>C ALL
24C   ?FSET 9
135    ?NCXQ            C=
060    184D                A*C
10E   A=C ALL
0F8   READ 3(X)
01D   ?NCXQ            C=
060    1807                A+C
10E   A=C ALL
0F8   READ 3(X)
0AE   A<>C ALL
0E8    WRIT 3(X)
3CC   ?KEY                       these 2 words written in blue may be deleted if you have a "newest" HP-41
360    ?C RTN                    simply press  ENTER^  ON  to stop an infinite loop
36E    ?A#C ALL
26F    JC -33h=-51d           If you have deleted  3CC  360 , replace this line by  27F  JC -31h=-49d
138    READ 4(L)
05E   C=0 MS
0AE   A<>C ALL
24C   ?FSET 9
135    ?NCXQ            C=
060    184D                A*C
0E8    WRIT 3(X)
3E0    RTN

-Now, flag F10 is unused and SIZE 010 is enough.
-The same results are returned in 13 seconds instead of 22s
-Registers  Z , T and Q are used
-The arguments must remain small ( for x = 1 execution time would be prohibitive )
-Press any key to stop the routine

4°)  Via Earth's Gravitational Potential

-The Earth's gravitational potential U may be computed by the formula ( spherical harmonics ):

U(R,L,B) = (GM/R) { 1 + Sumn=2,3,......,N  (a/R)n  Summ=0,1,2,.....,n  ( Cn,m cos m L + Sn,m sin m L ) [ k (2n+1) (n-m)! / (n+m)! ] Pn,m (sin B) }

with     R = distance from the center of the Earth ,  L = longitude,  B = geocentric latitude

GM = 3.986004415 1014  m3/s2 = geocentric gravitational constant                                         k = 1  if  m = 0
and        a   =  6378137 m =  equatorial radius of the Earth                                                                  k = 2  if  m # 0
Pn,m = (-1)m  Pnm  with  Pnm = associated Legendre function of the first kind.

-The dimensionless coefficients Cn,m and Sn,m are given in reference  up to degree and order N = 360
-"GRV3" takes  N = 5 but it's easy to use more terms if needed ( unfortunately not all terms however).

-Then, Earth's gravity vector  g = ( gR , gL , gB ) = Grad [ U(R,L,B) + (1/2)(omega)2 R2 cos2 B ]

where  omega = 7.292115 10 -5 rad/s = angular velocity of the Earth.

-Therefore, in spherical coordinates, the modulus   | g | = [ ( U/R )2 + ( 1/( R2 cos2B ) ) ( U/L )2 + ( 1/R2 ) ( U/B )2 ]1/2

-"GRV+" also uses the following recurrence relations:

Pn,m+1 - 2 m x ( 1 - x2 ) -1/2 Pn,m  =  [ m ( m - 1 ) - n ( n + 1 ) ]  Pn,m-1             ( x = sin B )
Pn,m+1 -  m x ( 1 - x2 ) -1/2 Pn,m  =  ( 1 - x2 )1/2 dPn,m/dx

with the starting values:   Pn,n  = ( 2 n - 1 ) !!  ( 1 - x2 )n/2  ,  Pn,n+1 = 0       where   ( 2 n - 1 ) !! = (2n-1)(2n-3)(2n-5)........5x3x1

Data Registers:    R00 & R17: temp                                     ( Registers R19 thru R54 are to be initialized before executing "GRV4" )

R01 = gR    R04 = partial sum   R07 = R          R10 = sin B    R13 = m             R16 = n(n+1)
R02 = gL    R05 = partial sum   R08 = L           R11 = a/R      R14 = Pn,m+1     R18 = counter
R03 = gB    R06 = partial sum   R09 = cos B     R12 = n         R15 = Pn,m     •  R19 thru R54 = Cn,m and Sn,m

•  R19 =  C2,0 = -0.484165339915E-03          •  R20 =  S2,0 =   0.000000000000E+00 = 0
•  R21 =  C2,1 =  -0.186987635955E-09         •  R22 =  S2,1 =  0.119528012031E-08
•  R23 =  C2,2 = 0.243920233891E-05            •  R24 =  S2,2 =  -0.140020057884E-05

•  R25 =  C3,0 = 0.957154667712E-06            •  R26 =  S3,0 =   0.000000000000E+00 = 0
•  R27 =  C3,1 = 0.202992864498E-05            •  R28 =  S3,1 =  0.248453317303E-06
•  R29 =  C3,2 = 0.904685216567E-06            •  R30 =  S3,2 = -0.618942656224E-06
•  R31 =  C3,3 = 0.721158450151E-06            •  R32 =  S3,3 =  0.141415394144E-05

•  R33 =  C4,0 = 0.539712296562E-06            •  R34 =  S4,0 =  0.000000000000E+00 = 0
•  R35 =  C4,1 = -0.536248859309E-06           •  R36 =  S4,1 = -0.473399003951E-06
•  R37 =  C4,2 = 0.350533469375E-06            •  R38 =  S4,2 =  0.662728713246E-06
•  R39 =  C4,3 = 0.990675816205E-06            •  R40 =  S4,3 = -0.200968156744E-06
•  R41 =  C4,4 = -0.188541606328E-06           •  R42 =  S4,4 =  0.308839948776E-06

•  R43 =  C5,0 = 0.686667576471E-07            •  R44 =  S5,0 =  0.000000000000E+00 = 0
•  R45 =  C5,1 = -0.621297740620E-07           •  R46 =  S5,1 = -0.941990411149E-07
•  R47 =  C5,2 = 0.652584722977E-06            •  R48 =  S5,2 = -0.323524051490E-06
•  R49 =  C5,3 = -0.452076694386E-06           •  R50 =  S5,3 = -0.214963972058E-06
•  R51 =  C5,4 = -0.295228720666E-06           •  R52 =  S5,4 =  0.498925339599E-07
•  R53 =  C5,5 = 0.174947930421E-06            •  R54 =  S5,5 = -0.669381580606E-06

Flags: /
Subroutines: /

-Line 04,  6378136.3  should be better
-Lines 19 and 21 must be changed according to the number of terms you want to use in the gravitational potential
-Line 150 is a three-byte GTO 01
-Line 151 may be replaced by 3986004.418 E8

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

( 262 bytes / SIZE 055 )

 STACK INPUTS OUTPUTS Z B (deg) / Y L (deg) / X R  (m) g ( m/s2)

where  R , L , B  are the geocentric coordinates.

R = distance from the center of the Earth
L = longitude,
B = geocentric latitude ( not geographic latitude )

-A program listed in "Terrestrial Geodesic Distance for the HP-41" ( paragraph 5-a )
may be used to calculate R & B from the altitude and the geographic latitude.

Example:      R = 6369806 m    L = -77.065556°    B = 38.733471°

38.733471   ENTER^
-77.065556    ENTER^
6369806     XEQ "GRV4"  >>>>   g = 9.800330872 m/s2            ---Execution time = 79s---

gR = R01 = -9.800278350  m/s2
gL = R02 =  0.000091744  m/s2
gB = R03 = -0.032085272  m/s2

•  If you want to use the spherical harmonics up to, say degree and order 10:

-Store the required coefficients in registers R19 to R144 ( in the order given in reference  )
-Replace line 19 by 10 and line 21 by 144

-With the same example, it yields ( in 4mn23s )     g = 9.800265047 m/s2

•  Likewise, with the spherical harmonics up to degree and order 12

-Store the coefficients in registers R19 to R194  ( 176 coefficients )
-Replace line 19 by 12 and line 21 by 194

-The HP-41 returns   g = 9.800287689 m/s2  ( in 6mn05s )

Notes:

-"GRV4" doesn't work if the latitude is exactly +/-90°
-The coefficients given in the reference below decrease very slowly and even with 176 coefficients,  the 4th decimal is still doubtful.
-EGM96 provides the spherical harmonics up to degree and order N = 360.
-The number of coefficients is  N2 + 3 N - 4
-So, it seems "difficult" to store all these numbers in a HP-41.

-With a SIZE 319, you can use the spherical harmonics up to degree and order 16.
-But you can also create a DATA file of 594 numbers ( store them in reverse order, delete the DSE 18 and replace  RCL IND 18 by  GETX )
and use the terms of order and degree 23.

-In fact, all the coefficients are needed at sea-level.
-But at high altitudes, say 1000 km, truncated series up to degree & order 12 provide a good accuracy.

Improvements:

-The component gL has actually a negligible effect on the modulus g
-So, the program may be simplified if we do not compute it!

-Furthermore, we can also save execution time if we multiply the coefficients C , S by the normalization factors.
-This may be done by the following routine:

LBL "CS-CS"          LBL 01            ST+ Y               ST* Y              /                          RTN               STO 02
STO 00                   RCL 02           RCL 02             +                    SQRT                   ISG 02            ISG 01
2                              INT                  INT                  RCL 01          ST* IND 00         GTO 01           CLX
STO 01                   X#0?                 -                      RCL 02           ISG 00                 RCL 02           GTO 01
E3                          SIGN                FACT              INT                ST* IND 00          INT                 END
/                              RCL 01             ST* Y              +                    ISG 00                   E3
STO 02                  ENTER^            +                     FACT             X=0?                     /                      ( 70 bytes )

-It takes the control number of the block of registers that contain Cn,m and Sn,m and it multiplies them by the normalization factors.
-This control number bbb.eee  must verify  eee = bbb + N2 + 3 N - 5  and   bbb > 002  ( bb = 15 now ).

-After that, we can use the simplified version hereunder:

Data Registers:    R00 = a/R                                                     ( Registers R15 thru R50 are to be initialized before executing "GRV5" )

R01 = gR    R03 = partial sum   R05 = R          R08 = sin B    R11 = n(n+1)
R02 = gB    R04 = partial sum   R06 = L           R09 = n         R12 = Pn,m+1           R14 = counter
R07 = cos B     R10 = m        R13 = Pn,m

•  R15 thru  •  R50 = Cn,m and Sn,m  in the same order as above

Flags: /
Subroutines: /

-Lines 18 and 20 must be changed according to the number of terms you want to use in the gravitational potential

 01  LBL "GRV5"        02  DEG   03  STO 05   04  6378137   05  X<>Y   06  /   07  STO 00   08  RDN   09  STO 06   10  CLX   11  STO 01   12  STO 02   13  SIGN   14  P-R   15  STO 07   16  X<>Y   17  STO 08   18  5   19  STO 09   20  50   21  STO 14   22  LBL 01   23  CLX   24  STO 03   25  STO 04   26  STO 12   27  RCL 09 28  ENTER^   29  STO 10                30  ST* Y   31  +   32  STO 11   33  LASTX   34  ST+ X   35  DSE X   36  FACT   37  RCL 09   38  DSE X   39  FACT   40  2   41  LASTX   42  Y^X   43  *   44  /   45  RCL 07   46  RCL 09   47  Y^X   48  *   49  LBL 02   50  STO 13   51  RCL 06   52  RCL 10   53  *   54  1 55  P-R   56  X<>Y   57  RCL IND 14        58  *   59  DSE 14   60  X<>Y   61  RCL IND 14   62  *   63  +   64  *   65  ST+ 03   66  LASTX   67  RCL 12   68  RCL 08   69  RCL 13   70  STO 12   71  *   72  RCL 10   73  *   74  RCL 07   75  /   76  STO T   77  -   78  *   79  ST+ 04   80  LASTX   81  R^ 82  -   83  RCL 10   84  DSE 14                85  X=0?   86  GTO 00   87  ENTER^   88  DSE X   89  ""    90  STO 10   91  *   92  RCL 11   93  -   94  /   95  GTO 02   96  LBL 00   97  RCL 00   98  ST* 01   99  ST* 02 100  RCL 09 101  RCL 03 102  ST* Y 103  + 104  ST+ 01 105  RCL 04 106  ST+ 02 107  DSE 09 108  RCL 09 109  DSE X 110  GTO 01 111  3986004.415 E8 112  RCL 05 113  X^2 114  / 115  ST* 02 116  CHS 117  ST* 01 118  ST+ 01 119  RCL 07 120  7.292115 E-5 121  X^2 122  RCL 05 123  * 124  RCL 07 125  * 126  * 127  ST+ 01 128  RCL 08 129  LASTX 130  * 131  ST- 02 132  RCL 01 133  RCL 02 134  R-P 135  END

( 200 bytes / SIZE 051 )

 STACK INPUTS OUTPUTS Z B (deg) / Y L (deg) / X R  (m) g ( m/s2)

where  R , L , B  are the geocentric coordinates.

R = distance from the center of the Earth
L = longitude,
B = geocentric latitude

-Before using "GRV5" , store the coefficients in R15 thru R50 and key in   15.050  XEQ "CS-CS"

Example:   With the same values:   R = 6369806 m    L = -77.065556°    B = 38.733471°

38.733471   ENTER^
-77.065556    ENTER^
6369806     XEQ "GRV5"  >>>>   g = 9.800330872 m/s2            ---Execution time = 47s---     ( instead of 77s )

•  If you want to use the spherical harmonics up to, say degree and order 10:

-Store the required coefficients in registers R15 to R140
-Replace line 18 by 10 and line 20 by 140
-Key in  15.140  XEQ "CS-CS"

-With the same example, it yields     g = 9.800265047 m/s2   ( in 2m30s instead of 4mn23s )

•  Likewise, with the spherical harmonics up to degree and order 12

-Store the coefficients in registers R15 to R190  ( 176 coefficients )
-Replace line 18 by 12 and line 20 by 190
-Key in  15.190  XEQ "CS-CS"

-The HP-41 returns   g = 9.800287689 m/s2  ( in 3m26s instead of 6m05s )

References:

-The new gravity model EGM2008 uses spherical harmonics up to order and degree 2159 + a few additional terms!
-Even more impossible for an HP-41...