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