/** * \file NormalGravity.hpp * \brief Header for GeographicLib::NormalGravity class * * Copyright (c) Charles Karney (2011-2020) and licensed * under the MIT/X11 License. For more information, see * https://geographiclib.sourceforge.io/ **********************************************************************/ #if !defined(GEOGRAPHICLIB_NORMALGRAVITY_HPP) #define GEOGRAPHICLIB_NORMALGRAVITY_HPP 1 #include #include namespace GeographicLib { /** * \brief The normal gravity of the earth * * "Normal" gravity refers to an idealization of the earth which is modeled * as an rotating ellipsoid. The eccentricity of the ellipsoid, the rotation * speed, and the distribution of mass within the ellipsoid are such that the * ellipsoid is a "level ellipoid", a surface of constant potential * (gravitational plus centrifugal). The acceleration due to gravity is * therefore perpendicular to the surface of the ellipsoid. * * Because the distribution of mass within the ellipsoid is unspecified, only * the potential exterior to the ellipsoid is well defined. In this class, * the mass is assumed to be to concentrated on a "focal disc" of radius, * (a2 − b2)1/2, where * \e a is the equatorial radius of the ellipsoid and \e b is its polar * semi-axis. In the case of an oblate ellipsoid, the mass is concentrated * on a "focal rod" of length 2(b2 − * a2)1/2. As a result the potential is well * defined everywhere. * * There is a closed solution to this problem which is implemented here. * Series "approximations" are only used to evaluate certain combinations of * elementary functions where use of the closed expression results in a loss * of accuracy for small arguments due to cancellation of the leading terms. * However these series include sufficient terms to give full machine * precision. * * Although the formulation used in this class applies to ellipsoids with * arbitrary flattening, in practice, its use should be limited to about * b/\e a ∈ [0.01, 100] or \e f ∈ [−99, 0.99]. * * Definitions: * - V0, the gravitational contribution to the normal * potential; * - Φ, the rotational contribution to the normal potential; * - \e U = V0 + Φ, the total potential; * - Γ = ∇V0, the acceleration due to * mass of the earth; * - f = ∇Φ, the centrifugal acceleration; * - γ = ∇\e U = Γ + f, the normal * acceleration; * - \e X, \e Y, \e Z, geocentric coordinates; * - \e x, \e y, \e z, local cartesian coordinates used to denote the east, * north and up directions. * * References: * - C. Somigliana, Teoria generale del campo gravitazionale dell'ellissoide * di rotazione, Mem. Soc. Astron. Ital, 4, 541--599 (1929). * - W. A. Heiskanen and H. Moritz, Physical Geodesy (Freeman, San * Francisco, 1967), Secs. 1-19, 2-7, 2-8 (2-9, 2-10), 6-2 (6-3). * - B. Hofmann-Wellenhof, H. Moritz, Physical Geodesy (Second edition, * Springer, 2006) https://doi.org/10.1007/978-3-211-33545-1 * - H. Moritz, Geodetic Reference System 1980, J. Geodesy 54(3), 395-405 * (1980) https://doi.org/10.1007/BF02521480 * * For more information on normal gravity see \ref normalgravity. * * Example of use: * \include example-NormalGravity.cpp **********************************************************************/ class GEOGRAPHICLIB_EXPORT NormalGravity { private: static const int maxit_ = 20; typedef Math::real real; friend class GravityModel; real _a, _GM, _omega, _f, _J2, _omega2, _aomega2; real _e2, _ep2, _b, _E, _U0, _gammae, _gammap, _Q0, _k, _fstar; Geocentric _earth; static real atanzz(real x, bool alt) { // This routine obeys the identity // atanzz(x, alt) = atanzz(-x/(1+x), !alt) // // Require x >= -1. Best to call with alt, s.t. x >= 0; this results in // a call to atan, instead of asin, or to asinh, instead of atanh. using std::sqrt; using std::abs; using std::atan; using std::asin; using std::asinh; using std::atanh; real z = sqrt(abs(x)); return x == 0 ? 1 : (alt ? (!(x < 0) ? asinh(z) : asin(z)) / sqrt(abs(x) / (1 + x)) : (!(x < 0) ? atan(z) : atanh(z)) / z); } static real atan7series(real x); static real atan5series(real x); static real Qf(real x, bool alt); static real Hf(real x, bool alt); static real QH3f(real x, bool alt); real Jn(int n) const; void Initialize(real a, real GM, real omega, real f_J2, bool geometricp); public: /** \name Setting up the normal gravity **********************************************************************/ ///@{ /** * Constructor for the normal gravity. * * @param[in] a equatorial radius (meters). * @param[in] GM mass constant of the ellipsoid * (meters3/seconds2); this is the product of \e G * the gravitational constant and \e M the mass of the earth (usually * including the mass of the earth's atmosphere). * @param[in] omega the angular velocity (rad s−1). * @param[in] f_J2 either the flattening of the ellipsoid \e f or the * the dynamical form factor \e J2. * @param[out] geometricp if true (the default), then \e f_J2 denotes the * flattening, else it denotes the dynamical form factor \e J2. * @exception if \e a is not positive or if the other parameters do not * obey the restrictions given below. * * The shape of the ellipsoid can be given in one of two ways: * - geometrically (\e geomtricp = true), the ellipsoid is defined by the * flattening \e f = (\e a − \e b) / \e a, where \e a and \e b are * the equatorial radius and the polar semi-axis. The parameters should * obey \e a > 0, \e f < 1. There are no restrictions on \e GM or * \e omega, in particular, \e GM need not be positive. * - physically (\e geometricp = false), the ellipsoid is defined by the * dynamical form factor J2 = (\e C − \e A) / * Ma2, where \e A and \e C are the equatorial and * polar moments of inertia and \e M is the mass of the earth. The * parameters should obey \e a > 0, \e GM > 0 and \e J2 < 1/3 * − (omega2a3/GM) * 8/(45π). There is no restriction on \e omega. **********************************************************************/ NormalGravity(real a, real GM, real omega, real f_J2, bool geometricp = true); /** * A default constructor for the normal gravity. This sets up an * uninitialized object and is used by GravityModel which constructs this * object before it has read in the parameters for the reference ellipsoid. **********************************************************************/ NormalGravity() : _a(-1) {} ///@} /** \name Compute the gravity **********************************************************************/ ///@{ /** * Evaluate the gravity on the surface of the ellipsoid. * * @param[in] lat the geographic latitude (degrees). * @return γ the acceleration due to gravity, positive downwards * (m s−2). * * Due to the axial symmetry of the ellipsoid, the result is independent of * the value of the longitude. This acceleration is perpendicular to the * surface of the ellipsoid. It includes the effects of the earth's * rotation. **********************************************************************/ Math::real SurfaceGravity(real lat) const; /** * Evaluate the gravity at an arbitrary point above (or below) the * ellipsoid. * * @param[in] lat the geographic latitude (degrees). * @param[in] h the height above the ellipsoid (meters). * @param[out] gammay the northerly component of the acceleration * (m s−2). * @param[out] gammaz the upward component of the acceleration * (m s−2); this is usually negative. * @return \e U the corresponding normal potential * (m2 s−2). * * Due to the axial symmetry of the ellipsoid, the result is independent of * the value of the longitude and the easterly component of the * acceleration vanishes, \e gammax = 0. The function includes the effects * of the earth's rotation. When \e h = 0, this function gives \e gammay = * 0 and the returned value matches that of NormalGravity::SurfaceGravity. **********************************************************************/ Math::real Gravity(real lat, real h, real& gammay, real& gammaz) const; /** * Evaluate the components of the acceleration due to gravity and the * centrifugal acceleration in geocentric coordinates. * * @param[in] X geocentric coordinate of point (meters). * @param[in] Y geocentric coordinate of point (meters). * @param[in] Z geocentric coordinate of point (meters). * @param[out] gammaX the \e X component of the acceleration * (m s−2). * @param[out] gammaY the \e Y component of the acceleration * (m s−2). * @param[out] gammaZ the \e Z component of the acceleration * (m s−2). * @return \e U = V0 + Φ the sum of the * gravitational and centrifugal potentials * (m2 s−2). * * The acceleration given by γ = ∇\e U = * ∇V0 + ∇Φ = Γ + f. **********************************************************************/ Math::real U(real X, real Y, real Z, real& gammaX, real& gammaY, real& gammaZ) const; /** * Evaluate the components of the acceleration due to the gravitational * force in geocentric coordinates. * * @param[in] X geocentric coordinate of point (meters). * @param[in] Y geocentric coordinate of point (meters). * @param[in] Z geocentric coordinate of point (meters). * @param[out] GammaX the \e X component of the acceleration due to the * gravitational force (m s−2). * @param[out] GammaY the \e Y component of the acceleration due to the * @param[out] GammaZ the \e Z component of the acceleration due to the * gravitational force (m s−2). * @return V0 the gravitational potential * (m2 s−2). * * This function excludes the centrifugal acceleration and is appropriate * to use for space applications. In terrestrial applications, the * function NormalGravity::U (which includes this effect) should usually be * used. **********************************************************************/ Math::real V0(real X, real Y, real Z, real& GammaX, real& GammaY, real& GammaZ) const; /** * Evaluate the centrifugal acceleration in geocentric coordinates. * * @param[in] X geocentric coordinate of point (meters). * @param[in] Y geocentric coordinate of point (meters). * @param[out] fX the \e X component of the centrifugal acceleration * (m s−2). * @param[out] fY the \e Y component of the centrifugal acceleration * (m s−2). * @return Φ the centrifugal potential (m2 * s−2). * * Φ is independent of \e Z, thus \e fZ = 0. This function * NormalGravity::U sums the results of NormalGravity::V0 and * NormalGravity::Phi. **********************************************************************/ Math::real Phi(real X, real Y, real& fX, real& fY) const; ///@} /** \name Inspector functions **********************************************************************/ ///@{ /** * @return true if the object has been initialized. **********************************************************************/ bool Init() const { return _a > 0; } /** * @return \e a the equatorial radius of the ellipsoid (meters). This is * the value used in the constructor. **********************************************************************/ Math::real EquatorialRadius() const { return Init() ? _a : Math::NaN(); } /** * @return \e GM the mass constant of the ellipsoid * (m3 s−2). This is the value used in the * constructor. **********************************************************************/ Math::real MassConstant() const { return Init() ? _GM : Math::NaN(); } /** * @return Jn the dynamical form factors of the * ellipsoid. * * If \e n = 2 (the default), this is the value of J2 * used in the constructor. Otherwise it is the zonal coefficient of the * Legendre harmonic sum of the normal gravitational potential. Note that * Jn = 0 if \e n is odd. In most gravity * applications, fully normalized Legendre functions are used and the * corresponding coefficient is Cn0 = * −Jn / sqrt(2 \e n + 1). **********************************************************************/ Math::real DynamicalFormFactor(int n = 2) const { return Init() ? ( n == 2 ? _J2 : Jn(n)) : Math::NaN(); } /** * @return ω the angular velocity of the ellipsoid (rad * s−1). This is the value used in the constructor. **********************************************************************/ Math::real AngularVelocity() const { return Init() ? _omega : Math::NaN(); } /** * @return f the flattening of the ellipsoid (\e a − \e b)/\e * a. **********************************************************************/ Math::real Flattening() const { return Init() ? _f : Math::NaN(); } /** * @return γe the normal gravity at equator (m * s−2). **********************************************************************/ Math::real EquatorialGravity() const { return Init() ? _gammae : Math::NaN(); } /** * @return γp the normal gravity at poles (m * s−2). **********************************************************************/ Math::real PolarGravity() const { return Init() ? _gammap : Math::NaN(); } /** * @return f* the gravity flattening (γp − * γe) / γe. **********************************************************************/ Math::real GravityFlattening() const { return Init() ? _fstar : Math::NaN(); } /** * @return U0 the constant normal potential for the * surface of the ellipsoid (m2 s−2). **********************************************************************/ Math::real SurfacePotential() const { return Init() ? _U0 : Math::NaN(); } /** * @return the Geocentric object used by this instance. **********************************************************************/ const Geocentric& Earth() const { return _earth; } /** * \deprecated An old name for EquatorialRadius(). **********************************************************************/ GEOGRAPHICLIB_DEPRECATED("Use EquatorialRadius()") Math::real MajorRadius() const { return EquatorialRadius(); } ///@} /** * A global instantiation of NormalGravity for the WGS84 ellipsoid. **********************************************************************/ static const NormalGravity& WGS84(); /** * A global instantiation of NormalGravity for the GRS80 ellipsoid. **********************************************************************/ static const NormalGravity& GRS80(); /** * Compute the flattening from the dynamical form factor. * * @param[in] a equatorial radius (meters). * @param[in] GM mass constant of the ellipsoid * (meters3/seconds2); this is the product of \e G * the gravitational constant and \e M the mass of the earth (usually * including the mass of the earth's atmosphere). * @param[in] omega the angular velocity (rad s−1). * @param[in] J2 the dynamical form factor. * @return \e f the flattening of the ellipsoid. * * This routine requires \e a > 0, \e GM > 0, \e J2 < 1/3 − * omega2a3/GM 8/(45π). A * NaN is returned if these conditions do not hold. The restriction to * positive \e GM is made because for negative \e GM two solutions are * possible. **********************************************************************/ static Math::real J2ToFlattening(real a, real GM, real omega, real J2); /** * Compute the dynamical form factor from the flattening. * * @param[in] a equatorial radius (meters). * @param[in] GM mass constant of the ellipsoid * (meters3/seconds2); this is the product of \e G * the gravitational constant and \e M the mass of the earth (usually * including the mass of the earth's atmosphere). * @param[in] omega the angular velocity (rad s−1). * @param[in] f the flattening of the ellipsoid. * @return \e J2 the dynamical form factor. * * This routine requires \e a > 0, \e GM ≠ 0, \e f < 1. The * values of these parameters are not checked. **********************************************************************/ static Math::real FlatteningToJ2(real a, real GM, real omega, real f); }; } // namespace GeographicLib #endif // GEOGRAPHICLIB_NORMALGRAVITY_HPP