GeographicLib  1.44
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-2014) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * http://geographiclib.sourceforge.net/
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  * surface of the ellipsoid is a surface of constant potential (gravitational
25  * plus centrifugal). The acceleration due to gravity is therefore
26  * perpendicular to the surface of the ellipsoid.
27  *
28  * There is a closed solution to this problem which is implemented here.
29  * Series "approximations" are only used to evaluate certain combinations of
30  * elementary functions where use of the closed expression results in a loss
31  * of accuracy for small arguments due to cancellation of the two leading
32  * terms. However these series include sufficient terms to give full machine
33  * precision.
34  *
35  * Definitions:
36  * - <i>V</i><sub>0</sub>, the gravitational contribution to the normal
37  * potential;
38  * - &Phi;, the rotational contribution to the normal potential;
39  * - \e U = <i>V</i><sub>0</sub> + &Phi;, the total
40  * potential;
41  * - <b>&Gamma;</b> = &nabla;<i>V</i><sub>0</sub>, the acceleration due to
42  * mass of the earth;
43  * - <b>f</b> = &nabla;&Phi;, the centrifugal acceleration;
44  * - <b>&gamma;</b> = &nabla;\e U = <b>&Gamma;</b> + <b>f</b>, the normal
45  * acceleration;
46  * - \e X, \e Y, \e Z, geocentric coordinates;
47  * - \e x, \e y, \e z, local cartesian coordinates used to denote the east,
48  * north and up directions.
49  *
50  * References:
51  * - W. A. Heiskanen and H. Moritz, Physical Geodesy (Freeman, San
52  * Francisco, 1967), Secs. 1-19, 2-7, 2-8 (2-9, 2-10), 6-2 (6-3).
53  * - H. Moritz, Geodetic Reference System 1980, J. Geodesy 54(3), 395-405
54  * (1980) https://dx.doi.org/10.1007/BF02521480
55  *
56  * Example of use:
57  * \include example-NormalGravity.cpp
58  **********************************************************************/
59 
61  private:
62  static const int maxit_ = 10;
63  typedef Math::real real;
64  friend class GravityModel;
65  real _a, _GM, _omega, _f, _J2, _omega2, _aomega2;
66  real _e2, _ep2, _b, _E, _U0, _gammae, _gammap, _q0, _m, _k, _fstar;
67  Geocentric _earth;
68  // (atan(y)-(y-y^3/3))/y^5 (y = sqrt(x)) = 1/5-y/7+y^2/9-y^3/11...
69  static real atan5(real x);
70  // (atan(y)-(y-y^3/3+y^5/5))/y^7 (y = sqrt(x)) = -1/7+x/9-x^2/11+x^3/13...
71  static real atan7(real x);
72  static real qf(real ep2);
73  static real dq(real ep2);
74  static real qpf(real ep2);
75  real Jn(int n) const;
76  public:
77 
78  /** \name Setting up the normal gravity
79  **********************************************************************/
80  ///@{
81  /**
82  * Constructor for the normal gravity.
83  *
84  * @param[in] a equatorial radius (meters).
85  * @param[in] GM mass constant of the ellipsoid
86  * (meters<sup>3</sup>/seconds<sup>2</sup>); this is the product of \e G
87  * the gravitational constant and \e M the mass of the earth (usually
88  * including the mass of the earth's atmosphere).
89  * @param[in] omega the angular velocity (rad s<sup>&minus;1</sup>).
90  * @param[in] f the flattening of the ellipsoid.
91  * @param[in] J2 the dynamical form factor.
92  * @exception if \e a is not positive or the other constants are
93  * inconsistent (see below).
94  *
95  * If \e omega is non-zero, then exactly one of \e f and \e J2 should be
96  * positive and this will be used to define the ellipsoid. The shape of
97  * the ellipsoid can be given in one of two ways:
98  * - geometrically, the ellipsoid is defined by the flattening \e f = (\e a
99  * &minus; \e b) / \e a, where \e a and \e b are the equatorial radius
100  * and the polar semi-axis.
101  * - physically, the ellipsoid is defined by the dynamical form factor
102  * <i>J</i><sub>2</sub> = (\e C &minus; \e A) / <i>Ma</i><sup>2</sup>,
103  * where \e A and \e C are the equatorial and polar moments of inertia
104  * and \e M is the mass of the earth.
105  * .
106  * If \e omega, \e f, and \e J2 are all zero, then the ellipsoid becomes a
107  * sphere.
108  **********************************************************************/
109  NormalGravity(real a, real GM, real omega, real f, real J2);
110 
111  /**
112  * A default constructor for the normal gravity. This sets up an
113  * uninitialized object and is used by GravityModel which constructs this
114  * object before it has read in the parameters for the reference ellipsoid.
115  **********************************************************************/
116  NormalGravity() : _a(-1) {}
117  ///@}
118 
119  /** \name Compute the gravity
120  **********************************************************************/
121  ///@{
122  /**
123  * Evaluate the gravity on the surface of the ellipsoid.
124  *
125  * @param[in] lat the geographic latitude (degrees).
126  * @return &gamma; the acceleration due to gravity, positive downwards
127  * (m s<sup>&minus;2</sup>).
128  *
129  * Due to the axial symmetry of the ellipsoid, the result is independent of
130  * the value of the longitude. This acceleration is perpendicular to the
131  * surface of the ellipsoid. It includes the effects of the earth's
132  * rotation.
133  **********************************************************************/
134  Math::real SurfaceGravity(real lat) const;
135 
136  /**
137  * Evaluate the gravity at an arbitrary point above (or below) the
138  * ellipsoid.
139  *
140  * @param[in] lat the geographic latitude (degrees).
141  * @param[in] h the height above the ellipsoid (meters).
142  * @param[out] gammay the northerly component of the acceleration
143  * (m s<sup>&minus;2</sup>).
144  * @param[out] gammaz the upward component of the acceleration
145  * (m s<sup>&minus;2</sup>); this is usually negative.
146  * @return \e U the corresponding normal potential.
147  *
148  * Due to the axial symmetry of the ellipsoid, the result is independent of
149  * the value of the longitude and the easterly component of the
150  * acceleration vanishes, \e gammax = 0. The function includes the effects
151  * of the earth's rotation. When \e h = 0, this function gives \e gammay =
152  * 0 and the returned value matches that of NormalGravity::SurfaceGravity.
153  **********************************************************************/
154  Math::real Gravity(real lat, real h, real& gammay, real& gammaz)
155  const;
156 
157  /**
158  * Evaluate the components of the acceleration due to gravity and the
159  * centrifugal acceleration in geocentric coordinates.
160  *
161  * @param[in] X geocentric coordinate of point (meters).
162  * @param[in] Y geocentric coordinate of point (meters).
163  * @param[in] Z geocentric coordinate of point (meters).
164  * @param[out] gammaX the \e X component of the acceleration
165  * (m s<sup>&minus;2</sup>).
166  * @param[out] gammaY the \e Y component of the acceleration
167  * (m s<sup>&minus;2</sup>).
168  * @param[out] gammaZ the \e Z component of the acceleration
169  * (m s<sup>&minus;2</sup>).
170  * @return \e U = <i>V</i><sub>0</sub> + &Phi; the sum of the
171  * gravitational and centrifugal potentials
172  * (m<sup>2</sup> s<sup>&minus;2</sup>).
173  *
174  * The acceleration given by <b>&gamma;</b> = &nabla;\e U =
175  * &nabla;<i>V</i><sub>0</sub> + &nabla;&Phi; = <b>&Gamma;</b> + <b>f</b>.
176  **********************************************************************/
177  Math::real U(real X, real Y, real Z,
178  real& gammaX, real& gammaY, real& gammaZ) const;
179 
180  /**
181  * Evaluate the components of the acceleration due to gravity alone in
182  * geocentric coordinates.
183  *
184  * @param[in] X geocentric coordinate of point (meters).
185  * @param[in] Y geocentric coordinate of point (meters).
186  * @param[in] Z geocentric coordinate of point (meters).
187  * @param[out] GammaX the \e X component of the acceleration due to gravity
188  * (m s<sup>&minus;2</sup>).
189  * @param[out] GammaY the \e Y component of the acceleration due to gravity
190  * (m s<sup>&minus;2</sup>).
191  * @param[out] GammaZ the \e Z component of the acceleration due to gravity
192  * (m s<sup>&minus;2</sup>).
193  * @return <i>V</i><sub>0</sub> the gravitational potential
194  * (m<sup>2</sup> s<sup>&minus;2</sup>).
195  *
196  * This function excludes the centrifugal acceleration and is appropriate
197  * to use for space applications. In terrestrial applications, the
198  * function NormalGravity::U (which includes this effect) should usually be
199  * used.
200  **********************************************************************/
201  Math::real V0(real X, real Y, real Z,
202  real& GammaX, real& GammaY, real& GammaZ) const;
203 
204  /**
205  * Evaluate the centrifugal acceleration in geocentric coordinates.
206  *
207  * @param[in] X geocentric coordinate of point (meters).
208  * @param[in] Y geocentric coordinate of point (meters).
209  * @param[out] fX the \e X component of the centrifugal acceleration
210  * (m s<sup>&minus;2</sup>).
211  * @param[out] fY the \e Y component of the centrifugal acceleration
212  * (m s<sup>&minus;2</sup>).
213  * @return &Phi; the centrifugal potential (m<sup>2</sup>
214  * s<sup>&minus;2</sup>).
215  *
216  * &Phi; is independent of \e Z, thus \e fZ = 0. This function
217  * NormalGravity::U sums the results of NormalGravity::V0 and
218  * NormalGravity::Phi.
219  **********************************************************************/
220  Math::real Phi(real X, real Y, real& fX, real& fY) const;
221  ///@}
222 
223  /** \name Inspector functions
224  **********************************************************************/
225  ///@{
226  /**
227  * @return true if the object has been initialized.
228  **********************************************************************/
229  bool Init() const { return _a > 0; }
230 
231  /**
232  * @return \e a the equatorial radius of the ellipsoid (meters). This is
233  * the value used in the constructor.
234  **********************************************************************/
236  { return Init() ? _a : Math::NaN(); }
237 
238  /**
239  * @return \e GM the mass constant of the ellipsoid
240  * (m<sup>3</sup> s<sup>&minus;2</sup>). This is the value used in the
241  * constructor.
242  **********************************************************************/
244  { return Init() ? _GM : Math::NaN(); }
245 
246  /**
247  * @return <i>J</i><sub><i>n</i></sub> the dynamical form factors of the
248  * ellipsoid.
249  *
250  * If \e n = 2 (the default), this is the value of <i>J</i><sub>2</sub>
251  * used in the constructor. Otherwise it is the zonal coefficient of the
252  * Legendre harmonic sum of the normal gravitational potential. Note that
253  * <i>J</i><sub><i>n</i></sub> = 0 if \e n is odd. In most gravity
254  * applications, fully normalized Legendre functions are used and the
255  * corresponding coefficient is <i>C</i><sub><i>n</i>0</sub> =
256  * &minus;<i>J</i><sub><i>n</i></sub> / sqrt(2 \e n + 1).
257  **********************************************************************/
259  { return Init() ? ( n == 2 ? _J2 : Jn(n)) : Math::NaN(); }
260 
261  /**
262  * @return &omega; the angular velocity of the ellipsoid (rad
263  * s<sup>&minus;1</sup>). This is the value used in the constructor.
264  **********************************************************************/
266  { return Init() ? _omega : Math::NaN(); }
267 
268  /**
269  * @return <i>f</i> the flattening of the ellipsoid (\e a &minus; \e b)/\e
270  * a.
271  **********************************************************************/
273  { return Init() ? _f : Math::NaN(); }
274 
275  /**
276  * @return &gamma;<sub>e</sub> the normal gravity at equator (m
277  * s<sup>&minus;2</sup>).
278  **********************************************************************/
280  { return Init() ? _gammae : Math::NaN(); }
281 
282  /**
283  * @return &gamma;<sub>p</sub> the normal gravity at poles (m
284  * s<sup>&minus;2</sup>).
285  **********************************************************************/
287  { return Init() ? _gammap : Math::NaN(); }
288 
289  /**
290  * @return <i>f*</i> the gravity flattening (&gamma;<sub>p</sub> &minus;
291  * &gamma;<sub>e</sub>) / &gamma;<sub>e</sub>.
292  **********************************************************************/
294  { return Init() ? _fstar : Math::NaN(); }
295 
296  /**
297  * @return <i>U</i><sub>0</sub> the constant normal potential for the
298  * surface of the ellipsoid (m<sup>2</sup> s<sup>&minus;2</sup>).
299  **********************************************************************/
301  { return Init() ? _U0 : Math::NaN(); }
302 
303  /**
304  * @return the Geocentric object used by this instance.
305  **********************************************************************/
306  const Geocentric& Earth() const { return _earth; }
307  ///@}
308 
309  /**
310  * A global instantiation of NormalGravity for the WGS84 ellipsoid.
311  **********************************************************************/
312  static const NormalGravity& WGS84();
313 
314  /**
315  * A global instantiation of NormalGravity for the GRS80 ellipsoid.
316  **********************************************************************/
317  static const NormalGravity& GRS80();
318 
319  /**
320  * Compute the flattening from the dynamical form factor.
321  *
322  * @param[in] a equatorial radius (meters).
323  * @param[in] GM mass constant of the ellipsoid
324  * (meters<sup>3</sup>/seconds<sup>2</sup>); this is the product of \e G
325  * the gravitational constant and \e M the mass of the earth (usually
326  * including the mass of the earth's atmosphere).
327  * @param[in] omega the angular velocity (rad s<sup>&minus;1</sup>).
328  * @param[in] J2 the dynamical form factor.
329  * @return \e f the flattening of the ellipsoid.
330  **********************************************************************/
331  static Math::real J2ToFlattening(real a, real GM, real omega, real J2);
332 
333  /**
334  * Compute the dynamical form factor from the flattening.
335  *
336  * @param[in] a equatorial radius (meters).
337  * @param[in] GM mass constant of the ellipsoid
338  * (meters<sup>3</sup>/seconds<sup>2</sup>); this is the product of \e G
339  * the gravitational constant and \e M the mass of the earth (usually
340  * including the mass of the earth's atmosphere).
341  * @param[in] omega the angular velocity (rad s<sup>&minus;1</sup>).
342  * @param[in] f the flattening of the ellipsoid.
343  * @return \e J2 the dynamical form factor.
344  **********************************************************************/
345  static Math::real FlatteningToJ2(real a, real GM, real omega, real f);
346  };
347 
348 } // namespace GeographicLib
349 
350 #endif // GEOGRAPHICLIB_NORMALGRAVITY_HPP
static T NaN()
Definition: Math.hpp:783
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:90
Math::real PolarGravity() const
GeographicLib::Math::real real
Definition: GeodSolve.cpp:32
The normal gravity of the earth.
Math::real MajorRadius() const
Math::real Flattening() const
Math::real EquatorialGravity() const
const Geocentric & Earth() const
Math::real AngularVelocity() const
Geocentric coordinates
Definition: Geocentric.hpp:67
Math::real SurfacePotential() const
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
Math::real DynamicalFormFactor(int n=2) const
Header for GeographicLib::Geocentric class.
Model of the earth's gravity field.
Header for GeographicLib::Constants class.
Math::real GravityFlattening() const
Math::real MassConstant() const