GeographicLib  1.50
GravityModel.cpp
Go to the documentation of this file.
1 /**
2  * \file GravityModel.cpp
3  * \brief Implementation for GeographicLib::GravityModel class
4  *
5  * Copyright (c) Charles Karney (2011-2019) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * https://geographiclib.sourceforge.io/
8  **********************************************************************/
9 
11 #include <fstream>
12 #include <limits>
16 
17 #if !defined(GEOGRAPHICLIB_DATA)
18 # if defined(_WIN32)
19 # define GEOGRAPHICLIB_DATA "C:/ProgramData/GeographicLib"
20 # else
21 # define GEOGRAPHICLIB_DATA "/usr/local/share/GeographicLib"
22 # endif
23 #endif
24 
25 #if !defined(GEOGRAPHICLIB_GRAVITY_DEFAULT_NAME)
26 # define GEOGRAPHICLIB_GRAVITY_DEFAULT_NAME "egm96"
27 #endif
28 
29 #if defined(_MSC_VER)
30 // Squelch warnings about unsafe use of getenv
31 # pragma warning (disable: 4996)
32 #endif
33 
34 namespace GeographicLib {
35 
36  using namespace std;
37 
38  GravityModel::GravityModel(const std::string& name, const std::string& path,
39  int Nmax, int Mmax)
40  : _name(name)
41  , _dir(path)
42  , _description("NONE")
43  , _date("UNKNOWN")
44  , _amodel(Math::NaN())
45  , _GMmodel(Math::NaN())
46  , _zeta0(0)
47  , _corrmult(1)
48  , _nmx(-1)
49  , _mmx(-1)
50  , _norm(SphericalHarmonic::FULL)
51  {
52  if (_dir.empty())
53  _dir = DefaultGravityPath();
54  bool truncate = Nmax >= 0 || Mmax >= 0;
55  if (truncate) {
56  if (Nmax >= 0 && Mmax < 0) Mmax = Nmax;
57  if (Nmax < 0) Nmax = numeric_limits<int>::max();
58  if (Mmax < 0) Mmax = numeric_limits<int>::max();
59  }
60  ReadMetadata(_name);
61  {
62  string coeff = _filename + ".cof";
63  ifstream coeffstr(coeff.c_str(), ios::binary);
64  if (!coeffstr.good())
65  throw GeographicErr("Error opening " + coeff);
66  char id[idlength_ + 1];
67  coeffstr.read(id, idlength_);
68  if (!coeffstr.good())
69  throw GeographicErr("No header in " + coeff);
70  id[idlength_] = '\0';
71  if (_id != string(id))
72  throw GeographicErr("ID mismatch: " + _id + " vs " + id);
73  int N, M;
74  if (truncate) { N = Nmax; M = Mmax; }
75  SphericalEngine::coeff::readcoeffs(coeffstr, N, M, _Cx, _Sx, truncate);
76  if (!(N >= 0 && M >= 0))
77  throw GeographicErr("Degree and order must be at least 0");
78  if (_Cx[0] != 0)
79  throw GeographicErr("The degree 0 term should be zero");
80  _Cx[0] = 1; // Include the 1/r term in the sum
81  _gravitational = SphericalHarmonic(_Cx, _Sx, N, N, M, _amodel, _norm);
82  if (truncate) { N = Nmax; M = Mmax; }
83  SphericalEngine::coeff::readcoeffs(coeffstr, N, M, _CC, _CS, truncate);
84  if (N < 0) {
85  N = M = 0;
86  _CC.resize(1, real(0));
87  }
88  _CC[0] += _zeta0 / _corrmult;
89  _correction = SphericalHarmonic(_CC, _CS, N, N, M, real(1), _norm);
90  int pos = int(coeffstr.tellg());
91  coeffstr.seekg(0, ios::end);
92  if (pos != coeffstr.tellg())
93  throw GeographicErr("Extra data in " + coeff);
94  }
95  int nmx = _gravitational.Coefficients().nmx();
96  _nmx = max(nmx, _correction.Coefficients().nmx());
97  _mmx = max(_gravitational.Coefficients().mmx(),
98  _correction.Coefficients().mmx());
99  // Adjust the normalization of the normal potential to match the model.
100  real mult = _earth._GM / _GMmodel;
101  real amult = Math::sq(_earth._a / _amodel);
102  // The 0th term in _zonal should be is 1 + _dzonal0. Instead set it to 1
103  // to give exact cancellation with the (0,0) term in the model and account
104  // for _dzonal0 separately.
105  _zonal.clear(); _zonal.push_back(1);
106  _dzonal0 = (_earth.MassConstant() - _GMmodel) / _GMmodel;
107  for (int n = 2; n <= nmx; n += 2) {
108  // Only include as many normal zonal terms as matter. Figuring the limit
109  // in this way works because the coefficients of the normal potential
110  // (which is smooth) decay much more rapidly that the corresponding
111  // coefficient of the model potential (which is bumpy). Typically this
112  // goes out to n = 18.
113  mult *= amult;
114  real
115  r = _Cx[n], // the model term
116  s = - mult * _earth.Jn(n) / sqrt(real(2 * n + 1)), // the normal term
117  t = r - s; // the difference
118  if (t == r) // the normal term is negligible
119  break;
120  _zonal.push_back(0); // index = n - 1; the odd terms are 0
121  _zonal.push_back(s);
122  }
123  int nmx1 = int(_zonal.size()) - 1;
124  _disturbing = SphericalHarmonic1(_Cx, _Sx,
125  _gravitational.Coefficients().N(),
126  nmx, _gravitational.Coefficients().mmx(),
127  _zonal,
128  _zonal, // This is not accessed!
129  nmx1, nmx1, 0,
130  _amodel,
132  }
133 
134  void GravityModel::ReadMetadata(const std::string& name) {
135  const char* spaces = " \t\n\v\f\r";
136  _filename = _dir + "/" + name + ".egm";
137  ifstream metastr(_filename.c_str());
138  if (!metastr.good())
139  throw GeographicErr("Cannot open " + _filename);
140  string line;
141  getline(metastr, line);
142  if (!(line.size() >= 6 && line.substr(0,5) == "EGMF-"))
143  throw GeographicErr(_filename + " does not contain EGMF-n signature");
144  string::size_type n = line.find_first_of(spaces, 5);
145  if (n != string::npos)
146  n -= 5;
147  string version(line, 5, n);
148  if (version != "1")
149  throw GeographicErr("Unknown version in " + _filename + ": " + version);
150  string key, val;
151  real a = Math::NaN(), GM = a, omega = a, f = a, J2 = a;
152  while (getline(metastr, line)) {
153  if (!Utility::ParseLine(line, key, val))
154  continue;
155  // Process key words
156  if (key == "Name")
157  _name = val;
158  else if (key == "Description")
159  _description = val;
160  else if (key == "ReleaseDate")
161  _date = val;
162  else if (key == "ModelRadius")
163  _amodel = Utility::val<real>(val);
164  else if (key == "ModelMass")
165  _GMmodel = Utility::val<real>(val);
166  else if (key == "AngularVelocity")
167  omega = Utility::val<real>(val);
168  else if (key == "ReferenceRadius")
169  a = Utility::val<real>(val);
170  else if (key == "ReferenceMass")
171  GM = Utility::val<real>(val);
172  else if (key == "Flattening")
173  f = Utility::fract<real>(val);
174  else if (key == "DynamicalFormFactor")
175  J2 = Utility::fract<real>(val);
176  else if (key == "HeightOffset")
177  _zeta0 = Utility::fract<real>(val);
178  else if (key == "CorrectionMultiplier")
179  _corrmult = Utility::fract<real>(val);
180  else if (key == "Normalization") {
181  if (val == "FULL" || val == "Full" || val == "full")
182  _norm = SphericalHarmonic::FULL;
183  else if (val == "SCHMIDT" || val == "Schmidt" || val == "schmidt")
185  else
186  throw GeographicErr("Unknown normalization " + val);
187  } else if (key == "ByteOrder") {
188  if (val == "Big" || val == "big")
189  throw GeographicErr("Only little-endian ordering is supported");
190  else if (!(val == "Little" || val == "little"))
191  throw GeographicErr("Unknown byte ordering " + val);
192  } else if (key == "ID")
193  _id = val;
194  // else unrecognized keywords are skipped
195  }
196  // Check values
197  if (!(Math::isfinite(_amodel) && _amodel > 0))
198  throw GeographicErr("Model radius must be positive");
199  if (!(Math::isfinite(_GMmodel) && _GMmodel > 0))
200  throw GeographicErr("Model mass constant must be positive");
201  if (!(Math::isfinite(_corrmult) && _corrmult > 0))
202  throw GeographicErr("Correction multiplier must be positive");
203  if (!(Math::isfinite(_zeta0)))
204  throw GeographicErr("Height offset must be finite");
205  if (int(_id.size()) != idlength_)
206  throw GeographicErr("Invalid ID");
207  if (Math::isfinite(f) && Math::isfinite(J2))
208  throw GeographicErr("Cannot specify both f and J2");
209  _earth = NormalGravity(a, GM, omega,
210  Math::isfinite(f) ? f : J2, Math::isfinite(f));
211  }
212 
213  Math::real GravityModel::InternalT(real X, real Y, real Z,
214  real& deltaX, real& deltaY, real& deltaZ,
215  bool gradp, bool correct) const {
216  // If correct, then produce the correct T = W - U. Otherwise, neglect the
217  // n = 0 term (which is proportial to the difference in the model and
218  // reference values of GM).
219  if (_dzonal0 == 0)
220  // No need to do the correction
221  correct = false;
222  real T, invR = correct ? 1 / Math::hypot(Math::hypot(X, Y), Z) : 1;
223  if (gradp) {
224  // initial values to suppress warnings
225  deltaX = deltaY = deltaZ = 0;
226  T = _disturbing(-1, X, Y, Z, deltaX, deltaY, deltaZ);
227  real f = _GMmodel / _amodel;
228  deltaX *= f;
229  deltaY *= f;
230  deltaZ *= f;
231  if (correct) {
232  invR = _GMmodel * _dzonal0 * invR * invR * invR;
233  deltaX += X * invR;
234  deltaY += Y * invR;
235  deltaZ += Z * invR;
236  }
237  } else
238  T = _disturbing(-1, X, Y, Z);
239  T = (T / _amodel - (correct ? _dzonal0 : 0) * invR) * _GMmodel;
240  return T;
241  }
242 
243  Math::real GravityModel::V(real X, real Y, real Z,
244  real& GX, real& GY, real& GZ) const {
245  real
246  Vres = _gravitational(X, Y, Z, GX, GY, GZ),
247  f = _GMmodel / _amodel;
248  Vres *= f;
249  GX *= f;
250  GY *= f;
251  GZ *= f;
252  return Vres;
253  }
254 
255  Math::real GravityModel::W(real X, real Y, real Z,
256  real& gX, real& gY, real& gZ) const {
257  real fX, fY,
258  Wres = V(X, Y, Z, gX, gY, gZ) + _earth.Phi(X, Y, fX, fY);
259  gX += fX;
260  gY += fY;
261  return Wres;
262  }
263 
264  void GravityModel::SphericalAnomaly(real lat, real lon, real h,
265  real& Dg01, real& xi, real& eta) const {
266  real X, Y, Z, M[Geocentric::dim2_];
267  _earth.Earth().IntForward(lat, lon, h, X, Y, Z, M);
268  real
269  deltax, deltay, deltaz,
270  T = InternalT(X, Y, Z, deltax, deltay, deltaz, true, false),
271  clam = M[3], slam = -M[0],
272  P = Math::hypot(X, Y),
273  R = Math::hypot(P, Z),
274  // psi is geocentric latitude
275  cpsi = R != 0 ? P / R : M[7],
276  spsi = R != 0 ? Z / R : M[8];
277  // Rotate cartesian into spherical coordinates
278  real MC[Geocentric::dim2_];
279  Geocentric::Rotation(spsi, cpsi, slam, clam, MC);
280  Geocentric::Unrotate(MC, deltax, deltay, deltaz, deltax, deltay, deltaz);
281  // H+M, Eq 2-151c
282  Dg01 = - deltaz - 2 * T / R;
283  real gammaX, gammaY, gammaZ;
284  _earth.U(X, Y, Z, gammaX, gammaY, gammaZ);
285  real gamma = Math::hypot( Math::hypot(gammaX, gammaY), gammaZ);
286  xi = -(deltay/gamma) / Math::degree();
287  eta = -(deltax/gamma) / Math::degree();
288  }
289 
290  Math::real GravityModel::GeoidHeight(real lat, real lon) const
291  {
292  real X, Y, Z;
293  _earth.Earth().IntForward(lat, lon, 0, X, Y, Z, NULL);
294  real
295  gamma0 = _earth.SurfaceGravity(lat),
296  dummy,
297  T = InternalT(X, Y, Z, dummy, dummy, dummy, false, false),
298  invR = 1 / Math::hypot(Math::hypot(X, Y), Z),
299  correction = _corrmult * _correction(invR * X, invR * Y, invR * Z);
300  // _zeta0 has been included in _correction
301  return T/gamma0 + correction;
302  }
303 
304  Math::real GravityModel::Gravity(real lat, real lon, real h,
305  real& gx, real& gy, real& gz) const {
306  real X, Y, Z, M[Geocentric::dim2_];
307  _earth.Earth().IntForward(lat, lon, h, X, Y, Z, M);
308  real Wres = W(X, Y, Z, gx, gy, gz);
309  Geocentric::Unrotate(M, gx, gy, gz, gx, gy, gz);
310  return Wres;
311  }
312  Math::real GravityModel::Disturbance(real lat, real lon, real h,
313  real& deltax, real& deltay,
314  real& deltaz) const {
315  real X, Y, Z, M[Geocentric::dim2_];
316  _earth.Earth().IntForward(lat, lon, h, X, Y, Z, M);
317  real Tres = InternalT(X, Y, Z, deltax, deltay, deltaz, true, true);
318  Geocentric::Unrotate(M, deltax, deltay, deltaz, deltax, deltay, deltaz);
319  return Tres;
320  }
321 
322  GravityCircle GravityModel::Circle(real lat, real h, unsigned caps) const {
323  if (h != 0)
324  // Disallow invoking GeoidHeight unless h is zero.
325  caps &= ~(CAP_GAMMA0 | CAP_C);
326  real X, Y, Z, M[Geocentric::dim2_];
327  _earth.Earth().IntForward(lat, 0, h, X, Y, Z, M);
328  // Y = 0, cphi = M[7], sphi = M[8];
329  real
330  invR = 1 / Math::hypot(X, Z),
331  gamma0 = (caps & CAP_GAMMA0 ?_earth.SurfaceGravity(lat)
332  : Math::NaN()),
333  fx, fy, fz, gamma;
334  if (caps & CAP_GAMMA) {
335  _earth.U(X, Y, Z, fx, fy, fz); // fy = 0
336  gamma = Math::hypot(fx, fz);
337  } else
338  gamma = Math::NaN();
339  _earth.Phi(X, Y, fx, fy);
340  return GravityCircle(GravityCircle::mask(caps),
341  _earth._a, _earth._f, lat, h, Z, X, M[7], M[8],
342  _amodel, _GMmodel, _dzonal0, _corrmult,
343  gamma0, gamma, fx,
344  caps & CAP_G ?
345  _gravitational.Circle(X, Z, true) :
346  CircularEngine(),
347  // N.B. If CAP_DELTA is set then CAP_T should be too.
348  caps & CAP_T ?
349  _disturbing.Circle(-1, X, Z, (caps&CAP_DELTA) != 0) :
350  CircularEngine(),
351  caps & CAP_C ?
352  _correction.Circle(invR * X, invR * Z, false) :
353  CircularEngine());
354  }
355 
357  string path;
358  char* gravitypath = getenv("GEOGRAPHICLIB_GRAVITY_PATH");
359  if (gravitypath)
360  path = string(gravitypath);
361  if (!path.empty())
362  return path;
363  char* datapath = getenv("GEOGRAPHICLIB_DATA");
364  if (datapath)
365  path = string(datapath);
366  return (!path.empty() ? path : string(GEOGRAPHICLIB_DATA)) + "/gravity";
367  }
368 
370  string name;
371  char* gravityname = getenv("GEOGRAPHICLIB_GRAVITY_NAME");
372  if (gravityname)
373  name = string(gravityname);
374  return !name.empty() ? name : string(GEOGRAPHICLIB_GRAVITY_DEFAULT_NAME);
375  }
376 
377 } // namespace GeographicLib
static void readcoeffs(std::istream &stream, int &N, int &M, std::vector< real > &C, std::vector< real > &S, bool truncate=false)
Math::real T(real X, real Y, real Z, real &deltaX, real &deltaY, real &deltaZ) const
static bool isfinite(T x)
Definition: Math.cpp:372
CircularEngine Circle(real tau, real p, real z, bool gradp) const
Header for GeographicLib::Utility class.
The normal gravity of the earth.
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:98
Header for GeographicLib::GravityModel class.
const SphericalEngine::coeff & Coefficients() const
#define GEOGRAPHICLIB_DATA
Math::real MassConstant() const
Math::real Gravity(real lat, real lon, real h, real &gx, real &gy, real &gz) const
static T sq(T x)
Definition: Math.hpp:201
CircularEngine Circle(real p, real z, bool gradp) const
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
static T NaN()
Definition: Math.cpp:389
Math::real Disturbance(real lat, real lon, real h, real &deltax, real &deltay, real &deltaz) const
static T degree()
Definition: Math.hpp:185
Spherical harmonic sums for a circle.
Math::real SurfaceGravity(real lat) const
static std::string DefaultGravityName()
GravityCircle Circle(real lat, real h, unsigned caps=ALL) const
Math::real V(real X, real Y, real Z, real &GX, real &GY, real &GZ) const
Exception handling for GeographicLib.
Definition: Constants.hpp:390
static std::string DefaultGravityPath()
Spherical harmonic series with a correction to the coefficients.
static T hypot(T x, T y)
Definition: Math.cpp:58
Spherical harmonic series.
Math::real W(real X, real Y, real Z, real &gX, real &gY, real &gZ) const
Header for GeographicLib::GravityCircle class.
Math::real Phi(real X, real Y, real &fX, real &fY) const
static bool ParseLine(const std::string &line, std::string &key, std::string &val)
Definition: Utility.cpp:22
Header for GeographicLib::SphericalEngine class.
#define GEOGRAPHICLIB_GRAVITY_DEFAULT_NAME
void SphericalAnomaly(real lat, real lon, real h, real &Dg01, real &xi, real &eta) const
const Geocentric & Earth() const
Math::real U(real X, real Y, real Z, real &gammaX, real &gammaY, real &gammaZ) const
Gravity on a circle of latitude.
Math::real GeoidHeight(real lat, real lon) const