00001 /** 00002 * \file MagneticModel.hpp 00003 * \brief Header 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 #if !defined(GEOGRAPHICLIB_MAGNETICMODEL_HPP) 00011 #define GEOGRAPHICLIB_MAGNETICMODEL_HPP 1 00012 00013 #include <GeographicLib/Constants.hpp> 00014 #include <GeographicLib/Geocentric.hpp> 00015 #include <GeographicLib/SphericalHarmonic.hpp> 00016 00017 #if defined(_MSC_VER) 00018 // Squelch warnings about dll vs vector 00019 # pragma warning (push) 00020 # pragma warning (disable: 4251) 00021 #endif 00022 00023 namespace GeographicLib { 00024 00025 class MagneticCircle; 00026 00027 /** 00028 * \brief Model of the earth's magnetic field 00029 * 00030 * Evaluate the earth's magnetic field according to a model. At present only 00031 * internal magnetic fields are handled. These are due to the earth's code 00032 * and crust; these vary slowly (over many years). Excluded are the effects 00033 * of currents in the ionosphere and magnetosphere which have daily and 00034 * annual variations. 00035 * 00036 * See \ref magnetic for details of how to install the magnetic model and the 00037 * data format. 00038 * 00039 * See 00040 * - General information: 00041 * - http://geomag.org/models/index.html 00042 * - WMM2010: 00043 * - http://ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml 00044 * - http://ngdc.noaa.gov/geomag/WMM/data/WMM2010/WMM2010COF.zip 00045 * - IGRF11: 00046 * - http://ngdc.noaa.gov/IAGA/vmod/igrf.html 00047 * - http://ngdc.noaa.gov/IAGA/vmod/igrf11coeffs.txt 00048 * - http://ngdc.noaa.gov/IAGA/vmod/geomag70_linux.tar.gz 00049 * - EMM2010: 00050 * - http://ngdc.noaa.gov/geomag/EMM/index.html 00051 * - http://ngdc.noaa.gov/geomag/EMM/data/geomag/EMM2010_Sph_Windows_Linux.zip 00052 * 00053 * Example of use: 00054 * \include example-MagneticModel.cpp 00055 * 00056 * <a href="MagneticField.1.html">MagneticField</a> is a command-line utility 00057 * providing access to the functionality of MagneticModel and MagneticCircle. 00058 **********************************************************************/ 00059 00060 class GEOGRAPHICLIB_EXPORT MagneticModel { 00061 private: 00062 typedef Math::real real; 00063 static const int idlength_ = 8; 00064 std::string _name, _dir, _description, _date, _filename, _id; 00065 real _t0, _dt0, _tmin, _tmax, _a, _hmin, _hmax; 00066 int _Nmodels; 00067 SphericalHarmonic::normalization _norm; 00068 Geocentric _earth; 00069 std::vector< std::vector<real> > _G; 00070 std::vector< std::vector<real> > _H; 00071 std::vector<SphericalHarmonic> _harm; 00072 void Field(real t, real lat, real lon, real h, bool diffp, 00073 real& Bx, real& By, real& Bz, 00074 real& Bxt, real& Byt, real& Bzt) const; 00075 void ReadMetadata(const std::string& name); 00076 MagneticModel(const MagneticModel&); // copy constructor not allowed 00077 MagneticModel& operator=(const MagneticModel&); // nor copy assignment 00078 public: 00079 00080 /** \name Setting up the magnetic model 00081 **********************************************************************/ 00082 ///@{ 00083 /** 00084 * Construct a magnetic model. 00085 * 00086 * @param[in] name the name of the model. 00087 * @param[in] path (optional) directory for data file. 00088 * @param[in] earth (optional) Geocentric object for converting 00089 * coordinates; default Geocentric::WGS84(). 00090 * @exception GeographicErr if the data file cannot be found, is 00091 * unreadable, or is corrupt. 00092 * @exception std::bad_alloc if the memory necessary for storing the model 00093 * can't be allocated. 00094 * 00095 * A filename is formed by appending ".wmm" (World Magnetic Model) to the 00096 * name. If \e path is specified (and is non-empty), then the file is 00097 * loaded from directory, \e path. Otherwise the path is given by the 00098 * DefaultMagneticPath(). 00099 * 00100 * This file contains the metadata which specifies the properties of the 00101 * model. The coefficients for the spherical harmonic sums are obtained 00102 * from a file obtained by appending ".cof" to metadata file (so the 00103 * filename ends in ".wwm.cof"). 00104 * 00105 * The model is not tied to a particular ellipsoidal model of the earth. 00106 * The final earth argument to the constructor specifies an ellipsoid to 00107 * allow geodetic coordinates to the transformed into the spherical 00108 * coordinates used in the spherical harmonic sum. 00109 **********************************************************************/ 00110 explicit MagneticModel(const std::string& name, 00111 const std::string& path = "", 00112 const Geocentric& earth = Geocentric::WGS84()); 00113 ///@} 00114 00115 /** \name Compute the magnetic field 00116 **********************************************************************/ 00117 ///@{ 00118 /** 00119 * Evaluate the components of the geomagnetic field. 00120 * 00121 * @param[in] t the time (years). 00122 * @param[in] lat latitude of the point (degrees). 00123 * @param[in] lon longitude of the point (degrees). 00124 * @param[in] h the height of the point above the ellipsoid (meters). 00125 * @param[out] Bx the easterly component of the magnetic field (nanotesla). 00126 * @param[out] By the northerly component of the magnetic field (nanotesla). 00127 * @param[out] Bz the vertical (up) component of the magnetic field 00128 * (nanotesla). 00129 **********************************************************************/ 00130 void operator()(real t, real lat, real lon, real h, 00131 real& Bx, real& By, real& Bz) const { 00132 real dummy; 00133 Field(t, lat, lon, h, false, Bx, By, Bz, dummy, dummy, dummy); 00134 } 00135 00136 /** 00137 * Evaluate the components of the geomagnetic field and their time 00138 * derivatives 00139 * 00140 * @param[in] t the time (years). 00141 * @param[in] lat latitude of the point (degrees). 00142 * @param[in] lon longitude of the point (degrees). 00143 * @param[in] h the height of the point above the ellipsoid (meters). 00144 * @param[out] Bx the easterly component of the magnetic field (nanotesla). 00145 * @param[out] By the northerly component of the magnetic field (nanotesla). 00146 * @param[out] Bz the vertical (up) component of the magnetic field 00147 * (nanotesla). 00148 * @param[out] Bxt the rate of change of \e Bx (nT/yr). 00149 * @param[out] Byt the rate of change of \e By (nT/yr). 00150 * @param[out] Bzt the rate of change of \e Bz (nT/yr). 00151 **********************************************************************/ 00152 void operator()(real t, real lat, real lon, real h, 00153 real& Bx, real& By, real& Bz, 00154 real& Bxt, real& Byt, real& Bzt) const { 00155 Field(t, lat, lon, h, true, Bx, By, Bz, Bxt, Byt, Bzt); 00156 } 00157 00158 /** 00159 * Create a MagneticCircle object to allow the geomagnetic field at many 00160 * points with constant \e lat, \e h, and \e t and varying \e lon to be 00161 * computed efficiently. 00162 * 00163 * @param[in] t the time (years). 00164 * @param[in] lat latitude of the point (degrees). 00165 * @param[in] h the height of the point above the ellipsoid (meters). 00166 * @exception std::bad_alloc if the memory necessary for creating a 00167 * MagneticCircle can't be allocated. 00168 * @return a MagneticCircle object whose MagneticCircle::operator()(real 00169 * lon) member function computes the field at particular values of \e 00170 * lon. 00171 * 00172 * If the field at several points on a circle of latitude need to be 00173 * calculated then creating a MagneticCircle and using its member functions 00174 * will be substantially faster, especially for high-degree models. 00175 **********************************************************************/ 00176 MagneticCircle Circle(real t, real lat, real h) const; 00177 00178 /** 00179 * Compute various quantities dependent on the magnetic field. 00180 * 00181 * @param[in] Bx the \e x (easterly) component of the magnetic field (nT). 00182 * @param[in] By the \e y (northerly) component of the magnetic field (nT). 00183 * @param[in] Bz the \e z (vertical, up positive) component of the magnetic 00184 * field (nT). 00185 * @param[out] H the horizontal magnetic field (nT). 00186 * @param[out] F the total magnetic field (nT). 00187 * @param[out] D the declination of the field (degrees east of north). 00188 * @param[out] I the inclination of the field (degrees down from 00189 * horizontal). 00190 **********************************************************************/ 00191 static void FieldComponents(real Bx, real By, real Bz, 00192 real& H, real& F, real& D, real& I) { 00193 real Ht, Ft, Dt, It; 00194 FieldComponents(Bx, By, Bz, real(0), real(1), real(0), 00195 H, F, D, I, Ht, Ft, Dt, It); 00196 } 00197 00198 /** 00199 * Compute various quantities dependent on the magnetic field and its rate 00200 * of change. 00201 * 00202 * @param[in] Bx the \e x (easterly) component of the magnetic field (nT). 00203 * @param[in] By the \e y (northerly) component of the magnetic field (nT). 00204 * @param[in] Bz the \e z (vertical, up positive) component of the magnetic 00205 * field (nT). 00206 * @param[in] Bxt the rate of change of \e Bx (nT/yr). 00207 * @param[in] Byt the rate of change of \e By (nT/yr). 00208 * @param[in] Bzt the rate of change of \e Bz (nT/yr). 00209 * @param[out] H the horizontal magnetic field (nT). 00210 * @param[out] F the total magnetic field (nT). 00211 * @param[out] D the declination of the field (degrees east of north). 00212 * @param[out] I the inclination of the field (degrees down from 00213 * horizontal). 00214 * @param[out] Ht the rate of change of \e H (nT/yr). 00215 * @param[out] Ft the rate of change of \e F (nT/yr). 00216 * @param[out] Dt the rate of change of \e D (degrees/yr). 00217 * @param[out] It the rate of change of \e I (degrees/yr). 00218 **********************************************************************/ 00219 static void FieldComponents(real Bx, real By, real Bz, 00220 real Bxt, real Byt, real Bzt, 00221 real& H, real& F, real& D, real& I, 00222 real& Ht, real& Ft, real& Dt, real& It); 00223 ///@} 00224 00225 /** \name Inspector functions 00226 **********************************************************************/ 00227 ///@{ 00228 /** 00229 * @return the description of the magnetic model, if available, from the 00230 * Description file in the data file; if absent, return "NONE". 00231 **********************************************************************/ 00232 const std::string& Description() const { return _description; } 00233 00234 /** 00235 * @return date of the model, if available, from the ReleaseDate field in 00236 * the data file; if absent, return "UNKNOWN". 00237 **********************************************************************/ 00238 const std::string& DateTime() const { return _date; } 00239 00240 /** 00241 * @return full file name used to load the magnetic model. 00242 **********************************************************************/ 00243 const std::string& MagneticFile() const { return _filename; } 00244 00245 /** 00246 * @return "name" used to load the magnetic model (from the first argument 00247 * of the constructor, but this may be overridden by the model file). 00248 **********************************************************************/ 00249 const std::string& MagneticModelName() const { return _name; } 00250 00251 /** 00252 * @return directory used to load the magnetic model. 00253 **********************************************************************/ 00254 const std::string& MagneticModelDirectory() const { return _dir; } 00255 00256 /** 00257 * @return the minimum height above the ellipsoid (in meters) for which 00258 * this MagneticModel should be used. 00259 * 00260 * Because the model will typically provide useful results 00261 * slightly outside the range of allowed heights, no check of \e t 00262 * argument is made by MagneticModel::operator()() or 00263 * MagneticModel::Circle. 00264 **********************************************************************/ 00265 Math::real MinHeight() const { return _hmin; } 00266 00267 /** 00268 * @return the maximum height above the ellipsoid (in meters) for which 00269 * this MagneticModel should be used. 00270 * 00271 * Because the model will typically provide useful results 00272 * slightly outside the range of allowed heights, no check of \e t 00273 * argument is made by MagneticModel::operator()() or 00274 * MagneticModel::Circle. 00275 **********************************************************************/ 00276 Math::real MaxHeight() const { return _hmax; } 00277 00278 /** 00279 * @return the minimum time (in years) for which this MagneticModel should 00280 * be used. 00281 * 00282 * Because the model will typically provide useful results 00283 * slightly outside the range of allowed times, no check of \e t 00284 * argument is made by MagneticModel::operator()() or 00285 * MagneticModel::Circle. 00286 **********************************************************************/ 00287 Math::real MinTime() const { return _tmin; } 00288 00289 /** 00290 * @return the maximum time (in years) for which this MagneticModel should 00291 * be used. 00292 * 00293 * Because the model will typically provide useful results 00294 * slightly outside the range of allowed times, no check of \e t 00295 * argument is made by MagneticModel::operator()() or 00296 * MagneticModel::Circle. 00297 **********************************************************************/ 00298 Math::real MaxTime() const { return _tmax; } 00299 00300 /** 00301 * @return \e a the equatorial radius of the ellipsoid (meters). This is 00302 * the value of \e a inherited from the Geocentric object used in the 00303 * constructor. 00304 **********************************************************************/ 00305 Math::real MajorRadius() const { return _earth.MajorRadius(); } 00306 00307 /** 00308 * @return \e f the flattening of the ellipsoid. This is the value 00309 * inherited from the Geocentric object used in the constructor. 00310 **********************************************************************/ 00311 Math::real Flattening() const { return _earth.Flattening(); } 00312 ///@} 00313 00314 /** 00315 * @return the default path for magnetic model data files. 00316 * 00317 * This is the value of the environment variable 00318 * GEOGRAPHICLIB_MAGNETIC_PATH, if set; otherwise, it is 00319 * $GEOGRAPHICLIB_DATA/magnetic if the environment variable 00320 * GEOGRAPHICLIB_DATA is set; otherwise, it is a compile-time default 00321 * (/usr/local/share/GeographicLib/magnetic on non-Windows systems and 00322 * C:/ProgramData/GeographicLib/magnetic on Windows systems). 00323 **********************************************************************/ 00324 static std::string DefaultMagneticPath(); 00325 00326 /** 00327 * @return the default name for the magnetic model. 00328 * 00329 * This is the value of the environment variable 00330 * GEOGRAPHICLIB_MAGNETIC_NAME, if set; otherwise, it is "wmm2010". The 00331 * MagneticModel class does not use this function; it is just provided as a 00332 * convenience for a calling program when constructing a MagneticModel 00333 * object. 00334 **********************************************************************/ 00335 static std::string DefaultMagneticName(); 00336 }; 00337 00338 } // namespace GeographicLib 00339 00340 #if defined(_MSC_VER) 00341 # pragma warning (pop) 00342 #endif 00343 00344 #endif // GEOGRAPHICLIB_MAGNETICMODEL_HPP