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

09212021, 07:04 AM
(This post was last modified: 01122023 01:45 PM by Gil.)
Post: #1




(49g 50g) Theoretical Earth gravity g = g(latitude, height), WGS84, GRS80/67
HP4950G
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.9318526464E3*SIN(lat)^2)/\v/(16.69437999014E3*SIN(lat)^2))(1.00139*SIN(lat)^2)*.0000030877*alt+7.2E13*alt^2' \>NUM "WGS 84" \>TAG \>> \>> Regards, Gil Campart 

09212021, 11:02 AM
Post: #2




RE: HP4950G : —>g gravity calculation = g(latitude, height) with WGS84
(09212021 07:04 AM)Gil Wrote: See, for example, for more details: It is amazing that formula has so many significant digits But, constants only matched 6significant digits from https://en.wikipedia.org/wiki/Gravity_of...tude_model 

09212021, 01:31 PM
Post: #3




RE: HP4950G : —>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 gformulae 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 nonintegers), 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). 

09242021, 12:02 AM
Post: #4




RE: HP4950G : —>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.9318526464E3*SIN(lat)^2)/\v/(16.69437999014E3*SIN(lat)^2))(1.00139*SIN(lat)^2)*.0000030877*alt+7.2E13*alt^2' \>NUM "WGS 84" \>TAG f STOF \>> \>> Regards, Gil 

09242021, 12:34 PM
Post: #5




RE: HP4950G : —>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 . 

09252021, 01:27 PM
(This post was last modified: 09252021 01:29 PM by Gil.)
Post: #6




RE: HP4950G : —>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(lat)^6+.0000000007*SIN(lat)^8)(1.00139*SIN(lat)^2)*.0000030877*alt+7.2E13*alt^2' \>NUM "Somigliana GRS 80" \>TAG '9.7803253359*((1+1.9318526464E3*SIN(lat)^2)/\v/(16.69437999014E3*SIN(lat)^2))(1.00139*SIN(lat)^2)*.0000030877*alt+7.2E13*alt^2' \>NUM "Somigliana WGS 84" \>TAG f STOF \>> \>> 

09262021, 09:25 PM
(This post was last modified: 09262021 09:31 PM by Gil.)
Post: #7




RE: HP4950G : —>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.indubioprogeo.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(lat)^6+.0000000007*SIN(lat)^8)' \>NUM lat alt 6378137 3.35281068118E3 3.44978600308E3 \> lat h a f m \<< '12/a*(1+f+m2*f*SIN(lat)^2)*h+3*(h/a)^2' * \>NUM \>> "Somigliana GRS 80" \>TAG '9.7803253359*((1+1.9318526464E3*SIN(lat)^2)/\v/(16.69437999014E3*SIN(lat)^2))' \>NUM lat alt 6378137 298.257223563 INV 3.44978650684E3 \> lat h a f m \<< '12/a*(1+f+m2*f*SIN(lat)^2)*h+3*(h/a)^2' * \>NUM \>> "Somigliana WGS 84" \>TAG f STOF \>> \>> Regards, Gil 

09282021, 03:19 PM
Post: #8




RE: HP4950G : —>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)*(1mm/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/é²)*(11/é²*ATAN(é²))1' .00268804118 q0 = '((1+3/é²)*ATAN(é²)3/é²)/2' .00007334625 é² = '(a*ab*b)/(b*b)' 6.73949674208E3 GM = 3.986004418E14 m = 'w*w*a*a*b/GM' 3.44978650683E3 w = .00007292115 a 6378137 a = 6378137 b = 'aa/298.257223563' 6356752.31425 The differences are now quite small, due to the roundings of the calculator. Regards, Gil 

10022021, 05:55 PM
(This post was last modified: 10032021 07:33 AM by Gil.)
Post: #9




HP4950G: 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 GSM84Values 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 

10062021, 10:14 PM
(This post was last modified: 10102021 02:50 PM by Gil.)
Post: #10




RE: HP4950G:Theorical Earth gravity g = g(latitude, height), WGS84, GRS80/67
RE: HP4950G:
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.3024E3 ß1.1 5.9E6 show inaccurateness Taylor new calc. Val ß.2 ~ 5.3022E3 ß1.2 ~ 5.8E6 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= aa×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 12digits 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.08269887351E3 : Published J2: .0010827 "CODE J2.CALC" "\\<< RCLF RAD f.STR OBJ\\> \\>NUM \\> f '(a^2(aa*f)^2)/a^2/3*(12/15*(w^2*a^2*(aa*f)/GM)*(\\v/((a^2(aa*f)^2)/(aa*f)^2)/(1/2*((1+3*(aa*f)^2/(a^2(aa*f)^2))*ATAN(\\v/(a^2(aa*f)^2)/(aa*f))3*(aa*f)/\\v/(a^2(aa*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(aa*f.SOL)^2)/a^2/3*(12/15*(w^2*a^2*(aa*f.SOL)/GM)*(\\v/((a^2(aa*f.SOL)^2)/(aa*f.SOL)^2)/(1/2*((1+3*(aa*f.SOL)^2/(a^2(aa*f.SOL)^2))*ATAN(\\v/(a^2(aa*f.SOL)^2)/(aa*f.SOL))3*(aa*f.SOL)/\\v/(a^2(aa*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.08262891487E3 : Published J2.STR: "0.0010826298213129219" "CODE J2.CALC" "\\<< RCLF RAD f.STR OBJ\\> \\>NUM \\> f '(a^2(aa*f)^2)/a^2/3*(12/15*(w^2*a^2*(aa*f)/GM)*(\\v/((a^2(aa*f)^2)/(aa*f)^2)/(1/2*((1+3*(aa*f)^2/(a^2(aa*f)^2))*ATAN(\\v/(a^2(aa*f)^2)/(aa*f))3*(aa*f)/\\v/(a^2(aa*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 '12/a*(1+f+m2*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/(1e\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 \<< 'af*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^2z\180^2' \>NUM 'd\180\180\178' STO 'r\180^2+z\180^2' \>NUM 'r\180\180\178' STO '\v/(a^2b^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/4D/2' \>NUM 0 MAX \v/  0 MAX \v/ 1 MIN ACOS '\Gb\180' STO '\v/(r\180\180\178E^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)*(1b\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)^21/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^2b^2)/a^2' \>NUM 'e\178' STO '(a^2b^2)/b^2' \>NUM '\233\178' STO 'w^2*a^2*b/GM' \>NUM 'm' STO 'e\178/3*(12/15*m*(\v/\233\178/q0))' \>NUM 'J2' STO \>> \>> Regards, Gil 

10112021, 09:08 AM
Post: #11




RE: HP4950G: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 

10172021, 04:22 PM
(This post was last modified: 10192021 07:42 PM by Gil.)
Post: #12




RE: HP4950G: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 selfexplaining NOTES at the end of each directory. 

10292021, 02:37 PM
(This post was last modified: 11042021 06:47 PM by Gil.)
Post: #13




HP4950G:Theorical Earth gravity g = g(latitude, height), WGS84, GRS80/67
Theorical Normal Earth Gravity g
for HP49HP50 Gravity version 16 Doublechecking 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, semifull 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^2b^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 = 12/a*(1+f+m2*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=aa*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 HP49HP50 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 2steps 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 HP49HP50 calculator into MAXIMA software command line. Alternatively to the 12 steps above, procede as follows with the steps 114 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 HP49HP50 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 HP49HP50 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 HP49HP50 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 HP49HP50 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 HP49HP50 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 HP49HP50 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 HP49HP50 calculator into MAXIMA software command line. My thanks here again to Charles Karney for his answers and help. Remarks are welcome. Gil Campart 

11052021, 09:31 AM
(This post was last modified: 11052021 12:20 PM by Gil.)
Post: #14




RE: HP4950G:Theorical Earth gravity g = g(latitude, height), WGS84, GRS80/67
HP4950G:Theorical Earth gravity g = g(latitude, height), WGS84, GRS80/67
Theorical Normal Earth Gravity g for HP49HP50 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 Doublechecking 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, semifull 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^2b^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 = 12/a*(1+f+m2*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=aa*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 HP49HP50 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 2steps 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 HP49HP50 calculator into MAXIMA software command line. Alternatively to the 12 steps above, procede as follows with the steps 114 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 HP49HP50 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 HP49HP50 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 HP49HP50 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 HP49HP50 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 HP49HP50 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 HP49HP50 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 HP49HP50 calculator into MAXIMA software command line. My thanks here again to Charles Karney for his answers and help. Remarks are welcome. Gil Campart 

11092021, 04:31 PM
Post: #15




RE: HP4950G:Theorical Earth gravity g = g(latitude, height), WGS84, GRS80/67
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 HP4950G): 9.755567739391774710519363610080945073955567777773578831194488706435195066202777754390572703114008923 (with HP50G: 9.7555677283) Li & Götze in Maxima (& with HP4950G): 9.815214937766250073333365482033892183981119210801868284213419088307958628663215851806608387826859838b0 (9.81521495521) 

11122021, 02:06 PM
Post: #16




RE: HP4950G:Theorical Earth gravity g = g(latitude, height), WGS84, GRS80/67
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. 

11152021, 03:42 PM
Post: #17




RE: HP4950G:Theorical Earth gravity g = g(latitude, height), WGS84, GRS80/67
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. 

01122023, 01:47 PM
Post: #18




RE: (49g 50g) Theoretical Earth gravity g = g(latitude, height), WGS84, GRS80/67
Gravity14b.hp
with .hp ending. 

« Next Oldest  Next Newest »

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