GeographicLib  1.21
GravityModel.cpp
Go to the documentation of this file.
00001 /**
00002  * \file GravityModel.cpp
00003  * \brief Implementation for GeographicLib::GravityModel class
00004  *
00005  * Copyright (c) Charles Karney (2011, 2012) <charles@karney.com> and licensed
00006  * under the MIT/X11 License.  For more information, see
00007  * http://geographiclib.sourceforge.net/
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 #define GEOGRAPHICLIB_GRAVITYMODEL_CPP \
00016   "$Id: 1897d0d53c7339ecdf20b1348637340e9f687f30 $"
00017 
00018 RCSID_DECL(GEOGRAPHICLIB_GRAVITYMODEL_CPP)
00019 RCSID_DECL(GEOGRAPHICLIB_GRAVITYMODEL_HPP)
00020 
00021 #if !defined(GEOGRAPHICLIB_DATA)
00022 #  if defined(_MSC_VER)
00023 #    define GEOGRAPHICLIB_DATA \
00024   "C:/Documents and Settings/All Users/Application Data/GeographicLib"
00025 #  else
00026 #    define GEOGRAPHICLIB_DATA "/usr/local/share/GeographicLib"
00027 #  endif
00028 #endif
00029 
00030 #if !defined(GRAVITY_DEFAULT_NAME)
00031 #  define GRAVITY_DEFAULT_NAME "egm96"
00032 #endif
00033 
00034 #if defined(_MSC_VER)
00035 // Squelch warnings about unsafe use of getenv
00036 #pragma warning (disable: 4996)
00037 #endif
00038 
00039 namespace GeographicLib {
00040 
00041   using namespace std;
00042 
00043   GravityModel::GravityModel(const std::string& name,const std::string& path)
00044     : _name(name)
00045     , _dir(path)
00046     , _description("NONE")
00047     , _date("UNKNOWN")
00048     , _amodel(Math::NaN<real>())
00049     , _GMmodel(Math::NaN<real>())
00050     , _zeta0(0)
00051     , _corrmult(1)
00052     , _norm(SphericalHarmonic::FULL)
00053   {
00054     if (_dir.empty())
00055       _dir = DefaultGravityPath();
00056     ReadMetadata(_name);
00057     {
00058       string coeff = _filename + ".cof";
00059       ifstream coeffstr(coeff.c_str(), ios::binary);
00060       if (!coeffstr.good())
00061         throw GeographicErr("Error opening " + coeff);
00062       char id[idlength_ + 1];
00063       coeffstr.read(id, idlength_);
00064       if (!coeffstr.good())
00065         throw GeographicErr("No header in " + coeff);
00066       id[idlength_] = '\0';
00067       if (_id != string(id))
00068         throw GeographicErr("ID mismatch: " + _id + " vs " + id);
00069       int N, M;
00070       SphericalEngine::coeff::readcoeffs(coeffstr, N, M, _C, _S);
00071       if (!(M < 0 || _C[0] == 0))
00072         throw GeographicErr("A degree 0 term should be zero");
00073       _C[0] = 1;                // Include the 1/r term in the sum
00074       _gravitational = SphericalHarmonic(_C, _S, N, N, M, _amodel, _norm);
00075       SphericalEngine::coeff::readcoeffs(coeffstr, N, M, _CC, _CS);
00076       if (N < 0) {
00077         N = M = 0;
00078         _CC.resize(1, real(0));
00079       }
00080       _CC[0] += _zeta0 / _corrmult;
00081       _correction = SphericalHarmonic(_CC, _CS, N, N, M, real(1), _norm);
00082       int pos = int(coeffstr.tellg());
00083       coeffstr.seekg(0, ios::end);
00084       if (pos != coeffstr.tellg())
00085         throw GeographicErr("Extra data in " + coeff);
00086     }
00087     int nmx = _gravitational.Coefficients().nmx();
00088     // Adjust the normalization of the normal potential to match the model.
00089     real mult = _earth._GM / _GMmodel;
00090     real amult = Math::sq(_earth._a / _amodel);
00091     // The 0th term in _zonal should be is 1 + _dzonal0.  Instead set it to 1
00092     // to give exact cancellation with the (0,0) term in the model and account
00093     // for _dzonal0 separately.
00094     _zonal.clear(); _zonal.push_back(1);
00095     _dzonal0 = (_earth.MassConstant() - _GMmodel) / _GMmodel;
00096     for (int n = 2; n <= nmx; n += 2) {
00097       // Only include as many normal zonal terms as matter.  Figuring the limit
00098       // in this way works because the coefficients of the normal potential
00099       // (which is smooth) decay much more rapidly that the corresponding
00100       // coefficient of the model potential (which is bumpy).  Typically this
00101       // goes out to n = 18.
00102       mult *= amult;
00103       real
00104         r = _C[n],                                         // the model term
00105         s = - mult * _earth.Jn(n) / sqrt(real(2 * n + 1)), // the normal term
00106         t = r - s;                                         // the difference
00107       if (t == r)               // the normal term is negligible
00108         break;
00109       _zonal.push_back(0);      // index = n - 1; the odd terms are 0
00110       _zonal.push_back(s);
00111     }
00112     int nmx1 = int(_zonal.size()) - 1;
00113     _disturbing = SphericalHarmonic1(_C, _S,
00114                                      _gravitational.Coefficients().N(),
00115                                      nmx, _gravitational.Coefficients().mmx(),
00116                                      _zonal,
00117                                      _zonal, // This is not accessed!
00118                                      nmx1, nmx1, 0,
00119                                      _amodel,
00120                                      SphericalHarmonic1::normalization(_norm));
00121   }
00122 
00123   void GravityModel::ReadMetadata(const std::string& name) {
00124     const char* spaces = " \t\n\v\f\r";
00125     _filename = _dir + "/" + name + ".egm";
00126     ifstream metastr(_filename.c_str());
00127     if (!metastr.good())
00128       throw GeographicErr("Cannot open " + _filename);
00129     string line;
00130     getline(metastr, line);
00131     if (!(line.size() >= 6 && line.substr(0,5) == "EGMF-"))
00132       throw GeographicErr(_filename + " does not contain EGMF-n signature");
00133     string::size_type n = line.find_first_of(spaces, 5);
00134     if (n != string::npos)
00135       n -= 5;
00136     string version = line.substr(5, n);
00137     if (version != "1")
00138       throw GeographicErr("Unknown version in " + _filename + ": " + version);
00139     string key, val;
00140     real a = Math::NaN<real>(), GM = a, omega = a, f = a, J2 = a;
00141     while (getline(metastr, line)) {
00142       if (!Utility::ParseLine(line, key, val))
00143         continue;
00144       // Process key words
00145       if (key == "Name")
00146         _name = val;
00147       else if (key == "Description")
00148         _description = val;
00149       else if (key == "ReleaseDate")
00150         _date = val;
00151       else if (key == "ModelRadius")
00152         _amodel = Utility::num<real>(val);
00153       else if (key == "ModelMass")
00154         _GMmodel = Utility::num<real>(val);
00155       else if (key == "AngularVelocity")
00156         omega = Utility::num<real>(val);
00157       else if (key == "ReferenceRadius")
00158         a = Utility::num<real>(val);
00159       else if (key == "ReferenceMass")
00160         GM = Utility::num<real>(val);
00161       else if (key == "Flattening")
00162         f = Utility::fract<real>(val);
00163       else if (key == "DynamicalFormFactor")
00164         J2 = Utility::fract<real>(val);
00165       else if (key == "HeightOffset")
00166         _zeta0 = Utility::fract<real>(val);
00167       else if (key == "CorrectionMultiplier")
00168         _corrmult = Utility::fract<real>(val);
00169       else if (key == "Normalization") {
00170         if (val == "FULL" || val == "Full" || val == "full")
00171           _norm = SphericalHarmonic::FULL;
00172         else if (val == "SCHMIDT" || val == "Schmidt" || val == "schmidt")
00173           _norm = SphericalHarmonic::SCHMIDT;
00174         else
00175           throw GeographicErr("Unknown normalization " + val);
00176       } else if (key == "ByteOrder") {
00177         if (val == "Big" || val == "big")
00178           throw GeographicErr("Only little-endian ordering is supported");
00179         else if (!(val == "Little" || val == "little"))
00180           throw GeographicErr("Unknown byte ordering " + val);
00181       } else if (key == "ID")
00182         _id = val;
00183       // else unrecognized keywords are skipped
00184     }
00185     // Check values
00186     if (!(Math::isfinite(_amodel) && _amodel > 0))
00187       throw GeographicErr("Model radius must be positive");
00188     if (!(Math::isfinite(_GMmodel) && _GMmodel > 0))
00189       throw GeographicErr("Model mass constant must be positive");
00190     if (!(Math::isfinite(_corrmult) && _corrmult > 0))
00191       throw GeographicErr("Correction multiplier must be positive");
00192     if (!(Math::isfinite(_zeta0)))
00193       throw GeographicErr("Height offset must be finite");
00194     if (int(_id.size()) != idlength_)
00195       throw GeographicErr("Invalid ID");
00196     _earth = NormalGravity(a, GM, omega, f, J2);
00197   }
00198 
00199   Math::real GravityModel::InternalT(real X, real Y, real Z,
00200                                      real& deltaX, real& deltaY, real& deltaZ,
00201                                      bool gradp, bool correct) const throw() {
00202     // If correct, then produce the correct T = W - U.  Otherwise, neglect the
00203     // n = 0 term (which is proportial to the difference in the model and
00204     // reference values of GM).
00205     if (_dzonal0 == 0)
00206       // No need to do the correction
00207       correct = false;
00208     real
00209       invR = correct ? 1 / Math::hypot(Math::hypot(X, Y), Z) : 1,
00210       T = (gradp
00211            ? _disturbing(-1, X, Y, Z, deltaX, deltaY, deltaZ)
00212            : _disturbing(-1, X, Y, Z));
00213     T = (T / _amodel - (correct ? _dzonal0 : 0) * invR) * _GMmodel;
00214     if (gradp) {
00215       real f = _GMmodel / _amodel;
00216       deltaX *= f;
00217       deltaY *= f;
00218       deltaZ *= f;
00219       if (correct) {
00220         invR = _GMmodel * _dzonal0 * invR * invR * invR;
00221         deltaX += X * invR;
00222         deltaY += Y * invR;
00223         deltaZ += Z * invR;
00224       }
00225     }
00226     return T;
00227   }
00228 
00229   Math::real GravityModel::V(real X, real Y, real Z,
00230                              real& GX, real& GY, real& GZ) const throw() {
00231     real
00232       Vres = _gravitational(X, Y, Z, GX, GY, GZ),
00233       f = _GMmodel / _amodel;
00234     Vres *= f;
00235     GX *= f;
00236     GY *= f;
00237     GZ *= f;
00238     return Vres;
00239   }
00240 
00241   Math::real GravityModel::W(real X, real Y, real Z,
00242                              real& gX, real& gY, real& gZ) const throw() {
00243     real fX, fY,
00244       Wres = V(X, Y, Z, gX, gY, gZ) + _earth.Phi(X, Y, fX, fY);
00245     gX += fX;
00246     gY += fY;
00247     return Wres;
00248   }
00249 
00250   void GravityModel::SphericalAnomaly(real lat, real lon, real h,
00251                                       real& Dg01, real& xi, real& eta)
00252     const throw() {
00253     real X, Y, Z, M[Geocentric::dim2_];
00254     _earth.Earth().IntForward(lat, lon, h, X, Y, Z, M);
00255     real
00256       deltax, deltay, deltaz,
00257       T = InternalT(X, Y, Z, deltax, deltay, deltaz, true, false),
00258       clam = M[3], slam = -M[0],
00259       P = Math::hypot(X, Y),
00260       R = Math::hypot(P, Z),
00261       // psi is geocentric latitude
00262       cpsi = R ? P / R : M[7],
00263       spsi = R ? Z / R : M[8];
00264     // Rotate cartesian into spherical coordinates
00265     real MC[Geocentric::dim2_];
00266     Geocentric::Rotation(spsi, cpsi, slam, clam, MC);
00267     Geocentric::Unrotate(MC, deltax, deltay, deltaz, deltax, deltay, deltaz);
00268     // H+M, Eq 2-151c
00269     Dg01 = - deltaz - 2 * T / R;
00270     real gammaX, gammaY, gammaZ;
00271     _earth.U(X, Y, Z, gammaX, gammaY, gammaZ);
00272     real gamma = Math::hypot( Math::hypot(gammaX, gammaY), gammaZ);
00273     xi  = -(deltay/gamma) / Math::degree<real>();
00274     eta = -(deltax/gamma) / Math::degree<real>();
00275   }
00276 
00277   Math::real GravityModel::GeoidHeight(real lat, real lon) const throw()
00278   {
00279     real X, Y, Z;
00280     _earth.Earth().IntForward(lat, lon, 0, X, Y, Z, NULL);
00281     real
00282       gamma0 = _earth.SurfaceGravity(lat),
00283       dummy,
00284       T = InternalT(X, Y, Z, dummy, dummy, dummy, false, false),
00285       invR = 1 / Math::hypot(Math::hypot(X, Y), Z),
00286       correction = _corrmult * _correction(invR * X, invR * Y, invR * Z);
00287     // _zeta0 has been included in _correction
00288     return T/gamma0 + correction;
00289   }
00290 
00291   Math::real GravityModel::Gravity(real lat, real lon, real h,
00292                                    real& gx, real& gy, real& gz) const throw() {
00293     real X, Y, Z, M[Geocentric::dim2_];
00294     _earth.Earth().IntForward(lat, lon, h, X, Y, Z, M);
00295     real Wres = W(X, Y, Z, gx, gy, gz);
00296     Geocentric::Unrotate(M, gx, gy, gz, gx, gy, gz);
00297     return Wres;
00298   }
00299   Math::real GravityModel::Disturbance(real lat, real lon, real h,
00300                                        real& deltax, real& deltay, real& deltaz)
00301     const throw() {
00302     real X, Y, Z, M[Geocentric::dim2_];
00303     _earth.Earth().IntForward(lat, lon, h, X, Y, Z, M);
00304     real Tres = InternalT(X, Y, Z, deltax, deltay, deltaz, true, true);
00305     Geocentric::Unrotate(M, deltax, deltay, deltaz, deltax, deltay, deltaz);
00306     return Tres;
00307   }
00308 
00309   GravityCircle GravityModel::Circle(real lat, real h, unsigned caps) const {
00310     if (h != 0)
00311       // Disallow invoking GeoidHeight unless h is zero.
00312       caps &= ~(CAP_GAMMA0 | CAP_C);
00313     real X, Y, Z, M[Geocentric::dim2_];
00314     _earth.Earth().IntForward(lat, 0, h, X, Y, Z, M);
00315     // Y = 0, cphi = M[7], sphi = M[8];
00316     real
00317       invR = 1 / Math::hypot(X, Z),
00318       gamma0 = (caps & CAP_GAMMA0 ?_earth.SurfaceGravity(lat)
00319                 : Math::NaN<real>()),
00320       fx, fy, fz, gamma;
00321     if (caps & CAP_GAMMA) {
00322       _earth.U(X, Y, Z, fx, fy, fz); // fy = 0
00323       gamma = Math::hypot(fx, fz);
00324     } else
00325       gamma = Math::NaN<real>();
00326     _earth.Phi(X, Y, fx, fy);
00327     return GravityCircle(GravityCircle::mask(caps),
00328                          _earth._a, _earth._f, lat, h, Z, X, M[7], M[8],
00329                          _amodel, _GMmodel, _dzonal0, _corrmult,
00330                          gamma0, gamma, fx,
00331                          caps & CAP_G ?
00332                          _gravitational.Circle(X, Z, true) :
00333                          CircularEngine(),
00334                          // N.B. If CAP_DELTA is set then CAP_T should be too.
00335                          caps & CAP_T ?
00336                          _disturbing.Circle(-1, X, Z, (caps & CAP_DELTA) != 0) :
00337                          CircularEngine(),
00338                          caps & CAP_C ?
00339                          _correction.Circle(invR * X, invR * Z, false) :
00340                          CircularEngine());
00341   }
00342 
00343   std::string GravityModel::DefaultGravityPath() {
00344     string path;
00345     char* gravitypath = getenv("GRAVITY_PATH");
00346     if (gravitypath)
00347       path = string(gravitypath);
00348     if (path.length())
00349       return path;
00350     char* datapath = getenv("GEOGRAPHICLIB_DATA");
00351     if (datapath)
00352       path = string(datapath);
00353     return (path.length() ? path : string(GEOGRAPHICLIB_DATA)) + "/gravity";
00354   }
00355 
00356   std::string GravityModel::DefaultGravityName() {
00357     string name;
00358     char* gravityname = getenv("GRAVITY_NAME");
00359     if (gravityname)
00360       name = string(gravityname);
00361     return name.length() ? name : string(GRAVITY_DEFAULT_NAME);
00362   }
00363 
00364 } // namespace GeographicLib