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