16 #if !defined(GEOGRAPHICLIB_DATA)
18 # define GEOGRAPHICLIB_DATA "C:/ProgramData/GeographicLib"
20 # define GEOGRAPHICLIB_DATA "/usr/local/share/GeographicLib"
24 #if !defined(GEOGRAPHICLIB_GRAVITY_DEFAULT_NAME)
25 # define GEOGRAPHICLIB_GRAVITY_DEFAULT_NAME "egm96"
30 # pragma warning (disable: 4996)
33 namespace GeographicLib {
37 GravityModel::GravityModel(
const std::string& name,
const std::string& path)
40 , _description(
"NONE")
42 , _amodel(
Math::NaN())
43 , _GMmodel(
Math::NaN())
52 string coeff = _filename +
".cof";
53 ifstream coeffstr(coeff.c_str(), ios::binary);
56 char id[idlength_ + 1];
57 coeffstr.read(
id, idlength_);
61 if (_id !=
string(
id))
65 if (!(M < 0 || _Cx[0] == 0))
72 _CC.resize(1, real(0));
74 _CC[0] += _zeta0 / _corrmult;
76 int pos = int(coeffstr.tellg());
77 coeffstr.seekg(0, ios::end);
78 if (pos != coeffstr.tellg())
83 real mult = _earth._GM / _GMmodel;
84 real amult =
Math::sq(_earth._a / _amodel);
88 _zonal.clear(); _zonal.push_back(1);
89 _dzonal0 = (_earth.
MassConstant() - _GMmodel) / _GMmodel;
90 for (
int n = 2; n <= nmx; n += 2) {
99 s = - mult * _earth.Jn(n) / sqrt(real(2 * n + 1)),
106 int nmx1 = int(_zonal.size()) - 1;
117 void GravityModel::ReadMetadata(
const std::string& name) {
118 const char* spaces =
" \t\n\v\f\r";
119 _filename = _dir +
"/" + name +
".egm";
120 ifstream metastr(_filename.c_str());
124 getline(metastr, line);
125 if (!(line.size() >= 6 && line.substr(0,5) ==
"EGMF-"))
126 throw GeographicErr(_filename +
" does not contain EGMF-n signature");
127 string::size_type n = line.find_first_of(spaces, 5);
128 if (n != string::npos)
130 string version = line.substr(5, n);
132 throw GeographicErr(
"Unknown version in " + _filename +
": " + version);
135 while (getline(metastr, line)) {
141 else if (key ==
"Description")
143 else if (key ==
"ReleaseDate")
145 else if (key ==
"ModelRadius")
146 _amodel = Utility::num<real>(val);
147 else if (key ==
"ModelMass")
148 _GMmodel = Utility::num<real>(val);
149 else if (key ==
"AngularVelocity")
150 omega = Utility::num<real>(val);
151 else if (key ==
"ReferenceRadius")
152 a = Utility::num<real>(val);
153 else if (key ==
"ReferenceMass")
154 GM = Utility::num<real>(val);
155 else if (key ==
"Flattening")
156 f = Utility::fract<real>(val);
157 else if (key ==
"DynamicalFormFactor")
158 J2 = Utility::fract<real>(val);
159 else if (key ==
"HeightOffset")
160 _zeta0 = Utility::fract<real>(val);
161 else if (key ==
"CorrectionMultiplier")
162 _corrmult = Utility::fract<real>(val);
163 else if (key ==
"Normalization") {
164 if (val ==
"FULL" || val ==
"Full" || val ==
"full")
166 else if (val ==
"SCHMIDT" || val ==
"Schmidt" || val ==
"schmidt")
170 }
else if (key ==
"ByteOrder") {
171 if (val ==
"Big" || val ==
"big")
172 throw GeographicErr(
"Only little-endian ordering is supported");
173 else if (!(val ==
"Little" || val ==
"little"))
174 throw GeographicErr(
"Unknown byte ordering " + val);
175 }
else if (key ==
"ID")
181 throw GeographicErr(
"Model radius must be positive");
183 throw GeographicErr(
"Model mass constant must be positive");
185 throw GeographicErr(
"Correction multiplier must be positive");
187 throw GeographicErr(
"Height offset must be finite");
188 if (
int(_id.size()) != idlength_)
189 throw GeographicErr(
"Invalid ID");
190 _earth = NormalGravity(a, GM, omega, f, J2);
195 bool gradp,
bool correct)
const {
205 deltaX = deltaY = deltaZ = 0;
206 T = _disturbing(-1, X, Y, Z, deltaX, deltaY, deltaZ);
207 real f = _GMmodel / _amodel;
212 invR = _GMmodel * _dzonal0 * invR * invR * invR;
218 T = _disturbing(-1, X, Y, Z);
219 T = (T / _amodel - (correct ? _dzonal0 : 0) * invR) * _GMmodel;
224 real& GX, real& GY, real& GZ)
const {
226 Vres = _gravitational(X, Y, Z, GX, GY, GZ),
227 f = _GMmodel / _amodel;
236 real& gX, real& gY, real& gZ)
const {
238 Wres =
V(X, Y, Z, gX, gY, gZ) + _earth.
Phi(X, Y, fX, fY);
245 real& Dg01, real& xi, real& eta)
247 real X, Y, Z, M[Geocentric::dim2_];
248 _earth.
Earth().IntForward(lat, lon, h, X, Y, Z, M);
250 deltax, deltay, deltaz,
251 T = InternalT(X, Y, Z, deltax, deltay, deltaz,
true,
false),
252 clam = M[3], slam = -M[0],
256 cpsi = R ? P / R : M[7],
257 spsi = R ? Z / R : M[8];
259 real MC[Geocentric::dim2_];
260 Geocentric::Rotation(spsi, cpsi, slam, clam, MC);
261 Geocentric::Unrotate(MC, deltax, deltay, deltaz, deltax, deltay, deltaz);
263 Dg01 = - deltaz - 2 * T / R;
264 real gammaX, gammaY, gammaZ;
265 _earth.
U(X, Y, Z, gammaX, gammaY, gammaZ);
274 _earth.
Earth().IntForward(lat, lon, 0, X, Y, Z, NULL);
278 T = InternalT(X, Y, Z, dummy, dummy, dummy,
false,
false),
280 correction = _corrmult * _correction(invR * X, invR * Y, invR * Z);
282 return T/gamma0 + correction;
286 real& gx, real& gy, real& gz)
const {
287 real X, Y, Z, M[Geocentric::dim2_];
288 _earth.
Earth().IntForward(lat, lon, h, X, Y, Z, M);
289 real Wres =
W(X, Y, Z, gx, gy, gz);
290 Geocentric::Unrotate(M, gx, gy, gz, gx, gy, gz);
294 real& deltax, real& deltay, real& deltaz)
296 real X, Y, Z, M[Geocentric::dim2_];
297 _earth.
Earth().IntForward(lat, lon, h, X, Y, Z, M);
298 real Tres = InternalT(X, Y, Z, deltax, deltay, deltaz,
true,
true);
299 Geocentric::Unrotate(M, deltax, deltay, deltaz, deltax, deltay, deltaz);
306 caps &= ~(CAP_GAMMA0 | CAP_C);
307 real X, Y, Z, M[Geocentric::dim2_];
308 _earth.
Earth().IntForward(lat, 0, h, X, Y, Z, M);
315 if (caps & CAP_GAMMA) {
316 _earth.
U(X, Y, Z, fx, fy, fz);
320 _earth.
Phi(X, Y, fx, fy);
322 _earth._a, _earth._f, lat, h, Z, X, M[7], M[8],
323 _amodel, _GMmodel, _dzonal0, _corrmult,
326 _gravitational.
Circle(X, Z,
true) :
330 _disturbing.
Circle(-1, X, Z, (caps & CAP_DELTA) != 0) :
333 _correction.
Circle(invR * X, invR * Z,
false) :
339 char* gravitypath = getenv(
"GEOGRAPHICLIB_GRAVITY_PATH");
341 path = string(gravitypath);
344 char* datapath = getenv(
"GEOGRAPHICLIB_DATA");
346 path = string(datapath);
352 char* gravityname = getenv(
"GEOGRAPHICLIB_GRAVITY_NAME");
354 name = string(gravityname);
Math::real SurfaceGravity(real lat) const
GeographicLib::Math::real real
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)
CircularEngine Circle(real p, real z, bool gradp) const
Math::real T(real X, real Y, real Z, real &deltaX, real &deltaY, real &deltaZ) const
Mathematical functions needed by GeographicLib.
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
GravityCircle Circle(real lat, real h, unsigned caps=ALL) const
friend class GravityCircle
const SphericalEngine::coeff & Coefficients() const
Spherical harmonic sums for a circle.
static std::string DefaultGravityName()
Exception handling for GeographicLib.
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)
Header for GeographicLib::SphericalEngine class.
#define GEOGRAPHICLIB_GRAVITY_DEFAULT_NAME
Gravity on a circle of latitude.
Math::real MassConstant() const