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

[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