GeographicLib  2.1.1
MagneticModel.hpp
Go to the documentation of this file.
1 /**
2  * \file MagneticModel.hpp
3  * \brief Header for GeographicLib::MagneticModel class
4  *
5  * Copyright (c) Charles Karney (2011-2022) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * https://geographiclib.sourceforge.io/
8  **********************************************************************/
9 
10 #if !defined(GEOGRAPHICLIB_MAGNETICMODEL_HPP)
11 #define GEOGRAPHICLIB_MAGNETICMODEL_HPP 1
12 
16 
17 #if defined(_MSC_VER)
18 // Squelch warnings about dll vs vector
19 # pragma warning (push)
20 # pragma warning (disable: 4251)
21 #endif
22 
23 namespace GeographicLib {
24 
25  class MagneticCircle;
26 
27  /**
28  * \brief Model of the earth's magnetic field
29  *
30  * Evaluate the earth's magnetic field according to a model. At present only
31  * internal magnetic fields are handled. These are due to the earth's code
32  * and crust; these vary slowly (over many years). Excluded are the effects
33  * of currents in the ionosphere and magnetosphere which have daily and
34  * annual variations.
35  *
36  * See \ref magnetic for details of how to install the magnetic models and
37  * the data format.
38  *
39  * See
40  * - General information:
41  * - http://geomag.org/models/index.html
42  * - WMM2010:
43  * - https://ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml
44  * - https://ngdc.noaa.gov/geomag/WMM/data/WMM2010/WMM2010COF.zip
45  * - WMM2015 (deprecated):
46  * - https://ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml
47  * - https://ngdc.noaa.gov/geomag/WMM/data/WMM2015/WMM2015COF.zip
48  * - WMM2015V2:
49  * - https://ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml
50  * - https://ngdc.noaa.gov/geomag/WMM/data/WMM2015/WMM2015v2COF.zip
51  * - WMM2020:
52  * - https://ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml
53  * - https://ngdc.noaa.gov/geomag/WMM/data/WMM2020/WMM2020COF.zip
54  * - IGRF11:
55  * - https://ngdc.noaa.gov/IAGA/vmod/igrf.html
56  * - https://ngdc.noaa.gov/IAGA/vmod/igrf11coeffs.txt
57  * - https://ngdc.noaa.gov/IAGA/vmod/geomag70_linux.tar.gz
58  * - EMM2010:
59  * - https://ngdc.noaa.gov/geomag/EMM/index.html
60  * - https://ngdc.noaa.gov/geomag/EMM/data/geomag/EMM2010_Sph_Windows_Linux.zip
61  * - EMM2015:
62  * - https://ngdc.noaa.gov/geomag/EMM/index.html
63  * - https://www.ngdc.noaa.gov/geomag/EMM/data/geomag/EMM2015_Sph_Linux.zip
64  * - EMM2017:
65  * - https://ngdc.noaa.gov/geomag/EMM/index.html
66  * - https://www.ngdc.noaa.gov/geomag/EMM/data/geomag/EMM2017_Sph_Linux.zip
67  *
68  * Example of use:
69  * \include example-MagneticModel.cpp
70  *
71  * <a href="MagneticField.1.html">MagneticField</a> is a command-line utility
72  * providing access to the functionality of MagneticModel and MagneticCircle.
73  **********************************************************************/
74 
76  private:
77  typedef Math::real real;
78  static const int idlength_ = 8;
79  std::string _name, _dir, _description, _date, _filename, _id;
80  real _t0, _dt0, _tmin, _tmax, _a, _hmin, _hmax;
81  int _nNmodels, _nNconstants, _nmx, _mmx;
83  Geocentric _earth;
84  std::vector< std::vector<real> > _gG;
85  std::vector< std::vector<real> > _hH;
86  std::vector<SphericalHarmonic> _harm;
87  void Field(real t, real lat, real lon, real h, bool diffp,
88  real& Bx, real& By, real& Bz,
89  real& Bxt, real& Byt, real& Bzt) const;
90  void ReadMetadata(const std::string& name);
91  // copy constructor not allowed
92  MagneticModel(const MagneticModel&) = delete;
93  // nor copy assignment
94  MagneticModel& operator=(const MagneticModel&) = delete;
95  public:
96 
97  /** \name Setting up the magnetic model
98  **********************************************************************/
99  ///@{
100  /**
101  * Construct a magnetic model.
102  *
103  * @param[in] name the name of the model.
104  * @param[in] path (optional) directory for data file.
105  * @param[in] earth (optional) Geocentric object for converting
106  * coordinates; default Geocentric::WGS84().
107  * @param[in] Nmax (optional) if non-negative, truncate the degree of the
108  * model this value.
109  * @param[in] Mmax (optional) if non-negative, truncate the order of the
110  * model this value.
111  * @exception GeographicErr if the data file cannot be found, is
112  * unreadable, or is corrupt, or if \e Mmax > \e Nmax.
113  * @exception std::bad_alloc if the memory necessary for storing the model
114  * can't be allocated.
115  *
116  * A filename is formed by appending ".wmm" (World Magnetic Model) to the
117  * name. If \e path is specified (and is non-empty), then the file is
118  * loaded from directory, \e path. Otherwise the path is given by the
119  * DefaultMagneticPath().
120  *
121  * This file contains the metadata which specifies the properties of the
122  * model. The coefficients for the spherical harmonic sums are obtained
123  * from a file obtained by appending ".cof" to metadata file (so the
124  * filename ends in ".wwm.cof").
125  *
126  * The model is not tied to a particular ellipsoidal model of the earth.
127  * The final earth argument to the constructor specifies an ellipsoid to
128  * allow geodetic coordinates to the transformed into the spherical
129  * coordinates used in the spherical harmonic sum.
130  *
131  * If \e Nmax &ge; 0 and \e Mmax < 0, then \e Mmax is set to \e Nmax.
132  * After the model is loaded, the maximum degree and order of the model can
133  * be found by the Degree() and Order() methods.
134  **********************************************************************/
135  explicit MagneticModel(const std::string& name,
136  const std::string& path = "",
137  const Geocentric& earth = Geocentric::WGS84(),
138  int Nmax = -1, int Mmax = -1);
139  ///@}
140 
141  /** \name Compute the magnetic field
142  **********************************************************************/
143  ///@{
144  /**
145  * Evaluate the components of the geomagnetic field.
146  *
147  * @param[in] t the time (fractional years).
148  * @param[in] lat latitude of the point (degrees).
149  * @param[in] lon longitude of the point (degrees).
150  * @param[in] h the height of the point above the ellipsoid (meters).
151  * @param[out] Bx the easterly component of the magnetic field (nanotesla).
152  * @param[out] By the northerly component of the magnetic field
153  * (nanotesla).
154  * @param[out] Bz the vertical (up) component of the magnetic field
155  * (nanotesla).
156  *
157  * Use Utility::fractionalyear to convert a date of the form yyyy-mm or
158  * yyyy-mm-dd into a fractional year.
159  **********************************************************************/
160  void operator()(real t, real lat, real lon, real h,
161  real& Bx, real& By, real& Bz) const {
162  real dummy;
163  Field(t, lat, lon, h, false, Bx, By, Bz, dummy, dummy, dummy);
164  }
165 
166  /**
167  * Evaluate the components of the geomagnetic field and their time
168  * derivatives
169  *
170  * @param[in] t the time (fractional years).
171  * @param[in] lat latitude of the point (degrees).
172  * @param[in] lon longitude of the point (degrees).
173  * @param[in] h the height of the point above the ellipsoid (meters).
174  * @param[out] Bx the easterly component of the magnetic field (nanotesla).
175  * @param[out] By the northerly component of the magnetic field
176  * (nanotesla).
177  * @param[out] Bz the vertical (up) component of the magnetic field
178  * (nanotesla).
179  * @param[out] Bxt the rate of change of \e Bx (nT/yr).
180  * @param[out] Byt the rate of change of \e By (nT/yr).
181  * @param[out] Bzt the rate of change of \e Bz (nT/yr).
182  *
183  * Use Utility::fractionalyear to convert a date of the form yyyy-mm or
184  * yyyy-mm-dd into a fractional year.
185  **********************************************************************/
186  void operator()(real t, real lat, real lon, real h,
187  real& Bx, real& By, real& Bz,
188  real& Bxt, real& Byt, real& Bzt) const {
189  Field(t, lat, lon, h, true, Bx, By, Bz, Bxt, Byt, Bzt);
190  }
191 
192  /**
193  * Create a MagneticCircle object to allow the geomagnetic field at many
194  * points with constant \e lat, \e h, and \e t and varying \e lon to be
195  * computed efficiently.
196  *
197  * @param[in] t the time (fractional years).
198  * @param[in] lat latitude of the point (degrees).
199  * @param[in] h the height of the point above the ellipsoid (meters).
200  * @exception std::bad_alloc if the memory necessary for creating a
201  * MagneticCircle can't be allocated.
202  * @return a MagneticCircle object whose MagneticCircle::operator()(real
203  * lon) member function computes the field at particular values of \e
204  * lon.
205  *
206  * If the field at several points on a circle of latitude need to be
207  * calculated then creating a MagneticCircle and using its member functions
208  * will be substantially faster, especially for high-degree models.
209  *
210  * Use Utility::fractionalyear to convert a date of the form yyyy-mm or
211  * yyyy-mm-dd into a fractional year.
212  **********************************************************************/
213  MagneticCircle Circle(real t, real lat, real h) const;
214 
215  /**
216  * Compute the magnetic field in geocentric coordinate.
217  *
218  * @param[in] t the time (fractional years).
219  * @param[in] X geocentric coordinate (meters).
220  * @param[in] Y geocentric coordinate (meters).
221  * @param[in] Z geocentric coordinate (meters).
222  * @param[out] BX the \e X component of the magnetic field (nT).
223  * @param[out] BY the \e Y component of the magnetic field (nT).
224  * @param[out] BZ the \e Z component of the magnetic field (nT).
225  * @param[out] BXt the rate of change of \e BX (nT/yr).
226  * @param[out] BYt the rate of change of \e BY (nT/yr).
227  * @param[out] BZt the rate of change of \e BZ (nT/yr).
228  *
229  * Use Utility::fractionalyear to convert a date of the form yyyy-mm or
230  * yyyy-mm-dd into a fractional year.
231  **********************************************************************/
232  void FieldGeocentric(real t, real X, real Y, real Z,
233  real& BX, real& BY, real& BZ,
234  real& BXt, real& BYt, real& BZt) const;
235 
236  /**
237  * Compute various quantities dependent on the magnetic field.
238  *
239  * @param[in] Bx the \e x (easterly) component of the magnetic field (nT).
240  * @param[in] By the \e y (northerly) component of the magnetic field (nT).
241  * @param[in] Bz the \e z (vertical, up positive) component of the magnetic
242  * field (nT).
243  * @param[out] H the horizontal magnetic field (nT).
244  * @param[out] F the total magnetic field (nT).
245  * @param[out] D the declination of the field (degrees east of north).
246  * @param[out] I the inclination of the field (degrees down from
247  * horizontal).
248  **********************************************************************/
249  static void FieldComponents(real Bx, real By, real Bz,
250  real& H, real& F, real& D, real& I) {
251  real Ht, Ft, Dt, It;
252  FieldComponents(Bx, By, Bz, real(0), real(1), real(0),
253  H, F, D, I, Ht, Ft, Dt, It);
254  }
255 
256  /**
257  * Compute various quantities dependent on the magnetic field and its rate
258  * of change.
259  *
260  * @param[in] Bx the \e x (easterly) component of the magnetic field (nT).
261  * @param[in] By the \e y (northerly) component of the magnetic field (nT).
262  * @param[in] Bz the \e z (vertical, up positive) component of the magnetic
263  * field (nT).
264  * @param[in] Bxt the rate of change of \e Bx (nT/yr).
265  * @param[in] Byt the rate of change of \e By (nT/yr).
266  * @param[in] Bzt the rate of change of \e Bz (nT/yr).
267  * @param[out] H the horizontal magnetic field (nT).
268  * @param[out] F the total magnetic field (nT).
269  * @param[out] D the declination of the field (degrees east of north).
270  * @param[out] I the inclination of the field (degrees down from
271  * horizontal).
272  * @param[out] Ht the rate of change of \e H (nT/yr).
273  * @param[out] Ft the rate of change of \e F (nT/yr).
274  * @param[out] Dt the rate of change of \e D (degrees/yr).
275  * @param[out] It the rate of change of \e I (degrees/yr).
276  **********************************************************************/
277  static void FieldComponents(real Bx, real By, real Bz,
278  real Bxt, real Byt, real Bzt,
279  real& H, real& F, real& D, real& I,
280  real& Ht, real& Ft, real& Dt, real& It);
281  ///@}
282 
283  /** \name Inspector functions
284  **********************************************************************/
285  ///@{
286  /**
287  * @return the description of the magnetic model, if available, from the
288  * Description file in the data file; if absent, return "NONE".
289  **********************************************************************/
290  const std::string& Description() const { return _description; }
291 
292  /**
293  * @return date of the model, if available, from the ReleaseDate field in
294  * the data file; if absent, return "UNKNOWN".
295  **********************************************************************/
296  const std::string& DateTime() const { return _date; }
297 
298  /**
299  * @return full file name used to load the magnetic model.
300  **********************************************************************/
301  const std::string& MagneticFile() const { return _filename; }
302 
303  /**
304  * @return "name" used to load the magnetic model (from the first argument
305  * of the constructor, but this may be overridden by the model file).
306  **********************************************************************/
307  const std::string& MagneticModelName() const { return _name; }
308 
309  /**
310  * @return directory used to load the magnetic model.
311  **********************************************************************/
312  const std::string& MagneticModelDirectory() const { return _dir; }
313 
314  /**
315  * @return the minimum height above the ellipsoid (in meters) for which
316  * this MagneticModel should be used.
317  *
318  * Because the model will typically provide useful results
319  * slightly outside the range of allowed heights, no check of \e t
320  * argument is made by MagneticModel::operator()() or
321  * MagneticModel::Circle.
322  **********************************************************************/
323  Math::real MinHeight() const { return _hmin; }
324 
325  /**
326  * @return the maximum height above the ellipsoid (in meters) for which
327  * this MagneticModel should be used.
328  *
329  * Because the model will typically provide useful results
330  * slightly outside the range of allowed heights, no check of \e t
331  * argument is made by MagneticModel::operator()() or
332  * MagneticModel::Circle.
333  **********************************************************************/
334  Math::real MaxHeight() const { return _hmax; }
335 
336  /**
337  * @return the minimum time (in years) for which this MagneticModel should
338  * be used.
339  *
340  * Because the model will typically provide useful results
341  * slightly outside the range of allowed times, no check of \e t
342  * argument is made by MagneticModel::operator()() or
343  * MagneticModel::Circle.
344  **********************************************************************/
345  Math::real MinTime() const { return _tmin; }
346 
347  /**
348  * @return the maximum time (in years) for which this MagneticModel should
349  * be used.
350  *
351  * Because the model will typically provide useful results
352  * slightly outside the range of allowed times, no check of \e t
353  * argument is made by MagneticModel::operator()() or
354  * MagneticModel::Circle.
355  **********************************************************************/
356  Math::real MaxTime() const { return _tmax; }
357 
358  /**
359  * @return \e a the equatorial radius of the ellipsoid (meters). This is
360  * the value of \e a inherited from the Geocentric object used in the
361  * constructor.
362  **********************************************************************/
363  Math::real EquatorialRadius() const { return _earth.EquatorialRadius(); }
364 
365  /**
366  * @return \e f the flattening of the ellipsoid. This is the value
367  * inherited from the Geocentric object used in the constructor.
368  **********************************************************************/
369  Math::real Flattening() const { return _earth.Flattening(); }
370 
371  /**
372  * @return \e Nmax the maximum degree of the components of the model.
373  **********************************************************************/
374  int Degree() const { return _nmx; }
375 
376  /**
377  * @return \e Mmax the maximum order of the components of the model.
378  **********************************************************************/
379  int Order() const { return _mmx; }
380  ///@}
381 
382  /**
383  * @return the default path for magnetic model data files.
384  *
385  * This is the value of the environment variable
386  * GEOGRAPHICLIB_MAGNETIC_PATH, if set; otherwise, it is
387  * $GEOGRAPHICLIB_DATA/magnetic if the environment variable
388  * GEOGRAPHICLIB_DATA is set; otherwise, it is a compile-time default
389  * (/usr/local/share/GeographicLib/magnetic on non-Windows systems and
390  * C:/ProgramData/GeographicLib/magnetic on Windows systems).
391  **********************************************************************/
392  static std::string DefaultMagneticPath();
393 
394  /**
395  * @return the default name for the magnetic model.
396  *
397  * This is the value of the environment variable
398  * GEOGRAPHICLIB_MAGNETIC_NAME, if set; otherwise, it is "wmm2020". The
399  * MagneticModel class does not use this function; it is just provided as a
400  * convenience for a calling program when constructing a MagneticModel
401  * object.
402  **********************************************************************/
403  static std::string DefaultMagneticName();
404  };
405 
406 } // namespace GeographicLib
407 
408 #if defined(_MSC_VER)
409 # pragma warning (pop)
410 #endif
411 
412 #endif // GEOGRAPHICLIB_MAGNETICMODEL_HPP
Header for GeographicLib::Constants class.
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:67
Header for GeographicLib::Geocentric class.
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
Header for GeographicLib::SphericalHarmonic class.
Geocentric coordinates
Definition: Geocentric.hpp:67
Math::real Flattening() const
Definition: Geocentric.hpp:255
Math::real EquatorialRadius() const
Definition: Geocentric.hpp:248
static const Geocentric & WGS84()
Definition: Geocentric.cpp:31
Geomagnetic field on a circle of latitude.
Model of the earth's magnetic field.
const std::string & DateTime() const
const std::string & Description() const
Math::real EquatorialRadius() const
void operator()(real t, real lat, real lon, real h, real &Bx, real &By, real &Bz) const
const std::string & MagneticModelDirectory() const
Math::real MaxHeight() const
const std::string & MagneticFile() const
Math::real Flattening() const
Math::real MinHeight() const
const std::string & MagneticModelName() const
void operator()(real t, real lat, real lon, real h, real &Bx, real &By, real &Bz, real &Bxt, real &Byt, real &Bzt) const
static void FieldComponents(real Bx, real By, real Bz, real &H, real &F, real &D, real &I)
Namespace for GeographicLib.
Definition: Accumulator.cpp:12