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 [1] ):
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 [2] 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 [2] )
-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:
[1] http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf
[2] http://cddis.nasa.gov/926/egm96/egm96.html
[3] http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/index.html
-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...