GeographicLib  1.35
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) <charles@karney.com> and licensed under
6  * the MIT/X11 License. For more information, see
7  * http://geographiclib.sourceforge.net/
8  **********************************************************************/
9 
11 
12 namespace GeographicLib {
13 
14  using namespace std;
15 
16  NormalGravity::NormalGravity(real a, real GM, real omega, real f, real J2)
17  : _a(a)
18  , _GM(GM)
19  , _omega(omega)
20  , _f(f)
21  , _J2(J2)
22  , _omega2(Math::sq(_omega))
23  , _aomega2(Math::sq(_omega * _a))
24  {
25  if (!(Math::isfinite(_a) && _a > 0))
26  throw GeographicErr("Major radius is not positive");
27  if (!(Math::isfinite(_GM) && _GM > 0))
28  throw GeographicErr("Gravitational constants is not positive");
29  bool flatp = _f > 0 && Math::isfinite(_f);
30  if (_J2 > 0 && Math::isfinite(_J2) && flatp)
31  throw GeographicErr("Cannot specify both f and J2");
32  if (!(_J2 > 0 && Math::isfinite(_J2)) && !flatp)
33  throw GeographicErr("Must specify one of f and J2");
34  if (!(Math::isfinite(_omega) && _omega != 0))
35  throw GeographicErr("Angular velocity is not non-zero");
36  real K = 2 * _aomega2 * _a / (15 * _GM);
37  if (flatp) {
38  _e2 = _f * (2 - _f);
39  _ep2 = _e2 / (1 - _e2);
40  _q0 = qf(_ep2);
41  _J2 = _e2 * ( 1 - K * sqrt(_e2) / _q0) / 3; // H+M, Eq 2-90
42  } else {
43  _e2 = 3 * _J2; // See Moritz (1980), p 398.
44  for (int j = 0; j < maxit_; ++j) {
45  real e2a = _e2;
46  real q0 = qf(_e2 / (1 - _e2));
47  _e2 = 3 * _J2 + K * _e2 * sqrt(_e2) / q0;
48  if (_e2 == e2a)
49  break;
50  }
51  _f = _e2 / (1 + sqrt(1 - _e2));
52  _ep2 = _e2 / (1 - _e2);
53  _q0 = qf(_ep2);
54  }
55  _earth = Geocentric(_a, _f);
56  _b = _a * (1 - _f);
57  _E = a * sqrt(_e2); // H+M, Eq 2-54
58  _U0 = _GM / _E * atan(sqrt(_ep2)) + _aomega2 / 3; // H+M, Eq 2-61
59  // The approximate ratio of the centrifugal acceleration (at the equator)
60  // to gravity.
61  _m = _aomega2 * _b / _GM; // H+M, Eq 2-70
62  real
63  Q = _m * sqrt(_ep2) * qpf(_ep2) / (3 * _q0),
64  G = (1 - _m - Q / 2);
65  _gammae = _GM / (_a * _b) * G; // H+M, Eq 2-73
66  _gammap = _GM / (_a * _a) * (1 + Q); // H+M, Eq 2-74
67  // k = b * gammap / (a * gammae) - 1
68  _k = (_m + 3 * Q / 2 - _e2 * (1 + Q)) / G;
69  // f* = (gammap - gammae) / gammae
70  _fstar = (_m + 3 * Q / 2 - _f * (1 + Q)) / G;
71  }
72 
73  const NormalGravity
74  NormalGravity::WGS84(Constants::WGS84_a<real>(), Constants::WGS84_GM<real>(),
75  Constants::WGS84_omega<real>(),
76  Constants::WGS84_f<real>(), 0);
77 
78  const NormalGravity
79  NormalGravity::GRS80(Constants::GRS80_a<real>(), Constants::GRS80_GM<real>(),
80  Constants::GRS80_omega<real>(),
81  0, Constants::GRS80_J2<real>());
82 
83  Math::real NormalGravity::qf(real ep2) throw() {
84  // Compute
85  //
86  // ((1 + 3/e'^2) * atan(e') - 3/e')/2
87  //
88  // See H+M, Eq 2-57, with E/u = e'. This suffers from two levels of
89  // cancelation. The e'^-1 and e'^1 terms drop out, so that the leading
90  // term is O(e'^3).
91  real ep = sqrt(ep2);
92  if (abs(ep2) > real(0.5)) // Use the closed expression
93  return ((1 + 3 / ep2) * atan(ep) - 3 / ep)/2;
94  else {
95  real ep2n = 1, q = 0; // The series expansion H+M, Eq 2-86
96  for (int n = 1; ; ++n) {
97  ep2n *= -ep2;
98  real
99  t = (ep2n * n) / ((2 * n + 1) * (2 * n + 3)),
100  qn = q + t;
101  if (qn == q)
102  break;
103  q = qn;
104  }
105  q *= -2 * ep;
106  return q;
107  }
108  }
109 
110  Math::real NormalGravity::qpf(real ep2) throw() {
111  // Compute
112  //
113  // 3*(1 + 1/e'^2) * (1 - atan(e')/e') - 1
114  //
115  // See H+M, Eq 2-67, with E/u = e'. This suffers from two levels of
116  // cancelation. The e'^-2 and e'^0 terms drop out, so that the leading
117  // term is O(e'^2).
118  if (abs(ep2) > real(0.5)) { // Use the closed expression
119  real ep = sqrt(ep2);
120  return 3 * (1 + 1 / ep2) * (1 - atan(ep) / ep) - 1;
121  } else {
122  real ep2n = 1, qp = 0; // The series expansion H+M, Eq 2-101c
123  for (int n = 1; ; ++n) {
124  ep2n *= -ep2;
125  real
126  t = ep2n / ((2 * n + 1) * (2 * n + 3)),
127  qpn = qp + t;
128  if (qpn == qp)
129  break;
130  qp = qpn;
131  }
132  qp *= -6;
133  return qp;
134  }
135  }
136 
137  Math::real NormalGravity::Jn(int n) const throw() {
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 
149  Math::real NormalGravity::SurfaceGravity(real lat) const throw() {
150  real
151  phi = lat * Math::degree<real>(),
152  sphi2 = abs(lat) == 90 ? 1 : Math::sq(sin(phi));
153  // H+M, Eq 2-78
154  return _gammae * (1 + _k * sphi2) / sqrt(1 - _e2 * sphi2);
155  }
156 
157  Math::real NormalGravity::V0(real X, real Y, real Z,
158  real& GammaX, real& GammaY, real& GammaZ)
159  const throw() {
160  // See H+M, Sec 6-2
161  real
162  p = Math::hypot(X, Y),
163  clam = p ? X/p : 1,
164  slam = p ? Y/p : 0,
165  r = Math::hypot(p, Z),
166  Q = Math::sq(r) - Math::sq(_E),
167  t2 = Math::sq(2 * _E * Z),
168  disc = sqrt(Math::sq(Q) + t2),
169  // This is H+M, Eq 6-8a, but generalized to deal with Q negative
170  // accurately.
171  u = sqrt((Q >= 0 ? (Q + disc) : t2 / (disc - Q)) / 2),
172  uE = Math::hypot(u, _E),
173  // H+M, Eq 6-8b
174  sbet = Z * uE,
175  cbet = p * u,
176  s = Math::hypot(cbet, sbet);
177  cbet = s ? cbet/s : 0;
178  sbet = s ? sbet/s : 1;
179  real
180  invw = uE / Math::hypot(u, _E * sbet), // H+M, Eq 2-63
181  ep = _E/u,
182  ep2 = Math::sq(ep),
183  q = qf(ep2) / _q0,
184  qp = qpf(ep2) / _q0,
185  // H+M, Eqs 2-62 + 6-9, but omitting last (rotational) term .
186  Vres = (_GM / _E * atan(_E / u)
187  + _aomega2 * q * (Math::sq(sbet) - 1/real(3)) / 2),
188  // H+M, Eq 6-10
189  gamu = - invw * (_GM
190  + (_aomega2 * _E * qp
191  * (Math::sq(sbet) - 1/real(3)) / 2)) / Math::sq(uE),
192  gamb = _aomega2 * q * sbet * cbet * invw / uE,
193  t = u * invw / uE;
194  // H+M, Eq 6-12
195  GammaX = t * cbet * gamu - invw * sbet * gamb;
196  GammaY = GammaX * slam;
197  GammaX *= clam;
198  GammaZ = invw * sbet * gamu + t * cbet * gamb;
199  return Vres;
200  }
201 
202  Math::real NormalGravity::Phi(real X, real Y, real& fX, real& fY)
203  const throw() {
204  fX = _omega2 * X;
205  fY = _omega2 * Y;
206  // N.B. fZ = 0;
207  return _omega2 * (Math::sq(X) + Math::sq(Y)) / 2;
208  }
209 
210  Math::real NormalGravity::U(real X, real Y, real Z,
211  real& gammaX, real& gammaY, real& gammaZ)
212  const throw() {
213  real fX, fY;
214  real Ures = V0(X, Y, Z, gammaX, gammaY, gammaZ) + Phi(X, Y, fX, fY);
215  gammaX += fX;
216  gammaY += fY;
217  return Ures;
218  }
219 
221  real& gammay, real& gammaz)
222  const throw() {
223  real X, Y, Z;
224  real M[Geocentric::dim2_];
225  _earth.IntForward(lat, 0, h, X, Y, Z, M);
226  real gammaX, gammaY, gammaZ,
227  Ures = U(X, Y, Z, gammaX, gammaY, gammaZ);
228  // gammax = M[0] * gammaX + M[3] * gammaY + M[6] * gammaZ;
229  gammay = M[1] * gammaX + M[4] * gammaY + M[7] * gammaZ;
230  gammaz = M[2] * gammaX + M[5] * gammaY + M[8] * gammaZ;
231  return Ures;
232  }
233 
234 } // namespace GeographicLib
Math::real SurfaceGravity(real lat) const
GeographicLib::Math::real real
Definition: GeodSolve.cpp:40
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:435
static const NormalGravity GRS80
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:73
Geocentric coordinates
Definition: Geocentric.hpp:61
static T hypot(T x, T y)
Definition: Math.hpp:165
static T sq(T x)
Definition: Math.hpp:153
Header for GeographicLib::NormalGravity class.
Exception handling for GeographicLib.
Definition: Constants.hpp:320
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
Math::real V0(real X, real Y, real Z, real &GammaX, real &GammaY, real &GammaZ) const