GeographicLib 2.5
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-2022) <karney@alum.mit.edu> 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
16namespace 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 * https://archive.org/details/HeiskanenMoritz1967PhysicalGeodesy
69 * - B. Hofmann-Wellenhof, H. Moritz, Physical Geodesy (Second edition,
70 * Springer, 2006). https://doi.org/10.1007/978-3-211-33545-1
71 * - H. Moritz, Geodetic Reference System 1980, J. Geodesy 54(3), 395-405
72 * (1980). https://doi.org/10.1007/BF02521480
73 *
74 * For more information on normal gravity see \ref normalgravity.
75 *
76 * Example of use:
77 * \include example-NormalGravity.cpp
78 **********************************************************************/
79
81 private:
82 static const int maxit_ = 20;
83 typedef Math::real real;
84 friend class GravityModel;
85 real _a, _gGM, _omega, _f, _jJ2, _omega2, _aomega2;
86 real _e2, _ep2, _b, _eE, _uU0, _gammae, _gammap, _qQ0, _k, _fstar;
87 Geocentric _earth;
88 static real atanzz(real x, bool alt) {
89 // This routine obeys the identity
90 // atanzz(x, alt) = atanzz(-x/(1+x), !alt)
91 //
92 // Require x >= -1. Best to call with alt, s.t. x >= 0; this results in
93 // a call to atan, instead of asin, or to asinh, instead of atanh.
94 using std::sqrt; using std::fabs; using std::atan; using std::asin;
95 using std::asinh; using std::atanh;
96 real z = sqrt(fabs(x));
97 return x == 0 ? 1 :
98 (alt ?
99 (!(x < 0) ? asinh(z) : asin(z)) / sqrt(fabs(x) / (1 + x)) :
100 (!(x < 0) ? atan(z) : atanh(z)) / z);
101 }
102 static real atan7series(real x);
103 static real atan5series(real x);
104 static real Qf(real x, bool alt);
105 static real Hf(real x, bool alt);
106 static real QH3f(real x, bool alt);
107 real Jn(int n) const;
108 void Initialize(real a, real GM, real omega, real f_J2, bool geometricp);
109 public:
110
111 /** \name Setting up the normal gravity
112 **********************************************************************/
113 ///@{
114 /**
115 * Constructor for the normal gravity.
116 *
117 * @param[in] a equatorial radius (meters).
118 * @param[in] GM mass constant of the ellipsoid
119 * (meters<sup>3</sup>/seconds<sup>2</sup>); this is the product of \e G
120 * the gravitational constant and \e M the mass of the earth (usually
121 * including the mass of the earth's atmosphere).
122 * @param[in] omega the angular velocity (rad s<sup>&minus;1</sup>).
123 * @param[in] f_J2 either the flattening of the ellipsoid \e f or the
124 * the dynamical form factor \e J2.
125 * @param[out] geometricp if true (the default), then \e f_J2 denotes the
126 * flattening, else it denotes the dynamical form factor \e J2.
127 * @exception if \e a is not positive or if the other parameters do not
128 * obey the restrictions given below.
129 *
130 * The shape of the ellipsoid can be given in one of two ways:
131 * - geometrically (\e geomtricp = true), the ellipsoid is defined by the
132 * flattening \e f = (\e a &minus; \e b) / \e a, where \e a and \e b are
133 * the equatorial radius and the polar semi-axis. The parameters should
134 * obey \e a &gt; 0, \e f &lt; 1. There are no restrictions on \e GM or
135 * \e omega, in particular, \e GM need not be positive.
136 * - physically (\e geometricp = false), the ellipsoid is defined by the
137 * dynamical form factor <i>J</i><sub>2</sub> = (\e C &minus; \e A) /
138 * <i>Ma</i><sup>2</sup>, where \e A and \e C are the equatorial and
139 * polar moments of inertia and \e M is the mass of the earth. The
140 * parameters should obey \e a &gt; 0, \e GM &gt; 0 and \e J2 &lt; 1/3
141 * &minus; (<i>omega</i><sup>2</sup><i>a</i><sup>3</sup>/<i>GM</i>)
142 * 8/(45&pi;). There is no restriction on \e omega.
143 **********************************************************************/
144 NormalGravity(real a, real GM, real omega, real f_J2,
145 bool geometricp = true);
146
147 /**
148 * A default constructor for the normal gravity. This sets up an
149 * uninitialized object and is used by GravityModel which constructs this
150 * object before it has read in the parameters for the reference ellipsoid.
151 **********************************************************************/
152 NormalGravity() : _a(-1) {}
153 ///@}
154
155 /** \name Compute the gravity
156 **********************************************************************/
157 ///@{
158 /**
159 * Evaluate the gravity on the surface of the ellipsoid.
160 *
161 * @param[in] lat the geographic latitude (degrees).
162 * @return &gamma; the acceleration due to gravity, positive downwards
163 * (m s<sup>&minus;2</sup>).
164 *
165 * Due to the axial symmetry of the ellipsoid, the result is independent of
166 * the value of the longitude. This acceleration is perpendicular to the
167 * surface of the ellipsoid. It includes the effects of the earth's
168 * rotation.
169 **********************************************************************/
170 Math::real SurfaceGravity(real lat) const;
171
172 /**
173 * Evaluate the gravity at an arbitrary point above (or below) the
174 * ellipsoid.
175 *
176 * @param[in] lat the geographic latitude (degrees).
177 * @param[in] h the height above the ellipsoid (meters).
178 * @param[out] gammay the northerly component of the acceleration
179 * (m s<sup>&minus;2</sup>).
180 * @param[out] gammaz the upward component of the acceleration
181 * (m s<sup>&minus;2</sup>); this is usually negative.
182 * @return \e U the corresponding normal potential
183 * (m<sup>2</sup> s<sup>&minus;2</sup>).
184 *
185 * Due to the axial symmetry of the ellipsoid, the result is independent of
186 * the value of the longitude and the easterly component of the
187 * acceleration vanishes, \e gammax = 0. The function includes the effects
188 * of the earth's rotation. When \e h = 0, this function gives \e gammay =
189 * 0 and the returned value matches that of NormalGravity::SurfaceGravity.
190 **********************************************************************/
191 Math::real Gravity(real lat, real h, real& gammay, real& gammaz)
192 const;
193
194 /**
195 * Evaluate the components of the acceleration due to gravity and the
196 * centrifugal acceleration in geocentric coordinates.
197 *
198 * @param[in] X geocentric coordinate of point (meters).
199 * @param[in] Y geocentric coordinate of point (meters).
200 * @param[in] Z geocentric coordinate of point (meters).
201 * @param[out] gammaX the \e X component of the acceleration
202 * (m s<sup>&minus;2</sup>).
203 * @param[out] gammaY the \e Y component of the acceleration
204 * (m s<sup>&minus;2</sup>).
205 * @param[out] gammaZ the \e Z component of the acceleration
206 * (m s<sup>&minus;2</sup>).
207 * @return \e U = <i>V</i><sub>0</sub> + &Phi; the sum of the
208 * gravitational and centrifugal potentials
209 * (m<sup>2</sup> s<sup>&minus;2</sup>).
210 *
211 * The acceleration given by <b>&gamma;</b> = &nabla;\e U =
212 * &nabla;<i>V</i><sub>0</sub> + &nabla;&Phi; = <b>&Gamma;</b> + <b>f</b>.
213 **********************************************************************/
214 Math::real U(real X, real Y, real Z,
215 real& gammaX, real& gammaY, real& gammaZ) const;
216
217 /**
218 * Evaluate the components of the acceleration due to the gravitational
219 * force in geocentric coordinates.
220 *
221 * @param[in] X geocentric coordinate of point (meters).
222 * @param[in] Y geocentric coordinate of point (meters).
223 * @param[in] Z geocentric coordinate of point (meters).
224 * @param[out] GammaX the \e X component of the acceleration due to the
225 * gravitational force (m s<sup>&minus;2</sup>).
226 * @param[out] GammaY the \e Y component of the acceleration due to the
227 * @param[out] GammaZ the \e Z component of the acceleration due to the
228 * gravitational force (m s<sup>&minus;2</sup>).
229 * @return <i>V</i><sub>0</sub> the gravitational potential
230 * (m<sup>2</sup> s<sup>&minus;2</sup>).
231 *
232 * This function excludes the centrifugal acceleration and is appropriate
233 * to use for space applications. In terrestrial applications, the
234 * function NormalGravity::U (which includes this effect) should usually be
235 * used.
236 **********************************************************************/
237 Math::real V0(real X, real Y, real Z,
238 real& GammaX, real& GammaY, real& GammaZ) const;
239
240 /**
241 * Evaluate the centrifugal acceleration in geocentric coordinates.
242 *
243 * @param[in] X geocentric coordinate of point (meters).
244 * @param[in] Y geocentric coordinate of point (meters).
245 * @param[out] fX the \e X component of the centrifugal acceleration
246 * (m s<sup>&minus;2</sup>).
247 * @param[out] fY the \e Y component of the centrifugal acceleration
248 * (m s<sup>&minus;2</sup>).
249 * @return &Phi; the centrifugal potential (m<sup>2</sup>
250 * s<sup>&minus;2</sup>).
251 *
252 * &Phi; is independent of \e Z, thus \e fZ = 0. This function
253 * NormalGravity::U sums the results of NormalGravity::V0 and
254 * NormalGravity::Phi.
255 **********************************************************************/
256 Math::real Phi(real X, real Y, real& fX, real& fY) const;
257 ///@}
258
259 /** \name Inspector functions
260 **********************************************************************/
261 ///@{
262 /**
263 * @return true if the object has been initialized.
264 **********************************************************************/
265 bool Init() const { return _a > 0; }
266
267 /**
268 * @return \e a the equatorial radius of the ellipsoid (meters). This is
269 * the value used in the constructor.
270 **********************************************************************/
272 { return Init() ? _a : Math::NaN(); }
273
274 /**
275 * @return \e GM the mass constant of the ellipsoid
276 * (m<sup>3</sup> s<sup>&minus;2</sup>). This is the value used in the
277 * constructor.
278 **********************************************************************/
280 { return Init() ? _gGM : Math::NaN(); }
281
282 /**
283 * @return <i>J</i><sub><i>n</i></sub> the dynamical form factors of the
284 * ellipsoid.
285 *
286 * If \e n = 2 (the default), this is the value of <i>J</i><sub>2</sub>
287 * used in the constructor. Otherwise it is the zonal coefficient of the
288 * Legendre harmonic sum of the normal gravitational potential. Note that
289 * <i>J</i><sub><i>n</i></sub> = 0 if \e n is odd. In most gravity
290 * applications, fully normalized Legendre functions are used and the
291 * corresponding coefficient is <i>C</i><sub><i>n</i>0</sub> =
292 * &minus;<i>J</i><sub><i>n</i></sub> / sqrt(2 \e n + 1).
293 **********************************************************************/
295 { return Init() ? ( n == 2 ? _jJ2 : Jn(n)) : Math::NaN(); }
296
297 /**
298 * @return &omega; the angular velocity of the ellipsoid (rad
299 * s<sup>&minus;1</sup>). This is the value used in the constructor.
300 **********************************************************************/
302 { return Init() ? _omega : Math::NaN(); }
303
304 /**
305 * @return <i>f</i> the flattening of the ellipsoid (\e a &minus; \e b)/\e
306 * a.
307 **********************************************************************/
309 { return Init() ? _f : Math::NaN(); }
310
311 /**
312 * @return &gamma;<sub>e</sub> the normal gravity at equator (m
313 * s<sup>&minus;2</sup>).
314 **********************************************************************/
316 { return Init() ? _gammae : Math::NaN(); }
317
318 /**
319 * @return &gamma;<sub>p</sub> the normal gravity at poles (m
320 * s<sup>&minus;2</sup>).
321 **********************************************************************/
323 { return Init() ? _gammap : Math::NaN(); }
324
325 /**
326 * @return <i>f*</i> the gravity flattening (&gamma;<sub>p</sub> &minus;
327 * &gamma;<sub>e</sub>) / &gamma;<sub>e</sub>.
328 **********************************************************************/
330 { return Init() ? _fstar : Math::NaN(); }
331
332 /**
333 * @return <i>U</i><sub>0</sub> the constant normal potential for the
334 * surface of the ellipsoid (m<sup>2</sup> s<sup>&minus;2</sup>).
335 **********************************************************************/
337 { return Init() ? _uU0 : Math::NaN(); }
338
339 /**
340 * @return the Geocentric object used by this instance.
341 **********************************************************************/
342 const Geocentric& Earth() const { return _earth; }
343 ///@}
344
345 /**
346 * A global instantiation of NormalGravity for the WGS84 ellipsoid.
347 **********************************************************************/
348 static const NormalGravity& WGS84();
349
350 /**
351 * A global instantiation of NormalGravity for the GRS80 ellipsoid.
352 **********************************************************************/
353 static const NormalGravity& GRS80();
354
355 /**
356 * Compute the flattening from the dynamical form factor.
357 *
358 * @param[in] a equatorial radius (meters).
359 * @param[in] GM mass constant of the ellipsoid
360 * (meters<sup>3</sup>/seconds<sup>2</sup>); this is the product of \e G
361 * the gravitational constant and \e M the mass of the earth (usually
362 * including the mass of the earth's atmosphere).
363 * @param[in] omega the angular velocity (rad s<sup>&minus;1</sup>).
364 * @param[in] J2 the dynamical form factor.
365 * @return \e f the flattening of the ellipsoid.
366 *
367 * This routine requires \e a &gt; 0, \e GM &gt; 0, \e J2 &lt; 1/3 &minus;
368 * <i>omega</i><sup>2</sup><i>a</i><sup>3</sup>/<i>GM</i> 8/(45&pi;). A
369 * NaN is returned if these conditions do not hold. The restriction to
370 * positive \e GM is made because for negative \e GM two solutions are
371 * possible.
372 **********************************************************************/
373 static Math::real J2ToFlattening(real a, real GM, real omega, real J2);
374
375 /**
376 * Compute the dynamical form factor from the flattening.
377 *
378 * @param[in] a equatorial radius (meters).
379 * @param[in] GM mass constant of the ellipsoid
380 * (meters<sup>3</sup>/seconds<sup>2</sup>); this is the product of \e G
381 * the gravitational constant and \e M the mass of the earth (usually
382 * including the mass of the earth's atmosphere).
383 * @param[in] omega the angular velocity (rad s<sup>&minus;1</sup>).
384 * @param[in] f the flattening of the ellipsoid.
385 * @return \e J2 the dynamical form factor.
386 *
387 * This routine requires \e a &gt; 0, \e GM &ne; 0, \e f &lt; 1. The
388 * values of these parameters are not checked.
389 **********************************************************************/
390 static Math::real FlatteningToJ2(real a, real GM, real omega, real f);
391 };
392
393} // namespace GeographicLib
394
395#endif // GEOGRAPHICLIB_NORMALGRAVITY_HPP
Header for GeographicLib::Constants class.
#define GEOGRAPHICLIB_EXPORT
Definition Constants.hpp:67
Header for GeographicLib::Geocentric class.
GeographicLib::Math::real real
Definition GeodSolve.cpp:28
Geocentric coordinates
Model of the earth's gravity field.
The normal gravity of the earth.
Math::real EquatorialRadius() const
const Geocentric & Earth() const
Math::real PolarGravity() const
Math::real DynamicalFormFactor(int n=2) const
Math::real AngularVelocity() const
Math::real MassConstant() const
Math::real GravityFlattening() const
Math::real SurfacePotential() const
Math::real EquatorialGravity() const
Namespace for GeographicLib.