GeographicLib  1.43
NormalGravity.cpp
Go to the documentation of this file.
1 /**
2  * \file NormalGravity.cpp
3  * \brief Implementation for GeographicLib::NormalGravity class
4  *
5  * Copyright (c) Charles Karney (2011-2014) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * http://geographiclib.sourceforge.net/
8  **********************************************************************/
9 
11 
12 #if defined(_MSC_VER)
13 // Squelch warnings about constant conditional expressions
14 # pragma warning (disable: 4127)
15 #endif
16 
17 namespace GeographicLib {
18 
19  using namespace std;
20 
21  NormalGravity::NormalGravity(real a, real GM, real omega, real f, real J2)
22  : _a(a)
23  , _GM(GM)
24  , _omega(omega)
25  , _f(f)
26  , _J2(J2)
27  , _omega2(Math::sq(_omega))
28  , _aomega2(Math::sq(_omega * _a))
29  {
30  if (!(Math::isfinite(_a) && _a > 0))
31  throw GeographicErr("Major radius is not positive");
32  if (!(Math::isfinite(_GM) && _GM > 0))
33  throw GeographicErr("Gravitational constants is not positive");
34  if (!(_omega == 0 && _f == 0 && _J2 == 0)) {
35  bool flatp = _f > 0 && Math::isfinite(_f);
36  if (_J2 > 0 && Math::isfinite(_J2) && flatp)
37  throw GeographicErr("Cannot specify both f and J2");
38  if (!(_J2 > 0 && Math::isfinite(_J2)) && !flatp)
39  throw GeographicErr("Must specify one of f and J2");
40  if (!(Math::isfinite(_omega) && _omega != 0))
41  throw GeographicErr("Angular velocity is not non-zero");
42  if (flatp)
43  _J2 = FlatteningToJ2(a, GM, omega, f);
44  else
45  _f = J2ToFlattening(a, GM, omega, J2);
46  } // else we have a sphere: omega = f = J2 = 0
47  _e2 = _f * (2 - _f);
48  _ep2 = _e2 / (1 - _e2);
49  _q0 = qf(_ep2);
50  _earth = Geocentric(_a, _f);
51  _b = _a * (1 - _f);
52  _E = a * sqrt(_e2); // H+M, Eq 2-54
53  _U0 = _GM / _E * atan(sqrt(_ep2)) + _aomega2 / 3; // H+M, Eq 2-61
54  // The approximate ratio of the centrifugal acceleration (at the equator)
55  // to gravity.
56  _m = _aomega2 * _b / _GM; // H+M, Eq 2-70
57  real
58  Q = _m * sqrt(_ep2) * qpf(_ep2) / (3 * _q0),
59  G = (1 - _m - Q / 2);
60  _gammae = _GM / (_a * _b) * G; // H+M, Eq 2-73
61  _gammap = _GM / (_a * _a) * (1 + Q); // H+M, Eq 2-74
62  // k = b * gammap / (a * gammae) - 1
63  _k = (_m + 3 * Q / 2 - _e2 * (1 + Q)) / G;
64  // f* = (gammap - gammae) / gammae
65  _fstar = (_m + 3 * Q / 2 - _f * (1 + Q)) / G;
66  }
67 
69  static const NormalGravity wgs84(Constants::WGS84_a(),
72  Constants::WGS84_f(), 0);
73  return wgs84;
74  }
75 
77  static const NormalGravity grs80(Constants::GRS80_a(),
80  0, Constants::GRS80_J2());
81  return grs80;
82  }
83 
84  // (atan(y)-(y-y^3/3+y^5/5))/y^7 (y = sqrt(x)) = -1/7+x/9-x^2/11+x^3/13...
85  Math::real NormalGravity::atan7(real x) {
86  if (abs(x) >= real(0.5)) {
87  real y = sqrt(abs(x)), x2 = Math::sq(x);
88  return ((x > 0 ? atan(y) : Math::atanh(y))- y * (1 - x / 3 + x2 / 5)) /
89  (x * x2 * y);
90  } else {
91  real xn = -1, q = 0;
92  for (int n = 7; ; n += 2) {
93  real qn = q + xn / n;
94  if (qn == q)
95  break;
96  q = qn;
97  xn *= -x;
98  }
99  return q;
100  }
101  }
102 
103  // (atan(y)-(y-y^3/3))/y^5 (y = sqrt(x)) = 1/5-x/7+x^2/9-x^3/11...
104  Math::real NormalGravity::atan5(real x)
105  { return 1/real(5) + x * atan7(x); }
106 
107  Math::real NormalGravity::qf(real ep2) {
108  // Compute
109  //
110  // ((1 + 3/e'^2) * atan(e') - 3/e')/2
111  //
112  // See H+M, Eq 2-57, with E/u = e'. This suffers from two levels of
113  // cancelation. The e'^-1 and e'^1 terms drop out, so that the leading
114  // term is O(e'^3). Substitute atan(e') = e' - e'^3/3 + e'^5*atan5(e'^2)
115  return sqrt(ep2) * ep2 * (3 * (3 + ep2) * atan5(ep2) - 1) / 6;
116  }
117 
118  Math::real NormalGravity::dq(real ep2) {
119  // Compute d qf(ep2) / d ep2 and substitute
120  // atan(e') = e' - e'^3/3 + e'^5/5 + e'^7*atan7(e'^2)
121  return sqrt(ep2) * (5 - 3 * (1 + ep2) * (1 + 5 * ep2 * atan7(ep2))) /
122  (10 * (1 + ep2));
123  }
124 
125  Math::real NormalGravity::qpf(real ep2) {
126  // Compute
127  //
128  // 3*(1 + 1/e'^2) * (1 - atan(e')/e') - 1
129  //
130  // See H+M, Eq 2-67, with E/u = e'. This suffers from two levels of
131  // cancelation. The e'^-2 and e'^0 terms drop out, so that the leading
132  // term is O(e'^2).
133  return ep2 * (1 - 3 * (1 + ep2) * atan5(ep2));
134  }
135 
136  Math::real NormalGravity::Jn(int n) const {
137  // Note Jn(0) = -1; Jn(2) = _J2; Jn(odd) = 0
138  if (n & 1 || n < 0)
139  return 0;
140  n /= 2;
141  real e2n = 1; // Perhaps this should just be e2n = pow(-_e2, n);
142  for (int j = n; j--;)
143  e2n *= -_e2;
144  return // H+M, Eq 2-92
145  -3 * e2n * ((1 - n) + 5 * n * _J2 / _e2) / ((2 * n + 1) * (2 * n + 3));
146  }
147 
149  real
150  phi = lat * Math::degree(),
151  sphi2 = abs(lat) == 90 ? 1 : Math::sq(sin(phi));
152  // H+M, Eq 2-78
153  return _gammae * (1 + _k * sphi2) / sqrt(1 - _e2 * sphi2);
154  }
155 
156  Math::real NormalGravity::V0(real X, real Y, real Z,
157  real& GammaX, real& GammaY, real& GammaZ)
158  const {
159  // See H+M, Sec 6-2
160  real
161  p = Math::hypot(X, Y),
162  clam = p ? X/p : 1,
163  slam = p ? Y/p : 0,
164  r = Math::hypot(p, Z),
165  Q = Math::sq(r) - Math::sq(_E),
166  t2 = Math::sq(2 * _E * Z),
167  disc = sqrt(Math::sq(Q) + t2),
168  // This is H+M, Eq 6-8a, but generalized to deal with Q negative
169  // accurately.
170  u = sqrt((Q >= 0 ? (Q + disc) : t2 / (disc - Q)) / 2),
171  uE = Math::hypot(u, _E),
172  // H+M, Eq 6-8b
173  sbet = Z * uE,
174  cbet = p * u,
175  s = Math::hypot(cbet, sbet);
176  cbet = s ? cbet/s : 0;
177  sbet = s ? sbet/s : 1;
178  real
179  invw = uE / Math::hypot(u, _E * sbet), // H+M, Eq 2-63
180  ep = _E/u,
181  ep2 = Math::sq(ep),
182  q = qf(ep2) / _q0,
183  qp = qpf(ep2) / _q0,
184  // H+M, Eqs 2-62 + 6-9, but omitting last (rotational) term .
185  Vres = (_GM / _E * atan(_E / u)
186  + _aomega2 * q * (Math::sq(sbet) - 1/real(3)) / 2),
187  // H+M, Eq 6-10
188  gamu = - invw * (_GM
189  + (_aomega2 * _E * qp
190  * (Math::sq(sbet) - 1/real(3)) / 2)) / Math::sq(uE),
191  gamb = _aomega2 * q * sbet * cbet * invw / uE,
192  t = u * invw / uE;
193  // H+M, Eq 6-12
194  GammaX = t * cbet * gamu - invw * sbet * gamb;
195  GammaY = GammaX * slam;
196  GammaX *= clam;
197  GammaZ = invw * sbet * gamu + t * cbet * gamb;
198  return Vres;
199  }
200 
201  Math::real NormalGravity::Phi(real X, real Y, real& fX, real& fY)
202  const {
203  fX = _omega2 * X;
204  fY = _omega2 * Y;
205  // N.B. fZ = 0;
206  return _omega2 * (Math::sq(X) + Math::sq(Y)) / 2;
207  }
208 
209  Math::real NormalGravity::U(real X, real Y, real Z,
210  real& gammaX, real& gammaY, real& gammaZ)
211  const {
212  real fX, fY;
213  real Ures = V0(X, Y, Z, gammaX, gammaY, gammaZ) + Phi(X, Y, fX, fY);
214  gammaX += fX;
215  gammaY += fY;
216  return Ures;
217  }
218 
220  real& gammay, real& gammaz)
221  const {
222  real X, Y, Z;
223  real M[Geocentric::dim2_];
224  _earth.IntForward(lat, 0, h, X, Y, Z, M);
225  real gammaX, gammaY, gammaZ,
226  Ures = U(X, Y, Z, gammaX, gammaY, gammaZ);
227  // gammax = M[0] * gammaX + M[3] * gammaY + M[6] * gammaZ;
228  gammay = M[1] * gammaX + M[4] * gammaY + M[7] * gammaZ;
229  gammaz = M[2] * gammaX + M[5] * gammaY + M[8] * gammaZ;
230  return Ures;
231  }
232 
234  real omega, real J2) {
235  real
236  K = 2 * Math::sq(a * omega) * a / (15 * GM),
237  e2 = 3 * J2; // See Moritz (1980), p 398.
238  // Solve using Newton's method
239  for (int j = 0; j < maxit_ || GEOGRAPHICLIB_PANIC; ++j) {
240  real e2a = e2,
241  ep2 = e2 / (1 - e2),
242  q0 = qf(ep2),
243  dq0 = dq(ep2) / Math::sq(1 - e2),
244  h = e2 * (1 - sqrt(e2) * K / q0) - 3 * J2,
245  dh = 1 - sqrt(e2) * K * (3 * q0 - 2 * e2 * dq0) / (2 * Math::sq(q0)),
246  de2 = - h / dh;
247  e2 = e2a + de2;
248  if (e2 == e2a)
249  break;
250  }
251  return e2 / (1 + sqrt(1 - e2));
252  }
253 
255  real omega, real f) {
256  real
257  K = 2 * Math::sq(a * omega) * a / (15 * GM),
258  e2 = f * (2 - f),
259  q0 = qf(e2 / (1 - e2));
260  return e2 * (1 - K * sqrt(e2) / q0) / 3; // H+M, Eq 2-90
261  }
262 
263 } // namespace GeographicLib
Math::real SurfaceGravity(real lat) const
GeographicLib::Math::real real
Definition: GeodSolve.cpp:32
Math::real Gravity(real lat, real h, real &gammay, real &gammaz) const
The normal gravity of the earth.
static bool isfinite(T x)
Definition: Math.hpp:614
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:102
static T atanh(T x)
Definition: Math.hpp:340
Geocentric coordinates
Definition: Geocentric.hpp:67
static T hypot(T x, T y)
Definition: Math.hpp:255
static T sq(T x)
Definition: Math.hpp:244
static Math::real FlatteningToJ2(real a, real GM, real omega, real f)
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
static Math::real J2ToFlattening(real a, real GM, real omega, real J2)
static T degree()
Definition: Math.hpp:228
static const NormalGravity & GRS80()
Header for GeographicLib::NormalGravity class.
Exception handling for GeographicLib.
Definition: Constants.hpp:382
Math::real U(real X, real Y, real Z, real &gammaX, real &gammaY, real &gammaZ) const
static const NormalGravity & WGS84()
Math::real Phi(real X, real Y, real &fX, real &fY) const
#define GEOGRAPHICLIB_PANIC
Definition: Math.hpp:87
Math::real V0(real X, real Y, real Z, real &GammaX, real &GammaY, real &GammaZ) const