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

[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