00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include <GeographicLib/GravityModel.hpp>
00011 #include <fstream>
00012 #include <GeographicLib/SphericalEngine.hpp>
00013 #include <GeographicLib/GravityCircle.hpp>
00014 #include <GeographicLib/Utility.hpp>
00015
00016 #if !defined(GEOGRAPHICLIB_DATA)
00017 # if defined(_WIN32)
00018 # define GEOGRAPHICLIB_DATA "C:/ProgramData/GeographicLib"
00019 # else
00020 # define GEOGRAPHICLIB_DATA "/usr/local/share/GeographicLib"
00021 # endif
00022 #endif
00023
00024 #if !defined(GEOGRAPHICLIB_GRAVITY_DEFAULT_NAME)
00025 # define GEOGRAPHICLIB_GRAVITY_DEFAULT_NAME "egm96"
00026 #endif
00027
00028 #if defined(_MSC_VER)
00029
00030 # pragma warning (disable: 4996)
00031 #endif
00032
00033 namespace GeographicLib {
00034
00035 using namespace std;
00036
00037 GravityModel::GravityModel(const std::string& name,const std::string& path)
00038 : _name(name)
00039 , _dir(path)
00040 , _description("NONE")
00041 , _date("UNKNOWN")
00042 , _amodel(Math::NaN())
00043 , _GMmodel(Math::NaN())
00044 , _zeta0(0)
00045 , _corrmult(1)
00046 , _norm(SphericalHarmonic::FULL)
00047 {
00048 if (_dir.empty())
00049 _dir = DefaultGravityPath();
00050 ReadMetadata(_name);
00051 {
00052 string coeff = _filename + ".cof";
00053 ifstream coeffstr(coeff.c_str(), ios::binary);
00054 if (!coeffstr.good())
00055 throw GeographicErr("Error opening " + coeff);
00056 char id[idlength_ + 1];
00057 coeffstr.read(id, idlength_);
00058 if (!coeffstr.good())
00059 throw GeographicErr("No header in " + coeff);
00060 id[idlength_] = '\0';
00061 if (_id != string(id))
00062 throw GeographicErr("ID mismatch: " + _id + " vs " + id);
00063 int N, M;
00064 SphericalEngine::coeff::readcoeffs(coeffstr, N, M, _Cx, _Sx);
00065 if (!(M < 0 || _Cx[0] == 0))
00066 throw GeographicErr("A degree 0 term should be zero");
00067 _Cx[0] = 1;
00068 _gravitational = SphericalHarmonic(_Cx, _Sx, N, N, M, _amodel, _norm);
00069 SphericalEngine::coeff::readcoeffs(coeffstr, N, M, _CC, _CS);
00070 if (N < 0) {
00071 N = M = 0;
00072 _CC.resize(1, real(0));
00073 }
00074 _CC[0] += _zeta0 / _corrmult;
00075 _correction = SphericalHarmonic(_CC, _CS, N, N, M, real(1), _norm);
00076 int pos = int(coeffstr.tellg());
00077 coeffstr.seekg(0, ios::end);
00078 if (pos != coeffstr.tellg())
00079 throw GeographicErr("Extra data in " + coeff);
00080 }
00081 int nmx = _gravitational.Coefficients().nmx();
00082
00083 real mult = _earth._GM / _GMmodel;
00084 real amult = Math::sq(_earth._a / _amodel);
00085
00086
00087
00088 _zonal.clear(); _zonal.push_back(1);
00089 _dzonal0 = (_earth.MassConstant() - _GMmodel) / _GMmodel;
00090 for (int n = 2; n <= nmx; n += 2) {
00091
00092
00093
00094
00095
00096 mult *= amult;
00097 real
00098 r = _Cx[n],
00099 s = - mult * _earth.Jn(n) / sqrt(real(2 * n + 1)),
00100 t = r - s;
00101 if (t == r)
00102 break;
00103 _zonal.push_back(0);
00104 _zonal.push_back(s);
00105 }
00106 int nmx1 = int(_zonal.size()) - 1;
00107 _disturbing = SphericalHarmonic1(_Cx, _Sx,
00108 _gravitational.Coefficients().N(),
00109 nmx, _gravitational.Coefficients().mmx(),
00110 _zonal,
00111 _zonal,
00112 nmx1, nmx1, 0,
00113 _amodel,
00114 SphericalHarmonic1::normalization(_norm));
00115 }
00116
00117 void GravityModel::ReadMetadata(const std::string& name) {
00118 const char* spaces = " \t\n\v\f\r";
00119 _filename = _dir + "/" + name + ".egm";
00120 ifstream metastr(_filename.c_str());
00121 if (!metastr.good())
00122 throw GeographicErr("Cannot open " + _filename);
00123 string line;
00124 getline(metastr, line);
00125 if (!(line.size() >= 6 && line.substr(0,5) == "EGMF-"))
00126 throw GeographicErr(_filename + " does not contain EGMF-n signature");
00127 string::size_type n = line.find_first_of(spaces, 5);
00128 if (n != string::npos)
00129 n -= 5;
00130 string version = line.substr(5, n);
00131 if (version != "1")
00132 throw GeographicErr("Unknown version in " + _filename + ": " + version);
00133 string key, val;
00134 real a = Math::NaN(), GM = a, omega = a, f = a, J2 = a;
00135 while (getline(metastr, line)) {
00136 if (!Utility::ParseLine(line, key, val))
00137 continue;
00138
00139 if (key == "Name")
00140 _name = val;
00141 else if (key == "Description")
00142 _description = val;
00143 else if (key == "ReleaseDate")
00144 _date = val;
00145 else if (key == "ModelRadius")
00146 _amodel = Utility::num<real>(val);
00147 else if (key == "ModelMass")
00148 _GMmodel = Utility::num<real>(val);
00149 else if (key == "AngularVelocity")
00150 omega = Utility::num<real>(val);
00151 else if (key == "ReferenceRadius")
00152 a = Utility::num<real>(val);
00153 else if (key == "ReferenceMass")
00154 GM = Utility::num<real>(val);
00155 else if (key == "Flattening")
00156 f = Utility::fract<real>(val);
00157 else if (key == "DynamicalFormFactor")
00158 J2 = Utility::fract<real>(val);
00159 else if (key == "HeightOffset")
00160 _zeta0 = Utility::fract<real>(val);
00161 else if (key == "CorrectionMultiplier")
00162 _corrmult = Utility::fract<real>(val);
00163 else if (key == "Normalization") {
00164 if (val == "FULL" || val == "Full" || val == "full")
00165 _norm = SphericalHarmonic::FULL;
00166 else if (val == "SCHMIDT" || val == "Schmidt" || val == "schmidt")
00167 _norm = SphericalHarmonic::SCHMIDT;
00168 else
00169 throw GeographicErr("Unknown normalization " + val);
00170 } else if (key == "ByteOrder") {
00171 if (val == "Big" || val == "big")
00172 throw GeographicErr("Only little-endian ordering is supported");
00173 else if (!(val == "Little" || val == "little"))
00174 throw GeographicErr("Unknown byte ordering " + val);
00175 } else if (key == "ID")
00176 _id = val;
00177
00178 }
00179
00180 if (!(Math::isfinite(_amodel) && _amodel > 0))
00181 throw GeographicErr("Model radius must be positive");
00182 if (!(Math::isfinite(_GMmodel) && _GMmodel > 0))
00183 throw GeographicErr("Model mass constant must be positive");
00184 if (!(Math::isfinite(_corrmult) && _corrmult > 0))
00185 throw GeographicErr("Correction multiplier must be positive");
00186 if (!(Math::isfinite(_zeta0)))
00187 throw GeographicErr("Height offset must be finite");
00188 if (int(_id.size()) != idlength_)
00189 throw GeographicErr("Invalid ID");
00190 _earth = NormalGravity(a, GM, omega, f, J2);
00191 }
00192
00193 Math::real GravityModel::InternalT(real X, real Y, real Z,
00194 real& deltaX, real& deltaY, real& deltaZ,
00195 bool gradp, bool correct) const {
00196
00197
00198
00199 if (_dzonal0 == 0)
00200
00201 correct = false;
00202 real T, invR = correct ? 1 / Math::hypot(Math::hypot(X, Y), Z) : 1;
00203 if (gradp) {
00204
00205 deltaX = deltaY = deltaZ = 0;
00206 T = _disturbing(-1, X, Y, Z, deltaX, deltaY, deltaZ);
00207 real f = _GMmodel / _amodel;
00208 deltaX *= f;
00209 deltaY *= f;
00210 deltaZ *= f;
00211 if (correct) {
00212 invR = _GMmodel * _dzonal0 * invR * invR * invR;
00213 deltaX += X * invR;
00214 deltaY += Y * invR;
00215 deltaZ += Z * invR;
00216 }
00217 } else
00218 T = _disturbing(-1, X, Y, Z);
00219 T = (T / _amodel - (correct ? _dzonal0 : 0) * invR) * _GMmodel;
00220 return T;
00221 }
00222
00223 Math::real GravityModel::V(real X, real Y, real Z,
00224 real& GX, real& GY, real& GZ) const {
00225 real
00226 Vres = _gravitational(X, Y, Z, GX, GY, GZ),
00227 f = _GMmodel / _amodel;
00228 Vres *= f;
00229 GX *= f;
00230 GY *= f;
00231 GZ *= f;
00232 return Vres;
00233 }
00234
00235 Math::real GravityModel::W(real X, real Y, real Z,
00236 real& gX, real& gY, real& gZ) const {
00237 real fX, fY,
00238 Wres = V(X, Y, Z, gX, gY, gZ) + _earth.Phi(X, Y, fX, fY);
00239 gX += fX;
00240 gY += fY;
00241 return Wres;
00242 }
00243
00244 void GravityModel::SphericalAnomaly(real lat, real lon, real h,
00245 real& Dg01, real& xi, real& eta)
00246 const {
00247 real X, Y, Z, M[Geocentric::dim2_];
00248 _earth.Earth().IntForward(lat, lon, h, X, Y, Z, M);
00249 real
00250 deltax, deltay, deltaz,
00251 T = InternalT(X, Y, Z, deltax, deltay, deltaz, true, false),
00252 clam = M[3], slam = -M[0],
00253 P = Math::hypot(X, Y),
00254 R = Math::hypot(P, Z),
00255
00256 cpsi = R ? P / R : M[7],
00257 spsi = R ? Z / R : M[8];
00258
00259 real MC[Geocentric::dim2_];
00260 Geocentric::Rotation(spsi, cpsi, slam, clam, MC);
00261 Geocentric::Unrotate(MC, deltax, deltay, deltaz, deltax, deltay, deltaz);
00262
00263 Dg01 = - deltaz - 2 * T / R;
00264 real gammaX, gammaY, gammaZ;
00265 _earth.U(X, Y, Z, gammaX, gammaY, gammaZ);
00266 real gamma = Math::hypot( Math::hypot(gammaX, gammaY), gammaZ);
00267 xi = -(deltay/gamma) / Math::degree();
00268 eta = -(deltax/gamma) / Math::degree();
00269 }
00270
00271 Math::real GravityModel::GeoidHeight(real lat, real lon) const
00272 {
00273 real X, Y, Z;
00274 _earth.Earth().IntForward(lat, lon, 0, X, Y, Z, NULL);
00275 real
00276 gamma0 = _earth.SurfaceGravity(lat),
00277 dummy,
00278 T = InternalT(X, Y, Z, dummy, dummy, dummy, false, false),
00279 invR = 1 / Math::hypot(Math::hypot(X, Y), Z),
00280 correction = _corrmult * _correction(invR * X, invR * Y, invR * Z);
00281
00282 return T/gamma0 + correction;
00283 }
00284
00285 Math::real GravityModel::Gravity(real lat, real lon, real h,
00286 real& gx, real& gy, real& gz) const {
00287 real X, Y, Z, M[Geocentric::dim2_];
00288 _earth.Earth().IntForward(lat, lon, h, X, Y, Z, M);
00289 real Wres = W(X, Y, Z, gx, gy, gz);
00290 Geocentric::Unrotate(M, gx, gy, gz, gx, gy, gz);
00291 return Wres;
00292 }
00293 Math::real GravityModel::Disturbance(real lat, real lon, real h,
00294 real& deltax, real& deltay, real& deltaz)
00295 const {
00296 real X, Y, Z, M[Geocentric::dim2_];
00297 _earth.Earth().IntForward(lat, lon, h, X, Y, Z, M);
00298 real Tres = InternalT(X, Y, Z, deltax, deltay, deltaz, true, true);
00299 Geocentric::Unrotate(M, deltax, deltay, deltaz, deltax, deltay, deltaz);
00300 return Tres;
00301 }
00302
00303 GravityCircle GravityModel::Circle(real lat, real h, unsigned caps) const {
00304 if (h != 0)
00305
00306 caps &= ~(CAP_GAMMA0 | CAP_C);
00307 real X, Y, Z, M[Geocentric::dim2_];
00308 _earth.Earth().IntForward(lat, 0, h, X, Y, Z, M);
00309
00310 real
00311 invR = 1 / Math::hypot(X, Z),
00312 gamma0 = (caps & CAP_GAMMA0 ?_earth.SurfaceGravity(lat)
00313 : Math::NaN()),
00314 fx, fy, fz, gamma;
00315 if (caps & CAP_GAMMA) {
00316 _earth.U(X, Y, Z, fx, fy, fz);
00317 gamma = Math::hypot(fx, fz);
00318 } else
00319 gamma = Math::NaN();
00320 _earth.Phi(X, Y, fx, fy);
00321 return GravityCircle(GravityCircle::mask(caps),
00322 _earth._a, _earth._f, lat, h, Z, X, M[7], M[8],
00323 _amodel, _GMmodel, _dzonal0, _corrmult,
00324 gamma0, gamma, fx,
00325 caps & CAP_G ?
00326 _gravitational.Circle(X, Z, true) :
00327 CircularEngine(),
00328
00329 caps & CAP_T ?
00330 _disturbing.Circle(-1, X, Z, (caps & CAP_DELTA) != 0) :
00331 CircularEngine(),
00332 caps & CAP_C ?
00333 _correction.Circle(invR * X, invR * Z, false) :
00334 CircularEngine());
00335 }
00336
00337 std::string GravityModel::DefaultGravityPath() {
00338 string path;
00339 char* gravitypath = getenv("GEOGRAPHICLIB_GRAVITY_PATH");
00340 if (gravitypath)
00341 path = string(gravitypath);
00342 if (!path.empty())
00343 return path;
00344 char* datapath = getenv("GEOGRAPHICLIB_DATA");
00345 if (datapath)
00346 path = string(datapath);
00347 return (!path.empty() ? path : string(GEOGRAPHICLIB_DATA)) + "/gravity";
00348 }
00349
00350 std::string GravityModel::DefaultGravityName() {
00351 string name;
00352 char* gravityname = getenv("GEOGRAPHICLIB_GRAVITY_NAME");
00353 if (gravityname)
00354 name = string(gravityname);
00355 return !name.empty() ? name : string(GEOGRAPHICLIB_GRAVITY_DEFAULT_NAME);
00356 }
00357
00358 }