HP4950G:Theorical Earth gravity g = g(latitude, height), WGS84, GRS80/67

09212021, 07:04 AM
(This post was last modified: 10032021 07:38 AM by Gil.)
Post: #1




HP4950G:Theorical 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. 

« Next Oldest  Next Newest »

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