GeographicLib  1.21
GravityCircle.hpp
Go to the documentation of this file.
00001 /**
00002  * \file GravityCircle.hpp
00003  * \brief Header for GeographicLib::GravityCircle 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_GRAVITYCIRCLE_HPP)
00011 #define GEOGRAPHICLIB_GRAVITYCIRCLE_HPP \
00012   "$Id: 6ae0869411185a48c9f55016e6d1fb14e69aaac1 $"
00013 
00014 #include <string>
00015 #include <vector>
00016 #include <GeographicLib/Constants.hpp>
00017 #include <GeographicLib/CircularEngine.hpp>
00018 #include <GeographicLib/GravityModel.hpp>
00019 
00020 namespace GeographicLib {
00021 
00022   /**
00023    * \brief Gravity on a circle of latitude
00024    *
00025    * Evaluate the earth's gravity field on a circle of constant height and
00026    * latitude.  This uses a CircleEngine to pre-evaluate the inner sum of the
00027    * spherical harmonic sum, allowing the values of the field at several
00028    * different longitudes to be evaluated rapidly.
00029    *
00030    * Use GravityModel::Circle to create a GravityCircle object.  (The
00031    * constructor for this class is private.)
00032    *
00033    * See \ref gravityparallel for an example of using GravityCircle (together
00034    * with OpenMP) to speed up the computation of geoid heights.
00035    *
00036    * Example of use:
00037    * \include example-GravityCircle.cpp
00038    *
00039    * <a href="Gravity.1.html">Gravity</a> is a command-line utility providing
00040    * access to the functionality of GravityModel and GravityCircle.
00041    **********************************************************************/
00042 
00043   class GEOGRAPHIC_EXPORT GravityCircle {
00044   private:
00045     typedef Math::real real;
00046     enum mask {
00047       NONE                 = GravityModel::NONE,
00048       GRAVITY              = GravityModel::GRAVITY,
00049       DISTURBANCE          = GravityModel::DISTURBANCE,
00050       DISTURBING_POTENTIAL = GravityModel::DISTURBING_POTENTIAL,
00051       GEOID_HEIGHT         = GravityModel::GEOID_HEIGHT,
00052       SPHERICAL_ANOMALY    = GravityModel::SPHERICAL_ANOMALY,
00053       ALL                  = GravityModel::ALL,
00054     };
00055 
00056     unsigned _caps;
00057     real _a, _f, _lat, _h, _Z, _P, _invR, _cpsi, _spsi,
00058       _cphi, _sphi, _amodel, _GMmodel, _dzonal0,
00059       _corrmult, _gamma0, _gamma, _frot;
00060     CircularEngine _gravitational, _disturbing, _correction;
00061 
00062     GravityCircle(mask caps, real a, real f, real lat, real h,
00063                   real Z, real P, real cphi, real sphi,
00064                   real amodel, real GMmodel, real dzonal0, real corrmult,
00065                   real gamma0, real gamma, real frot,
00066                   const CircularEngine& gravitational,
00067                   const CircularEngine& disturbing,
00068                   const CircularEngine& correction)
00069       : _caps(caps)
00070       , _a(a)
00071       , _f(f)
00072       , _lat(lat)
00073       , _h(h)
00074       , _Z(Z)
00075       , _P(P)
00076       , _invR(Math::hypot(_P, _Z))
00077       , _cpsi(_P * _invR)
00078       , _spsi(_Z * _invR)
00079       , _cphi(cphi)
00080       , _sphi(sphi)
00081       , _amodel(amodel)
00082       , _GMmodel(GMmodel)
00083       , _dzonal0(dzonal0)
00084       , _corrmult(corrmult)
00085       , _gamma0(gamma0)
00086       , _gamma(gamma)
00087       , _frot(frot)
00088       , _gravitational(gravitational)
00089       , _disturbing(disturbing)
00090       , _correction(correction)
00091     {}
00092 
00093     friend class GravityModel; // GravityModel calls the private constructor
00094     Math::real W(real clam, real slam,
00095                  real& gX, real& gY, real& gZ) const throw();
00096     Math::real V(real clam, real slam,
00097                  real& gX, real& gY, real& gZ) const throw();
00098     Math::real InternalT(real clam, real slam,
00099                          real& deltaX, real& deltaY, real& deltaZ,
00100                          bool gradp, bool correct) const throw();
00101   public:
00102     /**
00103      * A default constructor for the normal gravity.  This sets up an
00104      * uninitialized object which can be later replaced by the
00105      * GravityModel::Circle.
00106      **********************************************************************/
00107     GravityCircle() : _a(-1) {}
00108 
00109     /** \name Compute the gravitational field
00110      **********************************************************************/
00111     ///@{
00112     /**
00113      * Evaluate the gravity.
00114      *
00115      * @param[in] lon the geographic longitude (degrees).
00116      * @param[out] gx the easterly component of the acceleration
00117      *   (m s<sup>-2</sup>).
00118      * @param[out] gy the northerly component of the acceleration
00119      *   (m s<sup>-2</sup>).
00120      * @param[out] gz the upward component of the acceleration
00121      *   (m s<sup>-2</sup>); this is usually negative.
00122      * @return \e W the sum of the gravitational and centrifugal potentials.
00123      *
00124      * The function includes the effects of the earth's rotation.
00125      **********************************************************************/
00126     Math::real Gravity(real lon, real& gx, real& gy, real& gz) const throw();
00127 
00128     /**
00129      * Evaluate the gravity disturbance vector.
00130      *
00131      * @param[in] lon the geographic longitude (degrees).
00132      * @param[out] deltax the easterly component of the disturbance vector
00133      *   (m s<sup>-2</sup>).
00134      * @param[out] deltay the northerly component of the disturbance vector
00135      *   (m s<sup>-2</sup>).
00136      * @param[out] deltaz the upward component of the disturbance vector
00137      *   (m s<sup>-2</sup>).
00138      * @return \e T the corresponding disturbing potential.
00139      **********************************************************************/
00140     Math::real Disturbance(real lon, real& deltax, real& deltay, real& deltaz)
00141       const throw();
00142 
00143     /**
00144      * Evaluate the geoid height.
00145      *
00146      * @param[in] lon the geographic longitude (degrees).
00147      * @return \e N the height of the geoid above the reference ellipsoid
00148      *   (meters).
00149      *
00150      * Some approximations are made in computing the geoid height so that the
00151      * results of the NGA codes are reproduced accurately.  Details are given
00152      * in \ref gravitygeoid.
00153      **********************************************************************/
00154     Math::real GeoidHeight(real lon) const throw();
00155 
00156     /**
00157      * Evaluate the components of the gravity anomaly vector using the
00158      * spherical approximation.
00159      *
00160      * @param[in] lon the geographic longitude (degrees).
00161      * @param[out] Dg01 the gravity anomaly (m s<sup>-2</sup>).
00162      * @param[out] xi the northerly component of the deflection of the vertical
00163      *  (degrees).
00164      * @param[out] eta the easterly component of the deflection of the vertical
00165      *  (degrees).
00166      *
00167      * The spherical approximation (see Heiskanen and Moritz, Sec 2-14) is used
00168      * so that the results of the NGA codes are reproduced accurately.
00169      * approximations used here.  Details are given in \ref gravitygeoid.
00170      **********************************************************************/
00171     void SphericalAnomaly(real lon, real& Dg01, real& xi, real& eta)
00172       const throw();
00173 
00174     /**
00175      * Evaluate the components of the acceleration due to gravity and the
00176      * centrifugal acceleration in geocentric coordinates.
00177      *
00178      * @param[in] lon the geographic longitude (degrees).
00179      * @param[out] gX the \e X component of the acceleration
00180      *   (m s<sup>-2</sup>).
00181      * @param[out] gY the \e Y component of the acceleration
00182      *   (m s<sup>-2</sup>).
00183      * @param[out] gZ the \e Z component of the acceleration
00184      *   (m s<sup>-2</sup>).
00185      * @return \e W = \e V + \e Phi the sum of the gravitational and
00186      *   centrifugal potentials (m<sup>2</sup> s<sup>-2</sup>).
00187      **********************************************************************/
00188     Math::real W(real lon, real& gX, real& gY, real& gZ) const throw() {
00189       real clam, slam;
00190       CircularEngine::cossin(lon, clam, slam);
00191       return W(clam, slam, gX, gY, gZ);
00192     }
00193 
00194     /**
00195      * Evaluate the components of the acceleration due to gravity in geocentric
00196      * coordinates.
00197      *
00198      * @param[in] lon the geographic longitude (degrees).
00199      * @param[out] GX the \e X component of the acceleration
00200      *   (m s<sup>-2</sup>).
00201      * @param[out] GY the \e Y component of the acceleration
00202      *   (m s<sup>-2</sup>).
00203      * @param[out] GZ the \e Z component of the acceleration
00204      *   (m s<sup>-2</sup>).
00205      * @return \e V = \e W - \e Phi the gravitational potential
00206      *   (m<sup>2</sup> s<sup>-2</sup>).
00207      **********************************************************************/
00208     Math::real V(real lon, real& GX, real& GY, real& GZ) const throw() {
00209       real clam, slam;
00210       CircularEngine::cossin(lon, clam, slam);
00211       return V(clam, slam, GX, GY, GZ);
00212     }
00213 
00214 
00215     /**
00216      * Evaluate the components of the gravity disturbance in geocentric
00217      * coordinates.
00218      *
00219      * @param[in] lon the geographic longitude (degrees).
00220      * @param[out] deltaX the \e X component of the gravity disturbance
00221      *   (m s<sup>-2</sup>).
00222      * @param[out] deltaY the \e Y component of the gravity disturbance
00223      *   (m s<sup>-2</sup>).
00224      * @param[out] deltaZ the \e Z component of the gravity disturbance
00225      *   (m s<sup>-2</sup>).
00226      * @return \e T = \e W - \e U the disturbing potential (also called the
00227      *   anomalous potential) (m<sup>2</sup> s<sup>-2</sup>).
00228      **********************************************************************/
00229     Math::real T(real lon, real& deltaX, real& deltaY, real& deltaZ)
00230       const throw() {
00231       real clam, slam;
00232       CircularEngine::cossin(lon, clam, slam);
00233       return InternalT(clam, slam, deltaX, deltaY, deltaZ, true, true);
00234     }
00235 
00236     /**
00237      * Evaluate disturbing potential in geocentric coordinates.
00238      *
00239      * @param[in] lon the geographic longitude (degrees).
00240      * @return \e T = \e W - \e U the disturbing potential (also called the
00241      *   anomalous potential) (m<sup>2</sup> s<sup>-2</sup>).
00242      **********************************************************************/
00243     Math::real T(real lon) const throw() {
00244       real clam, slam, dummy;
00245       CircularEngine::cossin(lon, clam, slam);
00246       return InternalT(clam, slam, dummy, dummy, dummy, false, true);
00247     }
00248 
00249     ///@}
00250 
00251     /** \name Inspector functions
00252      **********************************************************************/
00253     ///@{
00254     /**
00255      * @return true if the object has been initialized.
00256      **********************************************************************/
00257     bool Init() const throw() { return _a > 0; }
00258 
00259     /**
00260      * @return \e a the equatorial radius of the ellipsoid (meters).  This is
00261      *   the value inherited from the GravityModel object used in the
00262      *   constructor.
00263      **********************************************************************/
00264     Math::real MajorRadius() const throw()
00265     { return Init() ? _a : Math::NaN<real>(); }
00266 
00267     /**
00268      * @return \e f the flattening of the ellipsoid.  This is the value
00269      *   inherited from the GravityModel object used in the constructor.
00270      **********************************************************************/
00271     Math::real Flattening() const throw()
00272     { return Init() ? _f : Math::NaN<real>(); }
00273 
00274     /**
00275      * @return the latitude of the circle (degrees).
00276      **********************************************************************/
00277     Math::real Latitude() const throw()
00278     { return Init() ? _lat : Math::NaN<real>(); }
00279 
00280     /**
00281      * @return the height of the circle (meters).
00282      **********************************************************************/
00283     Math::real Height() const throw()
00284     { return Init() ? _h : Math::NaN<real>(); }
00285 
00286     /**
00287      * @return \e caps the computational capabilities that this object was
00288      *   constructed with.
00289      **********************************************************************/
00290     unsigned Capabilities() const throw() { return _caps; }
00291 
00292     /**
00293      * @param[in] testcaps a set of bitor'ed GeodesicLine::mask values.
00294      * @return true if the GeodesicLine object has all these capabilities.
00295      **********************************************************************/
00296     bool Capabilities(unsigned testcaps) const throw() {
00297       return (_caps & testcaps) == testcaps;
00298     }
00299     ///@}
00300   };
00301 
00302 } // namespace GeographicLib
00303 
00304 #endif  // GEOGRAPHICLIB_GRAVITYCIRCLE_HPP