(49g 50g) Theoretical Earth gravity g = g(latitude, height), WGS84, GRS80/67

[Gil]

HP49-50G
Gravity g calculation in function of two arguments:
- latitude (in d.mmss), in stack level 2
- height/altitude (in m), in last stack level 1,

According to equations & parameters:
- International gravity 1967
- WGS 84.

See, for example, for more details:
https://en.m.wikipedia.org/wiki/Theoretical_gravity

The code is:

\<< "https://en.m.wikipedia.org/wiki/Theoretical_gravity
Version 1

2 Arg
. 
lat
 [in D.mmss]
. 
alt
 [in m]
" DROP "alt [m]" \->TAG SWAP "lat [D.mmss]" \->TAG SWAP DUP2 \-> lat alt
\<< lat HMS\-> 'lat' STO DEG '9.780327*(1+.0053024*SIN(lat)^2-.0000058*SIN(lat*2)^2)-.000003086*alt' \->NUM "
Int.grav 1967
" \->TAG '9.7803253359*((1+1.9318526464E-3*SIN(lat)^2)/\v/(1-6.69437999014E-3*SIN(lat)^2))-(1-.00139*SIN(lat)^2)*.0000030877*alt+7.2E-13*alt^2' \->NUM "
WGS 84
" \->TAG
\>>
\>>

Regards,
Gil Campart

[Albert Chan]

(09-21-2021 07:04 AM)Gil Wrote:  See, for example, for more details:
https://en.m.wikipedia.org/wiki/Theoretical_gravity

It is amazing that formula has so many significant digits



But, constants only matched 6-significant digits from https://en.wikipedia.org/wiki/Gravity_of...tude_model

[Gil]

You are perfectly correct — and in a way I am wrong.

The truth is the following : the formula in Wikipedia showed some wrong digits at the end, as well as the Chinese page.

To check, I used the a and f exact values of WGS84 and cut the decimals points and added zeros accordingly to operate with integers
(see my program:
\<< "Version 5.2

1 single Arg like
\[]'7/3*'
\[] or 
{
'7/3*' 
300}

for 
300
 digits
" DROP DUP TYPE 5 ==
IF
THEN OBJ\-> DROP
ELSE 100 "Put above 
200
 if
you want by default

200
 digits & not 100" DROP
END SWAP DUP UNROT 0 0 0 RCLF \-> digit x1 x2 x21 num f
\<< RAD STD -3 CF -105 CF digit \->STR "." "" SREPL DROP OBJ\-> 'digit' STO x1 \->STR "." "" SREPL 0 ==
IF
THEN DROP
ELSE OBJ\-> 'x2' STO x2 x1 / \->NUM LOG DUP FP 0 \=/
IF
THEN DROP "Instead of decimals (ab.c),

Try 
fractions

('abc/10') !" DOERR
END \->STR "." "" SREPL DROP OBJ\-> ALOG 'x21' STO
IF x21 1 >
THEN x2 x21 /
ELSE
IF x21 1 <
THEN x2 x21 INV *
ELSE x2
END
END \->STR "." "" SREPL DROP OBJ\->
END DUP EXPAND DUP2 SAME DROPN DUP \->NUM DUP 'num' STO num ABS 100000000000 > num FP 0 \=/ OR
IF
THEN OVER 10 digit ^ * PROPFRAC PROPFRAC -105 SF DUP TYPE 9 ==
IF
THEN OBJ\-> 3 DROPN
END
ELSE DROP
END f STOF
\>>
\>>

and get then the most accurate value for b.

I did the same when calculating the constants k and e².

So that the digits shown on the English and Chinese page for the WGS84 g-formulae are now all "correct", though meaningless, as you noticed.

The problem was : suppose I have an "effective")result
123456789012.6

Should I cut it into
123456789012 (the first 12 digits are correct)

Or prefer to have something printed incorrectly
123456789013
But that is nearer of the true value (and "better" for further calculations [on my calculator]).

I chose the second solution and decided to give the most complete values of k and e², leaving the choice for the reader to cut where it's the most "convenient" for him.

As I am limited on the digits of the values entered on my HP (for non-integers), you will see that the
constant k and e² in my programs do approximate correctly the theorical values (with many digits) of k and e².

But I could not decently write the initial values, for checking purposes,
to be 123456789013.

In fact, I cut the final digits of the k and e² values in Wikipedia, being sure that the first cut digit from the left was less than 5.


1234567890123456.
Could be cut to
1234567890123 × 10³ (because after the last digit 3 there is a digit < 5, here 4).

But not to
12345678901234 × 10², as
12345678901235 × 10² would be better (but not "nice looking", as the digit after 123 is a 4 and not a 5 as shown).

[Gil]

Just added the degrees calculation mode (below in bold) :

\<< "https://en.m.wikipedia.org/wiki/Theoretical_gravity
Version 1.b

2 Arg
. 
lat
 [in D.mmss]
. 
alt
 [in m]
" DROP "alt [m]" \->TAG SWAP "lat [D.mmss]" \->TAG SWAP DUP2 RCLF \-> lat alt f
\<< DEG lat HMS\-> 'lat' STO DEG '9.780327*(1+.0053024*SIN(lat)^2-.0000058*SIN(lat*2)^2)-.000003086*alt' \->NUM "
Int.grav 1967
" \->TAG '9.7803253359*((1+1.9318526464E-3*SIN(lat)^2)/\v/(1-6.69437999014E-3*SIN(lat)^2))-(1-.00139*SIN(lat)^2)*.0000030877*alt+7.2E-13*alt^2' \->NUM "
WGS 84
" \->TAG f STOF
\>>
\>>

Regards,
Gil

[Gil]

Please note for the value calculation of k and e² (e squared):

Most Internet sites and official papers or books show inaccurate "ending" digits, as they use too few digits for a truncate value of the axe b (Pole radius).


For instance they might give
b = 6 356 752.341
or b=6 356 752.3
or b=6 356 752.3412
or b=6 356 752.34125 asf,

instead of more digits, something like
b=6356752.31424517949756...

in order to be able to precess to a more appropriate calculation of k and e² (e squared) after roundings.

Regards,
Gil
.

[Gil]

Version 2b

Just changed some minor details, the main ones being the results labels.

Regards,
Gil
.
\<< "Version 2b 2021.09.25

2 Arg
. 
lat
 [in D.mmss]
. 
alt
 [in m]

https://en.m.wikipedia.org/wiki/Theoretical_gravity
https://eu.docworkspace.com/d/sIEG9949c5qq2igY
(GSR 80 by H Moritz)

\[]g(90,0)=9.8321849378
but form\->9.83218493787
\[]\GD calc poss. for alt
" DROP "alt [m]" \->TAG SWAP "lat D.mmss" \->TAG DUP2 RCLF \-> alt lat f
\<< DEG lat HMS\-> 'lat' STO DEG '9.780327*(1+.0053024*SIN(lat)^2-.0000058*SIN(lat*2)^2)-.000003086*alt' \->NUM "
Lambert GRS 80
" \->TAG '9.7803267715*(1+.0052790414*SIN(lat)^2+.0000232718*SIN(lat)^4+.0000001262*SIN(l​at)^6+.0000000007*SIN(lat)^8)-(1-.00139*SIN(lat)^2)*.0000030877*alt+7.2E-13*alt^2' \->NUM "
Somigliana GRS 80
" \->TAG '9.7803253359*((1+1.9318526464E-3*SIN(lat)^2)/\v/(1-6.69437999014E-3*SIN(lat)^2))-(1-.00139*SIN(lat)^2)*.0000030877*alt+7.2E-13*alt^2' \->NUM "
Somigliana WGS 84
" \->TAG f STOF
\>>
\>>

[Gil]

Version 3
Improved the details for the exact "correction" when taking into account the height/altitude for the calculation of the Earth gravity g.

Now the full parameters/constants a, f and m appear explicitly in the formulae with their corresponding 12 significant digits values according to Somigliana GRS 80 and Somigliana WGS 84.

I let here, however, a simplified version for Lambert GRS 80 regarding the height "correction".

See below in bold the main changes in the code:
\<< "Version 3 2021.09.26

2 Arg
. 
lat
 [in D.mmss]
. 
alt
 [in m]

https://en.m.wikipedia.org/wiki/Theoretical_gravity
https://eu.docworkspace.com/d/sIEG9949c5qq2igY
(GRS 80 by H Moritz)
http://www.in-dubio-pro-geo.de/index.php...lip/ngrav1

\[]gWGS84(lat 90, alt 0)
effect =9.8321849378
but form\->9.83218493787
" DROP "alt [m]" \->TAG SWAP "lat D.mmss" \->TAG DUP2 RCLF \-> alt lat f
\<< DEG lat HMS\-> 'lat' STO DEG '9.780327*(1+.0053024*SIN(lat)^2-.0000058*SIN(lat*2)^2)-.000003086*alt' \->NUM "
Lambert GRS 80
" \->TAG '9.7803267715*(1+.0052790414*SIN(lat)^2+.0000232718*SIN(lat)^4+.0000001262*SIN(l​at)^6+.0000000007*SIN(lat)^8)' \->NUM lat alt 6378137 3.35281068118E-3 3.44978600308E-3 \-> lat h a f m
\<< '1-2/a*(1+f+m-2*f*SIN(lat)^2)*h+3*(h/a)^2' * \->NUM
\>> "
Somigliana GRS 80
" \->TAG '9.7803253359*((1+1.9318526464E-3*SIN(lat)^2)/\v/(1-6.69437999014E-3*SIN(lat)^2))' \->NUM lat alt 6378137 298.257223563 INV 3.44978650684E-3 \-> lat h a f m
\<< '1-2/a*(1+f+m-2*f*SIN(lat)^2)*h+3*(h/a)^2' * \->NUM
\>> "
Somigliana WGS 84
" \->TAG f STOF
\>>
\>>

Regards,
Gil

[Gil]

To check,
the paper variables of Geodetic Reference System 1980, by H.Moritz, should be taken, instead of the simplifications given in Wikipedia "Theorical Gravity", Chinese page.

So

j.e =
'GM/(a*b)*(1-m-m/6*ƒé²*(q0´/q0))'
9.78032533482, from the above formulae
9.7803253359 official
—> almost equal value

j.p =
'GM/(a*a)*(1+m/3*ƒé²*(q0´/q0))'
9.83218494001, from the above formulae
9.8321849378, official
—> almost equal value!

With
ƒé² = sqrt (é²) =e'

And:

q0´ =
'3*(1+1/é²)*(1-1/ƒé²*ATAN(ƒé²))-1' .00268804118

q0 =
'((1+3/é²)*ATAN(ƒé²)-3/ƒé²)/2'
.00007334625

é² =
'(a*a-b*b)/(b*b)'
6.73949674208E-3

GM =
3.986004418E14

m =
'w*w*a*a*b/GM'
3.44978650683E-3

w =
.00007292115
a 6378137

a =
6378137

b =
'a-a/298.257223563'
6356752.31425

The differences are now quite small, due to the roundings of the calculator.

Regards,
Gil

[Gil]

Version 6e
(Theoretical) GRAVITY of Earth g
= g(latitude [D.mmss] ; height [m]).

With latitude [D.mmss] in stack level 2
and height [m]) in stack level 1.

Returns 4 results:

-g GRS67, according to International Gravity equation;
-g GRS80, according to Somigliana's equation;
-g WGS84, according to Somigliana's equation;
- g FREE, according to closed form,
Li & GÖtze: 'Tutorial Ellips,geoid,gravity'.

Main change 1

-You can choose your "FREE" ellipsoid.
-How?
Go in FREE directory (inside main file GRAVITY Dir) and
change/save any value of the four (normally fixed) variables GM, a, f or w.
-But don't suppress any of them (modify, yes ; delete, no).
-Then EXECUTE in that FREE directory
—>gFREE,
with, as usual, latitude [D.mmss] in stack 2
and height [m] in stack level 1.
-Everything is then calculated inside this FREE Dir automatically in a CLOSED form: no need like in SOMIGLIANA to compute/have the intermediary g official values at Equator & at Pole.
-The result will appear with the label/tag "CLOSED FORM FREE"
or "CLOSED FORM WGS84" (if all the 4 values GM a f w to be found in FREE Dir are the same as GM a f w official GSM84-Values located in G84EP Dir).
-The final, CLOSED result is quite accurate.
-However generally somewhat less accurate than the Somigliana's equation for GRS80 & WGS84.

Main change 2

-When executing —>g
(in the main GRAVITY Dir),
the program —>gFREE (to be found in FREE Dir & discussed above) will be launched automatically, with the corresponding label /tag "CLOSED FORM FREE"
or "CLOSED FORM WGS84" added to the final g result.

Main change 3

-Besides the FREE Dir, a g84EP Dir was added inside the main GRAVITY Dir.
-In the NAME of the Dir g84EP, 84 stands for WGS84, E for Equator and P for Pole.
-The four variables GM a f w in that file G84EP Dir belong to the official WGS84 model and therefore should not be deleted or even modified.
-The intermediary variables/equations inside that G84EP Dir are commonly to be found in the literature.
- They show a different, instructive way of calculating g WGS84 at Equator and at Pole when having the four, official WGS84 fixed Values GM, a, f and w.
-In that g84EP Dir, the final calculated results for g WGS84 at Equator and at Pole, though quite accurate, are — unfortunately — not perfectly in adequation with the official WGS84 g values at Equator and at Pole.

Main change 4

Most explanations/references are now given inside NOTES inside the adequate directories:

GRAVITY DIR: NOTE1 NOTE2 NOTE3 NOTE4

GRAVITY FREE Dir: NOTE

GRAVITY g84EP: NOTE1 NOTE2.

But the version number of the whole program is soon to appear at the beginning of the main program —>g.

Summary & Conclusion

Some doubts remain regarding the best equations to be used when refering to (old dated) GRS67.

For GRS80 or WGS84, Somigliana's equations give here most accurate results (to almost full calculator digits capacities).

CLOSED FORM equation is best fit for theorical gravity, when modelling an ellipsoid, changing one or several of its four "fixed" parameters GM, a, f and w.

Numerical Examples
Executing in main file GRAVITY Dir
—>g

with latitude [D.mmss] in stack level 2
and height [m] in stack level 1

will result in the following outputs:

:alt [m]: 0 :lat D.mmss: 0 :

Int Grav GRS 67
: 9.780318 :

So
m
igliana GRS 80
: 9.7803267715 :
So
m
igliana WGS 84
: 9.7803253359 :
Closed Form
 
WGS 84
: 9.78032532324

:alt [m]: 1000 :lat D.mmss: 0 :

Int Grav GRS 67
: 9.777232 :

So
m
igliana GRS 80
: 9.77723980166 :
So
m
igliana WGS 84
: 9.77723836651 :
Closed Form
 
WGS 84
: 9.77723826177

:alt [m]: 0 :lat D.mmss: 90 :

Int Grav GRS 67
: 9.83217715816 :

So
m
igliana GRS 80
: 9.83218636846 :
So
m
igliana WGS 84
: 9.83218493787 :
Closed Form
 
WGS 84
: 9.83218496308

:alt [m]: 1000 :lat D.mmss: 90 :

Int Grav GRS 67
: 9.82909115816 :

So
m
igliana GRS 80
: 9.82910370419 :
So
m
igliana WGS 84
: 9.82910227404 :
Closed Form
 
WGS 84
: 9.82910231835

:alt [m]: 0 :lat D.mmss: 45 :

Int Grav GRS 67
: 9.80618987521 :

So
m
igliana GRS 80
: 9.80619920255 :
So
m
igliana WGS 84
: 9.80619776931 :
Closed Form
 
WGS 84
: 9.80619777838

:alt [m]: 1000 :lat D.mmss: 45 :

Int Grav GRS 67
: 9.80310387521 :

So
m
igliana GRS 80
: 9.80311437628 :
So
m
igliana WGS 84
: 9.80311294349 :
Closed Form
 
WGS 84
: 9.80311291004

Last example practice

-Go into FREE Dir.
-Change the value of
WGS84 1/298.257223563
into 1/298.25722356
(cut the final, right digit 3 in the expression),

' 1/298.25722356' ENTER
'f' STO

Then
0 0 —>gFREE will return the following results:

alt [in m]: 0 :lat [in D.mmss]: 0 :

Closed Form 

FREE
: 9.78032534903

and 90 0 —>gFREE will return the following results:

alt [in m]: 0 :lat [in D.mmss]: 90 :

Closed Form 

FREE
: 9.83218491167

Regards,
Gil Campart

[Gil]

RE: HP49-50G:
Theorical Earth gravity g = g(latitude, height), WGS84, GRS80/67

Last Version 9k
Added the following NOTE5 in main directory GRAVITY:

"When you see xx xx.1

xx.2
 in 
GRS67 GRS80
WGS84
 Dir [*] ,
prog. uses always xx

Instead of having
xx=xx.1
 (offic. Val)

You can choose in [*]
xx=xx.2
 (altern.Val)
do: xx.2
 'xx' STO"

Version 9i included here
Added the following NOTE in GRS67 directory:

"Official Val
ß.1 5.302
4
E-3
ß1.1 -5.
9
E-6
show inaccurateness

Taylor new calc. Val
ß.2 ~ 5.302
2
E-3
ß1.2 ~ -5.
8
E-6
give better results".

Previous version 9e

In comparison to previous, published version 9d (no more present in this post), changed minor details in directories GRS67, GRS80 and WGS84.

Use
Just enter your latitude (in D.MMSS) and ellipsoid height
(in m) and launch —>g.

All the calculated values of g (by —>g) will appear on the right of the main program —>g (pages 1 & 2), with Chebyshev's approximation formulae of orders 2 labelled 2C) and 8 (labelled 8C) and
Somigliana's equations (labelled S)), according to 4 ellipsoid models (GRS67, GRS80, WGS84 and your own* ellipsoid).

* To have/use your own ellipsoid, just go in FREE Dir (page 3 of main directory) and modify the four parameters GM, a, f and w (not b, as b is always calculated automatically by b= a-a×f inside —>gFREE program). Then, enter as usual your lat and height on the stack. Now, again as usual, go in main directory and launch —>G.

** To see the official (fixed) parameters used, enter the respective directories GRS67, GRS80 or WGS84. Try however not to change the given values here (inside those directories) — unless, of course, you discover an error. Some of the formulae are given to compare the resulting outputs with the fixed, official values (the latter, as previously mentioned, in principle not to be modified in the directories GRS67, GRS80 and WGS84); example e² with e².STR.

Derived results (included the b Pole radius variable that is not to be changed) for the closed form of your FREE Dir (with the four parameters GM, a, f and w that you chose): all of them calculated once & automatically (by —>gFREE Program) and saved in FREE Dir.

*** Note that often in text books or Internet the given values of e², k, b and m are shown as truncated, possibly with wrong, final digits, because of the truncation of several intermediary results.

To avoid such a "problem" (and get the maximun 12-digits accurateness of the calculator), the above mentioned variables were reckoned (like integers) with no intermediary truncations and shown/saved in form of strings with many digits (for example e². STR).
To get the corresponding numerical, "correct" value, execute OBJ—> command (as in the main program —>g), for instance e². STR OBJ—>.

To understand & compare
J2 derived variable in GRS67 :

Go to Main Menu
PAGE3
Enter GRS67 Dir
PAGE4
Press J2.CALC
Result:
FromEQ & Roundings:
1.08269887351E-3 :
Published J2: .0010827


"CODE J2.CALC"
"\\<< RCLF RAD f.STR OBJ\\-> \\->NUM \\-> f '(a^2-(a-a*f)^2)/a^2/3*(1-2/15*(w^2*a^2*(a-a*f)/GM)*(\\v/((a^2-(a-a*f)^2)/(a-a*f)^2)/(1/2*((1+3*(a-a*f)^2/(a^2-(a-a*f)^2))*ATAN(\\v/(a^2-(a-a*f)^2)/(a-a*f))-3*(a-a*f)/\\v/(a^2-(a-a*f)^2)))))' \\->NUM \"FromEQ & Roundings\" \\->TAG SWAP STOF J2 \"Published J2\" \\->TAG
\\>>"

To understand & compare
f derived variable in GRS80 :

Go to Main Menu
PAGE3
Enter GRS80 Dir
Press f.CALC :
Result :
From ROOT Approxim:
"'1/298.257166516" :
Published f.STR:
"'1/298.2572221008827112431628366'"

Code f.CALC"
"\\<< RCLF RAD 'J2=(a^2-(a-a*f.SOL)^2)/a^2/3*(1-2/15*(w^2*a^2*(a-a*f.SOL)/GM)*(\\v/((a^2-(a-a*f.SOL)^2)/(a-a*f.SOL)^2)/(1/2*((1+3*(a-a*f.SOL)^2/(a^2-(a-a*f.SOL)^2))*ATAN(\\v/(a^2-(a-a*f.SOL)^2)/(a-a*f.SOL))-3*(a-a*f.SOL)/\\v/(a^2-(a-a*f.SOL)^2)))))' 'f.SOL' .03 ROOT 'f.SOL' PURGE INV \\->STR \"'1/\" SWAP + \"From ROOT Approxim\" \\->TAG SWAP STOF f.STR \"Published f.STR\" \\->TAG

To understand & compare
J2 derived variable in WGS84 :

Go to Main Menu
PAGE3
Enter WGS84 Dir
PAGE4
Press J2.CALC
Result:
FromEQ & Roundings:
& Roundings: 1.08262891487E-3 :
Published J2.STR: "0.0010826298213129219"

"CODE J2.CALC"
"\\<< RCLF RAD f.STR OBJ\\-> \\->NUM \\-> f '(a^2-(a-a*f)^2)/a^2/3*(1-2/15*(w^2*a^2*(a-a*f)/GM)*(\\v/((a^2-(a-a*f)^2)/(a-a*f)^2)/(1/2*((1+3*(a-a*f)^2/(a^2-(a-a*f)^2))*ATAN(\\v/(a^2-(a-a*f)^2)/(a-a*f))-3*(a-a*f)/\\v/(a^2-(a-a*f)^2)))))' \\->NUM \"FromEQ & Roundings\" \\->TAG SWAP STOF J2 \"Published J2\" \\->TAG
\\>>"

Code —>g
(in main directory) :

\<< "
Version 9d 
2021.10.08



2 Arg
. 
lat
 [in D.mmss]
. 
alt
 [in

Check param in 
GRS67
GRS80 WGS84 FREE
 Dir

Possib change of 4 par

GM a f w
: in 
FREE 
Dir

See 
NOTES
 at Last Page
" DROP STD "alt [m]" \->TAG SWAP "lat D.mmss" \->TAG DUP2 RCLF UNROT DUP2 SWAP FREE \->gFREE 5 DROPN UPDIR UNROT DUP2 5 ROLL { GRS67 GRS80 WGS84 } \-> alt lat fg dir
\<< -105 SF DEG lat HMS\-> 'lat' STO DEG 1 3
FOR i dir i GET EVAL 'gE*(1+\Gb*SIN(lat)^2+\Gb1*SIN(2*lat)^2)' \->NUM alt f.STR OBJ\-> \->NUM m.STR OBJ\-> \-> h f m '1-2/a*(1+f+m-2*f*SIN(lat)^2)*h+3*(h/a)^2' \->NUM DUP 'corr' STO * 'gE*(1+\Ga*SIN(lat)^2+\Ga1*SIN(lat)^4+\Ga2*SIN(lat)^6+\Ga3*SIN(lat)^8)' \->NUM corr * k.STR OBJ\-> e\178.STR OBJ\-> \-> k e\178 'gE*((1+k*SIN(lat)^2)/\v/(1-e\178*SIN(lat)^2))' \->NUM corr * 'corr' PURGE UPDIR
NEXT UNROT DROP2 'gS84' STO 'gS80' STO 'g8C80' STO 'g2C80' STO 'gS67' STO 'g8C67' STO 'g2C67' STO gS80 "
So
m
igliana GRS 80
" \->TAG gS84 "
So
m
igliana WGS 84
" \->TAG fg STOF "
Closed Form
" WGS84 GM w a f.STR OBJ\-> * * * \->NUM UPDIR FREE GM w a f * * * \->NUM == " 
WGS 84
" " 
FREE
" IFTE + gFREE UPDIR DUP 'gFREE' STO SWAP \->TAG
\>>
\>>

Code —>gNEW
(in FREE Dir):


\<< "2 Arg
. 
lat
 [in D.mmss]
. 
alt
 [in m]


GM w a f 
can be modif!
" DROP "alt [in m]" \->TAG SWAP "lat [in D.mmss]" \->TAG DUP2 DEG HMS\-> 180 \pi / / \->NUM \-> h lat
\<< 'a-f*a' \->NUM 'b' STO RAD 'ATAN(b/a*TAN(lat))' \->NUM '\Gb' STO 'b*SIN(\Gb)+h*SIN(lat)' \->NUM 'z\180' STO 'a*COS(\Gb)+h*COS(lat)' \->NUM 'r\180' STO 'r\180^2-z\180^2' \->NUM 'd\180\180\178' STO 'r\180^2+z\180^2' \->NUM 'r\180\180\178' STO '\v/(a^2-b^2)' \->NUM 'E' STO 'd\180\180\178/E^2' \->NUM 'D' STO 'r\180\180\178/E^2' \->NUM 'R' STO '1/2+R/2' \->NUM '1/4+R^2/4-D/2' \->NUM 0 MAX \v/ - 0 MAX \v/ 1 MIN ACOS '\Gb\180' STO '\v/(r\180\180\178-E^2*COS(\Gb\180)^2)' \->NUM 'b\180' STO '1/2*((1+3*b^2/E^2)*ATAN(E/b)-3*b/E)' \->NUM 'q0' STO '3*(1+b\180^2/E^2)*(1-b\180/E*ATAN(E/b\180))-1' \->NUM 'q\180' STO '\v/((b\180^2+E^2*SIN(\Gb\180)^2)/(b\180^2+E^2))' \->NUM 'W' STO '1/W*(GM/(b\180^2+E^2)+w^2*a^2*E*q\180/((b\180^2+E^2)*q0)*(1/2*SIN(\Gb\180)^2-1/6)-w^2*b\180*COS(\Gb\180)^2)' \->NUM DUP 'gFREE' STO "
Closed Form 
" a f GM w * * * \->NUM UPDIR WGS84 a f.STR OBJ\-> GM w * * * \->NUM == "
WGS84
" "
FREE
" IFTE + \->TAG UPDIR FREE '(a^2-b^2)/a^2' \->NUM 'e\178' STO '(a^2-b^2)/b^2' \->NUM '\233\178' STO 'w^2*a^2*b/GM' \->NUM 'm' STO 'e\178/3*(1-2/15*m*(\v/\233\178/q0))' \->NUM 'J2' STO
\>>
\>>

Regards,
Gil

[Gil]

Version 10

In directories GRS67, GRS80, WGS84 and FREE
was added the corresponding ellipsoid (GRS67, GRS80, WGS84 and FREE)
Earth mean theorical g calculation (gmu).

Regards,
Gil

[Gil]

New version 12b

Minor changes ("nice to have"),
above all relative to gE (g at Equator), gE.STR, gP (g at Pole), gP.STR, k and k.STR values in GRS80 and WGS84 directory.

I'm particularly indebted to Charles Karney for his patience and insight regarding some of the formulae and results in relation to the above and to the topic in general.

Version 11f

With latitude and height h, calculates automatically
g2C67 g8C67 gS67 g2C80 g8C80 gS80 gS84 gSFR g±SFR gClFR.

g 2C67: g
2C67: 
g
ravity result
g
2
C67: 
2
nd order appro
g2
C
67: 
C
hebyshev
g2C
67
: GRS
67

g
8
C80: 
8
th order appro
g8C
80
: GRS
80 ±S S
modif

g
S
84 : 
S
omigliana
gS
84 
: WGS
84 
g
µ
:g
Mean

g
FR(EE)
:g
your FREE 
ellips
Cl : closed form for g calculation, without having to introduce gE / GP (self calculated)

As previously, you can calculate your own FREE ellipsoid, changing the values of GM w f a in FREE Dir (no need to look for g at Equator or at Pole).
Two methods automatically calculated : closed form or modified Somigliana's.

J2.CAL: is to get a close value of J2 in GRS67 WGS84 and FREE.
It is not used, however, to further calculations (seemingly to many rounding intermediate errors).

F.CAL: is to get a close value of f in WGS80.
It is not used, however, to further calculations (seemingly to many rounding intermediate errors).
Only f and f.STR values (more accurate) are used instead for further calculations.

Many self-explaining NOTES at the end of each directory.

[Gil]

Theorical Normal Earth Gravity g
for HP49-HP50

Gravity version 16
Double-checking and correction of the string variables .STR when 100 digits are shown after 150 digits precision calculation with MAXIMA free software.
Note:
“x1 x2... x99 7(501)“ will be registered as "x1 x2... x99 8“;
“x1 x2... x99 7(499)“ will be registered as "x1 x2... x99 7“.

Gravity version 15 (then incorrectly labelled 15g)
Many changes, among others regarding the flags status (to be always set back at its initial stage), the variable J2 (and references to it) and, above all, the MAXIMA directory, completely reformulated. For the latter:
- full automatisation (with, inside MAXIMA dir, single subroutine GRS67 or GRS80 or WGS84 or FREE) for the string to be copied in one step into MAXIMA software;
- in case of calculator lack of space/memory, semi-full automatisation (with, instead of subroutine FREE in MAXIMA directory, subroutine FR3ST in MAXIMA directory) of the several strings to be copied into MAXIMA software.

Gravity Version 14j
One more Internet reference was added
under NOTE1 (at the end of the main directory).

Gravity Version 14i
Three minor changes depending on your use.

Previous subroutine
HdQ—>STR in MAXIMA directory
published under version 14h
could not write properly if flag -3 happened to be set (in "numerical" mode) instead of cleared (i.e. in "function" mode). The required instruction (-3 CF) was added; the variable z was also put in simple quotation mark ('z') to prevent a possible, prematurate (unwished) evaluation should z "appear" (already exist) in an upper ("parent") directory.

Here is, with the 3 mentioned changes,
the new, corrected code
(as it stood in former, published version j) for
HdQ—>STR (to be found in MAXIMA directory):
\<< RCLF -3 CF FORM 'z' H 'z' Q / MAXIMA \->STR "'" "" SREPL DROP "z" "(sqrt(a^2-b^2)/b)" SREPL DROP "ATAN" "atan" SREPL DROP SWAP STOF "HdivQ:" SWAP +
\>>

Gravity Version 14h
Effective calculation of gE and gP in FREE directory.

Besides many minor changes, a very important correction was made in this version h relative to previous published versions:
In the FREE directory (when possibly setting your own parameters GM, a, w and f), the gE (g at Equator, sea level) and gP (g at Pole, sea level) were then unduly taken as fixed in the closed formulae by Li and Götze.

Please note, by the way, that always was used gE for, most common in litterature, ga, and gP for gb.

Content
Use (A below)
Output (B below)
Directories (C, D & E below)

A
For a quick, normal use

- Just introduce, in stack level 2, the latitude in D.mmss (not decimal degrees!)
- Then, in stack level 1, the height h in m (even on sea level, do not forget to put here its value 0)
- And press, in main directory GRAVITY, —>g.

B
What is calculated and saved in the main directory GRAVITY

- g according to GRS67 model, 3 values:
* saved under g2C67: according to simplified Chebyshev form, order 2;
* saved under g8C67: according to simplified Chebyshev form, order 8;
* saved under gS67: according to Somigliana's equation, with the formula using the corresponding GRS67 calculated parameters k and e².

-g according to GRS80 model, 3 values:
* saved under g2C80: according to simplified Chebyshev form, order 2;
* saved under g8C80: according to simplified Chebyshev form, order 8;
* saved under gS80: according to Somigliana's equation, with the formula using the corresponding GRS80 calculated parameters k and e².

- g according to WGS84 model, 1 value:
* only one equation/value saved under gS84: according to Somigliana's equation, with the formula using the corresponding WGS84 calculated parameters k and e².

-g in your own, FR(EE) defined ellipsoid, in which the four parameters GM, f, w and a (but not b!) were previously modified (according to your wish) and saved in FREE directory, with 3 values:
* saved under gSFR: according to Somigliana's equation, with the formula using the corresponding calculated parameters k and e²;
* saved under gSaFR: according to Somigliana's alternative equation, without using the parameters k and e²; note that, without the rounding errors, gSaFR should theorically give exactly the same result as gSFR;
* saved under gClFR: according to exact Closed formulation by Li and Götze, without having to calculate previously the values at sea level of g at Equator and at Pole, formulae that integrates also directly the height h.

Note

height factor =
1-2/a*(1+f+m-2*f*SIN(lat)^2)*h+3*(h/a)^2.
height factor = g(height=h) / g(height =0).
It is the only approximated formulae used.
As already mentioned, last result gClFR is calculated in one closed formulae, including already the height h, without having to use this approximation height_factor.

The 6 (sub)directories
- GRS67 (C1 below, key GRS67 values)
- GRS80 (C2 below, key GRS80 values )
- WGS84 (C3 below, key WGS84 values)
- FREE (C4 below, your "own" ellipsoid)
- FORM (D below, recurring formulae)
- MAXIMA (E below, high precision checking)

C1
Entering GRS67 directory (inside the main directory)

- The first variables (that appear at page 1) GM, a, w and f are fixed variables, in principle not be modified.

- You will see also xxx.STR variables.
They are STRings that represent the exact values of xxx calculated with full 100 digits in MAXIMA program.
Example: b.STR.
Press b and you will see its numerical value.
But 'b' RCL will show how it was calculated (b=a-a*f).
(Note: f.STR is here in GRS67 just the result of f —>STR , with, of course, no 100 digits.)
For more details on high precision calculation, see also letter E below relative to MAXIMA directory.

- You will see also (in directory GRS67):

* gE: normally it corresponds to the numerical value of the variable gE.1, always at Equator
In the program or subroutines only gE (or gP) is used (and never forms like "xxx.1", "xxx.2" or "xxx.3").
* gE.1: is supposed to be the best resulting value (to 12 digits) to be employed for gE.
* gE.2 and gE.3 are alternative values found in the litterature.
Therefore, if you see, in a paper,
gE (GRS67) = 9.780318,
and want to check the value at lat=45/h=1000 according to that paper,
do here, previously, in this directory GRS67:
9.780318 ENTER (or gE.3 ENTER) 'gE' STO,
before launching, in the main directory, 45 0 —>g.
The best solution/result should be, however, to leave gE in the directory set = gE.1
(do for that: gE.1 ENTER 'gE' STO).

Next to gE.3 (in directory GRS67):

* gE.STR : it is the string representation of the best value of gE calculated (normally to 100 digits precidion by MAXIMA software, see further letter E) at Equator, according to exact formulae given here by program gE.CAL;
* gE.CAL: calculates the theorical "exact" value of gE at sea level, but, due to rounding errors, the resulting value shows discrepancies with the one calculated with the very same formulae in MAXIMA software (the latter appearing under gE.STR as mentioned above).

Other variables (in directory GRS67):

* Press for instance J2.CAL to better see the distinction between the given results for J2.
* Beta.i and alpha.i (i=2, 8) are the Chebyshev's coefficients used internally to calculated g.
* NOTE variables : some observations or explanations.

C2
Entering WGS84 directory (inside the main directory)

The same observations as for GRS67 directory apply here.

Just a point to be noted:
Fixed parameter w: the correct value is the one under the label w1 (3.986004418E14), but you may find in the litterature the value w2 (3.986004418E14). It's up to you to adopt the value to be used for w (not recommanded, though) doing w2 ENTER 'w' STO (to get [back] the initial correct value of w, do: w1 ENTER 'w' STO).

C3
Entering GRS80 directory (inside the main directory)

More or less the same observations as for GRS67 directory apply here.

Note at page 1 of this directory that, in this GRS67 Earth model, it is the form factor J2 that is considered fixed, and no more f, the flattening variable.
Consequently, f becomes a derived variable/value.
To calculate f, press f.CAL.
The exact formulae used here, as always, give — with the internal solver of the HP50G — somewhat inaccurate, calculated results due to rounding errors.
To view the exact 100 digits, press f.STR.
For full 100 (or more) digits exact calculation in MAXIMA program, see further down (letter E).

C4
Entering FREE directory (inside the main directory)

You may completely ignore this section.
Only useful if you want to use your "own" ellipsoid.
If necessary, change here in this directory named FREE the four values of GM, a, f and w.
Do not modify anything else (neither b, automatically calculated by b=a- a*f).
gE and gP (at h=0) will also be calculated automatically.
J2.1: the exact formulation to calculate J2.
J2.2: a simplified (unnecessary?) version (less accurate) to calculate J2.
k and e2 (and k.STR, e2.STR and f.STR) automatically calculated to be used internally, in particular in Somigliana's equation.

When you are in this FREE directory, you can launch directly (from here) your request for g FREE:
lat (in d.mmss) in stack level 2
height h (in m) in stack level 1
—>gFREE.

Two outputs of your FREE ellipsoid (saved under GSaFR, Somigliana alternative formulae for your FREE ellipsoid; and gClFR, Closed formulae by Li and Götze for your FREE ellipsoid) are then given.

The simpler, however,
is to use only the main program —>g in the main directory.

Please note that, if the four fixed parameters GM, a, w and f in FREE direcory are exactly the same as in WGS84 directory, the results will appear with an ending tag "WGS84" instead of "FREE".

D
Directory FORM (inside the main directory)

Place/Directory where some recurrent formulae are located.

E
Directory MAXIMA (inside the main directory)

You may completely ignore this section.

Should you need more free space (memory), the whole directly may even be suppressed.
How?
In the main directory GRAVITY, write 'MAXIMA' (with simple quotation marks), then PGDIR.

Directory MAXIMA is named after the eponymous MAXIMA program, a free, powerful software, enabling for example to calculate "complicate" algebraic formulae with variables and with as many exact digits as requested.

With the help of the instructions inside this HP MAXIMA directory together with the MAXIMA eponymous software, you can check some key results given in the different directories of this HP gravity program under the labels xxx.STR.
Alternatively, you can calculate, with very high accuracy, values of g(lat, h), gMean, J2, for any ellipsoid — and f, when J2 is fixed like in GRS80 instead of f like in WGS84.

How to use HP subroutines in MAXIMA directory to calculate — to 100 or more digits precision — gE, gP, gMEAN and other key values (at sea level & not considering the latitude) with the free MAXIMA software?

Step 1
Enter the MAXIMA directory.

Step 2
Then, with no argument, press GRS67
or press GRS80
or press WGS84.

Step 3
Copy the resulting string in stack from HP49-HP50 calculator into MAXIMA free software command line.

Furthermore, thanks to "concatenating" subroutine FREE in HP MAXIMA directory, it is also possible to reproduce easily, in MAXIMA free software, the succession of Li & Götze formulae and get the exact value (to 100 or more digits) of the theorical normal Earth gravity value g(lat, height), i.e. including the latitude & height h, following the 2-steps instructions below:

Step 1
* Enter the MAXIMA directory.
* Write the latitude in D.mmss.
* ENTER
* Write the height h in m above the ellipsoid.
* ENTER
* Then press FREE.

Step 2
Copy the resulting string in stack from HP49-HP50 calculator into MAXIMA software command line.

Alternatively to the 1-2 steps above, procede as follows with the steps 1-14 should you lack of space/memory in your HP calculator (using 3 times the subroutine FR3ST, instead of a single time the subroutine FREE):

Step 1
* Enter the MAXIMA directory.
* 0 0 (2 times the number 0)
* Then press FR3ST.

Step 2
Copy the resulting string in stack from HP49-HP50 calculator into MAXIMA software command line.

Step 3
* Back into Maxima directory
* Press gE (at page 2).

Step 4
Copy the resulting string in stack from HP49-HP50 calculator into MAXIMA software command line.

Step 5
* Back into Maxima directory
* 90 0 (90 & 0)
* Then press FR3ST.

Step 6
Copy the resulting string in stack from HP49-HP50 calculator into MAXIMA software command line.

Step 7
* Back into Maxima directory
* Press gP (at page 2).

Step 8
Copy the resulting string in stack from HP49-HP50 calculator into MAXIMA software command line.

Step 9
* Write the latitude in D.mmss.
* ENTER
* Write the height h in m above the ellipsoid.
* ENTER
* Then press FR3ST.

Step 10
Copy the resulting string in stack from HP49-HP50 calculator into MAXIMA software command line.

Step 11
* Back into Maxima directory
* Press k (at page 2).

Step 12
Copy the resulting string in stack from HP49-HP50 calculator into MAXIMA software command line.

Step 13
* Back into Maxima directory
* Press gMU—>STR (next to k at page 2).

Step 14
Copy the resulting string in stack from HP49-HP50 calculator into MAXIMA software command line.

My thanks here again to Charles Karney for his answers and help.

Remarks are welcome.

Gil Campart

[Gil]

HP49-50G:Theorical Earth gravity g = g(latitude, height), WGS84, GRS80/67

Theorical Normal Earth Gravity g
for HP49-HP50

Gravity version 16c
During checking process of variables accuracy of type .STR in version 16b, extraneous variables were introduced in the directories. Now they were duly suppressed.

Gravity version 16b
Precision for MAXIMA software set to 100 digits by default, instead of 150.

Gravity version 16
Double-checking and correction of the string variables .STR when 100 digits are shown after 150 digits precision calculation with MAXIMA free software.

Note 1:
“x1 x2... x99 7(501)“ will be registered as "x1 x2... x99 8“;
“x1 x2... x99 7(499)“ will be registered as "x1 x2... x99 7“.

Note 2:
Suppose gE & gP are shown and registered with exactly 100 correct digits, then further calculation with them like
k = b*gP/(a*gE)-1 will be "wrong" in their last digits. Therefore, when a 100 digits precision is required, play it safe, increasing the fpprecision value (for example, setting fpprec: 150 under gP in MAXIMA directory).

Gravity version 15 (then incorrectly labelled 15g)
Many changes, among others regarding the flags status (to be always set back at its initial stage), the variable J2 (and references to it) and, above all, the MAXIMA directory, completely reformulated. For the latter:
- full automatisation (with, inside MAXIMA dir, single subroutine GRS67 or GRS80 or WGS84 or FREE) for the string to be copied in one step into MAXIMA software;
- in case of calculator lack of space/memory, semi-full automatisation (with, instead of subroutine FREE in MAXIMA directory, subroutine FR3ST in MAXIMA directory) of the several strings to be copied into MAXIMA software.

Gravity Version 14j
One more Internet reference was added
under NOTE1 (at the end of the main directory).

Gravity Version 14i
Three minor changes depending on your use.

Previous subroutine
HdQ—>STR in MAXIMA directory
published under version 14h
could not write properly if flag -3 happened to be set (in "numerical" mode) instead of cleared (i.e. in "function" mode). The required instruction (-3 CF) was added; the variable z was also put in simple quotation mark ('z') to prevent a possible, prematurate (unwished) evaluation should z "appear" (already exist) in an upper ("parent") directory.

Here is, with the 3 mentioned changes,
the new, corrected code
(as it stood in former, published version j) for
HdQ—>STR (to be found in MAXIMA directory):
\<< RCLF -3 CF FORM 'z' H 'z' Q / MAXIMA \->STR "'" "" SREPL DROP "z" "(sqrt(a^2-b^2)/b)" SREPL DROP "ATAN" "atan" SREPL DROP SWAP STOF "HdivQ:" SWAP +
\>>

Gravity Version 14h
Effective calculation of gE and gP in FREE directory.

Besides many minor changes, a very important correction was made in this version h relative to previous published versions:
In the FREE directory (when possibly setting your own parameters GM, a, w and f), the gE (g at Equator, sea level) and gP (g at Pole, sea level) were then unduly taken as fixed in the closed formulae by Li and Götze.

Please note, by the way, that always was used gE for, most common in litterature, ga, and gP for gb.

Content
Use (A below)
Output (B below)
Directories (C, D & E below)

A
For a quick, normal use

- Just introduce, in stack level 2, the latitude in D.mmss (not decimal degrees!)
- Then, in stack level 1, the height h in m (even on sea level, do not forget to put here its value 0)
- And press, in main directory GRAVITY, —>g.

B
What is calculated and saved in the main directory GRAVITY

- g according to GRS67 model, 3 values:
* saved under g2C67: according to simplified Chebyshev form, order 2;
* saved under g8C67: according to simplified Chebyshev form, order 8;
* saved under gS67: according to Somigliana's equation, with the formula using the corresponding GRS67 calculated parameters k and e².

-g according to GRS80 model, 3 values:
* saved under g2C80: according to simplified Chebyshev form, order 2;
* saved under g8C80: according to simplified Chebyshev form, order 8;
* saved under gS80: according to Somigliana's equation, with the formula using the corresponding GRS80 calculated parameters k and e².

- g according to WGS84 model, 1 value:
* only one equation/value saved under gS84: according to Somigliana's equation, with the formula using the corresponding WGS84 calculated parameters k and e².

-g in your own, FR(EE) defined ellipsoid, in which the four parameters GM, f, w and a (but not b!) were previously modified (according to your wish) and saved in FREE directory, with 3 values:
* saved under gSFR: according to Somigliana's equation, with the formula using the corresponding calculated parameters k and e²;
* saved under gSaFR: according to Somigliana's alternative equation, without using the parameters k and e²; note that, without the rounding errors, gSaFR should theorically give exactly the same result as gSFR;
* saved under gClFR: according to exact Closed formulation by Li and Götze, without having to calculate previously the values at sea level of g at Equator and at Pole, formulae that integrates also directly the height h.

Note

height factor =
1-2/a*(1+f+m-2*f*SIN(lat)^2)*h+3*(h/a)^2.
height factor = g(height=h) / g(height =0).
It is the only approximated formulae used.
As already mentioned, last result gClFR is calculated in one closed formulae, including already the height h, without having to use this approximation height_factor.

The 6 (sub)directories
- GRS67 (C1 below, key GRS67 values)
- GRS80 (C2 below, key GRS80 values )
- WGS84 (C3 below, key WGS84 values)
- FREE (C4 below, your "own" ellipsoid)
- FORM (D below, recurring formulae)
- MAXIMA (E below, high precision checking)

C1
Entering GRS67 directory (inside the main directory)

- The first variables (that appear at page 1) GM, a, w and f are fixed variables, in principle not be modified.

- You will see also xxx.STR variables.
They are STRings that represent the exact values of xxx calculated with full 100 digits in MAXIMA program.
Example: b.STR.
Press b and you will see its numerical value.
But 'b' RCL will show how it was calculated (b=a-a*f).
(Note: f.STR is here in GRS67 just the result of f —>STR , with, of course, no 100 digits.)
For more details on high precision calculation, see also letter E below relative to MAXIMA directory.

- You will see also (in directory GRS67):

* gE: normally it corresponds to the numerical value of the variable gE.1, always at Equator
In the program or subroutines only gE (or gP) is used (and never forms like "xxx.1", "xxx.2" or "xxx.3").
* gE.1: is supposed to be the best resulting value (to 12 digits) to be employed for gE.
* gE.2 and gE.3 are alternative values found in the litterature.
Therefore, if you see, in a paper,
gE (GRS67) = 9.780318,
and want to check the value at lat=45/h=1000 according to that paper,
do here, previously, in this directory GRS67:
9.780318 ENTER (or gE.3 ENTER) 'gE' STO,
before launching, in the main directory, 45 0 —>g.
The best solution/result should be, however, to leave gE in the directory set = gE.1
(do for that: gE.1 ENTER 'gE' STO).

Next to gE.3 (in directory GRS67):

* gE.STR : it is the string representation of the best value of gE calculated (normally to 100 digits precidion by MAXIMA software, see further letter E) at Equator, according to exact formulae given here by program gE.CAL;
* gE.CAL: calculates the theorical "exact" value of gE at sea level, but, due to rounding errors, the resulting value shows discrepancies with the one calculated with the very same formulae in MAXIMA software (the latter appearing under gE.STR as mentioned above).

Other variables (in directory GRS67):

* Press for instance J2.CAL to better see the distinction between the given results for J2.
* Beta.i and alpha.i (i=2, 8) are the Chebyshev's coefficients used internally to calculated g.
* NOTE variables : some observations or explanations.

C2
Entering WGS84 directory (inside the main directory)

The same observations as for GRS67 directory apply here.

Just a point to be noted:
Fixed parameter w: the correct value is the one under the label w1 (3.986004418E14), but you may find in the litterature the value w2 (3.986004418E14). It's up to you to adopt the value to be used for w (not recommanded, though) doing w2 ENTER 'w' STO (to get [back] the initial correct value of w, do: w1 ENTER 'w' STO).

C3
Entering GRS80 directory (inside the main directory)

More or less the same observations as for GRS67 directory apply here.

Note at page 1 of this directory that, in this GRS67 Earth model, it is the form factor J2 that is considered fixed, and no more f, the flattening variable.
Consequently, f becomes a derived variable/value.
To calculate f, press f.CAL.
The exact formulae used here, as always, give — with the internal solver of the HP50G — somewhat inaccurate, calculated results due to rounding errors.
To view the exact 100 digits, press f.STR.
For full 100 (or more) digits exact calculation in MAXIMA program, see further down (letter E).

C4
Entering FREE directory (inside the main directory)

You may completely ignore this section.
Only useful if you want to use your "own" ellipsoid.
If necessary, change here in this directory named FREE the four values of GM, a, f and w.
Do not modify anything else (neither b, automatically calculated by b=a- a*f).
gE and gP (at h=0) will also be calculated automatically.
J2.1: the exact formulation to calculate J2.
J2.2: a simplified (unnecessary?) version (less accurate) to calculate J2.
k and e2 (and k.STR, e2.STR and f.STR) automatically calculated to be used internally, in particular in Somigliana's equation.

When you are in this FREE directory, you can launch directly (from here) your request for g FREE:
lat (in d.mmss) in stack level 2
height h (in m) in stack level 1
—>gFREE.

Two outputs of your FREE ellipsoid (saved under GSaFR, Somigliana alternative formulae for your FREE ellipsoid; and gClFR, Closed formulae by Li and Götze for your FREE ellipsoid) are then given.

The simpler, however,
is to use only the main program —>g in the main directory.

Please note that, if the four fixed parameters GM, a, w and f in FREE direcory are exactly the same as in WGS84 directory, the results will appear with an ending tag "WGS84" instead of "FREE".

D
Directory FORM (inside the main directory)

Place/Directory where some recurrent formulae are located.

E
Directory MAXIMA (inside the main directory)

You may completely ignore this section.

Should you need more free space (memory), the whole directly may even be suppressed.
How?
In the main directory GRAVITY, write 'MAXIMA' (with simple quotation marks), then PGDIR.

Directory MAXIMA is named after the eponymous MAXIMA program, a free, powerful software, enabling for example to calculate "complicate" algebraic formulae with variables and with as many exact digits as requested.

With the help of the instructions inside this HP MAXIMA directory together with the MAXIMA eponymous software, you can check some key results given in the different directories of this HP gravity program under the labels xxx.STR.
Alternatively, you can calculate, with very high accuracy, values of g(lat, h), gMean, J2, for any ellipsoid — and f, when J2 is fixed like in GRS80 instead of f like in WGS84.

How to use HP subroutines in MAXIMA directory to calculate — to 100 or more digits precision — gE, gP, gMEAN and other key values (at sea level & not considering the latitude) with the free MAXIMA software?

Step 1
Enter the MAXIMA directory.

Step 2
Then, with no argument, press GRS67
or press GRS80
or press WGS84.

Step 3
Copy the resulting string in stack from HP49-HP50 calculator into MAXIMA free software command line.

Furthermore, thanks to "concatenating" subroutine FREE in HP MAXIMA directory, it is also possible to reproduce easily, in MAXIMA free software, the succession of Li & Götze formulae and get the exact value (to 100 or more digits) of the theorical normal Earth gravity value g(lat, height), i.e. including the latitude & height h, following the 2-steps instructions below:

Step 1
* Enter the MAXIMA directory.
* Write the latitude in D.mmss.
* ENTER
* Write the height h in m above the ellipsoid.
* ENTER
* Then press FREE.

Step 2
Copy the resulting string in stack from HP49-HP50 calculator into MAXIMA software command line.

Alternatively to the 1-2 steps above, procede as follows with the steps 1-14 should you lack of space/memory in your HP calculator (using 3 times the subroutine FR3ST, instead of a single time the subroutine FREE):

Step 1
* Enter the MAXIMA directory.
* 0 0 (2 times the number 0)
* Then press FR3ST.

Step 2
Copy the resulting string in stack from HP49-HP50 calculator into MAXIMA software command line.

Step 3
* Back into Maxima directory
* Press gE (at page 2).

Step 4
Copy the resulting string in stack from HP49-HP50 calculator into MAXIMA software command line.

Step 5
* Back into Maxima directory
* 90 0 (90 & 0)
* Then press FR3ST.

Step 6
Copy the resulting string in stack from HP49-HP50 calculator into MAXIMA software command line.

Step 7
* Back into Maxima directory
* Press gP (at page 2).

Step 8
Copy the resulting string in stack from HP49-HP50 calculator into MAXIMA software command line.

Step 9
* Write the latitude in D.mmss.
* ENTER
* Write the height h in m above the ellipsoid.
* ENTER
* Then press FR3ST.

Step 10
Copy the resulting string in stack from HP49-HP50 calculator into MAXIMA software command line.

Step 11
* Back into Maxima directory
* Press k (at page 2).

Step 12
Copy the resulting string in stack from HP49-HP50 calculator into MAXIMA software command line.

Step 13
* Back into Maxima directory
* Press gMU—>STR (next to k at page 2).

Step 14
Copy the resulting string in stack from HP49-HP50 calculator into MAXIMA software command line.

My thanks here again to Charles Karney for his answers and help.

Remarks are welcome.

Gil Campart

[Gil]

Note, please, that calculations here with Li & Götze closed formulae do not match (could not match) the calculations by Charles Karney with full theory according to Heiskanen & Moritz.

WGS84
According to Karney:
"gammay and gammaz are the
northerly and up components of gamma (up is in the direction of the normal to the WGS84 ellipsoid). gamma is the magnitude of the gravity.

As expected these results differ from L+G's results for
gamma_u (because this is neither the up component of the gravity nor its magnitude).

lat h gammay gammaz gamma
20 10000 -0.00005231628077462272 -9.75556773962467229537
9.75556773976495081667

lat h gammay gammaz gamma
80 5000 -0.00001391201648793046 -9.81521493778272672782
9.81521493779258612489."

To be compared with Li & Götze in Maxima
(& with HP49-50G):
9.755567739391774710519363610080945073955567777773578831194488706435195066202777​754390572703114008923
(with HP50G: 9.7555677283)

Li & Götze in Maxima (& with HP49-50G):
9.815214937766250073333365482033892183981119210801868284213419088307958628663215​851806608387826859838b0
(9.81521495521)

[Gil]

Abandoned everywhere the closed notation for g calculation by Li and Götze in favor of more commun and clearer notation by Heiskanen and Moritz.

Changed the name CLOSED form, as what was calculated by Li and Götze when at height h above ellipsoid — with then abs(lat)≠90 or 0 — was not the normal g component, but guFR(EE).

For completeness, added, beside FREE variable guFR, the gbetaFR component, following the terminology of Heiskanen and Moritz.

The Maxima string (one huge set/string of instructions in HP) for checking — with the high digits precision — in MAXIMA free software the results was adapted consequently.

[Gil]

Version 17e

Here was unified the use of gravities guFR and gbetaFR with the ending FR (like FREE ellipsoid), instead of sometimes only gu and gbeta.

So that gMAGNIT is now correctly calculated.

[Gil]

Gravity14b.hp

with .hp ending.