GeographicLib
1.21
|
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