GeographicLib  1.35
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) <charles@karney.com> and licensed under
6  * 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) http://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  static Math::real qf(real ep2) throw();
69  static Math::real qpf(real ep2) throw();
70  Math::real Jn(int n) const throw();
71  public:
72 
73  /** \name Setting up the normal gravity
74  **********************************************************************/
75  ///@{
76  /**
77  * Constructor for the normal gravity.
78  *
79  * @param[in] a equatorial radius (meters).
80  * @param[in] GM mass constant of the ellipsoid
81  * (meters<sup>3</sup>/seconds<sup>2</sup>); this is the product of \e G
82  * the gravitational constant and \e M the mass of the earth (usually
83  * including the mass of the earth's atmosphere).
84  * @param[in] omega the angular velocity (rad s<sup>&minus;1</sup>).
85  * @param[in] f the flattening of the ellipsoid.
86  * @param[in] J2 dynamical form factor.
87  * @exception if \e a is not positive or the other constants are
88  * inconsistent (see below).
89  *
90  * Exactly one of \e f and \e J2 should be positive and this will be used
91  * to define the ellipsoid. The shape of the ellipsoid can be given in one
92  * of two ways:
93  * - geometrically, the ellipsoid is defined by the flattening \e f = (\e a
94  * &minus; \e b) / \e a, where \e a and \e b are the equatorial radius
95  * and the polar semi-axis.
96  * - physically, the ellipsoid is defined by the dynamical form factor
97  * <i>J</i><sub>2</sub> = (\e C &minus; \e A) / <i>Ma</i><sup>2</sup>,
98  * where \e A and \e C are the equatorial and polar moments of inertia
99  * and \e M is the mass of the earth.
100  **********************************************************************/
101  NormalGravity(real a, real GM, real omega, real f, real J2);
102 
103  /**
104  * A default constructor for the normal gravity. This sets up an
105  * uninitialized object and is used by GravityModel which constructs this
106  * object before it has read in the parameters for the reference ellipsoid.
107  **********************************************************************/
108  NormalGravity() : _a(-1) {}
109  ///@}
110 
111  /** \name Compute the gravity
112  **********************************************************************/
113  ///@{
114  /**
115  * Evaluate the gravity on the surface of the ellipsoid.
116  *
117  * @param[in] lat the geographic latitude (degrees).
118  * @return &gamma; the acceleration due to gravity, positive downwards
119  * (m s<sup>&minus;2</sup>).
120  *
121  * Due to the axial symmetry of the ellipsoid, the result is independent of
122  * the value of the longitude. This acceleration is perpendicular to the
123  * surface of the ellipsoid. It includes the effects of the earth's
124  * rotation.
125  **********************************************************************/
126  Math::real SurfaceGravity(real lat) const throw();
127 
128  /**
129  * Evaluate the gravity at an arbitrary point above (or below) the
130  * ellipsoid.
131  *
132  * @param[in] lat the geographic latitude (degrees).
133  * @param[in] h the height above the ellipsoid (meters).
134  * @param[out] gammay the northerly component of the acceleration
135  * (m s<sup>&minus;2</sup>).
136  * @param[out] gammaz the upward component of the acceleration
137  * (m s<sup>&minus;2</sup>); this is usually negative.
138  * @return \e U the corresponding normal potential.
139  *
140  * Due to the axial symmetry of the ellipsoid, the result is independent of
141  * the value of the longitude and the easterly component of the
142  * acceleration vanishes, \e gammax = 0. The function includes the effects
143  * of the earth's rotation. When \e h = 0, this function gives \e gammay =
144  * 0 and the returned value matches that of NormalGravity::SurfaceGravity.
145  **********************************************************************/
146  Math::real Gravity(real lat, real h, real& gammay, real& gammaz)
147  const throw();
148 
149  /**
150  * Evaluate the components of the acceleration due to gravity and the
151  * centrifugal acceleration in geocentric coordinates.
152  *
153  * @param[in] X geocentric coordinate of point (meters).
154  * @param[in] Y geocentric coordinate of point (meters).
155  * @param[in] Z geocentric coordinate of point (meters).
156  * @param[out] gammaX the \e X component of the acceleration
157  * (m s<sup>&minus;2</sup>).
158  * @param[out] gammaY the \e Y component of the acceleration
159  * (m s<sup>&minus;2</sup>).
160  * @param[out] gammaZ the \e Z component of the acceleration
161  * (m s<sup>&minus;2</sup>).
162  * @return \e U = <i>V</i><sub>0</sub> + &Phi; the sum of the
163  * gravitational and centrifugal potentials
164  * (m<sup>2</sup> s<sup>&minus;2</sup>).
165  *
166  * The acceleration given by <b>&gamma;</b> = &nabla;\e U =
167  * &nabla;<i>V</i><sub>0</sub> + &nabla;&Phi; = <b>&Gamma;</b> + <b>f</b>.
168  **********************************************************************/
169  Math::real U(real X, real Y, real Z,
170  real& gammaX, real& gammaY, real& gammaZ) const throw();
171 
172  /**
173  * Evaluate the components of the acceleration due to gravity alone in
174  * geocentric coordinates.
175  *
176  * @param[in] X geocentric coordinate of point (meters).
177  * @param[in] Y geocentric coordinate of point (meters).
178  * @param[in] Z geocentric coordinate of point (meters).
179  * @param[out] GammaX the \e X component of the acceleration due to gravity
180  * (m s<sup>&minus;2</sup>).
181  * @param[out] GammaY the \e Y component of the acceleration due to gravity
182  * (m s<sup>&minus;2</sup>).
183  * @param[out] GammaZ the \e Z component of the acceleration due to gravity
184  * (m s<sup>&minus;2</sup>).
185  * @return <i>V</i><sub>0</sub> the gravitational potential
186  * (m<sup>2</sup> s<sup>&minus;2</sup>).
187  *
188  * This function excludes the centrifugal acceleration and is appropriate
189  * to use for space applications. In terrestrial applications, the
190  * function NormalGravity::U (which includes this effect) should usually be
191  * used.
192  **********************************************************************/
193  Math::real V0(real X, real Y, real Z,
194  real& GammaX, real& GammaY, real& GammaZ) const throw();
195 
196  /**
197  * Evaluate the centrifugal acceleration in geocentric coordinates.
198  *
199  * @param[in] X geocentric coordinate of point (meters).
200  * @param[in] Y geocentric coordinate of point (meters).
201  * @param[out] fX the \e X component of the centrifugal acceleration
202  * (m s<sup>&minus;2</sup>).
203  * @param[out] fY the \e Y component of the centrifugal acceleration
204  * (m s<sup>&minus;2</sup>).
205  * @return &Phi; the centrifugal potential (m<sup>2</sup>
206  * s<sup>&minus;2</sup>).
207  *
208  * &Phi; is independent of \e Z, thus \e fZ = 0. This function
209  * NormalGravity::U sums the results of NormalGravity::V0 and
210  * NormalGravity::Phi.
211  **********************************************************************/
212  Math::real Phi(real X, real Y, real& fX, real& fY) const throw();
213  ///@}
214 
215  /** \name Inspector functions
216  **********************************************************************/
217  ///@{
218  /**
219  * @return true if the object has been initialized.
220  **********************************************************************/
221  bool Init() const throw() { return _a > 0; }
222 
223  /**
224  * @return \e a the equatorial radius of the ellipsoid (meters). This is
225  * the value used in the constructor.
226  **********************************************************************/
227  Math::real MajorRadius() const throw()
228  { return Init() ? _a : Math::NaN<real>(); }
229 
230  /**
231  * @return \e GM the mass constant of the ellipsoid
232  * (m<sup>3</sup> s<sup>&minus;2</sup>). This is the value used in the
233  * constructor.
234  **********************************************************************/
235  Math::real MassConstant() const throw()
236  { return Init() ? _GM : Math::NaN<real>(); }
237 
238  /**
239  * @return \e J<sub>n</sub> the dynamical form factors of the ellipsoid.
240  *
241  * If \e n = 2 (the default), this is the value of <i>J</i><sub>2</sub>
242  * used in the constructor. Otherwise it is the zonal coefficient of the
243  * Legendre harmonic sum of the normal gravitational potential. Note that
244  * \e J<sub>n</sub> = 0 if \e n is odd. In most gravity applications,
245  * fully normalized Legendre functions are used and the corresponding
246  * coefficient is <i>C</i><sub><i>n</i>0</sub> = &minus;\e J<sub>n</sub> /
247  * sqrt(2 \e n + 1).
248  **********************************************************************/
249  Math::real DynamicalFormFactor(int n = 2) const throw()
250  { return Init() ? ( n == 2 ? _J2 : Jn(n)) : Math::NaN<real>(); }
251 
252  /**
253  * @return &omega; the angular velocity of the ellipsoid (rad
254  * s<sup>&minus;1</sup>). This is the value used in the constructor.
255  **********************************************************************/
256  Math::real AngularVelocity() const throw()
257  { return Init() ? _omega : Math::NaN<real>(); }
258 
259  /**
260  * @return <i>f</i> the flattening of the ellipsoid (\e a &minus; \e b)/\e
261  * a.
262  **********************************************************************/
263  Math::real Flattening() const throw()
264  { return Init() ? _f : Math::NaN<real>(); }
265 
266  /**
267  * @return &gamma;<sub>e</sub> the normal gravity at equator (m
268  * s<sup>&minus;2</sup>).
269  **********************************************************************/
271  { return Init() ? _gammae : Math::NaN<real>(); }
272 
273  /**
274  * @return &gamma;<sub>p</sub> the normal gravity at poles (m
275  * s<sup>&minus;2</sup>).
276  **********************************************************************/
277  Math::real PolarGravity() const throw()
278  { return Init() ? _gammap : Math::NaN<real>(); }
279 
280  /**
281  * @return <i>f*</i> the gravity flattening (&gamma;<sub>p</sub> &minus;
282  * &gamma;<sub>e</sub>) / &gamma;<sub>e</sub>.
283  **********************************************************************/
285  { return Init() ? _fstar : Math::NaN<real>(); }
286 
287  /**
288  * @return <i>U</i><sub>0</sub> the constant normal potential for the
289  * surface of the ellipsoid (m<sup>2</sup> s<sup>&minus;2</sup>).
290  **********************************************************************/
291  Math::real SurfacePotential() const throw()
292  { return Init() ? _U0 : Math::NaN<real>(); }
293 
294  /**
295  * @return the Geocentric object used by this instance.
296  **********************************************************************/
297  const Geocentric& Earth() const throw() { return _earth; }
298  ///@}
299 
300  /**
301  * A global instantiation of NormalGravity for the WGS84 ellipsoid.
302  **********************************************************************/
303  static const NormalGravity WGS84;
304 
305  /**
306  * A global instantiation of NormalGravity for the GRS80 ellipsoid.
307  **********************************************************************/
308  static const NormalGravity GRS80;
309  };
310 
311 } // namespace GeographicLib
312 
313 #endif // GEOGRAPHICLIB_NORMALGRAVITY_HPP
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:52
Math::real PolarGravity() const
GeographicLib::Math::real real
Definition: GeodSolve.cpp:40
The normal gravity of the earth.
static const NormalGravity GRS80
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:61
Math::real SurfacePotential() const
Math::real DynamicalFormFactor(int n=2) const
Header for GeographicLib::Geocentric class.
Model of the earth's gravity field.
Header for GeographicLib::Constants class.
static const NormalGravity WGS84
Math::real GravityFlattening() const
Math::real MassConstant() const