GeographicLib  1.51
NormalGravity.hpp
Go to the documentation of this file.
1 /**
2  * \file NormalGravity.hpp
3  * \brief Header for GeographicLib::NormalGravity class
4  *
5  * Copyright (c) Charles Karney (2011-2020) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * https://geographiclib.sourceforge.io/
8  **********************************************************************/
9 
10 #if !defined(GEOGRAPHICLIB_NORMALGRAVITY_HPP)
11 #define GEOGRAPHICLIB_NORMALGRAVITY_HPP 1
12 
15 
16 namespace GeographicLib {
17 
18  /**
19  * \brief The normal gravity of the earth
20  *
21  * "Normal" gravity refers to an idealization of the earth which is modeled
22  * as an rotating ellipsoid. The eccentricity of the ellipsoid, the rotation
23  * speed, and the distribution of mass within the ellipsoid are such that the
24  * ellipsoid is a "level ellipoid", a surface of constant potential
25  * (gravitational plus centrifugal). The acceleration due to gravity is
26  * therefore perpendicular to the surface of the ellipsoid.
27  *
28  * Because the distribution of mass within the ellipsoid is unspecified, only
29  * the potential exterior to the ellipsoid is well defined. In this class,
30  * the mass is assumed to be to concentrated on a "focal disc" of radius,
31  * (<i>a</i><sup>2</sup> &minus; <i>b</i><sup>2</sup>)<sup>1/2</sup>, where
32  * \e a is the equatorial radius of the ellipsoid and \e b is its polar
33  * semi-axis. In the case of an oblate ellipsoid, the mass is concentrated
34  * on a "focal rod" of length 2(<i>b</i><sup>2</sup> &minus;
35  * <i>a</i><sup>2</sup>)<sup>1/2</sup>. As a result the potential is well
36  * defined everywhere.
37  *
38  * There is a closed solution to this problem which is implemented here.
39  * Series "approximations" are only used to evaluate certain combinations of
40  * elementary functions where use of the closed expression results in a loss
41  * of accuracy for small arguments due to cancellation of the leading terms.
42  * However these series include sufficient terms to give full machine
43  * precision.
44  *
45  * Although the formulation used in this class applies to ellipsoids with
46  * arbitrary flattening, in practice, its use should be limited to about
47  * <i>b</i>/\e a &isin; [0.01, 100] or \e f &isin; [&minus;99, 0.99].
48  *
49  * Definitions:
50  * - <i>V</i><sub>0</sub>, the gravitational contribution to the normal
51  * potential;
52  * - &Phi;, the rotational contribution to the normal potential;
53  * - \e U = <i>V</i><sub>0</sub> + &Phi;, the total potential;
54  * - <b>&Gamma;</b> = &nabla;<i>V</i><sub>0</sub>, the acceleration due to
55  * mass of the earth;
56  * - <b>f</b> = &nabla;&Phi;, the centrifugal acceleration;
57  * - <b>&gamma;</b> = &nabla;\e U = <b>&Gamma;</b> + <b>f</b>, the normal
58  * acceleration;
59  * - \e X, \e Y, \e Z, geocentric coordinates;
60  * - \e x, \e y, \e z, local cartesian coordinates used to denote the east,
61  * north and up directions.
62  *
63  * References:
64  * - C. Somigliana, Teoria generale del campo gravitazionale dell'ellissoide
65  * di rotazione, Mem. Soc. Astron. Ital, <b>4</b>, 541--599 (1929).
66  * - W. A. Heiskanen and H. Moritz, Physical Geodesy (Freeman, San
67  * Francisco, 1967), Secs. 1-19, 2-7, 2-8 (2-9, 2-10), 6-2 (6-3).
68  * - B. Hofmann-Wellenhof, H. Moritz, Physical Geodesy (Second edition,
69  * Springer, 2006) https://doi.org/10.1007/978-3-211-33545-1
70  * - H. Moritz, Geodetic Reference System 1980, J. Geodesy 54(3), 395-405
71  * (1980) https://doi.org/10.1007/BF02521480
72  *
73  * For more information on normal gravity see \ref normalgravity.
74  *
75  * Example of use:
76  * \include example-NormalGravity.cpp
77  **********************************************************************/
78 
80  private:
81  static const int maxit_ = 20;
82  typedef Math::real real;
83  friend class GravityModel;
84  real _a, _GM, _omega, _f, _J2, _omega2, _aomega2;
85  real _e2, _ep2, _b, _E, _U0, _gammae, _gammap, _Q0, _k, _fstar;
86  Geocentric _earth;
87  static real atanzz(real x, bool alt) {
88  // This routine obeys the identity
89  // atanzz(x, alt) = atanzz(-x/(1+x), !alt)
90  //
91  // Require x >= -1. Best to call with alt, s.t. x >= 0; this results in
92  // a call to atan, instead of asin, or to asinh, instead of atanh.
93  using std::sqrt; using std::abs; using std::atan; using std::asin;
94  using std::asinh; using std::atanh;
95  real z = sqrt(abs(x));
96  return x == 0 ? 1 :
97  (alt ?
98  (!(x < 0) ? asinh(z) : asin(z)) / sqrt(abs(x) / (1 + x)) :
99  (!(x < 0) ? atan(z) : atanh(z)) / z);
100  }
101  static real atan7series(real x);
102  static real atan5series(real x);
103  static real Qf(real x, bool alt);
104  static real Hf(real x, bool alt);
105  static real QH3f(real x, bool alt);
106  real Jn(int n) const;
107  void Initialize(real a, real GM, real omega, real f_J2, bool geometricp);
108  public:
109 
110  /** \name Setting up the normal gravity
111  **********************************************************************/
112  ///@{
113  /**
114  * Constructor for the normal gravity.
115  *
116  * @param[in] a equatorial radius (meters).
117  * @param[in] GM mass constant of the ellipsoid
118  * (meters<sup>3</sup>/seconds<sup>2</sup>); this is the product of \e G
119  * the gravitational constant and \e M the mass of the earth (usually
120  * including the mass of the earth's atmosphere).
121  * @param[in] omega the angular velocity (rad s<sup>&minus;1</sup>).
122  * @param[in] f_J2 either the flattening of the ellipsoid \e f or the
123  * the dynamical form factor \e J2.
124  * @param[out] geometricp if true (the default), then \e f_J2 denotes the
125  * flattening, else it denotes the dynamical form factor \e J2.
126  * @exception if \e a is not positive or if the other parameters do not
127  * obey the restrictions given below.
128  *
129  * The shape of the ellipsoid can be given in one of two ways:
130  * - geometrically (\e geomtricp = true), the ellipsoid is defined by the
131  * flattening \e f = (\e a &minus; \e b) / \e a, where \e a and \e b are
132  * the equatorial radius and the polar semi-axis. The parameters should
133  * obey \e a &gt; 0, \e f &lt; 1. There are no restrictions on \e GM or
134  * \e omega, in particular, \e GM need not be positive.
135  * - physically (\e geometricp = false), the ellipsoid is defined by the
136  * dynamical form factor <i>J</i><sub>2</sub> = (\e C &minus; \e A) /
137  * <i>Ma</i><sup>2</sup>, where \e A and \e C are the equatorial and
138  * polar moments of inertia and \e M is the mass of the earth. The
139  * parameters should obey \e a &gt; 0, \e GM &gt; 0 and \e J2 &lt; 1/3
140  * &minus; (<i>omega</i><sup>2</sup><i>a</i><sup>3</sup>/<i>GM</i>)
141  * 8/(45&pi;). There is no restriction on \e omega.
142  **********************************************************************/
143  NormalGravity(real a, real GM, real omega, real f_J2,
144  bool geometricp = true);
145 
146  /**
147  * A default constructor for the normal gravity. This sets up an
148  * uninitialized object and is used by GravityModel which constructs this
149  * object before it has read in the parameters for the reference ellipsoid.
150  **********************************************************************/
151  NormalGravity() : _a(-1) {}
152  ///@}
153 
154  /** \name Compute the gravity
155  **********************************************************************/
156  ///@{
157  /**
158  * Evaluate the gravity on the surface of the ellipsoid.
159  *
160  * @param[in] lat the geographic latitude (degrees).
161  * @return &gamma; the acceleration due to gravity, positive downwards
162  * (m s<sup>&minus;2</sup>).
163  *
164  * Due to the axial symmetry of the ellipsoid, the result is independent of
165  * the value of the longitude. This acceleration is perpendicular to the
166  * surface of the ellipsoid. It includes the effects of the earth's
167  * rotation.
168  **********************************************************************/
169  Math::real SurfaceGravity(real lat) const;
170 
171  /**
172  * Evaluate the gravity at an arbitrary point above (or below) the
173  * ellipsoid.
174  *
175  * @param[in] lat the geographic latitude (degrees).
176  * @param[in] h the height above the ellipsoid (meters).
177  * @param[out] gammay the northerly component of the acceleration
178  * (m s<sup>&minus;2</sup>).
179  * @param[out] gammaz the upward component of the acceleration
180  * (m s<sup>&minus;2</sup>); this is usually negative.
181  * @return \e U the corresponding normal potential
182  * (m<sup>2</sup> s<sup>&minus;2</sup>).
183  *
184  * Due to the axial symmetry of the ellipsoid, the result is independent of
185  * the value of the longitude and the easterly component of the
186  * acceleration vanishes, \e gammax = 0. The function includes the effects
187  * of the earth's rotation. When \e h = 0, this function gives \e gammay =
188  * 0 and the returned value matches that of NormalGravity::SurfaceGravity.
189  **********************************************************************/
190  Math::real Gravity(real lat, real h, real& gammay, real& gammaz)
191  const;
192 
193  /**
194  * Evaluate the components of the acceleration due to gravity and the
195  * centrifugal acceleration in geocentric coordinates.
196  *
197  * @param[in] X geocentric coordinate of point (meters).
198  * @param[in] Y geocentric coordinate of point (meters).
199  * @param[in] Z geocentric coordinate of point (meters).
200  * @param[out] gammaX the \e X component of the acceleration
201  * (m s<sup>&minus;2</sup>).
202  * @param[out] gammaY the \e Y component of the acceleration
203  * (m s<sup>&minus;2</sup>).
204  * @param[out] gammaZ the \e Z component of the acceleration
205  * (m s<sup>&minus;2</sup>).
206  * @return \e U = <i>V</i><sub>0</sub> + &Phi; the sum of the
207  * gravitational and centrifugal potentials
208  * (m<sup>2</sup> s<sup>&minus;2</sup>).
209  *
210  * The acceleration given by <b>&gamma;</b> = &nabla;\e U =
211  * &nabla;<i>V</i><sub>0</sub> + &nabla;&Phi; = <b>&Gamma;</b> + <b>f</b>.
212  **********************************************************************/
213  Math::real U(real X, real Y, real Z,
214  real& gammaX, real& gammaY, real& gammaZ) const;
215 
216  /**
217  * Evaluate the components of the acceleration due to the gravitational
218  * force in geocentric coordinates.
219  *
220  * @param[in] X geocentric coordinate of point (meters).
221  * @param[in] Y geocentric coordinate of point (meters).
222  * @param[in] Z geocentric coordinate of point (meters).
223  * @param[out] GammaX the \e X component of the acceleration due to the
224  * gravitational force (m s<sup>&minus;2</sup>).
225  * @param[out] GammaY the \e Y component of the acceleration due to the
226  * @param[out] GammaZ the \e Z component of the acceleration due to the
227  * gravitational force (m s<sup>&minus;2</sup>).
228  * @return <i>V</i><sub>0</sub> the gravitational potential
229  * (m<sup>2</sup> s<sup>&minus;2</sup>).
230  *
231  * This function excludes the centrifugal acceleration and is appropriate
232  * to use for space applications. In terrestrial applications, the
233  * function NormalGravity::U (which includes this effect) should usually be
234  * used.
235  **********************************************************************/
236  Math::real V0(real X, real Y, real Z,
237  real& GammaX, real& GammaY, real& GammaZ) const;
238 
239  /**
240  * Evaluate the centrifugal acceleration in geocentric coordinates.
241  *
242  * @param[in] X geocentric coordinate of point (meters).
243  * @param[in] Y geocentric coordinate of point (meters).
244  * @param[out] fX the \e X component of the centrifugal acceleration
245  * (m s<sup>&minus;2</sup>).
246  * @param[out] fY the \e Y component of the centrifugal acceleration
247  * (m s<sup>&minus;2</sup>).
248  * @return &Phi; the centrifugal potential (m<sup>2</sup>
249  * s<sup>&minus;2</sup>).
250  *
251  * &Phi; is independent of \e Z, thus \e fZ = 0. This function
252  * NormalGravity::U sums the results of NormalGravity::V0 and
253  * NormalGravity::Phi.
254  **********************************************************************/
255  Math::real Phi(real X, real Y, real& fX, real& fY) const;
256  ///@}
257 
258  /** \name Inspector functions
259  **********************************************************************/
260  ///@{
261  /**
262  * @return true if the object has been initialized.
263  **********************************************************************/
264  bool Init() const { return _a > 0; }
265 
266  /**
267  * @return \e a the equatorial radius of the ellipsoid (meters). This is
268  * the value used in the constructor.
269  **********************************************************************/
271  { return Init() ? _a : Math::NaN(); }
272 
273  /**
274  * @return \e GM the mass constant of the ellipsoid
275  * (m<sup>3</sup> s<sup>&minus;2</sup>). This is the value used in the
276  * constructor.
277  **********************************************************************/
279  { return Init() ? _GM : Math::NaN(); }
280 
281  /**
282  * @return <i>J</i><sub><i>n</i></sub> the dynamical form factors of the
283  * ellipsoid.
284  *
285  * If \e n = 2 (the default), this is the value of <i>J</i><sub>2</sub>
286  * used in the constructor. Otherwise it is the zonal coefficient of the
287  * Legendre harmonic sum of the normal gravitational potential. Note that
288  * <i>J</i><sub><i>n</i></sub> = 0 if \e n is odd. In most gravity
289  * applications, fully normalized Legendre functions are used and the
290  * corresponding coefficient is <i>C</i><sub><i>n</i>0</sub> =
291  * &minus;<i>J</i><sub><i>n</i></sub> / sqrt(2 \e n + 1).
292  **********************************************************************/
294  { return Init() ? ( n == 2 ? _J2 : Jn(n)) : Math::NaN(); }
295 
296  /**
297  * @return &omega; the angular velocity of the ellipsoid (rad
298  * s<sup>&minus;1</sup>). This is the value used in the constructor.
299  **********************************************************************/
301  { return Init() ? _omega : Math::NaN(); }
302 
303  /**
304  * @return <i>f</i> the flattening of the ellipsoid (\e a &minus; \e b)/\e
305  * a.
306  **********************************************************************/
308  { return Init() ? _f : Math::NaN(); }
309 
310  /**
311  * @return &gamma;<sub>e</sub> the normal gravity at equator (m
312  * s<sup>&minus;2</sup>).
313  **********************************************************************/
315  { return Init() ? _gammae : Math::NaN(); }
316 
317  /**
318  * @return &gamma;<sub>p</sub> the normal gravity at poles (m
319  * s<sup>&minus;2</sup>).
320  **********************************************************************/
322  { return Init() ? _gammap : Math::NaN(); }
323 
324  /**
325  * @return <i>f*</i> the gravity flattening (&gamma;<sub>p</sub> &minus;
326  * &gamma;<sub>e</sub>) / &gamma;<sub>e</sub>.
327  **********************************************************************/
329  { return Init() ? _fstar : Math::NaN(); }
330 
331  /**
332  * @return <i>U</i><sub>0</sub> the constant normal potential for the
333  * surface of the ellipsoid (m<sup>2</sup> s<sup>&minus;2</sup>).
334  **********************************************************************/
336  { return Init() ? _U0 : Math::NaN(); }
337 
338  /**
339  * @return the Geocentric object used by this instance.
340  **********************************************************************/
341  const Geocentric& Earth() const { return _earth; }
342 
343  /**
344  * \deprecated An old name for EquatorialRadius().
345  **********************************************************************/
346  GEOGRAPHICLIB_DEPRECATED("Use EquatorialRadius()")
347  Math::real MajorRadius() const { return EquatorialRadius(); }
348  ///@}
349 
350  /**
351  * A global instantiation of NormalGravity for the WGS84 ellipsoid.
352  **********************************************************************/
353  static const NormalGravity& WGS84();
354 
355  /**
356  * A global instantiation of NormalGravity for the GRS80 ellipsoid.
357  **********************************************************************/
358  static const NormalGravity& GRS80();
359 
360  /**
361  * Compute the flattening from the dynamical form factor.
362  *
363  * @param[in] a equatorial radius (meters).
364  * @param[in] GM mass constant of the ellipsoid
365  * (meters<sup>3</sup>/seconds<sup>2</sup>); this is the product of \e G
366  * the gravitational constant and \e M the mass of the earth (usually
367  * including the mass of the earth's atmosphere).
368  * @param[in] omega the angular velocity (rad s<sup>&minus;1</sup>).
369  * @param[in] J2 the dynamical form factor.
370  * @return \e f the flattening of the ellipsoid.
371  *
372  * This routine requires \e a &gt; 0, \e GM &gt; 0, \e J2 &lt; 1/3 &minus;
373  * <i>omega</i><sup>2</sup><i>a</i><sup>3</sup>/<i>GM</i> 8/(45&pi;). A
374  * NaN is returned if these conditions do not hold. The restriction to
375  * positive \e GM is made because for negative \e GM two solutions are
376  * possible.
377  **********************************************************************/
378  static Math::real J2ToFlattening(real a, real GM, real omega, real J2);
379 
380  /**
381  * Compute the dynamical form factor from the flattening.
382  *
383  * @param[in] a equatorial radius (meters).
384  * @param[in] GM mass constant of the ellipsoid
385  * (meters<sup>3</sup>/seconds<sup>2</sup>); this is the product of \e G
386  * the gravitational constant and \e M the mass of the earth (usually
387  * including the mass of the earth's atmosphere).
388  * @param[in] omega the angular velocity (rad s<sup>&minus;1</sup>).
389  * @param[in] f the flattening of the ellipsoid.
390  * @return \e J2 the dynamical form factor.
391  *
392  * This routine requires \e a &gt; 0, \e GM &ne; 0, \e f &lt; 1. The
393  * values of these parameters are not checked.
394  **********************************************************************/
395  static Math::real FlatteningToJ2(real a, real GM, real omega, real f);
396  };
397 
398 } // namespace GeographicLib
399 
400 #endif // GEOGRAPHICLIB_NORMALGRAVITY_HPP
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:66
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
The normal gravity of the earth.
Math::real Flattening() const
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:76
Math::real SurfacePotential() const
Geocentric coordinates
Definition: Geocentric.hpp:67
Math::real MassConstant() const
Math::real GravityFlattening() const
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
static T NaN()
Definition: Math.cpp:268
Header for GeographicLib::Geocentric class.
Math::real EquatorialGravity() const
Model of the earth&#39;s gravity field.
Math::real AngularVelocity() const
Math::real EquatorialRadius() const
Header for GeographicLib::Constants class.
#define GEOGRAPHICLIB_DEPRECATED(msg)
Definition: Constants.hpp:81
Math::real DynamicalFormFactor(int n=2) const
Math::real PolarGravity() const
const Geocentric & Earth() const