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