Post Reply 
HP49-50G:Theorical Earth gravity g = g(latitude, height), WGS84, GRS80/67
09-21-2021, 07:04 AM (This post was last modified: 10-03-2021 07:38 AM by Gil.)
Post: #1
HP49-50G:Theorical Earth gravity g = g(latitude, height), WGS84, GRS80/67
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


Attached File(s)
.doc  _G.DOC (Size: 553 bytes / Downloads: 2)
Find all posts by this user
Quote this message in a reply
09-21-2021, 11:02 AM
Post: #2
RE: HP49-50G : —>g gravity calculation = g(latitude, height) with WGS84
(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 Big Grin

[Image: d50ca6cb8c65f6eab2ac8e500df93bbf0dd0c76d]

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

[Image: 3f0a721aede8ed3761d1b68c5bab94116cb0ed6f]
Find all posts by this user
Quote this message in a reply
09-21-2021, 01:31 PM
Post: #3
RE: HP49-50G : —>g gravity calculation = g(latitude, height) with WGS84
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).
Find all posts by this user
Quote this message in a reply
09-24-2021, 12:02 AM
Post: #4
RE: HP49-50G : —>g gravity calculation = g(latitude, height) with WGS84
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


Attached File(s)
.doc  _G.vet1b.DOC (Size: 572 bytes / Downloads: 0)
Find all posts by this user
Quote this message in a reply
09-24-2021, 12:34 PM
Post: #5
RE: HP49-50G : —>g gravity calculation = g(latitude, height) with WGS84
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
.
Find all posts by this user
Quote this message in a reply
09-25-2021, 01:27 PM (This post was last modified: 09-25-2021 01:29 PM by Gil.)
Post: #6
RE: HP49-50G : —>g gravity calculation = g(latitude, height) with WGS84
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
\>>
\>>


Attached File(s)
.doc  _G.ver2b.DOC (Size: 981 bytes / Downloads: 1)
Find all posts by this user
Quote this message in a reply
09-26-2021, 09:25 PM (This post was last modified: 09-26-2021 09:31 PM by Gil.)
Post: #7
RE: HP49-50G : —>g gravity calculation = g(latitude, height) with WGS84
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


Attached File(s)
.doc  _G.ver3.DOC (Size: 1.2 KB / Downloads: 1)
Find all posts by this user
Quote this message in a reply
09-28-2021, 03:19 PM
Post: #8
RE: HP49-50G : —>g gravity calculation = g(latitude, height) with WGS84
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
Find all posts by this user
Quote this message in a reply
10-02-2021, 05:55 PM (This post was last modified: 10-03-2021 07:33 AM by Gil.)
Post: #9
HP49-50G: Theorical Earth gravity g = g(latitude, height) WGS84 GRS80/67
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 :
Somigliana GRS 80: 9.7803267715 :Somigliana WGS 84: 9.7803253359 :Closed Form WGS 84: 9.78032532324

:alt [m]: 1000 :lat D.mmss: 0 :
Int Grav GRS 67: 9.777232 :
Somigliana GRS 80: 9.77723980166 :Somigliana WGS 84: 9.77723836651 :Closed Form WGS 84: 9.77723826177

:alt [m]: 0 :lat D.mmss: 90 :
Int Grav GRS 67: 9.83217715816 :
Somigliana GRS 80: 9.83218636846 :Somigliana WGS 84: 9.83218493787 :Closed Form WGS 84: 9.83218496308

:alt [m]: 1000 :lat D.mmss: 90 :
Int Grav GRS 67: 9.82909115816 :
Somigliana GRS 80: 9.82910370419 :Somigliana WGS 84: 9.82910227404 :Closed Form WGS 84: 9.82910231835

:alt [m]: 0 :lat D.mmss: 45 :
Int Grav GRS 67: 9.80618987521 :
Somigliana GRS 80: 9.80619920255 :Somigliana WGS 84: 9.80619776931 :Closed Form WGS 84: 9.80619777838

:alt [m]: 1000 :lat D.mmss: 45 :
Int Grav GRS 67: 9.80310387521 :
Somigliana GRS 80: 9.80311437628 :Somigliana 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


Attached File(s)
.doc  GRAVITY.Vers6e.Doc (Size: 6.08 KB / Downloads: 1)
Find all posts by this user
Quote this message in a reply
10-06-2021, 10:14 PM (This post was last modified: 10-10-2021 02:50 PM by Gil.)
Post: #10
RE: HP49-50G:Theorical Earth gravity g = g(latitude, height), WGS84, GRS80/67
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.3024E-3
ß1.1 -5.9E-6
show inaccurateness

Taylor new calc. Val
ß.2 ~ 5.3022E-3
ß1.2 ~ -5.8E-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 "Somigliana GRS 80" \->TAG gS84 "Somigliana 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


Attached File(s)
.doc  GRAVITY.Vers9i.Doc (Size: 11.27 KB / Downloads: 0)
.doc  GRAVITY.Vers9e.Doc (Size: 10.92 KB / Downloads: 0)
.doc  GRAVITY.Vers9k.Doc (Size: 11.7 KB / Downloads: 0)
Find all posts by this user
Quote this message in a reply
10-11-2021, 09:08 AM
Post: #11
RE: HP49-50G:Theorical Earth gravity g = g(latitude, height), WGS84, GRS80/67
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


Attached File(s)
.doc  GRAVITY.Vers10.Doc (Size: 12.18 KB / Downloads: 1)
Find all posts by this user
Quote this message in a reply
10-17-2021, 04:22 PM (This post was last modified: 10-19-2021 07:42 PM by Gil.)
Post: #12
RE: HP49-50G:Theorical Earth gravity g = g(latitude, height), WGS84, GRS80/67
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: g2C67: gravity result
g2C67: 2nd order appro
g2C67: Chebyshev
g2C67: GRS67
g8C80: 8th order appro
g8C80: GRS80 ±S Smodif
gS84 : Somigliana
gS84 : WGS84 gµ:gMean
gFR(EE):gyour 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.


Attached File(s)
.doc  Gravity.Ver11f.DOC (Size: 14.74 KB / Downloads: 3)
.doc  Gravity.Ver12b.DOC (Size: 14.56 KB / Downloads: 1)
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)