378 lines
14 KiB
C++
378 lines
14 KiB
C++
/**
|
|
* \file GravityModel.cpp
|
|
* \brief Implementation for GeographicLib::GravityModel class
|
|
*
|
|
* Copyright (c) Charles Karney (2011-2020) <charles@karney.com> and licensed
|
|
* under the MIT/X11 License. For more information, see
|
|
* https://geographiclib.sourceforge.io/
|
|
**********************************************************************/
|
|
|
|
#include <GeographicLib/GravityModel.hpp>
|
|
#include <fstream>
|
|
#include <limits>
|
|
#include <GeographicLib/SphericalEngine.hpp>
|
|
#include <GeographicLib/GravityCircle.hpp>
|
|
#include <GeographicLib/Utility.hpp>
|
|
|
|
#if !defined(GEOGRAPHICLIB_DATA)
|
|
# if defined(_WIN32)
|
|
# define GEOGRAPHICLIB_DATA "C:/ProgramData/GeographicLib"
|
|
# else
|
|
# define GEOGRAPHICLIB_DATA "/usr/local/share/GeographicLib"
|
|
# endif
|
|
#endif
|
|
|
|
#if !defined(GEOGRAPHICLIB_GRAVITY_DEFAULT_NAME)
|
|
# define GEOGRAPHICLIB_GRAVITY_DEFAULT_NAME "egm96"
|
|
#endif
|
|
|
|
#if defined(_MSC_VER)
|
|
// Squelch warnings about unsafe use of getenv
|
|
# pragma warning (disable: 4996)
|
|
#endif
|
|
|
|
namespace GeographicLib {
|
|
|
|
using namespace std;
|
|
|
|
GravityModel::GravityModel(const std::string& name, const std::string& path,
|
|
int Nmax, int Mmax)
|
|
: _name(name)
|
|
, _dir(path)
|
|
, _description("NONE")
|
|
, _date("UNKNOWN")
|
|
, _amodel(Math::NaN())
|
|
, _gGMmodel(Math::NaN())
|
|
, _zeta0(0)
|
|
, _corrmult(1)
|
|
, _nmx(-1)
|
|
, _mmx(-1)
|
|
, _norm(SphericalHarmonic::FULL)
|
|
{
|
|
if (_dir.empty())
|
|
_dir = DefaultGravityPath();
|
|
bool truncate = Nmax >= 0 || Mmax >= 0;
|
|
if (truncate) {
|
|
if (Nmax >= 0 && Mmax < 0) Mmax = Nmax;
|
|
if (Nmax < 0) Nmax = numeric_limits<int>::max();
|
|
if (Mmax < 0) Mmax = numeric_limits<int>::max();
|
|
}
|
|
ReadMetadata(_name);
|
|
{
|
|
string coeff = _filename + ".cof";
|
|
ifstream coeffstr(coeff.c_str(), ios::binary);
|
|
if (!coeffstr.good())
|
|
throw GeographicErr("Error opening " + coeff);
|
|
char id[idlength_ + 1];
|
|
coeffstr.read(id, idlength_);
|
|
if (!coeffstr.good())
|
|
throw GeographicErr("No header in " + coeff);
|
|
id[idlength_] = '\0';
|
|
if (_id != string(id))
|
|
throw GeographicErr("ID mismatch: " + _id + " vs " + id);
|
|
int N, M;
|
|
if (truncate) { N = Nmax; M = Mmax; }
|
|
SphericalEngine::coeff::readcoeffs(coeffstr, N, M, _cCx, _sSx, truncate);
|
|
if (!(N >= 0 && M >= 0))
|
|
throw GeographicErr("Degree and order must be at least 0");
|
|
if (_cCx[0] != 0)
|
|
throw GeographicErr("The degree 0 term should be zero");
|
|
_cCx[0] = 1; // Include the 1/r term in the sum
|
|
_gravitational = SphericalHarmonic(_cCx, _sSx, N, N, M, _amodel, _norm);
|
|
if (truncate) { N = Nmax; M = Mmax; }
|
|
SphericalEngine::coeff::readcoeffs(coeffstr, N, M, _cCC, _cCS, truncate);
|
|
if (N < 0) {
|
|
N = M = 0;
|
|
_cCC.resize(1, real(0));
|
|
}
|
|
_cCC[0] += _zeta0 / _corrmult;
|
|
_correction = SphericalHarmonic(_cCC, _cCS, N, N, M, real(1), _norm);
|
|
int pos = int(coeffstr.tellg());
|
|
coeffstr.seekg(0, ios::end);
|
|
if (pos != coeffstr.tellg())
|
|
throw GeographicErr("Extra data in " + coeff);
|
|
}
|
|
int nmx = _gravitational.Coefficients().nmx();
|
|
_nmx = max(nmx, _correction.Coefficients().nmx());
|
|
_mmx = max(_gravitational.Coefficients().mmx(),
|
|
_correction.Coefficients().mmx());
|
|
// Adjust the normalization of the normal potential to match the model.
|
|
real mult = _earth._gGM / _gGMmodel;
|
|
real amult = Math::sq(_earth._a / _amodel);
|
|
// The 0th term in _zonal should be is 1 + _dzonal0. Instead set it to 1
|
|
// to give exact cancellation with the (0,0) term in the model and account
|
|
// for _dzonal0 separately.
|
|
_zonal.clear(); _zonal.push_back(1);
|
|
_dzonal0 = (_earth.MassConstant() - _gGMmodel) / _gGMmodel;
|
|
for (int n = 2; n <= nmx; n += 2) {
|
|
// Only include as many normal zonal terms as matter. Figuring the limit
|
|
// in this way works because the coefficients of the normal potential
|
|
// (which is smooth) decay much more rapidly that the corresponding
|
|
// coefficient of the model potential (which is bumpy). Typically this
|
|
// goes out to n = 18.
|
|
mult *= amult;
|
|
real
|
|
r = _cCx[n], // the model term
|
|
s = - mult * _earth.Jn(n) / sqrt(real(2 * n + 1)), // the normal term
|
|
t = r - s; // the difference
|
|
if (t == r) // the normal term is negligible
|
|
break;
|
|
_zonal.push_back(0); // index = n - 1; the odd terms are 0
|
|
_zonal.push_back(s);
|
|
}
|
|
int nmx1 = int(_zonal.size()) - 1;
|
|
_disturbing = SphericalHarmonic1(_cCx, _sSx,
|
|
_gravitational.Coefficients().N(),
|
|
nmx, _gravitational.Coefficients().mmx(),
|
|
_zonal,
|
|
_zonal, // This is not accessed!
|
|
nmx1, nmx1, 0,
|
|
_amodel,
|
|
SphericalHarmonic1::normalization(_norm));
|
|
}
|
|
|
|
void GravityModel::ReadMetadata(const string& name) {
|
|
const char* spaces = " \t\n\v\f\r";
|
|
_filename = _dir + "/" + name + ".egm";
|
|
ifstream metastr(_filename.c_str());
|
|
if (!metastr.good())
|
|
throw GeographicErr("Cannot open " + _filename);
|
|
string line;
|
|
getline(metastr, line);
|
|
if (!(line.size() >= 6 && line.substr(0,5) == "EGMF-"))
|
|
throw GeographicErr(_filename + " does not contain EGMF-n signature");
|
|
string::size_type n = line.find_first_of(spaces, 5);
|
|
if (n != string::npos)
|
|
n -= 5;
|
|
string version(line, 5, n);
|
|
if (version != "1")
|
|
throw GeographicErr("Unknown version in " + _filename + ": " + version);
|
|
string key, val;
|
|
real a = Math::NaN(), GM = a, omega = a, f = a, J2 = a;
|
|
while (getline(metastr, line)) {
|
|
if (!Utility::ParseLine(line, key, val))
|
|
continue;
|
|
// Process key words
|
|
if (key == "Name")
|
|
_name = val;
|
|
else if (key == "Description")
|
|
_description = val;
|
|
else if (key == "ReleaseDate")
|
|
_date = val;
|
|
else if (key == "ModelRadius")
|
|
_amodel = Utility::val<real>(val);
|
|
else if (key == "ModelMass")
|
|
_gGMmodel = Utility::val<real>(val);
|
|
else if (key == "AngularVelocity")
|
|
omega = Utility::val<real>(val);
|
|
else if (key == "ReferenceRadius")
|
|
a = Utility::val<real>(val);
|
|
else if (key == "ReferenceMass")
|
|
GM = Utility::val<real>(val);
|
|
else if (key == "Flattening")
|
|
f = Utility::fract<real>(val);
|
|
else if (key == "DynamicalFormFactor")
|
|
J2 = Utility::fract<real>(val);
|
|
else if (key == "HeightOffset")
|
|
_zeta0 = Utility::fract<real>(val);
|
|
else if (key == "CorrectionMultiplier")
|
|
_corrmult = Utility::fract<real>(val);
|
|
else if (key == "Normalization") {
|
|
if (val == "FULL" || val == "Full" || val == "full")
|
|
_norm = SphericalHarmonic::FULL;
|
|
else if (val == "SCHMIDT" || val == "Schmidt" || val == "schmidt")
|
|
_norm = SphericalHarmonic::SCHMIDT;
|
|
else
|
|
throw GeographicErr("Unknown normalization " + val);
|
|
} else if (key == "ByteOrder") {
|
|
if (val == "Big" || val == "big")
|
|
throw GeographicErr("Only little-endian ordering is supported");
|
|
else if (!(val == "Little" || val == "little"))
|
|
throw GeographicErr("Unknown byte ordering " + val);
|
|
} else if (key == "ID")
|
|
_id = val;
|
|
// else unrecognized keywords are skipped
|
|
}
|
|
// Check values
|
|
if (!(isfinite(_amodel) && _amodel > 0))
|
|
throw GeographicErr("Model radius must be positive");
|
|
if (!(isfinite(_gGMmodel) && _gGMmodel > 0))
|
|
throw GeographicErr("Model mass constant must be positive");
|
|
if (!(isfinite(_corrmult) && _corrmult > 0))
|
|
throw GeographicErr("Correction multiplier must be positive");
|
|
if (!(isfinite(_zeta0)))
|
|
throw GeographicErr("Height offset must be finite");
|
|
if (int(_id.size()) != idlength_)
|
|
throw GeographicErr("Invalid ID");
|
|
if (isfinite(f) && isfinite(J2))
|
|
throw GeographicErr("Cannot specify both f and J2");
|
|
_earth = NormalGravity(a, GM, omega,
|
|
isfinite(f) ? f : J2, isfinite(f));
|
|
}
|
|
|
|
Math::real GravityModel::InternalT(real X, real Y, real Z,
|
|
real& deltaX, real& deltaY, real& deltaZ,
|
|
bool gradp, bool correct) const {
|
|
// If correct, then produce the correct T = W - U. Otherwise, neglect the
|
|
// n = 0 term (which is proportial to the difference in the model and
|
|
// reference values of GM).
|
|
if (_dzonal0 == 0)
|
|
// No need to do the correction
|
|
correct = false;
|
|
real T, invR = correct ? 1 / hypot(hypot(X, Y), Z) : 1;
|
|
if (gradp) {
|
|
// initial values to suppress warnings
|
|
deltaX = deltaY = deltaZ = 0;
|
|
T = _disturbing(-1, X, Y, Z, deltaX, deltaY, deltaZ);
|
|
real f = _gGMmodel / _amodel;
|
|
deltaX *= f;
|
|
deltaY *= f;
|
|
deltaZ *= f;
|
|
if (correct) {
|
|
invR = _gGMmodel * _dzonal0 * invR * invR * invR;
|
|
deltaX += X * invR;
|
|
deltaY += Y * invR;
|
|
deltaZ += Z * invR;
|
|
}
|
|
} else
|
|
T = _disturbing(-1, X, Y, Z);
|
|
T = (T / _amodel - (correct ? _dzonal0 : 0) * invR) * _gGMmodel;
|
|
return T;
|
|
}
|
|
|
|
Math::real GravityModel::V(real X, real Y, real Z,
|
|
real& GX, real& GY, real& GZ) const {
|
|
real
|
|
Vres = _gravitational(X, Y, Z, GX, GY, GZ),
|
|
f = _gGMmodel / _amodel;
|
|
Vres *= f;
|
|
GX *= f;
|
|
GY *= f;
|
|
GZ *= f;
|
|
return Vres;
|
|
}
|
|
|
|
Math::real GravityModel::W(real X, real Y, real Z,
|
|
real& gX, real& gY, real& gZ) const {
|
|
real fX, fY,
|
|
Wres = V(X, Y, Z, gX, gY, gZ) + _earth.Phi(X, Y, fX, fY);
|
|
gX += fX;
|
|
gY += fY;
|
|
return Wres;
|
|
}
|
|
|
|
void GravityModel::SphericalAnomaly(real lat, real lon, real h,
|
|
real& Dg01, real& xi, real& eta) const {
|
|
real X, Y, Z, M[Geocentric::dim2_];
|
|
_earth.Earth().IntForward(lat, lon, h, X, Y, Z, M);
|
|
real
|
|
deltax, deltay, deltaz,
|
|
T = InternalT(X, Y, Z, deltax, deltay, deltaz, true, false),
|
|
clam = M[3], slam = -M[0],
|
|
P = hypot(X, Y),
|
|
R = hypot(P, Z),
|
|
// psi is geocentric latitude
|
|
cpsi = R != 0 ? P / R : M[7],
|
|
spsi = R != 0 ? Z / R : M[8];
|
|
// Rotate cartesian into spherical coordinates
|
|
real MC[Geocentric::dim2_];
|
|
Geocentric::Rotation(spsi, cpsi, slam, clam, MC);
|
|
Geocentric::Unrotate(MC, deltax, deltay, deltaz, deltax, deltay, deltaz);
|
|
// H+M, Eq 2-151c
|
|
Dg01 = - deltaz - 2 * T / R;
|
|
real gammaX, gammaY, gammaZ;
|
|
_earth.U(X, Y, Z, gammaX, gammaY, gammaZ);
|
|
real gamma = hypot( hypot(gammaX, gammaY), gammaZ);
|
|
xi = -(deltay/gamma) / Math::degree();
|
|
eta = -(deltax/gamma) / Math::degree();
|
|
}
|
|
|
|
Math::real GravityModel::GeoidHeight(real lat, real lon) const
|
|
{
|
|
real X, Y, Z;
|
|
_earth.Earth().IntForward(lat, lon, 0, X, Y, Z, NULL);
|
|
real
|
|
gamma0 = _earth.SurfaceGravity(lat),
|
|
dummy,
|
|
T = InternalT(X, Y, Z, dummy, dummy, dummy, false, false),
|
|
invR = 1 / hypot(hypot(X, Y), Z),
|
|
correction = _corrmult * _correction(invR * X, invR * Y, invR * Z);
|
|
// _zeta0 has been included in _correction
|
|
return T/gamma0 + correction;
|
|
}
|
|
|
|
Math::real GravityModel::Gravity(real lat, real lon, real h,
|
|
real& gx, real& gy, real& gz) const {
|
|
real X, Y, Z, M[Geocentric::dim2_];
|
|
_earth.Earth().IntForward(lat, lon, h, X, Y, Z, M);
|
|
real Wres = W(X, Y, Z, gx, gy, gz);
|
|
Geocentric::Unrotate(M, gx, gy, gz, gx, gy, gz);
|
|
return Wres;
|
|
}
|
|
Math::real GravityModel::Disturbance(real lat, real lon, real h,
|
|
real& deltax, real& deltay,
|
|
real& deltaz) const {
|
|
real X, Y, Z, M[Geocentric::dim2_];
|
|
_earth.Earth().IntForward(lat, lon, h, X, Y, Z, M);
|
|
real Tres = InternalT(X, Y, Z, deltax, deltay, deltaz, true, true);
|
|
Geocentric::Unrotate(M, deltax, deltay, deltaz, deltax, deltay, deltaz);
|
|
return Tres;
|
|
}
|
|
|
|
GravityCircle GravityModel::Circle(real lat, real h, unsigned caps) const {
|
|
if (h != 0)
|
|
// Disallow invoking GeoidHeight unless h is zero.
|
|
caps &= ~(CAP_GAMMA0 | CAP_C);
|
|
real X, Y, Z, M[Geocentric::dim2_];
|
|
_earth.Earth().IntForward(lat, 0, h, X, Y, Z, M);
|
|
// Y = 0, cphi = M[7], sphi = M[8];
|
|
real
|
|
invR = 1 / hypot(X, Z),
|
|
gamma0 = (caps & CAP_GAMMA0 ?_earth.SurfaceGravity(lat)
|
|
: Math::NaN()),
|
|
fx, fy, fz, gamma;
|
|
if (caps & CAP_GAMMA) {
|
|
_earth.U(X, Y, Z, fx, fy, fz); // fy = 0
|
|
gamma = hypot(fx, fz);
|
|
} else
|
|
gamma = Math::NaN();
|
|
_earth.Phi(X, Y, fx, fy);
|
|
return GravityCircle(GravityCircle::mask(caps),
|
|
_earth._a, _earth._f, lat, h, Z, X, M[7], M[8],
|
|
_amodel, _gGMmodel, _dzonal0, _corrmult,
|
|
gamma0, gamma, fx,
|
|
caps & CAP_G ?
|
|
_gravitational.Circle(X, Z, true) :
|
|
CircularEngine(),
|
|
// N.B. If CAP_DELTA is set then CAP_T should be too.
|
|
caps & CAP_T ?
|
|
_disturbing.Circle(-1, X, Z, (caps&CAP_DELTA) != 0) :
|
|
CircularEngine(),
|
|
caps & CAP_C ?
|
|
_correction.Circle(invR * X, invR * Z, false) :
|
|
CircularEngine());
|
|
}
|
|
|
|
string GravityModel::DefaultGravityPath() {
|
|
string path;
|
|
char* gravitypath = getenv("GEOGRAPHICLIB_GRAVITY_PATH");
|
|
if (gravitypath)
|
|
path = string(gravitypath);
|
|
if (!path.empty())
|
|
return path;
|
|
char* datapath = getenv("GEOGRAPHICLIB_DATA");
|
|
if (datapath)
|
|
path = string(datapath);
|
|
return (!path.empty() ? path : string(GEOGRAPHICLIB_DATA)) + "/gravity";
|
|
}
|
|
|
|
string GravityModel::DefaultGravityName() {
|
|
string name;
|
|
char* gravityname = getenv("GEOGRAPHICLIB_GRAVITY_NAME");
|
|
if (gravityname)
|
|
name = string(gravityname);
|
|
return !name.empty() ? name : string(GEOGRAPHICLIB_GRAVITY_DEFAULT_NAME);
|
|
}
|
|
|
|
} // namespace GeographicLib
|