GeographicLib
1.21
|
00001 /** 00002 * \file NormalGravity.cpp 00003 * \brief Implementation for GeographicLib::NormalGravity class 00004 * 00005 * Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed under 00006 * the MIT/X11 License. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #include <GeographicLib/NormalGravity.hpp> 00011 00012 #define GEOGRAPHICLIB_NORMALGRAVITY_CPP \ 00013 "$Id: ccd3d6bd4bddefb0b40cbce0006d863f08d64db4 $" 00014 00015 RCSID_DECL(GEOGRAPHICLIB_NORMALGRAVITY_CPP) 00016 RCSID_DECL(GEOGRAPHICLIB_NORMALGRAVITY_HPP) 00017 00018 namespace GeographicLib { 00019 00020 using namespace std; 00021 00022 NormalGravity::NormalGravity(real a, real GM, real omega, real f, real J2) 00023 : _a(a) 00024 , _GM(GM) 00025 , _omega(omega) 00026 , _f(f) 00027 , _J2(J2) 00028 , _omega2(Math::sq(_omega)) 00029 , _aomega2(Math::sq(_omega * _a)) 00030 { 00031 if (!(Math::isfinite(_a) && _a > 0)) 00032 throw GeographicErr("Major radius is not positive"); 00033 if (!(Math::isfinite(_GM) && _GM > 0)) 00034 throw GeographicErr("Gravitational constants is not positive"); 00035 bool flatp = _f > 0 && Math::isfinite(_f); 00036 if (_J2 > 0 && Math::isfinite(_J2) && flatp) 00037 throw GeographicErr("Cannot specify both f and J2"); 00038 if (!(_J2 > 0 && Math::isfinite(_J2)) && !flatp) 00039 throw GeographicErr("Must specify one of f and J2"); 00040 if (!(Math::isfinite(_omega) && _omega != 0)) 00041 throw GeographicErr("Angular velocity is not non-zero"); 00042 real K = 2 * _aomega2 * _a / (15 * _GM); 00043 if (flatp) { 00044 _e2 = _f * (2 - _f); 00045 _ep2 = _e2 / (1 - _e2); 00046 _q0 = qf(_ep2); 00047 _J2 = _e2 * ( 1 - K * sqrt(_e2) / _q0) / 3; // H+M, Eq 2-90 00048 } else { 00049 _e2 = 3 * _J2; // See Moritz (1980), p 398. 00050 for (int j = 0; j < maxit_; ++j) { 00051 real e2a = _e2; 00052 real q0 = qf(_e2 / (1 - _e2)); 00053 _e2 = 3 * _J2 + K * _e2 * sqrt(_e2) / q0; 00054 if (_e2 == e2a) 00055 break; 00056 } 00057 _f = _e2 / (1 + sqrt(1 - _e2)); 00058 _ep2 = _e2 / (1 - _e2); 00059 _q0 = qf(_ep2); 00060 } 00061 _earth = Geocentric(_a, _f); 00062 _b = _a * (1 - _f); 00063 _E = a * sqrt(_e2); // H+M, Eq 2-54 00064 _U0 = _GM / _E * atan(sqrt(_ep2)) + _aomega2 / 3; // H+M, Eq 2-61 00065 // The approximate ratio of the centrifugal acceleration (at the equator) 00066 // to gravity. 00067 _m = _aomega2 * _b / _GM; // H+M, Eq 2-70 00068 real 00069 Q = _m * sqrt(_ep2) * qpf(_ep2) / (3 * _q0), 00070 G = (1 - _m - Q / 2); 00071 _gammae = _GM / (_a * _b) * G; // H+M, Eq 2-73 00072 _gammap = _GM / (_a * _a) * (1 + Q); // H+M, Eq 2-74 00073 // k = b * gammap / (a * gammae) - 1 00074 _k = (_m + 3 * Q / 2 - _e2 * (1 + Q)) / G; 00075 // f* = (gammap - gammae) / gammae 00076 _fstar = (_m + 3 * Q / 2 - _f * (1 + Q)) / G; 00077 } 00078 00079 const NormalGravity 00080 NormalGravity::WGS84(Constants::WGS84_a<real>(), Constants::WGS84_GM<real>(), 00081 Constants::WGS84_omega<real>(), 00082 Constants::WGS84_f<real>(), 0); 00083 00084 const NormalGravity 00085 NormalGravity::GRS80(Constants::GRS80_a<real>(), Constants::GRS80_GM<real>(), 00086 Constants::GRS80_omega<real>(), 00087 0, Constants::GRS80_J2<real>()); 00088 00089 Math::real NormalGravity::qf(real ep2) throw() { 00090 // Compute 00091 // 00092 // ((1 + 3/e'^2) * atan(e') - 3/e')/2 00093 // 00094 // See H+M, Eq 2-57, with E/u = e'. This suffers from two levels of 00095 // cancelation. The e'^-1 and e'^1 terms drop out, so that the leading 00096 // term is O(e'^3). 00097 real ep = sqrt(ep2); 00098 if (abs(ep2) > real(0.5)) // Use the closed expression 00099 return ((1 + 3 / ep2) * atan(ep) - 3 / ep)/2; 00100 else { 00101 real ep2n = 1, q = 0; // The series expansion H+M, Eq 2-86 00102 for (int n = 1; ; ++n) { 00103 ep2n *= -ep2; 00104 real 00105 t = (ep2n * n) / ((2 * n + 1) * (2 * n + 3)), 00106 qn = q + t; 00107 if (qn == q) 00108 break; 00109 q = qn; 00110 } 00111 q *= -2 * ep; 00112 return q; 00113 } 00114 } 00115 00116 Math::real NormalGravity::qpf(real ep2) throw() { 00117 // Compute 00118 // 00119 // 3*(1 + 1/e'^2) * (1 - atan(e')/e') - 1 00120 // 00121 // See H+M, Eq 2-67, with E/u = e'. This suffers from two levels of 00122 // cancelation. The e'^-2 and e'^0 terms drop out, so that the leading 00123 // term is O(e'^2). 00124 if (abs(ep2) > real(0.5)) { // Use the closed expression 00125 real ep = sqrt(ep2); 00126 return 3 * (1 + 1 / ep2) * (1 - atan(ep) / ep) - 1; 00127 } else { 00128 real ep2n = 1, qp = 0; // The series expansion H+M, Eq 2-101c 00129 for (int n = 1; ; ++n) { 00130 ep2n *= -ep2; 00131 real 00132 t = ep2n / ((2 * n + 1) * (2 * n + 3)), 00133 qpn = qp + t; 00134 if (qpn == qp) 00135 break; 00136 qp = qpn; 00137 } 00138 qp *= -6; 00139 return qp; 00140 } 00141 } 00142 00143 Math::real NormalGravity::Jn(int n) const throw() { 00144 // Note Jn(0) = -1; Jn(2) = _J2; Jn(odd) = 0 00145 if (n & 1 || n < 0) 00146 return 0; 00147 n /= 2; 00148 real e2n = 1; // Perhaps this should just be e2n = pow(-_e2, n); 00149 for (int j = n; j--;) 00150 e2n *= -_e2; 00151 return // H+M, Eq 2-92 00152 -3 * e2n * (1 - n + 5 * n * _J2 / _e2) / ((2 * n + 1) * (2 * n + 3)); 00153 } 00154 00155 Math::real NormalGravity::SurfaceGravity(real lat) const throw() { 00156 real 00157 phi = lat * Math::degree<real>(), 00158 sphi2 = abs(lat) == 90 ? 1 : Math::sq(sin(phi)); 00159 // H+M, Eq 2-78 00160 return _gammae * (1 + _k * sphi2) / sqrt(1 - _e2 * sphi2); 00161 } 00162 00163 Math::real NormalGravity::V0(real X, real Y, real Z, 00164 real& GammaX, real& GammaY, real& GammaZ) 00165 const throw() { 00166 // See H+M, Sec 6-2 00167 real 00168 p = Math::hypot(X, Y), 00169 clam = p ? X/p : 1, 00170 slam = p ? Y/p : 0, 00171 r = Math::hypot(p, Z), 00172 Q = Math::sq(r) - Math::sq(_E), 00173 t2 = Math::sq(2 * _E * Z), 00174 disc = sqrt(Math::sq(Q) + t2), 00175 // This is H+M, Eq 6-8a, but generalized to deal with Q negative 00176 // accurately. 00177 u = sqrt((Q >= 0 ? (Q + disc) : t2 / (disc - Q)) / 2), 00178 uE = Math::hypot(u, _E), 00179 // H+M, Eq 6-8b 00180 sbet = Z * uE, 00181 cbet = p * u, 00182 s = Math::hypot(cbet, sbet); 00183 cbet = s ? cbet/s : 0; 00184 sbet = s ? sbet/s : 1; 00185 real 00186 invw = uE / Math::hypot(u, _E * sbet), // H+M, Eq 2-63 00187 ep = _E/u, 00188 ep2 = Math::sq(ep), 00189 q = qf(ep2) / _q0, 00190 qp = qpf(ep2) / _q0, 00191 // H+M, Eqs 2-62 + 6-9, but omitting last (rotational) term . 00192 Vres = (_GM / _E * atan(_E / u) 00193 + _aomega2 * q * (Math::sq(sbet) - 1/real(3)) / 2), 00194 // H+M, Eq 6-10 00195 gamu = - invw * (_GM 00196 + (_aomega2 * _E * qp 00197 * (Math::sq(sbet) - 1/real(3)) / 2)) / Math::sq(uE), 00198 gamb = _aomega2 * q * sbet * cbet * invw / uE, 00199 t = u * invw / uE; 00200 // H+M, Eq 6-12 00201 GammaX = t * cbet * gamu - invw * sbet * gamb; 00202 GammaY = GammaX * slam; 00203 GammaX *= clam; 00204 GammaZ = invw * sbet * gamu + t * cbet * gamb; 00205 return Vres; 00206 } 00207 00208 Math::real NormalGravity::Phi(real X, real Y, real& fX, real& fY) 00209 const throw() { 00210 fX = _omega2 * X; 00211 fY = _omega2 * Y; 00212 // N.B. fZ = 0; 00213 return _omega2 * (Math::sq(X) + Math::sq(Y)) / 2; 00214 } 00215 00216 Math::real NormalGravity::U(real X, real Y, real Z, 00217 real& gammaX, real& gammaY, real& gammaZ) 00218 const throw() { 00219 real fX, fY; 00220 real Ures = V0(X, Y, Z, gammaX, gammaY, gammaZ) + Phi(X, Y, fX, fY); 00221 gammaX += fX; 00222 gammaY += fY; 00223 return Ures; 00224 } 00225 00226 Math::real NormalGravity::Gravity(real lat, real h, 00227 real& gammay, real& gammaz) 00228 const throw() { 00229 real X, Y, Z; 00230 real M[Geocentric::dim2_]; 00231 _earth.IntForward(lat, 0, h, X, Y, Z, M); 00232 real gammaX, gammaY, gammaZ, 00233 Ures = U(X, Y, Z, gammaX, gammaY, gammaZ); 00234 // gammax = M[0] * gammaX + M[3] * gammaY + M[6] * gammaZ; 00235 gammay = M[1] * gammaX + M[4] * gammaY + M[7] * gammaZ; 00236 gammaz = M[2] * gammaX + M[5] * gammaY + M[8] * gammaZ; 00237 return Ures; 00238 } 00239 00240 } // namespace GeographicLib