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) 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 (!(N >= 0 && M >= 0))
74 _CC.resize(1, real(0));
76 _CC[0] += _zeta0 / _corrmult;
78 int pos = int(coeffstr.tellg());
79 coeffstr.seekg(0, ios::end);
80 if (pos != coeffstr.tellg())
85 real mult = _earth._GM / _GMmodel;
86 real amult =
Math::sq(_earth._a / _amodel);
90 _zonal.clear(); _zonal.push_back(1);
91 _dzonal0 = (_earth.
MassConstant() - _GMmodel) / _GMmodel;
92 for (
int n = 2; n <= nmx; n += 2) {
101 s = - mult * _earth.Jn(n) / sqrt(real(2 * n + 1)),
108 int nmx1 = int(_zonal.size()) - 1;
119 void GravityModel::ReadMetadata(
const std::string& name) {
120 const char* spaces =
" \t\n\v\f\r";
121 _filename = _dir +
"/" + name +
".egm";
122 ifstream metastr(_filename.c_str());
126 getline(metastr, line);
127 if (!(line.size() >= 6 && line.substr(0,5) ==
"EGMF-"))
128 throw GeographicErr(_filename +
" does not contain EGMF-n signature");
129 string::size_type n = line.find_first_of(spaces, 5);
130 if (n != string::npos)
132 string version(line, 5, n);
134 throw GeographicErr(
"Unknown version in " + _filename +
": " + version);
136 real a =
Math::NaN(), GM = a, omega = a, f = a, J2 = a;
137 while (getline(metastr, line)) {
143 else if (key ==
"Description")
145 else if (key ==
"ReleaseDate")
147 else if (key ==
"ModelRadius")
148 _amodel = Utility::val<real>(val);
149 else if (key ==
"ModelMass")
150 _GMmodel = Utility::val<real>(val);
151 else if (key ==
"AngularVelocity")
152 omega = Utility::val<real>(val);
153 else if (key ==
"ReferenceRadius")
154 a = Utility::val<real>(val);
155 else if (key ==
"ReferenceMass")
156 GM = Utility::val<real>(val);
157 else if (key ==
"Flattening")
158 f = Utility::fract<real>(val);
159 else if (key ==
"DynamicalFormFactor")
160 J2 = Utility::fract<real>(val);
161 else if (key ==
"HeightOffset")
162 _zeta0 = Utility::fract<real>(val);
163 else if (key ==
"CorrectionMultiplier")
164 _corrmult = Utility::fract<real>(val);
165 else if (key ==
"Normalization") {
166 if (val ==
"FULL" || val ==
"Full" || val ==
"full")
168 else if (val ==
"SCHMIDT" || val ==
"Schmidt" || val ==
"schmidt")
172 }
else if (key ==
"ByteOrder") {
173 if (val ==
"Big" || val ==
"big")
174 throw GeographicErr(
"Only little-endian ordering is supported");
175 else if (!(val ==
"Little" || val ==
"little"))
177 }
else if (key ==
"ID")
187 throw GeographicErr(
"Correction multiplier must be positive");
190 if (
int(_id.size()) != idlength_)
198 Math::real GravityModel::InternalT(real X, real Y, real Z,
199 real& deltaX, real& deltaY, real& deltaZ,
200 bool gradp,
bool correct)
const {
210 deltaX = deltaY = deltaZ = 0;
211 T = _disturbing(-1, X, Y, Z, deltaX, deltaY, deltaZ);
212 real f = _GMmodel / _amodel;
217 invR = _GMmodel * _dzonal0 * invR * invR * invR;
223 T = _disturbing(-1, X, Y, Z);
224 T = (T / _amodel - (correct ? _dzonal0 : 0) * invR) * _GMmodel;
229 real& GX, real& GY, real& GZ)
const {
231 Vres = _gravitational(X, Y, Z, GX, GY, GZ),
232 f = _GMmodel / _amodel;
241 real& gX, real& gY, real& gZ)
const {
243 Wres =
V(X, Y, Z, gX, gY, gZ) + _earth.
Phi(X, Y, fX, fY);
250 real& Dg01, real& xi, real& eta)
252 real X, Y, Z, M[Geocentric::dim2_];
253 _earth.
Earth().IntForward(lat, lon, h, X, Y, Z, M);
255 deltax, deltay, deltaz,
256 T = InternalT(X, Y, Z, deltax, deltay, deltaz,
true,
false),
257 clam = M[3], slam = -M[0],
261 cpsi = R ? P / R : M[7],
262 spsi = R ? Z / R : M[8];
264 real MC[Geocentric::dim2_];
265 Geocentric::Rotation(spsi, cpsi, slam, clam, MC);
266 Geocentric::Unrotate(MC, deltax, deltay, deltaz, deltax, deltay, deltaz);
268 Dg01 = - deltaz - 2 * T / R;
269 real gammaX, gammaY, gammaZ;
270 _earth.
U(X, Y, Z, gammaX, gammaY, gammaZ);
279 _earth.
Earth().IntForward(lat, lon, 0, X, Y, Z, NULL);
283 T = InternalT(X, Y, Z, dummy, dummy, dummy,
false,
false),
285 correction = _corrmult * _correction(invR * X, invR * Y, invR * Z);
287 return T/gamma0 + correction;
291 real& gx, real& gy, real& gz)
const {
292 real X, Y, Z, M[Geocentric::dim2_];
293 _earth.
Earth().IntForward(lat, lon, h, X, Y, Z, M);
294 real Wres =
W(X, Y, Z, gx, gy, gz);
295 Geocentric::Unrotate(M, gx, gy, gz, gx, gy, gz);
299 real& deltax, real& deltay, real& deltaz)
301 real X, Y, Z, M[Geocentric::dim2_];
302 _earth.
Earth().IntForward(lat, lon, h, X, Y, Z, M);
303 real Tres = InternalT(X, Y, Z, deltax, deltay, deltaz,
true,
true);
304 Geocentric::Unrotate(M, deltax, deltay, deltaz, deltax, deltay, deltaz);
311 caps &= ~(CAP_GAMMA0 | CAP_C);
312 real X, Y, Z, M[Geocentric::dim2_];
313 _earth.
Earth().IntForward(lat, 0, h, X, Y, Z, M);
320 if (caps & CAP_GAMMA) {
321 _earth.
U(X, Y, Z, fx, fy, fz);
325 _earth.
Phi(X, Y, fx, fy);
327 _earth._a, _earth._f, lat, h, Z, X, M[7], M[8],
328 _amodel, _GMmodel, _dzonal0, _corrmult,
331 _gravitational.
Circle(X, Z,
true) :
335 _disturbing.
Circle(-1, X, Z, (caps & CAP_DELTA) != 0) :
338 _correction.
Circle(invR * X, invR * Z,
false) :
344 char* gravitypath = getenv(
"GEOGRAPHICLIB_GRAVITY_PATH");
346 path = string(gravitypath);
349 char* datapath = getenv(
"GEOGRAPHICLIB_DATA");
351 path = string(datapath);
357 char* gravityname = getenv(
"GEOGRAPHICLIB_GRAVITY_NAME");
359 name = string(gravityname);
Math::real T(real X, real Y, real Z, real &deltaX, real &deltaY, real &deltaZ) const
CircularEngine Circle(real tau, real p, real z, bool gradp) const
Header for GeographicLib::Utility class.
The normal gravity of the earth.
static bool isfinite(T x)
Mathematical functions needed by GeographicLib.
Header for GeographicLib::GravityModel class.
const SphericalEngine::coeff & Coefficients() const
#define GEOGRAPHICLIB_DATA
static void readcoeffs(std::istream &stream, int &N, int &M, std::vector< real > &C, std::vector< real > &S)
Math::real MassConstant() const
Math::real Gravity(real lat, real lon, real h, real &gx, real &gy, real &gz) const
CircularEngine Circle(real p, real z, bool gradp) const
friend class GravityCircle
Namespace for GeographicLib.
Math::real Disturbance(real lat, real lon, real h, real &deltax, real &deltay, real &deltaz) const
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.
static std::string DefaultGravityPath()
Spherical harmonic series with a correction to the coefficients.
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)
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