GeographicLib  1.45
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-2015) <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 constant 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  // H+M, Eq 2-61
54  _U0 = _GM / (_E ? _E / atan(sqrt(_ep2)) : _b) + _aomega2 / 3;
55  // The approximate ratio of the centrifugal acceleration (at the equator)
56  // to gravity.
57  _m = _aomega2 * _b / _GM; // H+M, Eq 2-70
58  real
59  Q = _m * (_q0 ? sqrt(_ep2) * qpf(_ep2) / (3 * _q0) : 1),
60  G = (1 - _m - Q / 2);
61  _gammae = _GM / (_a * _b) * G; // H+M, Eq 2-73
62  _gammap = _GM / (_a * _a) * (1 + Q); // H+M, Eq 2-74
63  // k = b * gammap / (a * gammae) - 1
64  _k = (_m + 3 * Q / 2 - _e2 * (1 + Q)) / G;
65  // f* = (gammap - gammae) / gammae
66  _fstar = (_m + 3 * Q / 2 - _f * (1 + Q)) / G;
67  }
68 
70  static const NormalGravity wgs84(Constants::WGS84_a(),
73  Constants::WGS84_f(), 0);
74  return wgs84;
75  }
76 
78  static const NormalGravity grs80(Constants::GRS80_a(),
81  0, Constants::GRS80_J2());
82  return grs80;
83  }
84 
85  // (atan(y)-(y-y^3/3+y^5/5))/y^7 (y = sqrt(x)) = -1/7+x/9-x^2/11+x^3/13...
86  Math::real NormalGravity::atan7(real x) {
87  if (abs(x) >= real(0.5)) {
88  real y = sqrt(abs(x)), x2 = Math::sq(x);
89  return ((x > 0 ? atan(y) : Math::atanh(y))- y * (1 - x / 3 + x2 / 5)) /
90  (x * x2 * y);
91  } else {
92  real xn = -1, q = 0;
93  for (int n = 7; ; n += 2) {
94  real qn = q + xn / n;
95  if (qn == q)
96  break;
97  q = qn;
98  xn *= -x;
99  }
100  return q;
101  }
102  }
103 
104  // (atan(y)-(y-y^3/3))/y^5 (y = sqrt(x)) = 1/5-x/7+x^2/9-x^3/11...
105  Math::real NormalGravity::atan5(real x)
106  { return 1/real(5) + x * atan7(x); }
107 
108  Math::real NormalGravity::qf(real ep2) {
109  // Compute
110  //
111  // ((1 + 3/e'^2) * atan(e') - 3/e')/2
112  //
113  // See H+M, Eq 2-57, with E/u = e'. This suffers from two levels of
114  // cancelation. The e'^-1 and e'^1 terms drop out, so that the leading
115  // term is O(e'^3). Substitute atan(e') = e' - e'^3/3 + e'^5*atan5(e'^2)
116  return sqrt(ep2) * ep2 * (3 * (3 + ep2) * atan5(ep2) - 1) / 6;
117  }
118 
119  Math::real NormalGravity::dq(real ep2) {
120  // Compute d qf(ep2) / d ep2 and substitute
121  // atan(e') = e' - e'^3/3 + e'^5/5 + e'^7*atan7(e'^2)
122  return sqrt(ep2) * (5 - 3 * (1 + ep2) * (1 + 5 * ep2 * atan7(ep2))) /
123  (10 * (1 + ep2));
124  }
125 
126  Math::real NormalGravity::qpf(real ep2) {
127  // Compute
128  //
129  // 3*(1 + 1/e'^2) * (1 - atan(e')/e') - 1
130  //
131  // See H+M, Eq 2-67, with E/u = e'. This suffers from two levels of
132  // cancelation. The e'^-2 and e'^0 terms drop out, so that the leading
133  // term is O(e'^2).
134  return ep2 * (1 - 3 * (1 + ep2) * atan5(ep2));
135  }
136 
137  Math::real NormalGravity::Jn(int n) const {
138  // Note Jn(0) = -1; Jn(2) = _J2; Jn(odd) = 0
139  if (n & 1 || n < 0)
140  return 0;
141  n /= 2;
142  real e2n = 1; // Perhaps this should just be e2n = pow(-_e2, n);
143  for (int j = n; j--;)
144  e2n *= -_e2;
145  return // H+M, Eq 2-92
146  -3 * e2n * ((1 - n) + 5 * n * _J2 / _e2) / ((2 * n + 1) * (2 * n + 3));
147  }
148 
150  real sphi = Math::sind(Math::LatFix(lat));
151  // H+M, Eq 2-78
152  return _gammae * (1 + _k * Math::sq(sphi)) / sqrt(1 - _e2 * Math::sq(sphi));
153  }
154 
155  Math::real NormalGravity::V0(real X, real Y, real Z,
156  real& GammaX, real& GammaY, real& GammaZ)
157  const {
158  // See H+M, Sec 6-2
159  real
160  p = Math::hypot(X, Y),
161  clam = p ? X/p : 1,
162  slam = p ? Y/p : 0,
163  r = Math::hypot(p, Z),
164  Q = Math::sq(r) - Math::sq(_E),
165  t2 = Math::sq(2 * _E * Z),
166  disc = sqrt(Math::sq(Q) + t2),
167  // This is H+M, Eq 6-8a, but generalized to deal with Q negative
168  // accurately.
169  u = sqrt((Q >= 0 ? (Q + disc) : t2 / (disc - Q)) / 2),
170  uE = Math::hypot(u, _E),
171  // H+M, Eq 6-8b
172  sbet = Z * uE,
173  cbet = p * u,
174  s = Math::hypot(cbet, sbet);
175  cbet = s ? cbet/s : 0;
176  sbet = s ? sbet/s : 1;
177  real
178  invw = uE / Math::hypot(u, _E * sbet), // H+M, Eq 2-63
179  ep = _E/u,
180  ep2 = Math::sq(ep),
181  q = _q0 ? qf(ep2) / _q0 : pow(_a / u, 3),
182  qp = qpf(ep2) / _q0,
183  // H+M, Eqs 2-62 + 6-9, but omitting last (rotational) term .
184  Vres = (_GM / (_E ? _E / atan(_E / u) : u)
185  + _aomega2 * q * (Math::sq(sbet) - 1/real(3)) / 2),
186  // H+M, Eq 6-10
187  gamu = - invw * (_GM
188  + (_aomega2 * (_q0 ? _E * qp : 3 * q * u)
189  * (Math::sq(sbet) - 1/real(3)) / 2)) / Math::sq(uE),
190  gamb = _aomega2 * q * sbet * cbet * invw / uE,
191  t = u * invw / uE;
192  // H+M, Eq 6-12
193  GammaX = t * cbet * gamu - invw * sbet * gamb;
194  GammaY = GammaX * slam;
195  GammaX *= clam;
196  GammaZ = invw * sbet * gamu + t * cbet * gamb;
197  return Vres;
198  }
199 
200  Math::real NormalGravity::Phi(real X, real Y, real& fX, real& fY)
201  const {
202  fX = _omega2 * X;
203  fY = _omega2 * Y;
204  // N.B. fZ = 0;
205  return _omega2 * (Math::sq(X) + Math::sq(Y)) / 2;
206  }
207 
208  Math::real NormalGravity::U(real X, real Y, real Z,
209  real& gammaX, real& gammaY, real& gammaZ)
210  const {
211  real fX, fY;
212  real Ures = V0(X, Y, Z, gammaX, gammaY, gammaZ) + Phi(X, Y, fX, fY);
213  gammaX += fX;
214  gammaY += fY;
215  return Ures;
216  }
217 
219  real& gammay, real& gammaz)
220  const {
221  real X, Y, Z;
222  real M[Geocentric::dim2_];
223  _earth.IntForward(lat, 0, h, X, Y, Z, M);
224  real gammaX, gammaY, gammaZ,
225  Ures = U(X, Y, Z, gammaX, gammaY, gammaZ);
226  // gammax = M[0] * gammaX + M[3] * gammaY + M[6] * gammaZ;
227  gammay = M[1] * gammaX + M[4] * gammaY + M[7] * gammaZ;
228  gammaz = M[2] * gammaX + M[5] * gammaY + M[8] * gammaZ;
229  return Ures;
230  }
231 
233  real omega, real J2) {
234  real
235  K = 2 * Math::sq(a * omega) * a / (15 * GM),
236  e2 = 3 * J2; // See Moritz (1980), p 398.
237  // Solve using Newton's method
238  for (int j = 0; j < maxit_ || GEOGRAPHICLIB_PANIC; ++j) {
239  real e2a = e2,
240  ep2 = e2 / (1 - e2),
241  q0 = qf(ep2),
242  dq0 = dq(ep2) / Math::sq(1 - e2),
243  h = e2 * (1 - sqrt(e2) * K / q0) - 3 * J2,
244  dh = 1 - sqrt(e2) * K * (3 * q0 - 2 * e2 * dq0) / (2 * Math::sq(q0)),
245  de2 = - h / dh;
246  e2 = e2a + de2;
247  if (e2 == e2a)
248  break;
249  }
250  return e2 / (1 + sqrt(1 - e2));
251  }
252 
254  real omega, real f) {
255  real
256  K = 2 * Math::sq(a * omega) * a / (15 * GM),
257  e2 = f * (2 - f),
258  q0 = qf(e2 / (1 - e2));
259  return e2 * (1 - K * sqrt(e2) / q0) / 3; // H+M, Eq 2-90
260  }
261 
262 } // 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:768
static T LatFix(T x)
Definition: Math.hpp:482
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:102
static T sind(T x)
Definition: Math.hpp:597
static T atanh(T x)
Definition: Math.hpp:342
Geocentric coordinates
Definition: Geocentric.hpp:67
static T hypot(T x, T y)
Definition: Math.hpp:257
static T sq(T x)
Definition: Math.hpp:246
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 const NormalGravity & GRS80()
Header for GeographicLib::NormalGravity class.
Exception handling for GeographicLib.
Definition: Constants.hpp:386
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