GeographicLib  1.21
NormalGravity.cpp
Go to the documentation of this file.
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