00001 /** 00002 * \file NormalGravity.hpp 00003 * \brief Header for GeographicLib::NormalGravity class 00004 * 00005 * Copyright (c) Charles Karney (2011-2014) <charles@karney.com> and licensed 00006 * under the MIT/X11 License. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #if !defined(GEOGRAPHICLIB_NORMALGRAVITY_HPP) 00011 #define GEOGRAPHICLIB_NORMALGRAVITY_HPP 1 00012 00013 #include <GeographicLib/Constants.hpp> 00014 #include <GeographicLib/Geocentric.hpp> 00015 00016 namespace GeographicLib { 00017 00018 /** 00019 * \brief The normal gravity of the earth 00020 * 00021 * "Normal" gravity refers to an idealization of the earth which is modeled 00022 * as an rotating ellipsoid. The eccentricity of the ellipsoid, the rotation 00023 * speed, and the distribution of mass within the ellipsoid are such that the 00024 * surface of the ellipsoid is a surface of constant potential (gravitational 00025 * plus centrifugal). The acceleration due to gravity is therefore 00026 * perpendicular to the surface of the ellipsoid. 00027 * 00028 * There is a closed solution to this problem which is implemented here. 00029 * Series "approximations" are only used to evaluate certain combinations of 00030 * elementary functions where use of the closed expression results in a loss 00031 * of accuracy for small arguments due to cancellation of the two leading 00032 * terms. However these series include sufficient terms to give full machine 00033 * precision. 00034 * 00035 * Definitions: 00036 * - <i>V</i><sub>0</sub>, the gravitational contribution to the normal 00037 * potential; 00038 * - Φ, the rotational contribution to the normal potential; 00039 * - \e U = <i>V</i><sub>0</sub> + Φ, the total 00040 * potential; 00041 * - <b>Γ</b> = ∇<i>V</i><sub>0</sub>, the acceleration due to 00042 * mass of the earth; 00043 * - <b>f</b> = ∇Φ, the centrifugal acceleration; 00044 * - <b>γ</b> = ∇\e U = <b>Γ</b> + <b>f</b>, the normal 00045 * acceleration; 00046 * - \e X, \e Y, \e Z, geocentric coordinates; 00047 * - \e x, \e y, \e z, local cartesian coordinates used to denote the east, 00048 * north and up directions. 00049 * 00050 * References: 00051 * - W. A. Heiskanen and H. Moritz, Physical Geodesy (Freeman, San 00052 * Francisco, 1967), Secs. 1-19, 2-7, 2-8 (2-9, 2-10), 6-2 (6-3). 00053 * - H. Moritz, Geodetic Reference System 1980, J. Geodesy 54(3), 395-405 00054 * (1980) http://dx.doi.org/10.1007/BF02521480 00055 * 00056 * Example of use: 00057 * \include example-NormalGravity.cpp 00058 **********************************************************************/ 00059 00060 class GEOGRAPHICLIB_EXPORT NormalGravity { 00061 private: 00062 static const int maxit_ = 10; 00063 typedef Math::real real; 00064 friend class GravityModel; 00065 real _a, _GM, _omega, _f, _J2, _omega2, _aomega2; 00066 real _e2, _ep2, _b, _E, _U0, _gammae, _gammap, _q0, _m, _k, _fstar; 00067 Geocentric _earth; 00068 // (atan(y)-(y-y^3/3))/y^5 (y = sqrt(x)) = 1/5-y/7+y^2/9-y^3/11... 00069 static real atan5(real x); 00070 // (atan(y)-(y-y^3/3+y^5/5))/y^7 (y = sqrt(x)) = -1/7+x/9-x^2/11+x^3/13... 00071 static real atan7(real x); 00072 static real qf(real ep2); 00073 static real dq(real ep2); 00074 static real qpf(real ep2); 00075 real Jn(int n) const; 00076 public: 00077 00078 /** \name Setting up the normal gravity 00079 **********************************************************************/ 00080 ///@{ 00081 /** 00082 * Constructor for the normal gravity. 00083 * 00084 * @param[in] a equatorial radius (meters). 00085 * @param[in] GM mass constant of the ellipsoid 00086 * (meters<sup>3</sup>/seconds<sup>2</sup>); this is the product of \e G 00087 * the gravitational constant and \e M the mass of the earth (usually 00088 * including the mass of the earth's atmosphere). 00089 * @param[in] omega the angular velocity (rad s<sup>−1</sup>). 00090 * @param[in] f the flattening of the ellipsoid. 00091 * @param[in] J2 the dynamical form factor. 00092 * @exception if \e a is not positive or the other constants are 00093 * inconsistent (see below). 00094 * 00095 * If \e omega is non-zero, then exactly one of \e f and \e J2 should be 00096 * positive and this will be used to define the ellipsoid. The shape of 00097 * the ellipsoid can be given in one of two ways: 00098 * - geometrically, the ellipsoid is defined by the flattening \e f = (\e a 00099 * − \e b) / \e a, where \e a and \e b are the equatorial radius 00100 * and the polar semi-axis. 00101 * - physically, the ellipsoid is defined by the dynamical form factor 00102 * <i>J</i><sub>2</sub> = (\e C − \e A) / <i>Ma</i><sup>2</sup>, 00103 * where \e A and \e C are the equatorial and polar moments of inertia 00104 * and \e M is the mass of the earth. 00105 * . 00106 * If \e omega, \e f, and \e J2 are all zero, then the ellipsoid becomes a 00107 * sphere. 00108 **********************************************************************/ 00109 NormalGravity(real a, real GM, real omega, real f, real J2); 00110 00111 /** 00112 * A default constructor for the normal gravity. This sets up an 00113 * uninitialized object and is used by GravityModel which constructs this 00114 * object before it has read in the parameters for the reference ellipsoid. 00115 **********************************************************************/ 00116 NormalGravity() : _a(-1) {} 00117 ///@} 00118 00119 /** \name Compute the gravity 00120 **********************************************************************/ 00121 ///@{ 00122 /** 00123 * Evaluate the gravity on the surface of the ellipsoid. 00124 * 00125 * @param[in] lat the geographic latitude (degrees). 00126 * @return γ the acceleration due to gravity, positive downwards 00127 * (m s<sup>−2</sup>). 00128 * 00129 * Due to the axial symmetry of the ellipsoid, the result is independent of 00130 * the value of the longitude. This acceleration is perpendicular to the 00131 * surface of the ellipsoid. It includes the effects of the earth's 00132 * rotation. 00133 **********************************************************************/ 00134 Math::real SurfaceGravity(real lat) const; 00135 00136 /** 00137 * Evaluate the gravity at an arbitrary point above (or below) the 00138 * ellipsoid. 00139 * 00140 * @param[in] lat the geographic latitude (degrees). 00141 * @param[in] h the height above the ellipsoid (meters). 00142 * @param[out] gammay the northerly component of the acceleration 00143 * (m s<sup>−2</sup>). 00144 * @param[out] gammaz the upward component of the acceleration 00145 * (m s<sup>−2</sup>); this is usually negative. 00146 * @return \e U the corresponding normal potential. 00147 * 00148 * Due to the axial symmetry of the ellipsoid, the result is independent of 00149 * the value of the longitude and the easterly component of the 00150 * acceleration vanishes, \e gammax = 0. The function includes the effects 00151 * of the earth's rotation. When \e h = 0, this function gives \e gammay = 00152 * 0 and the returned value matches that of NormalGravity::SurfaceGravity. 00153 **********************************************************************/ 00154 Math::real Gravity(real lat, real h, real& gammay, real& gammaz) 00155 const; 00156 00157 /** 00158 * Evaluate the components of the acceleration due to gravity and the 00159 * centrifugal acceleration in geocentric coordinates. 00160 * 00161 * @param[in] X geocentric coordinate of point (meters). 00162 * @param[in] Y geocentric coordinate of point (meters). 00163 * @param[in] Z geocentric coordinate of point (meters). 00164 * @param[out] gammaX the \e X component of the acceleration 00165 * (m s<sup>−2</sup>). 00166 * @param[out] gammaY the \e Y component of the acceleration 00167 * (m s<sup>−2</sup>). 00168 * @param[out] gammaZ the \e Z component of the acceleration 00169 * (m s<sup>−2</sup>). 00170 * @return \e U = <i>V</i><sub>0</sub> + Φ the sum of the 00171 * gravitational and centrifugal potentials 00172 * (m<sup>2</sup> s<sup>−2</sup>). 00173 * 00174 * The acceleration given by <b>γ</b> = ∇\e U = 00175 * ∇<i>V</i><sub>0</sub> + ∇Φ = <b>Γ</b> + <b>f</b>. 00176 **********************************************************************/ 00177 Math::real U(real X, real Y, real Z, 00178 real& gammaX, real& gammaY, real& gammaZ) const; 00179 00180 /** 00181 * Evaluate the components of the acceleration due to gravity alone in 00182 * geocentric coordinates. 00183 * 00184 * @param[in] X geocentric coordinate of point (meters). 00185 * @param[in] Y geocentric coordinate of point (meters). 00186 * @param[in] Z geocentric coordinate of point (meters). 00187 * @param[out] GammaX the \e X component of the acceleration due to gravity 00188 * (m s<sup>−2</sup>). 00189 * @param[out] GammaY the \e Y component of the acceleration due to gravity 00190 * (m s<sup>−2</sup>). 00191 * @param[out] GammaZ the \e Z component of the acceleration due to gravity 00192 * (m s<sup>−2</sup>). 00193 * @return <i>V</i><sub>0</sub> the gravitational potential 00194 * (m<sup>2</sup> s<sup>−2</sup>). 00195 * 00196 * This function excludes the centrifugal acceleration and is appropriate 00197 * to use for space applications. In terrestrial applications, the 00198 * function NormalGravity::U (which includes this effect) should usually be 00199 * used. 00200 **********************************************************************/ 00201 Math::real V0(real X, real Y, real Z, 00202 real& GammaX, real& GammaY, real& GammaZ) const; 00203 00204 /** 00205 * Evaluate the centrifugal acceleration in geocentric coordinates. 00206 * 00207 * @param[in] X geocentric coordinate of point (meters). 00208 * @param[in] Y geocentric coordinate of point (meters). 00209 * @param[out] fX the \e X component of the centrifugal acceleration 00210 * (m s<sup>−2</sup>). 00211 * @param[out] fY the \e Y component of the centrifugal acceleration 00212 * (m s<sup>−2</sup>). 00213 * @return Φ the centrifugal potential (m<sup>2</sup> 00214 * s<sup>−2</sup>). 00215 * 00216 * Φ is independent of \e Z, thus \e fZ = 0. This function 00217 * NormalGravity::U sums the results of NormalGravity::V0 and 00218 * NormalGravity::Phi. 00219 **********************************************************************/ 00220 Math::real Phi(real X, real Y, real& fX, real& fY) const; 00221 ///@} 00222 00223 /** \name Inspector functions 00224 **********************************************************************/ 00225 ///@{ 00226 /** 00227 * @return true if the object has been initialized. 00228 **********************************************************************/ 00229 bool Init() const { return _a > 0; } 00230 00231 /** 00232 * @return \e a the equatorial radius of the ellipsoid (meters). This is 00233 * the value used in the constructor. 00234 **********************************************************************/ 00235 Math::real MajorRadius() const 00236 { return Init() ? _a : Math::NaN(); } 00237 00238 /** 00239 * @return \e GM the mass constant of the ellipsoid 00240 * (m<sup>3</sup> s<sup>−2</sup>). This is the value used in the 00241 * constructor. 00242 **********************************************************************/ 00243 Math::real MassConstant() const 00244 { return Init() ? _GM : Math::NaN(); } 00245 00246 /** 00247 * @return <i>J</i><sub><i>n</i></sub> the dynamical form factors of the 00248 * ellipsoid. 00249 * 00250 * If \e n = 2 (the default), this is the value of <i>J</i><sub>2</sub> 00251 * used in the constructor. Otherwise it is the zonal coefficient of the 00252 * Legendre harmonic sum of the normal gravitational potential. Note that 00253 * <i>J</i><sub><i>n</i></sub> = 0 if \e n is odd. In most gravity 00254 * applications, fully normalized Legendre functions are used and the 00255 * corresponding coefficient is <i>C</i><sub><i>n</i>0</sub> = 00256 * −<i>J</i><sub><i>n</i></sub> / sqrt(2 \e n + 1). 00257 **********************************************************************/ 00258 Math::real DynamicalFormFactor(int n = 2) const 00259 { return Init() ? ( n == 2 ? _J2 : Jn(n)) : Math::NaN(); } 00260 00261 /** 00262 * @return ω the angular velocity of the ellipsoid (rad 00263 * s<sup>−1</sup>). This is the value used in the constructor. 00264 **********************************************************************/ 00265 Math::real AngularVelocity() const 00266 { return Init() ? _omega : Math::NaN(); } 00267 00268 /** 00269 * @return <i>f</i> the flattening of the ellipsoid (\e a − \e b)/\e 00270 * a. 00271 **********************************************************************/ 00272 Math::real Flattening() const 00273 { return Init() ? _f : Math::NaN(); } 00274 00275 /** 00276 * @return γ<sub>e</sub> the normal gravity at equator (m 00277 * s<sup>−2</sup>). 00278 **********************************************************************/ 00279 Math::real EquatorialGravity() const 00280 { return Init() ? _gammae : Math::NaN(); } 00281 00282 /** 00283 * @return γ<sub>p</sub> the normal gravity at poles (m 00284 * s<sup>−2</sup>). 00285 **********************************************************************/ 00286 Math::real PolarGravity() const 00287 { return Init() ? _gammap : Math::NaN(); } 00288 00289 /** 00290 * @return <i>f*</i> the gravity flattening (γ<sub>p</sub> − 00291 * γ<sub>e</sub>) / γ<sub>e</sub>. 00292 **********************************************************************/ 00293 Math::real GravityFlattening() const 00294 { return Init() ? _fstar : Math::NaN(); } 00295 00296 /** 00297 * @return <i>U</i><sub>0</sub> the constant normal potential for the 00298 * surface of the ellipsoid (m<sup>2</sup> s<sup>−2</sup>). 00299 **********************************************************************/ 00300 Math::real SurfacePotential() const 00301 { return Init() ? _U0 : Math::NaN(); } 00302 00303 /** 00304 * @return the Geocentric object used by this instance. 00305 **********************************************************************/ 00306 const Geocentric& Earth() const { return _earth; } 00307 ///@} 00308 00309 /** 00310 * A global instantiation of NormalGravity for the WGS84 ellipsoid. 00311 **********************************************************************/ 00312 static const NormalGravity& WGS84(); 00313 00314 /** 00315 * A global instantiation of NormalGravity for the GRS80 ellipsoid. 00316 **********************************************************************/ 00317 static const NormalGravity& GRS80(); 00318 00319 /** 00320 * Compute the flattening from the dynamical form factor. 00321 * 00322 * @param[in] a equatorial radius (meters). 00323 * @param[in] GM mass constant of the ellipsoid 00324 * (meters<sup>3</sup>/seconds<sup>2</sup>); this is the product of \e G 00325 * the gravitational constant and \e M the mass of the earth (usually 00326 * including the mass of the earth's atmosphere). 00327 * @param[in] omega the angular velocity (rad s<sup>−1</sup>). 00328 * @param[in] J2 the dynamical form factor. 00329 * @return \e f the flattening of the ellipsoid. 00330 **********************************************************************/ 00331 static Math::real J2ToFlattening(real a, real GM, real omega, real J2); 00332 00333 /** 00334 * Compute the dynamical form factor from the flattening. 00335 * 00336 * @param[in] a equatorial radius (meters). 00337 * @param[in] GM mass constant of the ellipsoid 00338 * (meters<sup>3</sup>/seconds<sup>2</sup>); this is the product of \e G 00339 * the gravitational constant and \e M the mass of the earth (usually 00340 * including the mass of the earth's atmosphere). 00341 * @param[in] omega the angular velocity (rad s<sup>−1</sup>). 00342 * @param[in] f the flattening of the ellipsoid. 00343 * @return \e J2 the dynamical form factor. 00344 **********************************************************************/ 00345 static Math::real FlatteningToJ2(real a, real GM, real omega, real f); 00346 }; 00347 00348 } // namespace GeographicLib 00349 00350 #endif // GEOGRAPHICLIB_NORMALGRAVITY_HPP