GeographicLib  1.21
GravityModel.hpp
Go to the documentation of this file.
00001 /**
00002  * \file GravityModel.hpp
00003  * \brief Header for GeographicLib::GravityModel 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_GRAVITYMODEL_HPP)
00011 #define GEOGRAPHICLIB_GRAVITYMODEL_HPP \
00012   "$Id: e1a573fb0148fa5bc408b2dbdb096d4cd3091bac $"
00013 
00014 #include <string>
00015 #include <sstream>
00016 #include <vector>
00017 #include <GeographicLib/Constants.hpp>
00018 #include <GeographicLib/NormalGravity.hpp>
00019 #include <GeographicLib/SphericalHarmonic.hpp>
00020 #include <GeographicLib/SphericalHarmonic1.hpp>
00021 
00022 #if defined(_MSC_VER)
00023 // Squelch warnings about dll vs vector
00024 #pragma warning (push)
00025 #pragma warning (disable: 4251)
00026 #endif
00027 
00028 namespace GeographicLib {
00029 
00030   class GravityCircle;
00031 
00032   /**
00033    * \brief Model of the earth's gravity field
00034    *
00035    * Evaluate the earth's gravity field according to a model.  The supported
00036    * models treat only the gravitational field exterior to the mass of the
00037    * earth.  When computing the field at points near (but above) the surface of
00038    * the earth a small correction can be applied to account for the mass of the
00039    * atomsphere above the point in question; see \ref gravityatmos.
00040    * Determining the geoid height entails correcting for the mass of the earth
00041    * above the geoid.  The egm96 and egm2008 include separate correction terms
00042    * to account for this mass.
00043    *
00044    * Definitions and terminology (from Heiskanen and Moritz, Sec 2-13):
00045    * - \e V = gravitational potential;
00046    * - \e Phi = rotational potential;
00047    * - \e W = \e V + \e Phi = \e T + \e U = total potential;
00048    * - <i>V</i><sub>0</sub> = normal gravitation potential;
00049    * - \e U = <i>V</i><sub>0</sub> + \e Phi = total normal potential;
00050    * - \e T = \e W - \e U = \e V - <i>V</i><sub>0</sub> = anomalous or
00051    *   disturbing potential;
00052    * - <b>g</b> = <b>grad</b> \e W = <b>gamma</b> + <b>delta</b>;
00053    * - <b>f</b> = <b>grad</b> \e Phi;
00054    * - <b>Gamma</b> = <b>grad</b> <i>V</i><sub>0</sub>;
00055    * - <b>gamma</b> = <b>grad</b> \e U;
00056    * - <b>delta</b> = <b>grad</b> \e T = gravity disturbance vector
00057    *   = <b>g</b><sub><i>P</i></sub> - <b>gamma</b><sub><i>P</i></sub>;
00058    * - delta \e g = gravity disturbance = \e g<sub><i>P</i></sub> - \e
00059    *   gamma<sub><i>P</i></sub>;
00060    * - Delta <b>g</b> = gravity anomaly vector =
00061    *   <b>g</b><sub><i>P</i></sub> - <b>gamma</b><sub><i>Q</i></sub>; here the
00062    *   line \e PQ is perpendicular to ellipsoid and the potential at \e P
00063    *   equals the normal potential at \e Q;
00064    * - Delta \e g = gravity anomaly = \e g<sub><i>P</i></sub> - \e
00065    *   gamma<sub><i>Q</i></sub>;
00066    * - (\e xi, \e eta) deflection of the vertical, the difference in
00067    *   directions of <b>g</b><sub><i>P</i></sub> and
00068    *   <b>gamma</b><sub><i>Q</i></sub>, \e xi = NS, \e eta = EW.
00069    * - \e X, \e Y, \e Z, geocentric coordinates;
00070    * - \e x, \e y, \e z, local cartesian coordinates used to denote the east,
00071    *   north and up directions.
00072    *
00073    * See \ref gravity for details of how to install the gravity model and the
00074    * data format.
00075    *
00076    * References:
00077    * - W. A. Heiskanen and H. Moritz, Physical Geodesy (Freeman, San
00078    *   Francisco, 1967).
00079    *
00080    * Example of use:
00081    * \include example-GravityModel.cpp
00082    *
00083    * <a href="Gravity.1.html">Gravity</a> is a command-line utility providing
00084    * access to the functionality of GravityModel and GravityCircle.
00085    **********************************************************************/
00086 
00087   class GEOGRAPHIC_EXPORT GravityModel {
00088   private:
00089     typedef Math::real real;
00090     friend class GravityCircle;
00091     static const int idlength_ = 8;
00092     std::string _name, _dir, _description, _date, _filename, _id;
00093     real _amodel, _GMmodel, _zeta0, _corrmult;
00094     SphericalHarmonic::normalization _norm;
00095     NormalGravity _earth;
00096     std::vector<real> _C, _S, _CC, _CS, _zonal;
00097     real _dzonal0;              // A left over contribution to _zonal.
00098     SphericalHarmonic _gravitational;
00099     SphericalHarmonic1 _disturbing;
00100     SphericalHarmonic _correction;
00101     void ReadMetadata(const std::string& name);
00102     Math::real InternalT(real X, real Y, real Z,
00103                          real& deltaX, real& deltaY, real& deltaZ,
00104                          bool gradp, bool correct) const throw();
00105     enum captype {
00106       CAP_NONE   = 0U,
00107       CAP_G      = 1U<<0,       // implies potentials W and V
00108       CAP_T      = 1U<<1,
00109       CAP_DELTA  = 1U<<2 | CAP_T, // delta implies T?
00110       CAP_C      = 1U<<3,
00111       CAP_GAMMA0 = 1U<<4,
00112       CAP_GAMMA  = 1U<<5,
00113       CAP_ALL    = 0x3FU,
00114     };
00115 
00116   public:
00117 
00118     /**
00119      * Bit masks for the capabilities to be given to the GravityCircle object
00120      * produced by Circle.
00121      **********************************************************************/
00122     enum mask {
00123       /**
00124        * No capabilities.
00125        * @hideinitializer
00126        **********************************************************************/
00127       NONE = 0U,
00128       /**
00129        * Allow calls to GravityCircle::Gravity, GravityCircle::W, and
00130        * GravityCircle::V.
00131        * @hideinitializer
00132        **********************************************************************/
00133       GRAVITY = CAP_G,
00134       /**
00135        * Allow calls to GravityCircle::Disturbance and GravityCircle::T.
00136        * @hideinitializer
00137        **********************************************************************/
00138       DISTURBANCE = CAP_DELTA,
00139       /**
00140        * Allow calls to GravityCircle::T(real lon) (i.e., computing the
00141        * disturbing potential and not the gravity disturbance vector).
00142        * @hideinitializer
00143        **********************************************************************/
00144       DISTURBING_POTENTIAL = CAP_T,
00145       /**
00146        * Allow calls to GravityCircle::SphericalAnomaly.
00147        * @hideinitializer
00148        **********************************************************************/
00149       SPHERICAL_ANOMALY = CAP_DELTA | CAP_GAMMA,
00150       /**
00151        * Allow calls to GravityCircle::GeoidHeight.
00152        * @hideinitializer
00153        **********************************************************************/
00154       GEOID_HEIGHT = CAP_T | CAP_C | CAP_GAMMA0,
00155       /**
00156        * All capabilities.
00157        * @hideinitializer
00158        **********************************************************************/
00159       ALL = CAP_ALL,
00160     };
00161     /** \name Setting up the gravity model
00162      **********************************************************************/
00163     ///@{
00164     /**
00165      * Construct a gravity model.
00166      *
00167      * @param[in] name the name of the model.
00168      * @param[in] path (optional) directory for data file.
00169      *
00170      * A filename is formed by appending ".egm" (World Gravity Model) to the
00171      * name.  If \e path is specified (and is non-empty), then the file is
00172      * loaded from directory, \e path.  Otherwise the path is given by
00173      * DefaultGravityPath().  This may throw an exception because the file does
00174      * not exist, is unreadable, or is corrupt.
00175      *
00176      * This file contains the metadata which specifies the properties of the
00177      * model.  The coefficients for the spherical harmonic sums are obtained
00178      * from a file obtained by appending ".cof" to metadata file (so the
00179      * filename ends in ".egm.cof").
00180      **********************************************************************/
00181     explicit GravityModel(const std::string& name,
00182                           const std::string& path = "");
00183     ///@}
00184 
00185     /** \name Compute gravity in geodetic coordinates
00186      **********************************************************************/
00187     ///@{
00188     /**
00189      * Evaluate the gravity at an arbitrary point above (or below) the
00190      * ellipsoid.
00191      *
00192      * @param[in] lat the geographic latitude (degrees).
00193      * @param[in] lon the geographic longitude (degrees).
00194      * @param[in] h the height above the ellipsoid (meters).
00195      * @param[out] gx the easterly component of the acceleration
00196      *   (m s<sup>-2</sup>).
00197      * @param[out] gy the northerly component of the acceleration
00198      *   (m s<sup>-2</sup>).
00199      * @param[out] gz the upward component of the acceleration
00200      *   (m s<sup>-2</sup>); this is usually negative.
00201      * @return \e W the sum of the gravitational and centrifugal potentials.
00202      *
00203      * The function includes the effects of the earth's rotation.
00204      **********************************************************************/
00205     Math::real Gravity(real lat, real lon, real h,
00206                        real& gx, real& gy, real& gz) const throw();
00207 
00208     /**
00209      * Evaluate the gravity disturbance vector at an arbitrary point above (or
00210      * below) the ellipsoid.
00211      *
00212      * @param[in] lat the geographic latitude (degrees).
00213      * @param[in] lon the geographic longitude (degrees).
00214      * @param[in] h the height above the ellipsoid (meters).
00215      * @param[out] deltax the easterly component of the disturbance vector
00216      *   (m s<sup>-2</sup>).
00217      * @param[out] deltay the northerly component of the disturbance vector
00218      *   (m s<sup>-2</sup>).
00219      * @param[out] deltaz the upward component of the disturbance vector
00220      *   (m s<sup>-2</sup>).
00221      * @return \e T the corresponding disturbing potential.
00222      **********************************************************************/
00223     Math::real Disturbance(real lat, real lon, real h,
00224                            real& deltax, real& deltay, real& deltaz)
00225       const throw();
00226 
00227     /**
00228      * Evaluate the geoid height.
00229      *
00230      * @param[in] lat the geographic latitude (degrees).
00231      * @param[in] lon the geographic longitude (degrees).
00232      * @return \e N the height of the geoid above the ReferenceEllipsoid()
00233      *   (meters).
00234      *
00235      * This calls NormalGravity::U for ReferenceEllipsoid().  Some
00236      * approximations are made in computing the geoid height so that the
00237      * results of the NGA codes are reproduced accurately.  Details are given
00238      * in \ref gravitygeoid.
00239      **********************************************************************/
00240     Math::real GeoidHeight(real lat, real lon) const throw();
00241 
00242     /**
00243      * Evaluate the components of the gravity anomaly vector using the
00244      * spherical approximation.
00245      *
00246      * @param[in] lat the geographic latitude (degrees).
00247      * @param[in] lon the geographic longitude (degrees).
00248      * @param[in] h the height above the ellipsoid (meters).
00249      * @param[out] Dg01 the gravity anomaly (m s<sup>-2</sup>).
00250      * @param[out] xi the northerly component of the deflection of the vertical
00251      *  (degrees).
00252      * @param[out] eta the easterly component of the deflection of the vertical
00253      *  (degrees).
00254      *
00255      * The spherical approximation (see Heiskanen and Moritz, Sec 2-14) is used
00256      * so that the results of the NGA codes are reproduced accurately.
00257      * approximations used here.  Details are given in \ref gravitygeoid.
00258      **********************************************************************/
00259     void SphericalAnomaly(real lat, real lon, real h,
00260                           real& Dg01, real& xi, real& eta) const throw();
00261     ///@}
00262 
00263     /** \name Compute gravity in geocentric coordinates
00264      **********************************************************************/
00265     ///@{
00266     /**
00267      * Evaluate the components of the acceleration due to gravity and the
00268      * centrifugal acceleration in geocentric coordinates.
00269      *
00270      * @param[in] X geocentric coordinate of point (meters).
00271      * @param[in] Y geocentric coordinate of point (meters).
00272      * @param[in] Z geocentric coordinate of point (meters).
00273      * @param[out] gX the \e X component of the acceleration
00274      *   (m s<sup>-2</sup>).
00275      * @param[out] gY the \e Y component of the acceleration
00276      *   (m s<sup>-2</sup>).
00277      * @param[out] gZ the \e Z component of the acceleration
00278      *   (m s<sup>-2</sup>).
00279      * @return \e W = \e V + \e Phi the sum of the gravitational and
00280      *   centrifugal potentials (m<sup>2</sup> s<sup>-2</sup>).
00281      *
00282      * This calls NormalGravity::U for  ReferenceEllipsoid().
00283      **********************************************************************/
00284     Math::real W(real X, real Y, real Z,
00285                  real& gX, real& gY, real& gZ) const throw();
00286 
00287     /**
00288      * Evaluate the components of the acceleration due to gravity in geocentric
00289      * coordinates.
00290      *
00291      * @param[in] X geocentric coordinate of point (meters).
00292      * @param[in] Y geocentric coordinate of point (meters).
00293      * @param[in] Z geocentric coordinate of point (meters).
00294      * @param[out] GX the \e X component of the acceleration
00295      *   (m s<sup>-2</sup>).
00296      * @param[out] GY the \e Y component of the acceleration
00297      *   (m s<sup>-2</sup>).
00298      * @param[out] GZ the \e Z component of the acceleration
00299      *   (m s<sup>-2</sup>).
00300      * @return \e V = \e W - \e Phi the gravitational potential
00301      *   (m<sup>2</sup> s<sup>-2</sup>).
00302      **********************************************************************/
00303     Math::real V(real X, real Y, real Z,
00304                  real& GX, real& GY, real& GZ) const throw();
00305 
00306     /**
00307      * Evaluate the components of the gravity disturbance in geocentric
00308      * coordinates.
00309      *
00310      * @param[in] X geocentric coordinate of point (meters).
00311      * @param[in] Y geocentric coordinate of point (meters).
00312      * @param[in] Z geocentric coordinate of point (meters).
00313      * @param[out] deltaX the \e X component of the gravity disturbance
00314      *   (m s<sup>-2</sup>).
00315      * @param[out] deltaY the \e Y component of the gravity disturbance
00316      *   (m s<sup>-2</sup>).
00317      * @param[out] deltaZ the \e Z component of the gravity disturbance
00318      *   (m s<sup>-2</sup>).
00319      * @return \e T = \e W - \e U the disturbing potential (also called the
00320      *   anomalous potential) (m<sup>2</sup> s<sup>-2</sup>).
00321      **********************************************************************/
00322     Math::real T(real X, real Y, real Z,
00323                  real& deltaX, real& deltaY, real& deltaZ) const throw()
00324     { return InternalT(X, Y, Z, deltaX, deltaY, deltaZ, true, true); }
00325 
00326     /**
00327      * Evaluate disturbing potential in geocentric coordinates.
00328      *
00329      * @param[in] X geocentric coordinate of point (meters).
00330      * @param[in] Y geocentric coordinate of point (meters).
00331      * @param[in] Z geocentric coordinate of point (meters).
00332      * @return \e T = \e W - \e U the disturbing potential (also called the
00333      *   anomalous potential) (m<sup>2</sup> s<sup>-2</sup>).
00334      **********************************************************************/
00335     Math::real T(real X, real Y, real Z) const throw() {
00336       real dummy;
00337       return InternalT(X, Y, Z, dummy, dummy, dummy, false, true);
00338     }
00339 
00340     /**
00341      * Evaluate the components of the acceleration due to normal gravity and the
00342      * centrifugal acceleration in geocentric coordinates.
00343      *
00344      * @param[in] X geocentric coordinate of point (meters).
00345      * @param[in] Y geocentric coordinate of point (meters).
00346      * @param[in] Z geocentric coordinate of point (meters).
00347      * @param[out] gammaX the \e X component of the normal acceleration
00348      *   (m s<sup>-2</sup>).
00349      * @param[out] gammaY the \e Y component of the normal acceleration
00350      *   (m s<sup>-2</sup>).
00351      * @param[out] gammaZ the \e Z component of the normal acceleration
00352      *   (m s<sup>-2</sup>).
00353      * @return \e U = <i>V</i><sub>0</sub> + \e Phi the sum of the
00354      *   normal gravitational and centrifugal potentials
00355      *   (m<sup>2</sup> s<sup>-2</sup>).
00356      *
00357      * This calls NormalGravity::U for  ReferenceEllipsoid().
00358      **********************************************************************/
00359     Math::real U(real X, real Y, real Z,
00360                  real& gammaX, real& gammaY, real& gammaZ) const throw()
00361     { return _earth.U(X, Y, Z, gammaX, gammaY, gammaZ); }
00362 
00363     /**
00364      * Evaluate the centrifugal acceleration in geocentric coordinates.
00365      *
00366      * @param[in] X geocentric coordinate of point (meters).
00367      * @param[in] Y geocentric coordinate of point (meters).
00368      * @param[out] fX the \e X component of the centrifugal acceleration
00369      *   (m s<sup>-2</sup>).
00370      * @param[out] fY the \e Y component of the centrifugal acceleration
00371      *   (m s<sup>-2</sup>).
00372      * @return \e Phi the centrifugal potential (m<sup>2</sup> s<sup>-2</sup>).
00373      *
00374      * This calls NormalGravity::Phi for  ReferenceEllipsoid().
00375      **********************************************************************/
00376     Math::real Phi(real X, real Y, real& fX, real& fY) const throw()
00377     { return _earth.Phi(X, Y, fX, fY); }
00378     ///@}
00379 
00380     /** \name Compute gravity on a circle of constant latitude
00381      **********************************************************************/
00382     ///@{
00383     /**
00384      * Create a GravityCircle object to allow the gravity field at many points
00385      * with constant \e lat and \e h and varying \e lon to be computed
00386      * efficiently.
00387      *
00388      * @param[in] lat latitude of the point (degrees).
00389      * @param[in] h the height of the point above the ellipsoid (meters).
00390      * @param[in] caps bitor'ed combination of GravityModel::mask values
00391      *   specifying the capabilities of the resulting GravityCircle object.
00392      * @return a GravityCircle object whose member functions computes the
00393      *   gravitational field at a particular values of \e lon.
00394      *
00395      * The GravityModel::mask values are
00396      * - \e caps |= GravityModel::GRAVITY
00397      * - \e caps |= GravityModel::DISTURBANCE
00398      * - \e caps |= GravityModel::DISTURBING_POTENTIAL
00399      * - \e caps |= GravityModel::SPHERICAL_ANOMALY
00400      * - \e caps |= GravityModel::GEOID_HEIGHT
00401      * .
00402      * The default value of \e caps is GravityModel::ALL which turns on all the
00403      * capabilities.  If an unsupported function is invoked, it will return
00404      * NaNs.  Note that GravityModel::GEOID_HEIGHT will only be honored if \e h
00405      * = 0.
00406      *
00407      * If the field at several points on a circle of latitude need to be
00408      * calculated then creating a GravityCircle object and using its member
00409      * functions will be substantially faster, especially for high-degree
00410      * models.  See \ref gravityparallel for an example of using GravityCircle
00411      * (together with OpenMP) to speed up the computation of geoid heights.
00412      **********************************************************************/
00413     GravityCircle Circle(real lat, real h, unsigned caps = ALL) const;
00414     ///@}
00415 
00416     /** \name Inspector functions
00417      **********************************************************************/
00418     ///@{
00419 
00420     /**
00421      * @return the NormalGravity object for the reference ellipsoid.
00422      **********************************************************************/
00423     const NormalGravity& ReferenceEllipsoid() const throw() { return _earth; }
00424 
00425     /**
00426      * @return the description of the gravity model, if available, in the data
00427      *   file; if absent, return "NONE".
00428      **********************************************************************/
00429     const std::string& Description() const throw() { return _description; }
00430 
00431     /**
00432      * @return date of the model; if absent, return "UNKNOWN".
00433      **********************************************************************/
00434     const std::string& DateTime() const throw() { return _date; }
00435 
00436     /**
00437      * @return full file name used to load the gravity model.
00438      **********************************************************************/
00439     const std::string& GravityFile() const throw() { return _filename; }
00440 
00441     /**
00442      * @return "name" used to load the gravity model (from the first argument
00443      *   of the constructor, but this may be overridden by the model file).
00444      **********************************************************************/
00445     const std::string& GravityModelName() const throw() { return _name; }
00446 
00447     /**
00448      * @return directory used to load the gravity model.
00449      **********************************************************************/
00450     const std::string& GravityModelDirectory() const throw() { return _dir; }
00451 
00452     /**
00453      * @return \e a the equatorial radius of the ellipsoid (meters).
00454      **********************************************************************/
00455     Math::real MajorRadius() const throw() { return _earth.MajorRadius(); }
00456 
00457     /**
00458      * @return \e GM the mass constant of the model
00459      *   (m<sup>3</sup> s<sup>-2</sup>); this is the product of \e G the
00460      *   gravitational constant and \e M the mass of the earth (usually
00461      *   including the mass of the earth's atmosphere).
00462      **********************************************************************/
00463     Math::real MassConstant() const throw() { return _GMmodel; }
00464 
00465     /**
00466      * @return \e GM the mass constant of the ReferenceEllipsoid()
00467      *   (m<sup>3</sup> s<sup>-2</sup>).
00468      **********************************************************************/
00469     Math::real ReferenceMassConstant() const throw()
00470     { return _earth.MassConstant(); }
00471 
00472     /**
00473      * @return \e omega the angular velocity of the model and the
00474      *   ReferenceEllipsoid() (rad s<sup>-1</sup>).
00475      **********************************************************************/
00476     Math::real AngularVelocity() const throw()
00477     { return _earth.AngularVelocity(); }
00478 
00479     /**
00480      * @return \e f the flattening of the ellipsoid.
00481      **********************************************************************/
00482     Math::real Flattening() const throw() { return _earth.Flattening(); }
00483     ///@}
00484 
00485     /**
00486      * @return the default path for gravity model data files.
00487      *
00488      * This is the value of the environment variable GRAVITY_PATH, if set;
00489      * otherwise, it is $GEOGRAPHICLIB_DATA/gravity if the environment variable
00490      * GEOGRAPHICLIB_DATA is set; otherwise, it is a compile-time default
00491      * (/usr/local/share/GeographicLib/gravity on non-Windows systems and
00492      * C:/Documents and Settings/All Users/Application
00493      * Data/GeographicLib/gravity on Windows systems).
00494      **********************************************************************/
00495     static std::string DefaultGravityPath();
00496 
00497     /**
00498      * @return the default name for the gravity model.
00499      *
00500      * This is the value of the environment variable GRAVITY_NAME, if set,
00501      * otherwise, it is "egm96".  The GravityModel class does not use
00502      * this function; it is just provided as a convenience for a calling
00503      * program when constructing a GravityModel object.
00504      **********************************************************************/
00505     static std::string DefaultGravityName();
00506   };
00507 
00508 } // namespace GeographicLib
00509 
00510 #if defined(_MSC_VER)
00511 #pragma warning (pop)
00512 #endif
00513 
00514 #endif  // GEOGRAPHICLIB_GRAVITYMODEL_HPP