1 /** 2 * \file GravityModel.hpp 3 * \brief Header for GeographicLib::GravityModel class 4 * 5 * Copyright (c) Charles Karney (2011-2021) <charles@karney.com> and licensed 6 * under the MIT/X11 License. For more information, see 7 * https://geographiclib.sourceforge.io/ 8 **********************************************************************/ 9 10 #if !defined(GEOGRAPHICLIB_GRAVITYMODEL_HPP) 11 #define GEOGRAPHICLIB_GRAVITYMODEL_HPP 1 12 13 #include <GeographicLib/Constants.hpp> 14 #include <GeographicLib/NormalGravity.hpp> 15 #include <GeographicLib/SphericalHarmonic.hpp> 16 #include <GeographicLib/SphericalHarmonic1.hpp> 17 18 #if defined(_MSC_VER) 19 // Squelch warnings about dll vs vector 20 # pragma warning (push) 21 # pragma warning (disable: 4251) 22 #endif 23 24 namespace GeographicLib { 25 26 class GravityCircle; 27 28 /** 29 * \brief Model of the earth's gravity field 30 * 31 * Evaluate the earth's gravity field according to a model. The supported 32 * models treat only the gravitational field exterior to the mass of the 33 * earth. When computing the field at points near (but above) the surface of 34 * the earth a small correction can be applied to account for the mass of the 35 * atmosphere above the point in question; see \ref gravityatmos. 36 * Determining the height of the geoid above the ellipsoid entails correcting 37 * for the mass of the earth above the geoid. The egm96 and egm2008 include 38 * separate correction terms to account for this mass. 39 * 40 * Definitions and terminology (from Heiskanen and Moritz, Sec 2-13): 41 * - \e V = gravitational potential; 42 * - Φ = rotational potential; 43 * - \e W = \e V + Φ = \e T + \e U = total potential; 44 * - <i>V</i><sub>0</sub> = normal gravitation potential; 45 * - \e U = <i>V</i><sub>0</sub> + Φ = total normal potential; 46 * - \e T = \e W − \e U = \e V − <i>V</i><sub>0</sub> = anomalous 47 * or disturbing potential; 48 * - <b>g</b> = ∇\e W = <b>γ</b> + <b>δ</b>; 49 * - <b>f</b> = ∇Φ; 50 * - <b>Γ</b> = ∇<i>V</i><sub>0</sub>; 51 * - <b>γ</b> = ∇\e U; 52 * - <b>δ</b> = ∇\e T = gravity disturbance vector 53 * = <b>g</b><sub><i>P</i></sub> − <b>γ</b><sub><i>P</i></sub>; 54 * - δ\e g = gravity disturbance = <i>g</i><sub><i>P</i></sub> − 55 * γ<sub><i>P</i></sub>; 56 * - Δ<b>g</b> = gravity anomaly vector = <b>g</b><sub><i>P</i></sub> 57 * − <b>γ</b><sub><i>Q</i></sub>; here the line \e PQ is 58 * perpendicular to ellipsoid and the potential at \e P equals the normal 59 * potential at \e Q; 60 * - Δ\e g = gravity anomaly = <i>g</i><sub><i>P</i></sub> − 61 * γ<sub><i>Q</i></sub>; 62 * - (ξ, η) deflection of the vertical, the difference in 63 * directions of <b>g</b><sub><i>P</i></sub> and 64 * <b>γ</b><sub><i>Q</i></sub>, ξ = NS, η = EW. 65 * - \e X, \e Y, \e Z, geocentric coordinates; 66 * - \e x, \e y, \e z, local cartesian coordinates used to denote the east, 67 * north and up directions. 68 * 69 * See \ref gravity for details of how to install the gravity models and the 70 * data format. 71 * 72 * References: 73 * - W. A. Heiskanen and H. Moritz, Physical Geodesy (Freeman, San 74 * Francisco, 1967). 75 * 76 * Example of use: 77 * \include example-GravityModel.cpp 78 * 79 * <a href="Gravity.1.html">Gravity</a> is a command-line utility providing 80 * access to the functionality of GravityModel and GravityCircle. 81 **********************************************************************/ 82 83 class GEOGRAPHICLIB_EXPORT GravityModel { 84 private: 85 typedef Math::real real; 86 friend class GravityCircle; 87 static const int idlength_ = 8; 88 std::string _name, _dir, _description, _date, _filename, _id; 89 real _amodel, _GMmodel, _zeta0, _corrmult; 90 int _nmx, _mmx; 91 SphericalHarmonic::normalization _norm; 92 NormalGravity _earth; 93 std::vector<real> _Cx, _Sx, _CC, _CS, _zonal; 94 real _dzonal0; // A left over contribution to _zonal. 95 SphericalHarmonic _gravitational; 96 SphericalHarmonic1 _disturbing; 97 SphericalHarmonic _correction; 98 void ReadMetadata(const std::string& name); 99 Math::real InternalT(real X, real Y, real Z, 100 real& deltaX, real& deltaY, real& deltaZ, 101 bool gradp, bool correct) const; 102 GravityModel(const GravityModel&) = delete; // copy constructor not allowed 103 // nor copy assignment 104 GravityModel& operator=(const GravityModel&) = delete; 105 106 enum captype { 107 CAP_NONE = 0U, 108 CAP_G = 1U<<0, // implies potentials W and V 109 CAP_T = 1U<<1, 110 CAP_DELTA = 1U<<2 | CAP_T, // delta implies T? 111 CAP_C = 1U<<3, 112 CAP_GAMMA0 = 1U<<4, 113 CAP_GAMMA = 1U<<5, 114 CAP_ALL = 0x3FU, 115 }; 116 117 public: 118 119 /** 120 * Bit masks for the capabilities to be given to the GravityCircle object 121 * produced by Circle. 122 **********************************************************************/ 123 enum mask { 124 /** 125 * No capabilities. 126 * @hideinitializer 127 **********************************************************************/ 128 NONE = 0U, 129 /** 130 * Allow calls to GravityCircle::Gravity, GravityCircle::W, and 131 * GravityCircle::V. 132 * @hideinitializer 133 **********************************************************************/ 134 GRAVITY = CAP_G, 135 /** 136 * Allow calls to GravityCircle::Disturbance and GravityCircle::T. 137 * @hideinitializer 138 **********************************************************************/ 139 DISTURBANCE = CAP_DELTA, 140 /** 141 * Allow calls to GravityCircle::T(real lon) (i.e., computing the 142 * disturbing potential and not the gravity disturbance vector). 143 * @hideinitializer 144 **********************************************************************/ 145 DISTURBING_POTENTIAL = CAP_T, 146 /** 147 * Allow calls to GravityCircle::SphericalAnomaly. 148 * @hideinitializer 149 **********************************************************************/ 150 SPHERICAL_ANOMALY = CAP_DELTA | CAP_GAMMA, 151 /** 152 * Allow calls to GravityCircle::GeoidHeight. 153 * @hideinitializer 154 **********************************************************************/ 155 GEOID_HEIGHT = CAP_T | CAP_C | CAP_GAMMA0, 156 /** 157 * All capabilities. 158 * @hideinitializer 159 **********************************************************************/ 160 ALL = CAP_ALL, 161 }; 162 /** \name Setting up the gravity model 163 **********************************************************************/ 164 ///@{ 165 /** 166 * Construct a gravity model. 167 * 168 * @param[in] name the name of the model. 169 * @param[in] path (optional) directory for data file. 170 * @param[in] Nmax (optional) if non-negative, truncate the degree of the 171 * model this value. 172 * @param[in] Mmax (optional) if non-negative, truncate the order of the 173 * model this value. 174 * @exception GeographicErr if the data file cannot be found, is 175 * unreadable, or is corrupt, or if \e Mmax > \e Nmax. 176 * @exception std::bad_alloc if the memory necessary for storing the model 177 * can't be allocated. 178 * 179 * A filename is formed by appending ".egm" (World Gravity Model) to the 180 * name. If \e path is specified (and is non-empty), then the file is 181 * loaded from directory, \e path. Otherwise the path is given by 182 * DefaultGravityPath(). 183 * 184 * This file contains the metadata which specifies the properties of the 185 * model. The coefficients for the spherical harmonic sums are obtained 186 * from a file obtained by appending ".cof" to metadata file (so the 187 * filename ends in ".egm.cof"). 188 * 189 * If \e Nmax ≥ 0 and \e Mmax < 0, then \e Mmax is set to \e Nmax. 190 * After the model is loaded, the maximum degree and order of the model can 191 * be found by the Degree() and Order() methods. 192 **********************************************************************/ 193 explicit GravityModel(const std::string& name, 194 const std::string& path = "", 195 int Nmax = -1, int Mmax = -1); 196 ///@} 197 198 /** \name Compute gravity in geodetic coordinates 199 **********************************************************************/ 200 ///@{ 201 /** 202 * Evaluate the gravity at an arbitrary point above (or below) the 203 * ellipsoid. 204 * 205 * @param[in] lat the geographic latitude (degrees). 206 * @param[in] lon the geographic longitude (degrees). 207 * @param[in] h the height above the ellipsoid (meters). 208 * @param[out] gx the easterly component of the acceleration 209 * (m s<sup>−2</sup>). 210 * @param[out] gy the northerly component of the acceleration 211 * (m s<sup>−2</sup>). 212 * @param[out] gz the upward component of the acceleration 213 * (m s<sup>−2</sup>); this is usually negative. 214 * @return \e W the sum of the gravitational and centrifugal potentials 215 * (m<sup>2</sup> s<sup>−2</sup>). 216 * 217 * The function includes the effects of the earth's rotation. 218 **********************************************************************/ 219 Math::real Gravity(real lat, real lon, real h, 220 real& gx, real& gy, real& gz) const; 221 222 /** 223 * Evaluate the gravity disturbance vector at an arbitrary point above (or 224 * below) the ellipsoid. 225 * 226 * @param[in] lat the geographic latitude (degrees). 227 * @param[in] lon the geographic longitude (degrees). 228 * @param[in] h the height above the ellipsoid (meters). 229 * @param[out] deltax the easterly component of the disturbance vector 230 * (m s<sup>−2</sup>). 231 * @param[out] deltay the northerly component of the disturbance vector 232 * (m s<sup>−2</sup>). 233 * @param[out] deltaz the upward component of the disturbance vector 234 * (m s<sup>−2</sup>). 235 * @return \e T the corresponding disturbing potential 236 * (m<sup>2</sup> s<sup>−2</sup>). 237 **********************************************************************/ 238 Math::real Disturbance(real lat, real lon, real h, 239 real& deltax, real& deltay, real& deltaz) 240 const; 241 242 /** 243 * Evaluate the geoid height. 244 * 245 * @param[in] lat the geographic latitude (degrees). 246 * @param[in] lon the geographic longitude (degrees). 247 * @return \e N the height of the geoid above the ReferenceEllipsoid() 248 * (meters). 249 * 250 * This calls NormalGravity::U for ReferenceEllipsoid(). Some 251 * approximations are made in computing the geoid height so that the 252 * results of the NGA codes are reproduced accurately. Details are given 253 * in \ref gravitygeoid. 254 **********************************************************************/ 255 Math::real GeoidHeight(real lat, real lon) const; 256 257 /** 258 * Evaluate the components of the gravity anomaly vector using the 259 * spherical approximation. 260 * 261 * @param[in] lat the geographic latitude (degrees). 262 * @param[in] lon the geographic longitude (degrees). 263 * @param[in] h the height above the ellipsoid (meters). 264 * @param[out] Dg01 the gravity anomaly (m s<sup>−2</sup>). 265 * @param[out] xi the northerly component of the deflection of the vertical 266 * (degrees). 267 * @param[out] eta the easterly component of the deflection of the vertical 268 * (degrees). 269 * 270 * The spherical approximation (see Heiskanen and Moritz, Sec 2-14) is used 271 * so that the results of the NGA codes are reproduced accurately. 272 * approximations used here. Details are given in \ref gravitygeoid. 273 **********************************************************************/ 274 void SphericalAnomaly(real lat, real lon, real h, 275 real& Dg01, real& xi, real& eta) const; 276 ///@} 277 278 /** \name Compute gravity in geocentric coordinates 279 **********************************************************************/ 280 ///@{ 281 /** 282 * Evaluate the components of the acceleration due to gravity and the 283 * centrifugal acceleration in geocentric coordinates. 284 * 285 * @param[in] X geocentric coordinate of point (meters). 286 * @param[in] Y geocentric coordinate of point (meters). 287 * @param[in] Z geocentric coordinate of point (meters). 288 * @param[out] gX the \e X component of the acceleration 289 * (m s<sup>−2</sup>). 290 * @param[out] gY the \e Y component of the acceleration 291 * (m s<sup>−2</sup>). 292 * @param[out] gZ the \e Z component of the acceleration 293 * (m s<sup>−2</sup>). 294 * @return \e W = \e V + Φ the sum of the gravitational and 295 * centrifugal potentials (m<sup>2</sup> s<sup>−2</sup>). 296 * 297 * This calls NormalGravity::U for ReferenceEllipsoid(). 298 **********************************************************************/ 299 Math::real W(real X, real Y, real Z, 300 real& gX, real& gY, real& gZ) const; 301 302 /** 303 * Evaluate the components of the acceleration due to gravity in geocentric 304 * coordinates. 305 * 306 * @param[in] X geocentric coordinate of point (meters). 307 * @param[in] Y geocentric coordinate of point (meters). 308 * @param[in] Z geocentric coordinate of point (meters). 309 * @param[out] GX the \e X component of the acceleration 310 * (m s<sup>−2</sup>). 311 * @param[out] GY the \e Y component of the acceleration 312 * (m s<sup>−2</sup>). 313 * @param[out] GZ the \e Z component of the acceleration 314 * (m s<sup>−2</sup>). 315 * @return \e V = \e W - Φ the gravitational potential 316 * (m<sup>2</sup> s<sup>−2</sup>). 317 **********************************************************************/ 318 Math::real V(real X, real Y, real Z, 319 real& GX, real& GY, real& GZ) const; 320 321 /** 322 * Evaluate the components of the gravity disturbance in geocentric 323 * coordinates. 324 * 325 * @param[in] X geocentric coordinate of point (meters). 326 * @param[in] Y geocentric coordinate of point (meters). 327 * @param[in] Z geocentric coordinate of point (meters). 328 * @param[out] deltaX the \e X component of the gravity disturbance 329 * (m s<sup>−2</sup>). 330 * @param[out] deltaY the \e Y component of the gravity disturbance 331 * (m s<sup>−2</sup>). 332 * @param[out] deltaZ the \e Z component of the gravity disturbance 333 * (m s<sup>−2</sup>). 334 * @return \e T = \e W - \e U the disturbing potential (also called the 335 * anomalous potential) (m<sup>2</sup> s<sup>−2</sup>). 336 **********************************************************************/ T(real X,real Y,real Z,real & deltaX,real & deltaY,real & deltaZ) const337 Math::real T(real X, real Y, real Z, 338 real& deltaX, real& deltaY, real& deltaZ) const 339 { return InternalT(X, Y, Z, deltaX, deltaY, deltaZ, true, true); } 340 341 /** 342 * Evaluate disturbing potential in geocentric coordinates. 343 * 344 * @param[in] X geocentric coordinate of point (meters). 345 * @param[in] Y geocentric coordinate of point (meters). 346 * @param[in] Z geocentric coordinate of point (meters). 347 * @return \e T = \e W - \e U the disturbing potential (also called the 348 * anomalous potential) (m<sup>2</sup> s<sup>−2</sup>). 349 **********************************************************************/ T(real X,real Y,real Z) const350 Math::real T(real X, real Y, real Z) const { 351 real dummy; 352 return InternalT(X, Y, Z, dummy, dummy, dummy, false, true); 353 } 354 355 /** 356 * Evaluate the components of the acceleration due to normal gravity and 357 * the centrifugal acceleration in geocentric coordinates. 358 * 359 * @param[in] X geocentric coordinate of point (meters). 360 * @param[in] Y geocentric coordinate of point (meters). 361 * @param[in] Z geocentric coordinate of point (meters). 362 * @param[out] gammaX the \e X component of the normal acceleration 363 * (m s<sup>−2</sup>). 364 * @param[out] gammaY the \e Y component of the normal acceleration 365 * (m s<sup>−2</sup>). 366 * @param[out] gammaZ the \e Z component of the normal acceleration 367 * (m s<sup>−2</sup>). 368 * @return \e U = <i>V</i><sub>0</sub> + Φ the sum of the 369 * normal gravitational and centrifugal potentials 370 * (m<sup>2</sup> s<sup>−2</sup>). 371 * 372 * This calls NormalGravity::U for ReferenceEllipsoid(). 373 **********************************************************************/ U(real X,real Y,real Z,real & gammaX,real & gammaY,real & gammaZ) const374 Math::real U(real X, real Y, real Z, 375 real& gammaX, real& gammaY, real& gammaZ) const 376 { return _earth.U(X, Y, Z, gammaX, gammaY, gammaZ); } 377 378 /** 379 * Evaluate the centrifugal acceleration in geocentric coordinates. 380 * 381 * @param[in] X geocentric coordinate of point (meters). 382 * @param[in] Y geocentric coordinate of point (meters). 383 * @param[out] fX the \e X component of the centrifugal acceleration 384 * (m s<sup>−2</sup>). 385 * @param[out] fY the \e Y component of the centrifugal acceleration 386 * (m s<sup>−2</sup>). 387 * @return Φ the centrifugal potential (m<sup>2</sup> 388 * s<sup>−2</sup>). 389 * 390 * This calls NormalGravity::Phi for ReferenceEllipsoid(). 391 **********************************************************************/ Phi(real X,real Y,real & fX,real & fY) const392 Math::real Phi(real X, real Y, real& fX, real& fY) const 393 { return _earth.Phi(X, Y, fX, fY); } 394 ///@} 395 396 /** \name Compute gravity on a circle of constant latitude 397 **********************************************************************/ 398 ///@{ 399 /** 400 * Create a GravityCircle object to allow the gravity field at many points 401 * with constant \e lat and \e h and varying \e lon to be computed 402 * efficiently. 403 * 404 * @param[in] lat latitude of the point (degrees). 405 * @param[in] h the height of the point above the ellipsoid (meters). 406 * @param[in] caps bitor'ed combination of GravityModel::mask values 407 * specifying the capabilities of the resulting GravityCircle object. 408 * @exception std::bad_alloc if the memory necessary for creating a 409 * GravityCircle can't be allocated. 410 * @return a GravityCircle object whose member functions computes the 411 * gravitational field at a particular values of \e lon. 412 * 413 * The GravityModel::mask values are 414 * - \e caps |= GravityModel::GRAVITY 415 * - \e caps |= GravityModel::DISTURBANCE 416 * - \e caps |= GravityModel::DISTURBING_POTENTIAL 417 * - \e caps |= GravityModel::SPHERICAL_ANOMALY 418 * - \e caps |= GravityModel::GEOID_HEIGHT 419 * . 420 * The default value of \e caps is GravityModel::ALL which turns on all the 421 * capabilities. If an unsupported function is invoked, it will return 422 * NaNs. Note that GravityModel::GEOID_HEIGHT will only be honored if \e h 423 * = 0. 424 * 425 * If the field at several points on a circle of latitude need to be 426 * calculated then creating a GravityCircle object and using its member 427 * functions will be substantially faster, especially for high-degree 428 * models. See \ref gravityparallel for an example of using GravityCircle 429 * (together with OpenMP) to speed up the computation of geoid heights. 430 **********************************************************************/ 431 GravityCircle Circle(real lat, real h, unsigned caps = ALL) const; 432 ///@} 433 434 /** \name Inspector functions 435 **********************************************************************/ 436 ///@{ 437 438 /** 439 * @return the NormalGravity object for the reference ellipsoid. 440 **********************************************************************/ ReferenceEllipsoid() const441 const NormalGravity& ReferenceEllipsoid() const { return _earth; } 442 443 /** 444 * @return the description of the gravity model, if available, in the data 445 * file; if absent, return "NONE". 446 **********************************************************************/ Description() const447 const std::string& Description() const { return _description; } 448 449 /** 450 * @return date of the model; if absent, return "UNKNOWN". 451 **********************************************************************/ DateTime() const452 const std::string& DateTime() const { return _date; } 453 454 /** 455 * @return full file name used to load the gravity model. 456 **********************************************************************/ GravityFile() const457 const std::string& GravityFile() const { return _filename; } 458 459 /** 460 * @return "name" used to load the gravity model (from the first argument 461 * of the constructor, but this may be overridden by the model file). 462 **********************************************************************/ GravityModelName() const463 const std::string& GravityModelName() const { return _name; } 464 465 /** 466 * @return directory used to load the gravity model. 467 **********************************************************************/ GravityModelDirectory() const468 const std::string& GravityModelDirectory() const { return _dir; } 469 470 /** 471 * @return \e a the equatorial radius of the ellipsoid (meters). 472 **********************************************************************/ EquatorialRadius() const473 Math::real EquatorialRadius() const { return _earth.EquatorialRadius(); } 474 475 /** 476 * @return \e GM the mass constant of the model (m<sup>3</sup> 477 * s<sup>−2</sup>); this is the product of \e G the gravitational 478 * constant and \e M the mass of the earth (usually including the mass of 479 * the earth's atmosphere). 480 **********************************************************************/ MassConstant() const481 Math::real MassConstant() const { return _GMmodel; } 482 483 /** 484 * @return \e GM the mass constant of the ReferenceEllipsoid() 485 * (m<sup>3</sup> s<sup>−2</sup>). 486 **********************************************************************/ ReferenceMassConstant() const487 Math::real ReferenceMassConstant() const 488 { return _earth.MassConstant(); } 489 490 /** 491 * @return ω the angular velocity of the model and the 492 * ReferenceEllipsoid() (rad s<sup>−1</sup>). 493 **********************************************************************/ AngularVelocity() const494 Math::real AngularVelocity() const 495 { return _earth.AngularVelocity(); } 496 497 /** 498 * @return \e f the flattening of the ellipsoid. 499 **********************************************************************/ Flattening() const500 Math::real Flattening() const { return _earth.Flattening(); } 501 502 /** 503 * @return \e Nmax the maximum degree of the components of the model. 504 **********************************************************************/ Degree() const505 int Degree() const { return _nmx; } 506 507 /** 508 * @return \e Mmax the maximum order of the components of the model. 509 **********************************************************************/ Order() const510 int Order() const { return _mmx; } 511 512 /** 513 * \deprecated An old name for EquatorialRadius(). 514 **********************************************************************/ 515 GEOGRAPHICLIB_DEPRECATED("Use EquatorialRadius()") MajorRadius() const516 Math::real MajorRadius() const { return EquatorialRadius(); } 517 ///@} 518 519 /** 520 * @return the default path for gravity model data files. 521 * 522 * This is the value of the environment variable 523 * GEOGRAPHICLIB_GRAVITY_PATH, if set; otherwise, it is 524 * $GEOGRAPHICLIB_DATA/gravity if the environment variable 525 * GEOGRAPHICLIB_DATA is set; otherwise, it is a compile-time default 526 * (/usr/local/share/GeographicLib/gravity on non-Windows systems and 527 * C:/ProgramData/GeographicLib/gravity on Windows systems). 528 **********************************************************************/ 529 static std::string DefaultGravityPath(); 530 531 /** 532 * @return the default name for the gravity model. 533 * 534 * This is the value of the environment variable 535 * GEOGRAPHICLIB_GRAVITY_NAME, if set; otherwise, it is "egm96". The 536 * GravityModel class does not use this function; it is just provided as a 537 * convenience for a calling program when constructing a GravityModel 538 * object. 539 **********************************************************************/ 540 static std::string DefaultGravityName(); 541 }; 542 543 } // namespace GeographicLib 544 545 #if defined(_MSC_VER) 546 # pragma warning (pop) 547 #endif 548 549 #endif // GEOGRAPHICLIB_GRAVITYMODEL_HPP 550