GeographicLib
1.49
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-2017) <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
13
#include <
GeographicLib/Constants.hpp
>
14
#include <
GeographicLib/Geocentric.hpp
>
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> − <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> −
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 ∈ [0.01, 100] or \e f ∈ [−99, 0.99].
48
*
49
* Definitions:
50
* - <i>V</i><sub>0</sub>, the gravitational contribution to the normal
51
* potential;
52
* - Φ, the rotational contribution to the normal potential;
53
* - \e U = <i>V</i><sub>0</sub> + Φ, the total potential;
54
* - <b>Γ</b> = ∇<i>V</i><sub>0</sub>, the acceleration due to
55
* mass of the earth;
56
* - <b>f</b> = ∇Φ, the centrifugal acceleration;
57
* - <b>γ</b> = ∇\e U = <b>Γ</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
79
class
GEOGRAPHICLIB_EXPORT
NormalGravity
{
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
real z = sqrt(abs(x));
95
return
x == 0 ? 1 :
96
(alt ?
97
(!(x < 0) ?
Math::asinh
(z) : asin(z)) / sqrt(abs(x) / (1 + x)) :
98
(!(x < 0) ? atan(z) :
Math
::atanh(z)) / z);
99
}
100
static
real
atan7series(
real
x);
101
static
real
atan5series(
real
x);
102
static
real
Qf(
real
x,
bool
alt);
103
static
real
Hf(
real
x,
bool
alt);
104
static
real
QH3f(
real
x,
bool
alt);
105
real
Jn(
int
n)
const
;
106
void
Initialize(
real
a,
real
GM,
real
omega,
real
f_J2,
bool
geometricp);
107
public
:
108
109
/** \name Setting up the normal gravity
110
**********************************************************************/
111
///@{
112
/**
113
* Constructor for the normal gravity.
114
*
115
* @param[in] a equatorial radius (meters).
116
* @param[in] GM mass constant of the ellipsoid
117
* (meters<sup>3</sup>/seconds<sup>2</sup>); this is the product of \e G
118
* the gravitational constant and \e M the mass of the earth (usually
119
* including the mass of the earth's atmosphere).
120
* @param[in] omega the angular velocity (rad s<sup>−1</sup>).
121
* @param[in] f_J2 either the flattening of the ellipsoid \e f or the
122
* the dynamical form factor \e J2.
123
* @param[out] geometricp if true (the default), then \e f_J2 denotes the
124
* flattening, else it denotes the dynamical form factor \e J2.
125
* @exception if \e a is not positive or if the other parameters do not
126
* obey the restrictions given below.
127
*
128
* The shape of the ellipsoid can be given in one of two ways:
129
* - geometrically (\e geomtricp = true), the ellipsoid is defined by the
130
* flattening \e f = (\e a − \e b) / \e a, where \e a and \e b are
131
* the equatorial radius and the polar semi-axis. The parameters should
132
* obey \e a > 0, \e f < 1. There are no restrictions on \e GM or
133
* \e omega, in particular, \e GM need not be positive.
134
* - physically (\e geometricp = false), the ellipsoid is defined by the
135
* dynamical form factor <i>J</i><sub>2</sub> = (\e C − \e A) /
136
* <i>Ma</i><sup>2</sup>, where \e A and \e C are the equatorial and
137
* polar moments of inertia and \e M is the mass of the earth. The
138
* parameters should obey \e a > 0, \e GM > 0 and \e J2 < 1/3
139
* − (<i>omega</i><sup>2</sup><i>a</i><sup>3</sup>/<i>GM</i>)
140
* 8/(45π). There is no restriction on \e omega.
141
**********************************************************************/
142
NormalGravity
(
real
a,
real
GM,
real
omega,
real
f_J2,
143
bool
geometricp =
true
);
144
/**
145
* \deprecated Old constructor for the normal gravity.
146
*
147
* @param[in] a equatorial radius (meters).
148
* @param[in] GM mass constant of the ellipsoid
149
* (meters<sup>3</sup>/seconds<sup>2</sup>); this is the product of \e G
150
* the gravitational constant and \e M the mass of the earth (usually
151
* including the mass of the earth's atmosphere).
152
* @param[in] omega the angular velocity (rad s<sup>−1</sup>).
153
* @param[in] f the flattening of the ellipsoid.
154
* @param[in] J2 the dynamical form factor.
155
* @exception if \e a is not positive or the other constants are
156
* inconsistent (see below).
157
*
158
* If \e omega is non-zero, then exactly one of \e f and \e J2 should be
159
* positive and this will be used to define the ellipsoid. The shape of
160
* the ellipsoid can be given in one of two ways:
161
* - geometrically, the ellipsoid is defined by the flattening \e f = (\e a
162
* − \e b) / \e a, where \e a and \e b are the equatorial radius
163
* and the polar semi-axis.
164
* - physically, the ellipsoid is defined by the dynamical form factor
165
* <i>J</i><sub>2</sub> = (\e C − \e A) / <i>Ma</i><sup>2</sup>,
166
* where \e A and \e C are the equatorial and polar moments of inertia
167
* and \e M is the mass of the earth.
168
* .
169
* If \e omega, \e f, and \e J2 are all zero, then the ellipsoid becomes a
170
* sphere.
171
**********************************************************************/
172
GEOGRAPHICLIB_DEPRECATED
(
"Use new NormalGravity constructor"
)
173
NormalGravity
(
real
a,
real
GM,
real
omega,
real
f,
real
J2);
174
175
/**
176
* A default constructor for the normal gravity. This sets up an
177
* uninitialized object and is used by GravityModel which constructs this
178
* object before it has read in the parameters for the reference ellipsoid.
179
**********************************************************************/
180
NormalGravity
() : _a(-1) {}
181
///@}
182
183
/** \name Compute the gravity
184
**********************************************************************/
185
///@{
186
/**
187
* Evaluate the gravity on the surface of the ellipsoid.
188
*
189
* @param[in] lat the geographic latitude (degrees).
190
* @return γ the acceleration due to gravity, positive downwards
191
* (m s<sup>−2</sup>).
192
*
193
* Due to the axial symmetry of the ellipsoid, the result is independent of
194
* the value of the longitude. This acceleration is perpendicular to the
195
* surface of the ellipsoid. It includes the effects of the earth's
196
* rotation.
197
**********************************************************************/
198
Math::real
SurfaceGravity(
real
lat)
const
;
199
200
/**
201
* Evaluate the gravity at an arbitrary point above (or below) the
202
* ellipsoid.
203
*
204
* @param[in] lat the geographic latitude (degrees).
205
* @param[in] h the height above the ellipsoid (meters).
206
* @param[out] gammay the northerly component of the acceleration
207
* (m s<sup>−2</sup>).
208
* @param[out] gammaz the upward component of the acceleration
209
* (m s<sup>−2</sup>); this is usually negative.
210
* @return \e U the corresponding normal potential
211
* (m<sup>2</sup> s<sup>−2</sup>).
212
*
213
* Due to the axial symmetry of the ellipsoid, the result is independent of
214
* the value of the longitude and the easterly component of the
215
* acceleration vanishes, \e gammax = 0. The function includes the effects
216
* of the earth's rotation. When \e h = 0, this function gives \e gammay =
217
* 0 and the returned value matches that of NormalGravity::SurfaceGravity.
218
**********************************************************************/
219
Math::real
Gravity(
real
lat,
real
h,
real
& gammay,
real
& gammaz)
220
const
;
221
222
/**
223
* Evaluate the components of the acceleration due to gravity and the
224
* centrifugal acceleration in geocentric coordinates.
225
*
226
* @param[in] X geocentric coordinate of point (meters).
227
* @param[in] Y geocentric coordinate of point (meters).
228
* @param[in] Z geocentric coordinate of point (meters).
229
* @param[out] gammaX the \e X component of the acceleration
230
* (m s<sup>−2</sup>).
231
* @param[out] gammaY the \e Y component of the acceleration
232
* (m s<sup>−2</sup>).
233
* @param[out] gammaZ the \e Z component of the acceleration
234
* (m s<sup>−2</sup>).
235
* @return \e U = <i>V</i><sub>0</sub> + Φ the sum of the
236
* gravitational and centrifugal potentials
237
* (m<sup>2</sup> s<sup>−2</sup>).
238
*
239
* The acceleration given by <b>γ</b> = ∇\e U =
240
* ∇<i>V</i><sub>0</sub> + ∇Φ = <b>Γ</b> + <b>f</b>.
241
**********************************************************************/
242
Math::real
U(
real
X,
real
Y,
real
Z,
243
real
& gammaX,
real
& gammaY,
real
& gammaZ)
const
;
244
245
/**
246
* Evaluate the components of the acceleration due to the gravitational
247
* force in geocentric coordinates.
248
*
249
* @param[in] X geocentric coordinate of point (meters).
250
* @param[in] Y geocentric coordinate of point (meters).
251
* @param[in] Z geocentric coordinate of point (meters).
252
* @param[out] GammaX the \e X component of the acceleration due to the
253
* gravitational force (m s<sup>−2</sup>).
254
* @param[out] GammaY the \e Y component of the acceleration due to the
255
* @param[out] GammaZ the \e Z component of the acceleration due to the
256
* gravitational force (m s<sup>−2</sup>).
257
* @return <i>V</i><sub>0</sub> the gravitational potential
258
* (m<sup>2</sup> s<sup>−2</sup>).
259
*
260
* This function excludes the centrifugal acceleration and is appropriate
261
* to use for space applications. In terrestrial applications, the
262
* function NormalGravity::U (which includes this effect) should usually be
263
* used.
264
**********************************************************************/
265
Math::real
V0(
real
X,
real
Y,
real
Z,
266
real
& GammaX,
real
& GammaY,
real
& GammaZ)
const
;
267
268
/**
269
* Evaluate the centrifugal acceleration in geocentric coordinates.
270
*
271
* @param[in] X geocentric coordinate of point (meters).
272
* @param[in] Y geocentric coordinate of point (meters).
273
* @param[out] fX the \e X component of the centrifugal acceleration
274
* (m s<sup>−2</sup>).
275
* @param[out] fY the \e Y component of the centrifugal acceleration
276
* (m s<sup>−2</sup>).
277
* @return Φ the centrifugal potential (m<sup>2</sup>
278
* s<sup>−2</sup>).
279
*
280
* Φ is independent of \e Z, thus \e fZ = 0. This function
281
* NormalGravity::U sums the results of NormalGravity::V0 and
282
* NormalGravity::Phi.
283
**********************************************************************/
284
Math::real
Phi(
real
X,
real
Y,
real
& fX,
real
& fY)
const
;
285
///@}
286
287
/** \name Inspector functions
288
**********************************************************************/
289
///@{
290
/**
291
* @return true if the object has been initialized.
292
**********************************************************************/
293
bool
Init
()
const
{
return
_a > 0; }
294
295
/**
296
* @return \e a the equatorial radius of the ellipsoid (meters). This is
297
* the value used in the constructor.
298
**********************************************************************/
299
Math::real
MajorRadius
()
const
300
{
return
Init() ? _a :
Math::NaN
(); }
301
302
/**
303
* @return \e GM the mass constant of the ellipsoid
304
* (m<sup>3</sup> s<sup>−2</sup>). This is the value used in the
305
* constructor.
306
**********************************************************************/
307
Math::real
MassConstant
()
const
308
{
return
Init() ? _GM :
Math::NaN
(); }
309
310
/**
311
* @return <i>J</i><sub><i>n</i></sub> the dynamical form factors of the
312
* ellipsoid.
313
*
314
* If \e n = 2 (the default), this is the value of <i>J</i><sub>2</sub>
315
* used in the constructor. Otherwise it is the zonal coefficient of the
316
* Legendre harmonic sum of the normal gravitational potential. Note that
317
* <i>J</i><sub><i>n</i></sub> = 0 if \e n is odd. In most gravity
318
* applications, fully normalized Legendre functions are used and the
319
* corresponding coefficient is <i>C</i><sub><i>n</i>0</sub> =
320
* −<i>J</i><sub><i>n</i></sub> / sqrt(2 \e n + 1).
321
**********************************************************************/
322
Math::real
DynamicalFormFactor
(
int
n = 2)
const
323
{
return
Init() ? ( n == 2 ? _J2 : Jn(n)) :
Math::NaN
(); }
324
325
/**
326
* @return ω the angular velocity of the ellipsoid (rad
327
* s<sup>−1</sup>). This is the value used in the constructor.
328
**********************************************************************/
329
Math::real
AngularVelocity
()
const
330
{
return
Init() ? _omega :
Math::NaN
(); }
331
332
/**
333
* @return <i>f</i> the flattening of the ellipsoid (\e a − \e b)/\e
334
* a.
335
**********************************************************************/
336
Math::real
Flattening
()
const
337
{
return
Init() ? _f :
Math::NaN
(); }
338
339
/**
340
* @return γ<sub>e</sub> the normal gravity at equator (m
341
* s<sup>−2</sup>).
342
**********************************************************************/
343
Math::real
EquatorialGravity
()
const
344
{
return
Init() ? _gammae :
Math::NaN
(); }
345
346
/**
347
* @return γ<sub>p</sub> the normal gravity at poles (m
348
* s<sup>−2</sup>).
349
**********************************************************************/
350
Math::real
PolarGravity
()
const
351
{
return
Init() ? _gammap :
Math::NaN
(); }
352
353
/**
354
* @return <i>f*</i> the gravity flattening (γ<sub>p</sub> −
355
* γ<sub>e</sub>) / γ<sub>e</sub>.
356
**********************************************************************/
357
Math::real
GravityFlattening
()
const
358
{
return
Init() ? _fstar :
Math::NaN
(); }
359
360
/**
361
* @return <i>U</i><sub>0</sub> the constant normal potential for the
362
* surface of the ellipsoid (m<sup>2</sup> s<sup>−2</sup>).
363
**********************************************************************/
364
Math::real
SurfacePotential
()
const
365
{
return
Init() ? _U0 :
Math::NaN
(); }
366
367
/**
368
* @return the Geocentric object used by this instance.
369
**********************************************************************/
370
const
Geocentric
&
Earth
()
const
{
return
_earth; }
371
///@}
372
373
/**
374
* A global instantiation of NormalGravity for the WGS84 ellipsoid.
375
**********************************************************************/
376
static
const
NormalGravity
& WGS84();
377
378
/**
379
* A global instantiation of NormalGravity for the GRS80 ellipsoid.
380
**********************************************************************/
381
static
const
NormalGravity
& GRS80();
382
383
/**
384
* Compute the flattening from the dynamical form factor.
385
*
386
* @param[in] a equatorial radius (meters).
387
* @param[in] GM mass constant of the ellipsoid
388
* (meters<sup>3</sup>/seconds<sup>2</sup>); this is the product of \e G
389
* the gravitational constant and \e M the mass of the earth (usually
390
* including the mass of the earth's atmosphere).
391
* @param[in] omega the angular velocity (rad s<sup>−1</sup>).
392
* @param[in] J2 the dynamical form factor.
393
* @return \e f the flattening of the ellipsoid.
394
*
395
* This routine requires \e a > 0, \e GM > 0, \e J2 < 1/3 −
396
* <i>omega</i><sup>2</sup><i>a</i><sup>3</sup>/<i>GM</i> 8/(45π). A
397
* NaN is returned if these conditions do not hold. The restriction to
398
* positive \e GM is made because for negative \e GM two solutions are
399
* possible.
400
**********************************************************************/
401
static
Math::real
J2ToFlattening(
real
a,
real
GM,
real
omega,
real
J2);
402
403
/**
404
* Compute the dynamical form factor from the flattening.
405
*
406
* @param[in] a equatorial radius (meters).
407
* @param[in] GM mass constant of the ellipsoid
408
* (meters<sup>3</sup>/seconds<sup>2</sup>); this is the product of \e G
409
* the gravitational constant and \e M the mass of the earth (usually
410
* including the mass of the earth's atmosphere).
411
* @param[in] omega the angular velocity (rad s<sup>−1</sup>).
412
* @param[in] f the flattening of the ellipsoid.
413
* @return \e J2 the dynamical form factor.
414
*
415
* This routine requires \e a > 0, \e GM ≠ 0, \e f < 1. The
416
* values of these parameters are not checked.
417
**********************************************************************/
418
static
Math::real
FlatteningToJ2(
real
a,
real
GM,
real
omega,
real
f);
419
};
420
421
}
// namespace GeographicLib
422
423
#endif // GEOGRAPHICLIB_NORMALGRAVITY_HPP
real
GeographicLib::Math::real real
Definition:
GeodSolve.cpp:31
GeographicLib::NormalGravity::PolarGravity
Math::real PolarGravity() const
Definition:
NormalGravity.hpp:350
GeographicLib::Math::NaN
static T NaN()
Definition:
Math.hpp:830
GeographicLib
Namespace for GeographicLib.
Definition:
Accumulator.cpp:12
GeographicLib::NormalGravity::GravityFlattening
Math::real GravityFlattening() const
Definition:
NormalGravity.hpp:357
GeographicLib::NormalGravity::MassConstant
Math::real MassConstant() const
Definition:
NormalGravity.hpp:307
GEOGRAPHICLIB_EXPORT
#define GEOGRAPHICLIB_EXPORT
Definition:
Constants.hpp:91
GeographicLib::NormalGravity::Flattening
Math::real Flattening() const
Definition:
NormalGravity.hpp:336
GeographicLib::Geocentric
Geocentric coordinates
Definition:
Geocentric.hpp:67
GeographicLib::NormalGravity::AngularVelocity
Math::real AngularVelocity() const
Definition:
NormalGravity.hpp:329
GeographicLib::Math::real
double real
Definition:
Math.hpp:129
GeographicLib::NormalGravity::SurfacePotential
Math::real SurfacePotential() const
Definition:
NormalGravity.hpp:364
GeographicLib::NormalGravity::Init
bool Init() const
Definition:
NormalGravity.hpp:293
GeographicLib::Math
Mathematical functions needed by GeographicLib.
Definition:
Math.hpp:102
Geocentric.hpp
Header for GeographicLib::Geocentric class.
GeographicLib::NormalGravity
The normal gravity of the earth.
Definition:
NormalGravity.hpp:79
GEOGRAPHICLIB_DEPRECATED
#define GEOGRAPHICLIB_DEPRECATED(msg)
Definition:
Constants.hpp:106
GeographicLib::Math::asinh
static T asinh(T x)
Definition:
Math.hpp:311
Constants.hpp
Header for GeographicLib::Constants class.
GeographicLib::NormalGravity::DynamicalFormFactor
Math::real DynamicalFormFactor(int n=2) const
Definition:
NormalGravity.hpp:322
GeographicLib::NormalGravity::Earth
const Geocentric & Earth() const
Definition:
NormalGravity.hpp:370
GeographicLib::NormalGravity::EquatorialGravity
Math::real EquatorialGravity() const
Definition:
NormalGravity.hpp:343
GeographicLib::GravityModel
Model of the earth's gravity field.
Definition:
GravityModel.hpp:83
GeographicLib::NormalGravity::MajorRadius
Math::real MajorRadius() const
Definition:
NormalGravity.hpp:299
include
GeographicLib
NormalGravity.hpp
Generated by
1.8.20