16 #if !defined(GEOGRAPHICLIB_DATA)
18 # define GEOGRAPHICLIB_DATA \
19 "C:/Documents and Settings/All Users/Application Data/GeographicLib"
21 # define GEOGRAPHICLIB_DATA "/usr/local/share/GeographicLib"
25 #if !defined(GRAVITY_DEFAULT_NAME)
26 # define GRAVITY_DEFAULT_NAME "egm96"
31 # pragma warning (disable: 4996)
34 namespace GeographicLib {
38 GravityModel::GravityModel(
const std::string& name,
const std::string& path)
41 , _description(
"NONE")
43 , _amodel(
Math::NaN<real>())
44 , _GMmodel(
Math::NaN<real>())
53 string coeff = _filename +
".cof";
54 ifstream coeffstr(coeff.c_str(), ios::binary);
57 char id[idlength_ + 1];
58 coeffstr.read(
id, idlength_);
62 if (_id !=
string(
id))
66 if (!(M < 0 || _Cx[0] == 0))
73 _CC.resize(1, real(0));
75 _CC[0] += _zeta0 / _corrmult;
77 int pos = int(coeffstr.tellg());
78 coeffstr.seekg(0, ios::end);
79 if (pos != coeffstr.tellg())
84 real mult = _earth._GM / _GMmodel;
85 real amult =
Math::sq(_earth._a / _amodel);
89 _zonal.clear(); _zonal.push_back(1);
90 _dzonal0 = (_earth.
MassConstant() - _GMmodel) / _GMmodel;
91 for (
int n = 2; n <= nmx; n += 2) {
100 s = - mult * _earth.Jn(n) / sqrt(real(2 * n + 1)),
107 int nmx1 = int(_zonal.size()) - 1;
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());
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)
131 string version = line.substr(5, n);
133 throw GeographicErr(
"Unknown version in " + _filename +
": " + version);
135 real a = Math::NaN<real>(), GM = a, omega = a, f = a, J2 = a;
136 while (getline(metastr, line)) {
142 else if (key ==
"Description")
144 else if (key ==
"ReleaseDate")
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")
167 else if (val ==
"SCHMIDT" || val ==
"Schmidt" || val ==
"schmidt")
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")
182 throw GeographicErr(
"Model radius must be positive");
184 throw GeographicErr(
"Model mass constant must be positive");
186 throw GeographicErr(
"Correction multiplier must be positive");
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);
196 bool gradp,
bool correct)
const throw() {
206 ? _disturbing(-1, X, Y, Z, deltaX, deltaY, deltaZ)
207 : _disturbing(-1, X, Y, Z));
208 T = (T / _amodel - (correct ? _dzonal0 : 0) * invR) * _GMmodel;
210 real f = _GMmodel / _amodel;
215 invR = _GMmodel * _dzonal0 * invR * invR * invR;
225 real& GX, real& GY, real& GZ)
const throw() {
227 Vres = _gravitational(X, Y, Z, GX, GY, GZ),
228 f = _GMmodel / _amodel;
237 real& gX, real& gY, real& gZ)
const throw() {
239 Wres = V(X, Y, Z, gX, gY, gZ) + _earth.Phi(X, Y, fX, fY);
246 real& Dg01, real& xi, real& eta)
248 real X, Y, Z, M[Geocentric::dim2_];
249 _earth.Earth().IntForward(lat, lon, h, X, Y, Z, M);
251 deltax, deltay, deltaz,
252 T = InternalT(X, Y, Z, deltax, deltay, deltaz,
true,
false),
253 clam = M[3], slam = -M[0],
257 cpsi = R ? P / R : M[7],
258 spsi = R ? Z / R : M[8];
260 real MC[Geocentric::dim2_];
261 Geocentric::Rotation(spsi, cpsi, slam, clam, MC);
262 Geocentric::Unrotate(MC, deltax, deltay, deltaz, deltax, deltay, deltaz);
264 Dg01 = - deltaz - 2 * T / R;
265 real gammaX, gammaY, gammaZ;
266 _earth.U(X, Y, Z, gammaX, gammaY, gammaZ);
268 xi = -(deltay/gamma) / Math::degree<real>();
269 eta = -(deltax/gamma) / Math::degree<real>();
275 _earth.Earth().IntForward(lat, lon, 0, X, Y, Z, NULL);
277 gamma0 = _earth.SurfaceGravity(lat),
279 T = InternalT(X, Y, Z, dummy, dummy, dummy,
false,
false),
281 correction = _corrmult * _correction(invR * X, invR * Y, invR * Z);
283 return T/gamma0 + correction;
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);
295 real& deltax, real& deltay, real& deltaz)
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);
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);
314 : Math::NaN<real>()),
316 if (caps & CAP_GAMMA) {
317 _earth.
U(X, Y, Z, fx, fy, fz);
320 gamma = Math::NaN<real>();
321 _earth.
Phi(X, Y, fx, fy);
323 _earth._a, _earth._f, lat, h, Z, X, M[7], M[8],
324 _amodel, _GMmodel, _dzonal0, _corrmult,
327 _gravitational.
Circle(X, Z,
true) :
331 _disturbing.
Circle(-1, X, Z, (caps & CAP_DELTA) != 0) :
334 _correction.
Circle(invR * X, invR * Z,
false) :
340 char* gravitypath = getenv(
"GRAVITY_PATH");
342 path = string(gravitypath);
345 char* datapath = getenv(
"GEOGRAPHICLIB_DATA");
347 path = string(datapath);
353 char* gravityname = getenv(
"GRAVITY_NAME");
355 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
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()
#define GRAVITY_DEFAULT_NAME
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.
Gravity on a circle of latitude.
Math::real MassConstant() const