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