(49g 50g) Theoretical Earth gravity g = g(latitude, height), WGS84, GRS80/67
|
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 |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 1 Guest(s)