00001
00002
00003
00004
00005
00006
00007
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 #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_MAGNETIC_DEFAULT_NAME)
00025 # define GEOGRAPHICLIB_MAGNETIC_DEFAULT_NAME "wmm2010"
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 MagneticModel::MagneticModel(const std::string& name,const std::string& path,
00038 const Geocentric& earth)
00039 : _name(name)
00040 , _dir(path)
00041 , _description("NONE")
00042 , _date("UNKNOWN")
00043 , _t0(Math::NaN())
00044 , _dt0(1)
00045 , _tmin(Math::NaN())
00046 , _tmax(Math::NaN())
00047 , _a(Math::NaN())
00048 , _hmin(Math::NaN())
00049 , _hmax(Math::NaN())
00050 , _Nmodels(1)
00051 , _norm(SphericalHarmonic::SCHMIDT)
00052 , _earth(earth)
00053 {
00054 if (_dir.empty())
00055 _dir = DefaultMagneticPath();
00056 ReadMetadata(_name);
00057 _G.resize(_Nmodels + 1);
00058 _H.resize(_Nmodels + 1);
00059 {
00060 string coeff = _filename + ".cof";
00061 ifstream coeffstr(coeff.c_str(), ios::binary);
00062 if (!coeffstr.good())
00063 throw GeographicErr("Error opening " + coeff);
00064 char id[idlength_ + 1];
00065 coeffstr.read(id, idlength_);
00066 if (!coeffstr.good())
00067 throw GeographicErr("No header in " + coeff);
00068 id[idlength_] = '\0';
00069 if (_id != string(id))
00070 throw GeographicErr("ID mismatch: " + _id + " vs " + id);
00071 for (int i = 0; i <= _Nmodels; ++i) {
00072 int N, M;
00073 SphericalEngine::coeff::readcoeffs(coeffstr, N, M, _G[i], _H[i]);
00074 if (!(M < 0 || _G[i][0] == 0))
00075 throw GeographicErr("A degree 0 term is not permitted");
00076 _harm.push_back(SphericalHarmonic(_G[i], _H[i], N, N, M, _a, _norm));
00077 }
00078 int pos = int(coeffstr.tellg());
00079 coeffstr.seekg(0, ios::end);
00080 if (pos != coeffstr.tellg())
00081 throw GeographicErr("Extra data in " + coeff);
00082 }
00083 }
00084
00085 void MagneticModel::ReadMetadata(const std::string& name) {
00086 const char* spaces = " \t\n\v\f\r";
00087 _filename = _dir + "/" + name + ".wmm";
00088 ifstream metastr(_filename.c_str());
00089 if (!metastr.good())
00090 throw GeographicErr("Cannot open " + _filename);
00091 string line;
00092 getline(metastr, line);
00093 if (!(line.size() >= 6 && line.substr(0,5) == "WMMF-"))
00094 throw GeographicErr(_filename + " does not contain WMMF-n signature");
00095 string::size_type n = line.find_first_of(spaces, 5);
00096 if (n != string::npos)
00097 n -= 5;
00098 string version = line.substr(5, n);
00099 if (version != "1")
00100 throw GeographicErr("Unknown version in " + _filename + ": " + version);
00101 string key, val;
00102 while (getline(metastr, line)) {
00103 if (!Utility::ParseLine(line, key, val))
00104 continue;
00105
00106 if (key == "Name")
00107 _name = val;
00108 else if (key == "Description")
00109 _description = val;
00110 else if (key == "ReleaseDate")
00111 _date = val;
00112 else if (key == "Radius")
00113 _a = Utility::num<real>(val);
00114 else if (key == "Type") {
00115 if (!(val == "Linear" || val == "linear"))
00116 throw GeographicErr("Only linear models are supported");
00117 } else if (key == "Epoch")
00118 _t0 = Utility::num<real>(val);
00119 else if (key == "DeltaEpoch")
00120 _dt0 = Utility::num<real>(val);
00121 else if (key == "NumModels")
00122 _Nmodels = Utility::num<int>(val);
00123 else if (key == "MinTime")
00124 _tmin = Utility::num<real>(val);
00125 else if (key == "MaxTime")
00126 _tmax = Utility::num<real>(val);
00127 else if (key == "MinHeight")
00128 _hmin = Utility::num<real>(val);
00129 else if (key == "MaxHeight")
00130 _hmax = Utility::num<real>(val);
00131 else if (key == "Normalization") {
00132 if (val == "FULL" || val == "Full" || val == "full")
00133 _norm = SphericalHarmonic::FULL;
00134 else if (val == "SCHMIDT" || val == "Schmidt" || val == "schmidt")
00135 _norm = SphericalHarmonic::SCHMIDT;
00136 else
00137 throw GeographicErr("Unknown normalization " + val);
00138 } else if (key == "ByteOrder") {
00139 if (val == "Big" || val == "big")
00140 throw GeographicErr("Only little-endian ordering is supported");
00141 else if (!(val == "Little" || val == "little"))
00142 throw GeographicErr("Unknown byte ordering " + val);
00143 } else if (key == "ID")
00144 _id = val;
00145
00146 }
00147
00148 if (!(Math::isfinite(_a) && _a > 0))
00149 throw GeographicErr("Reference radius must be positive");
00150 if (!(_t0 > 0))
00151 throw GeographicErr("Epoch time not defined");
00152 if (_tmin >= _tmax)
00153 throw GeographicErr("Min time exceeds max time");
00154 if (_hmin >= _hmax)
00155 throw GeographicErr("Min height exceeds max height");
00156 if (int(_id.size()) != idlength_)
00157 throw GeographicErr("Invalid ID");
00158 if (!(_dt0 > 0)) {
00159 if (_Nmodels > 1)
00160 throw GeographicErr("DeltaEpoch must be positive");
00161 else
00162 _dt0 = 1;
00163 }
00164 }
00165
00166 void MagneticModel::Field(real t, real lat, real lon, real h, bool diffp,
00167 real& Bx, real& By, real& Bz,
00168 real& Bxt, real& Byt, real& Bzt) const {
00169 t -= _t0;
00170 int n = max(min(int(floor(t / _dt0)), _Nmodels - 1), 0);
00171 bool interpolate = n + 1 < _Nmodels;
00172 t -= n * _dt0;
00173 real X, Y, Z;
00174 real M[Geocentric::dim2_];
00175 _earth.IntForward(lat, lon, h, X, Y, Z, M);
00176
00177
00178 real BX0 = 0, BY0 = 0, BZ0 = 0, BX1 = 0, BY1 = 0, BZ1 = 0;
00179 _harm[n](X, Y, Z, BX0, BY0, BZ0);
00180 _harm[n + 1](X, Y, Z, BX1, BY1, BZ1);
00181 if (interpolate) {
00182
00183 BX1 = (BX1 - BX0) / _dt0;
00184 BY1 = (BY1 - BY0) / _dt0;
00185 BZ1 = (BZ1 - BZ0) / _dt0;
00186 }
00187 BX0 += t * BX1;
00188 BY0 += t * BY1;
00189 BZ0 += t * BZ1;
00190 if (diffp) {
00191 Geocentric::Unrotate(M, BX1, BY1, BZ1, Bxt, Byt, Bzt);
00192 Bxt *= - _a;
00193 Byt *= - _a;
00194 Bzt *= - _a;
00195 }
00196 Geocentric::Unrotate(M, BX0, BY0, BZ0, Bx, By, Bz);
00197 Bx *= - _a;
00198 By *= - _a;
00199 Bz *= - _a;
00200 }
00201
00202 MagneticCircle MagneticModel::Circle(real t, real lat, real h) const {
00203 real t1 = t - _t0;
00204 int n = max(min(int(floor(t1 / _dt0)), _Nmodels - 1), 0);
00205 bool interpolate = n + 1 < _Nmodels;
00206 t1 -= n * _dt0;
00207 real X, Y, Z, M[Geocentric::dim2_];
00208 _earth.IntForward(lat, 0, h, X, Y, Z, M);
00209
00210
00211 return MagneticCircle(_a, _earth._f, lat, h, t,
00212 M[7], M[8], t1, _dt0, interpolate,
00213 _harm[n].Circle(X, Z, true),
00214 _harm[n + 1].Circle(X, Z, true));
00215 }
00216
00217 void MagneticModel::FieldComponents(real Bx, real By, real Bz,
00218 real Bxt, real Byt, real Bzt,
00219 real& H, real& F, real& D, real& I,
00220 real& Ht, real& Ft,
00221 real& Dt, real& It) {
00222 H = Math::hypot(Bx, By);
00223 Ht = H ? (Bx * Bxt + By * Byt) / H : Math::hypot(Bxt, Byt);
00224 D = (0 - (H ? atan2(-Bx, By) : atan2(-Bxt, Byt))) / Math::degree();
00225 Dt = (H ? (By * Bxt - Bx * Byt) / Math::sq(H) : 0) / Math::degree();
00226 F = Math::hypot(H, Bz);
00227 Ft = F ? (H * Ht + Bz * Bzt) / F : Math::hypot(Ht, Bzt);
00228 I = (F ? atan2(-Bz, H) : atan2(-Bzt, Ht)) / Math::degree();
00229 It = (F ? (Bz * Ht - H * Bzt) / Math::sq(F) : 0) / Math::degree();
00230 }
00231
00232 std::string MagneticModel::DefaultMagneticPath() {
00233 string path;
00234 char* magneticpath = getenv("GEOGRAPHICLIB_MAGNETIC_PATH");
00235 if (magneticpath)
00236 path = string(magneticpath);
00237 if (!path.empty())
00238 return path;
00239 char* datapath = getenv("GEOGRAPHICLIB_DATA");
00240 if (datapath)
00241 path = string(datapath);
00242 return (!path.empty() ? path : string(GEOGRAPHICLIB_DATA)) + "/magnetic";
00243 }
00244
00245 std::string MagneticModel::DefaultMagneticName() {
00246 string name;
00247 char* magneticname = getenv("GEOGRAPHICLIB_MAGNETIC_NAME");
00248 if (magneticname)
00249 name = string(magneticname);
00250 return !name.empty() ? name : string(GEOGRAPHICLIB_MAGNETIC_DEFAULT_NAME);
00251 }
00252
00253 }