GeographicLib  1.21
NormalGravity.hpp
Go to the documentation of this file.
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