GeographicLib  1.21
MagneticModel.cpp
Go to the documentation of this file.
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