GeographicLib
1.21
|
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