GeographicLib  2.1.2
MagneticModel.cpp
Go to the documentation of this file.
1 /**
2  * \file MagneticModel.cpp
3  * \brief Implementation for GeographicLib::MagneticModel class
4  *
5  * Copyright (c) Charles Karney (2011-2021) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * https://geographiclib.sourceforge.io/
8  **********************************************************************/
9 
11 #include <fstream>
15 
16 #if !defined(GEOGRAPHICLIB_DATA)
17 # if defined(_WIN32)
18 # define GEOGRAPHICLIB_DATA "C:/ProgramData/GeographicLib"
19 # else
20 # define GEOGRAPHICLIB_DATA "/usr/local/share/GeographicLib"
21 # endif
22 #endif
23 
24 #if !defined(GEOGRAPHICLIB_MAGNETIC_DEFAULT_NAME)
25 # define GEOGRAPHICLIB_MAGNETIC_DEFAULT_NAME "wmm2020"
26 #endif
27 
28 #if defined(_MSC_VER)
29 // Squelch warnings about unsafe use of getenv
30 # pragma warning (disable: 4996)
31 #endif
32 
33 namespace GeographicLib {
34 
35  using namespace std;
36 
37  MagneticModel::MagneticModel(const std::string& name, const std::string& path,
38  const Geocentric& earth, int Nmax, int Mmax)
39  : _name(name)
40  , _dir(path)
41  , _description("NONE")
42  , _date("UNKNOWN")
43  , _t0(Math::NaN())
44  , _dt0(1)
45  , _tmin(Math::NaN())
46  , _tmax(Math::NaN())
47  , _a(Math::NaN())
48  , _hmin(Math::NaN())
49  , _hmax(Math::NaN())
50  , _nNmodels(1)
51  , _nNconstants(0)
52  , _nmx(-1)
53  , _mmx(-1)
54  , _norm(SphericalHarmonic::SCHMIDT)
55  , _earth(earth)
56  {
57  if (_dir.empty())
58  _dir = DefaultMagneticPath();
59  bool truncate = Nmax >= 0 || Mmax >= 0;
60  if (truncate) {
61  if (Nmax >= 0 && Mmax < 0) Mmax = Nmax;
62  if (Nmax < 0) Nmax = numeric_limits<int>::max();
63  if (Mmax < 0) Mmax = numeric_limits<int>::max();
64  }
65  ReadMetadata(_name);
66  _gG.resize(_nNmodels + 1 + _nNconstants);
67  _hH.resize(_nNmodels + 1 + _nNconstants);
68  {
69  string coeff = _filename + ".cof";
70  ifstream coeffstr(coeff.c_str(), ios::binary);
71  if (!coeffstr.good())
72  throw GeographicErr("Error opening " + coeff);
73  char id[idlength_ + 1];
74  coeffstr.read(id, idlength_);
75  if (!coeffstr.good())
76  throw GeographicErr("No header in " + coeff);
77  id[idlength_] = '\0';
78  if (_id != string(id))
79  throw GeographicErr("ID mismatch: " + _id + " vs " + id);
80  for (int i = 0; i < _nNmodels + 1 + _nNconstants; ++i) {
81  int N, M;
82  if (truncate) { N = Nmax; M = Mmax; }
83  SphericalEngine::coeff::readcoeffs(coeffstr, N, M, _gG[i], _hH[i],
84  truncate);
85  if (!(M < 0 || _gG[i][0] == 0))
86  throw GeographicErr("A degree 0 term is not permitted");
87  _harm.push_back(SphericalHarmonic(_gG[i], _hH[i], N, N, M, _a, _norm));
88  _nmx = max(_nmx, _harm.back().Coefficients().nmx());
89  _mmx = max(_mmx, _harm.back().Coefficients().mmx());
90  }
91  int pos = int(coeffstr.tellg());
92  coeffstr.seekg(0, ios::end);
93  if (pos != coeffstr.tellg())
94  throw GeographicErr("Extra data in " + coeff);
95  }
96  }
97 
98  void MagneticModel::ReadMetadata(const string& name) {
99  const char* spaces = " \t\n\v\f\r";
100  _filename = _dir + "/" + name + ".wmm";
101  ifstream metastr(_filename.c_str());
102  if (!metastr.good())
103  throw GeographicErr("Cannot open " + _filename);
104  string line;
105  getline(metastr, line);
106  if (!(line.size() >= 6 && line.substr(0,5) == "WMMF-"))
107  throw GeographicErr(_filename + " does not contain WMMF-n signature");
108  string::size_type n = line.find_first_of(spaces, 5);
109  if (n != string::npos)
110  n -= 5;
111  string version(line, 5, n);
112  if (!(version == "1" || version == "2"))
113  throw GeographicErr("Unknown version in " + _filename + ": " + version);
114  string key, val;
115  while (getline(metastr, line)) {
116  if (!Utility::ParseLine(line, key, val))
117  continue;
118  // Process key words
119  if (key == "Name")
120  _name = val;
121  else if (key == "Description")
122  _description = val;
123  else if (key == "ReleaseDate")
124  _date = val;
125  else if (key == "Radius")
126  _a = Utility::val<real>(val);
127  else if (key == "Type") {
128  if (!(val == "Linear" || val == "linear"))
129  throw GeographicErr("Only linear models are supported");
130  } else if (key == "Epoch")
131  _t0 = Utility::val<real>(val);
132  else if (key == "DeltaEpoch")
133  _dt0 = Utility::val<real>(val);
134  else if (key == "NumModels")
135  _nNmodels = Utility::val<int>(val);
136  else if (key == "NumConstants")
137  _nNconstants = Utility::val<int>(val);
138  else if (key == "MinTime")
139  _tmin = Utility::val<real>(val);
140  else if (key == "MaxTime")
141  _tmax = Utility::val<real>(val);
142  else if (key == "MinHeight")
143  _hmin = Utility::val<real>(val);
144  else if (key == "MaxHeight")
145  _hmax = Utility::val<real>(val);
146  else if (key == "Normalization") {
147  if (val == "FULL" || val == "Full" || val == "full")
148  _norm = SphericalHarmonic::FULL;
149  else if (val == "SCHMIDT" || val == "Schmidt" || val == "schmidt")
151  else
152  throw GeographicErr("Unknown normalization " + val);
153  } else if (key == "ByteOrder") {
154  if (val == "Big" || val == "big")
155  throw GeographicErr("Only little-endian ordering is supported");
156  else if (!(val == "Little" || val == "little"))
157  throw GeographicErr("Unknown byte ordering " + val);
158  } else if (key == "ID")
159  _id = val;
160  // else unrecognized keywords are skipped
161  }
162  // Check values
163  if (!(isfinite(_a) && _a > 0))
164  throw GeographicErr("Reference radius must be positive");
165  if (!(_t0 > 0))
166  throw GeographicErr("Epoch time not defined");
167  if (_tmin >= _tmax)
168  throw GeographicErr("Min time exceeds max time");
169  if (_hmin >= _hmax)
170  throw GeographicErr("Min height exceeds max height");
171  if (int(_id.size()) != idlength_)
172  throw GeographicErr("Invalid ID");
173  if (_nNmodels < 1)
174  throw GeographicErr("NumModels must be positive");
175  if (!(_nNconstants == 0 || _nNconstants == 1))
176  throw GeographicErr("NumConstants must be 0 or 1");
177  if (!(_dt0 > 0)) {
178  if (_nNmodels > 1)
179  throw GeographicErr("DeltaEpoch must be positive");
180  else
181  _dt0 = 1;
182  }
183  }
184 
185  void MagneticModel::FieldGeocentric(real t, real X, real Y, real Z,
186  real& BX, real& BY, real& BZ,
187  real& BXt, real& BYt, real& BZt) const {
188  t -= _t0;
189  int n = max(min(int(floor(t / _dt0)), _nNmodels - 1), 0);
190  bool interpolate = n + 1 < _nNmodels;
191  t -= n * _dt0;
192  // Components in geocentric basis
193  // initial values to suppress warning
194  real BXc = 0, BYc = 0, BZc = 0;
195  _harm[n](X, Y, Z, BX, BY, BZ);
196  _harm[n + 1](X, Y, Z, BXt, BYt, BZt);
197  if (_nNconstants)
198  _harm[_nNmodels + 1](X, Y, Z, BXc, BYc, BZc);
199  if (interpolate) {
200  // Convert to a time derivative
201  BXt = (BXt - BX) / _dt0;
202  BYt = (BYt - BY) / _dt0;
203  BZt = (BZt - BZ) / _dt0;
204  }
205  BX += t * BXt + BXc;
206  BY += t * BYt + BYc;
207  BZ += t * BZt + BZc;
208 
209  BXt = BXt * - _a;
210  BYt = BYt * - _a;
211  BZt = BZt * - _a;
212 
213  BX *= - _a;
214  BY *= - _a;
215  BZ *= - _a;
216  }
217 
218  void MagneticModel::Field(real t, real lat, real lon, real h, bool diffp,
219  real& Bx, real& By, real& Bz,
220  real& Bxt, real& Byt, real& Bzt) const {
221  real X, Y, Z;
222  real M[Geocentric::dim2_];
223  _earth.IntForward(lat, lon, h, X, Y, Z, M);
224  // Components in geocentric basis
225  // initial values to suppress warning
226  real BX = 0, BY = 0, BZ = 0, BXt = 0, BYt = 0, BZt = 0;
227  FieldGeocentric(t, X, Y, Z, BX, BY, BZ, BXt, BYt, BZt);
228  if (diffp)
229  Geocentric::Unrotate(M, BXt, BYt, BZt, Bxt, Byt, Bzt);
230  Geocentric::Unrotate(M, BX, BY, BZ, Bx, By, Bz);
231  }
232 
233  MagneticCircle MagneticModel::Circle(real t, real lat, real h) const {
234  real t1 = t - _t0;
235  int n = max(min(int(floor(t1 / _dt0)), _nNmodels - 1), 0);
236  bool interpolate = n + 1 < _nNmodels;
237  t1 -= n * _dt0;
238  real X, Y, Z, M[Geocentric::dim2_];
239  _earth.IntForward(lat, 0, h, X, Y, Z, M);
240  // Y = 0, cphi = M[7], sphi = M[8];
241 
242  return (_nNconstants == 0 ?
243  MagneticCircle(_a, _earth._f, lat, h, t,
244  M[7], M[8], t1, _dt0, interpolate,
245  _harm[n].Circle(X, Z, true),
246  _harm[n + 1].Circle(X, Z, true)) :
247  MagneticCircle(_a, _earth._f, lat, h, t,
248  M[7], M[8], t1, _dt0, interpolate,
249  _harm[n].Circle(X, Z, true),
250  _harm[n + 1].Circle(X, Z, true),
251  _harm[_nNmodels + 1].Circle(X, Z, true)));
252  }
253 
254  void MagneticModel::FieldComponents(real Bx, real By, real Bz,
255  real Bxt, real Byt, real Bzt,
256  real& H, real& F, real& D, real& I,
257  real& Ht, real& Ft,
258  real& Dt, real& It) {
259  H = hypot(Bx, By);
260  Ht = H != 0 ? (Bx * Bxt + By * Byt) / H : hypot(Bxt, Byt);
261  D = H != 0 ? Math::atan2d(Bx, By) : Math::atan2d(Bxt, Byt);
262  Dt = (H != 0 ? (By * Bxt - Bx * Byt) / Math::sq(H) : 0) / Math::degree();
263  F = hypot(H, Bz);
264  Ft = F != 0 ? (H * Ht + Bz * Bzt) / F : hypot(Ht, Bzt);
265  I = F != 0 ? Math::atan2d(-Bz, H) : Math::atan2d(-Bzt, Ht);
266  It = (F != 0 ? (Bz * Ht - H * Bzt) / Math::sq(F) : 0) / Math::degree();
267  }
268 
270  string path;
271  char* magneticpath = getenv("GEOGRAPHICLIB_MAGNETIC_PATH");
272  if (magneticpath)
273  path = string(magneticpath);
274  if (!path.empty())
275  return path;
276  char* datapath = getenv("GEOGRAPHICLIB_DATA");
277  if (datapath)
278  path = string(datapath);
279  return (!path.empty() ? path : string(GEOGRAPHICLIB_DATA)) + "/magnetic";
280  }
281 
283  string name;
284  char* magneticname = getenv("GEOGRAPHICLIB_MAGNETIC_NAME");
285  if (magneticname)
286  name = string(magneticname);
287  return !name.empty() ? name : string(GEOGRAPHICLIB_MAGNETIC_DEFAULT_NAME);
288  }
289 
290 } // namespace GeographicLib
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
Header for GeographicLib::MagneticCircle class.
#define GEOGRAPHICLIB_DATA
#define GEOGRAPHICLIB_MAGNETIC_DEFAULT_NAME
Header for GeographicLib::MagneticModel class.
Header for GeographicLib::SphericalEngine class.
Header for GeographicLib::Utility class.
Geocentric coordinates
Definition: Geocentric.hpp:67
Exception handling for GeographicLib.
Definition: Constants.hpp:316
Geomagnetic field on a circle of latitude.
void FieldGeocentric(real t, real X, real Y, real Z, real &BX, real &BY, real &BZ, real &BXt, real &BYt, real &BZt) const
MagneticCircle Circle(real t, real lat, real h) const
static std::string DefaultMagneticPath()
static std::string DefaultMagneticName()
static void FieldComponents(real Bx, real By, real Bz, real &H, real &F, real &D, real &I)
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:76
static T degree()
Definition: Math.hpp:200
static T atan2d(T y, T x)
Definition: Math.cpp:183
static T sq(T x)
Definition: Math.hpp:212
static void readcoeffs(std::istream &stream, int &N, int &M, std::vector< real > &C, std::vector< real > &S, bool truncate=false)
Spherical harmonic series.
static bool ParseLine(const std::string &line, std::string &key, std::string &value, char equals='\0', char comment='#')
Definition: Utility.cpp:170
Namespace for GeographicLib.
Definition: Accumulator.cpp:12