GeographicLib
1.21
|
00001 /** 00002 * \file MagneticModel.cpp 00003 * \brief Implementation for GeographicLib::MagneticModel class 00004 * 00005 * Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed under 00006 * the MIT/X11 License. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #include <GeographicLib/MagneticModel.hpp> 00011 #include <fstream> 00012 #include <GeographicLib/SphericalEngine.hpp> 00013 #include <GeographicLib/MagneticCircle.hpp> 00014 #include <GeographicLib/Utility.hpp> 00015 00016 #define GEOGRAPHICLIB_MAGNETICMODEL_CPP \ 00017 "$Id: b0287ac014f10e4c6656b67f21c764432a47559a $" 00018 00019 RCSID_DECL(GEOGRAPHICLIB_MAGNETICMODEL_CPP) 00020 RCSID_DECL(GEOGRAPHICLIB_MAGNETICMODEL_HPP) 00021 00022 #if !defined(GEOGRAPHICLIB_DATA) 00023 # if defined(_MSC_VER) 00024 # define GEOGRAPHICLIB_DATA \ 00025 "C:/Documents and Settings/All Users/Application Data/GeographicLib" 00026 # else 00027 # define GEOGRAPHICLIB_DATA "/usr/local/share/GeographicLib" 00028 # endif 00029 #endif 00030 00031 #if !defined(MAGNETIC_DEFAULT_NAME) 00032 # define MAGNETIC_DEFAULT_NAME "wmm2010" 00033 #endif 00034 00035 #if defined(_MSC_VER) 00036 // Squelch warnings about unsafe use of getenv 00037 #pragma warning (disable: 4996) 00038 #endif 00039 00040 namespace GeographicLib { 00041 00042 using namespace std; 00043 00044 MagneticModel::MagneticModel(const std::string& name,const std::string& path, 00045 const Geocentric& earth) 00046 : _name(name) 00047 , _dir(path) 00048 , _description("NONE") 00049 , _date("UNKNOWN") 00050 , _t0(Math::NaN<real>()) 00051 , _dt0(1) 00052 , _tmin(Math::NaN<real>()) 00053 , _tmax(Math::NaN<real>()) 00054 , _a(Math::NaN<real>()) 00055 , _hmin(Math::NaN<real>()) 00056 , _hmax(Math::NaN<real>()) 00057 , _Nmodels(1) 00058 , _norm(SphericalHarmonic::SCHMIDT) 00059 , _earth(earth) 00060 { 00061 if (_dir.empty()) 00062 _dir = DefaultMagneticPath(); 00063 ReadMetadata(_name); 00064 _G.resize(_Nmodels + 1); 00065 _H.resize(_Nmodels + 1); 00066 { 00067 string coeff = _filename + ".cof"; 00068 ifstream coeffstr(coeff.c_str(), ios::binary); 00069 if (!coeffstr.good()) 00070 throw GeographicErr("Error opening " + coeff); 00071 char id[idlength_ + 1]; 00072 coeffstr.read(id, idlength_); 00073 if (!coeffstr.good()) 00074 throw GeographicErr("No header in " + coeff); 00075 id[idlength_] = '\0'; 00076 if (_id != string(id)) 00077 throw GeographicErr("ID mismatch: " + _id + " vs " + id); 00078 for (int i = 0; i <= _Nmodels; ++i) { 00079 int N, M; 00080 SphericalEngine::coeff::readcoeffs(coeffstr, N, M, _G[i], _H[i]); 00081 if (!(M < 0 || _G[i][0] == 0)) 00082 throw GeographicErr("A degree 0 term is not permitted"); 00083 _harm.push_back(SphericalHarmonic(_G[i], _H[i], N, N, M, _a, _norm)); 00084 } 00085 int pos = int(coeffstr.tellg()); 00086 coeffstr.seekg(0, ios::end); 00087 if (pos != coeffstr.tellg()) 00088 throw GeographicErr("Extra data in " + coeff); 00089 } 00090 } 00091 00092 void MagneticModel::ReadMetadata(const std::string& name) { 00093 const char* spaces = " \t\n\v\f\r"; 00094 _filename = _dir + "/" + name + ".wmm"; 00095 ifstream metastr(_filename.c_str()); 00096 if (!metastr.good()) 00097 throw GeographicErr("Cannot open " + _filename); 00098 string line; 00099 getline(metastr, line); 00100 if (!(line.size() >= 6 && line.substr(0,5) == "WMMF-")) 00101 throw GeographicErr(_filename + " does not contain WMMF-n signature"); 00102 string::size_type n = line.find_first_of(spaces, 5); 00103 if (n != string::npos) 00104 n -= 5; 00105 string version = line.substr(5, n); 00106 if (version != "1") 00107 throw GeographicErr("Unknown version in " + _filename + ": " + version); 00108 string key, val; 00109 while (getline(metastr, line)) { 00110 if (!Utility::ParseLine(line, key, val)) 00111 continue; 00112 // Process key words 00113 if (key == "Name") 00114 _name = val; 00115 else if (key == "Description") 00116 _description = val; 00117 else if (key == "ReleaseDate") 00118 _date = val; 00119 else if (key == "Radius") 00120 _a = Utility::num<real>(val); 00121 else if (key == "Type") { 00122 if (!(val == "Linear" || val == "linear")) 00123 throw GeographicErr("Only linear models are supported"); 00124 } else if (key == "Epoch") 00125 _t0 = Utility::num<real>(val); 00126 else if (key == "DeltaEpoch") 00127 _dt0 = Utility::num<real>(val); 00128 else if (key == "NumModels") 00129 _Nmodels = Utility::num<int>(val); 00130 else if (key == "MinTime") 00131 _tmin = Utility::num<real>(val); 00132 else if (key == "MaxTime") 00133 _tmax = Utility::num<real>(val); 00134 else if (key == "MinHeight") 00135 _hmin = Utility::num<real>(val); 00136 else if (key == "MaxHeight") 00137 _hmax = Utility::num<real>(val); 00138 else if (key == "Normalization") { 00139 if (val == "FULL" || val == "Full" || val == "full") 00140 _norm = SphericalHarmonic::FULL; 00141 else if (val == "SCHMIDT" || val == "Schmidt" || val == "schmidt") 00142 _norm = SphericalHarmonic::SCHMIDT; 00143 else 00144 throw GeographicErr("Unknown normalization " + val); 00145 } else if (key == "ByteOrder") { 00146 if (val == "Big" || val == "big") 00147 throw GeographicErr("Only little-endian ordering is supported"); 00148 else if (!(val == "Little" || val == "little")) 00149 throw GeographicErr("Unknown byte ordering " + val); 00150 } else if (key == "ID") 00151 _id = val; 00152 // else unrecognized keywords are skipped 00153 } 00154 // Check values 00155 if (!(Math::isfinite(_a) && _a > 0)) 00156 throw GeographicErr("Reference radius must be positive"); 00157 if (!(_t0 > 0)) 00158 throw GeographicErr("Epoch time not defined"); 00159 if (_tmin >= _tmax) 00160 throw GeographicErr("Min time exceeds max time"); 00161 if (_hmin >= _hmax) 00162 throw GeographicErr("Min height exceeds max height"); 00163 if (int(_id.size()) != idlength_) 00164 throw GeographicErr("Invalid ID"); 00165 if (!(_dt0 > 0)) { 00166 if (_Nmodels > 1) 00167 throw GeographicErr("DeltaEpoch must be positive"); 00168 else 00169 _dt0 = 1; 00170 } 00171 } 00172 00173 void MagneticModel::Field(real t, real lat, real lon, real h, bool diffp, 00174 real& Bx, real& By, real& Bz, 00175 real& Bxt, real& Byt, real& Bzt) const throw() { 00176 t -= _t0; 00177 int n = max(min(int(floor(t / _dt0)), _Nmodels - 1), 0); 00178 bool interpolate = n + 1 < _Nmodels; 00179 t -= n * _dt0; 00180 real X, Y, Z; 00181 real M[Geocentric::dim2_]; 00182 _earth.IntForward(lat, lon, h, X, Y, Z, M); 00183 real BX0, BY0, BZ0, BX1, BY1, BZ1; // Components in geocentric basis 00184 _harm[n](X, Y, Z, BX0, BY0, BZ0); 00185 _harm[n + 1](X, Y, Z, BX1, BY1, BZ1); 00186 if (interpolate) { 00187 // Convert to a time derivative 00188 BX1 = (BX1 - BX0) / _dt0; 00189 BY1 = (BY1 - BY0) / _dt0; 00190 BZ1 = (BZ1 - BZ0) / _dt0; 00191 } 00192 BX0 += t * BX1; 00193 BY0 += t * BY1; 00194 BZ0 += t * BZ1; 00195 if (diffp) { 00196 Geocentric::Unrotate(M, BX1, BY1, BZ1, Bxt, Byt, Bzt); 00197 Bxt *= - _a; 00198 Byt *= - _a; 00199 Bzt *= - _a; 00200 } 00201 Geocentric::Unrotate(M, BX0, BY0, BZ0, Bx, By, Bz); 00202 Bx *= - _a; 00203 By *= - _a; 00204 Bz *= - _a; 00205 } 00206 00207 MagneticCircle MagneticModel::Circle(real t, real lat, real h) const { 00208 real t1 = t - _t0; 00209 int n = max(min(int(floor(t1 / _dt0)), _Nmodels - 1), 0); 00210 bool interpolate = n + 1 < _Nmodels; 00211 t1 -= n * _dt0; 00212 real X, Y, Z, M[Geocentric::dim2_]; 00213 _earth.IntForward(lat, 0, h, X, Y, Z, M); 00214 // Y = 0, cphi = M[7], sphi = M[8]; 00215 00216 return MagneticCircle(_a, _earth._f, lat, h, t, 00217 M[7], M[8], t1, _dt0, interpolate, 00218 _harm[n].Circle(X, Z, true), 00219 _harm[n + 1].Circle(X, Z, true)); 00220 } 00221 00222 void MagneticModel::FieldComponents(real Bx, real By, real Bz, 00223 real Bxt, real Byt, real Bzt, 00224 real& H, real& F, real& D, real& I, 00225 real& Ht, real& Ft, 00226 real& Dt, real& It) throw() { 00227 H = Math::hypot(Bx, By); 00228 Ht = H ? (Bx * Bxt + By * Byt) / H : Math::hypot(Bxt, Byt); 00229 D = (0 - (H ? atan2(-Bx, By) : atan2(-Bxt, Byt))) / Math::degree<real>(); 00230 Dt = (H ? (By * Bxt - Bx * Byt) / Math::sq(H) : 0) / Math::degree<real>(); 00231 F = Math::hypot(H, Bz); 00232 Ft = F ? (H * Ht + Bz * Bzt) / F : Math::hypot(Ht, Bzt); 00233 I = (F ? atan2(-Bz, H) : atan2(-Bzt, Ht)) / Math::degree<real>(); 00234 It = (F ? (Bz * Ht - H * Bzt) / Math::sq(F) : 0) / Math::degree<real>(); 00235 } 00236 00237 std::string MagneticModel::DefaultMagneticPath() { 00238 string path; 00239 char* magneticpath = getenv("MAGNETIC_PATH"); 00240 if (magneticpath) 00241 path = string(magneticpath); 00242 if (path.length()) 00243 return path; 00244 char* datapath = getenv("GEOGRAPHICLIB_DATA"); 00245 if (datapath) 00246 path = string(datapath); 00247 return (path.length() ? path : string(GEOGRAPHICLIB_DATA)) + "/magnetic"; 00248 } 00249 00250 std::string MagneticModel::DefaultMagneticName() { 00251 string name; 00252 char* magneticname = getenv("MAGNETIC_NAME"); 00253 if (magneticname) 00254 name = string(magneticname); 00255 return name.length() ? name : string(MAGNETIC_DEFAULT_NAME); 00256 } 00257 00258 } // namespace GeographicLib