GeographicLib  1.50
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-2019) <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 "wmm2015v2"
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  , _Nmodels(1)
51  , _Nconstants(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  _G.resize(_Nmodels + 1 + _Nconstants);
67  _H.resize(_Nmodels + 1 + _Nconstants);
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 < _Nmodels + 1 + _Nconstants; ++i) {
81  int N, M;
82  if (truncate) { N = Nmax; M = Mmax; }
83  SphericalEngine::coeff::readcoeffs(coeffstr, N, M, _G[i], _H[i],
84  truncate);
85  if (!(M < 0 || _G[i][0] == 0))
86  throw GeographicErr("A degree 0 term is not permitted");
87  _harm.push_back(SphericalHarmonic(_G[i], _H[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 std::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  _Nmodels = Utility::val<int>(val);
136  else if (key == "NumConstants")
137  _Nconstants = 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 (!(Math::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 (_Nmodels < 1)
174  throw GeographicErr("NumModels must be positive");
175  if (!(_Nconstants == 0 || _Nconstants == 1))
176  throw GeographicErr("NumConstants must be 0 or 1");
177  if (!(_dt0 > 0)) {
178  if (_Nmodels > 1)
179  throw GeographicErr("DeltaEpoch must be positive");
180  else
181  _dt0 = 1;
182  }
183  }
184 
185  void MagneticModel::Field(real t, real lat, real lon, real h, bool diffp,
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)), _Nmodels - 1), 0);
190  bool interpolate = n + 1 < _Nmodels;
191  t -= n * _dt0;
192  real X, Y, Z;
193  real M[Geocentric::dim2_];
194  _earth.IntForward(lat, lon, h, X, Y, Z, M);
195  // Components in geocentric basis
196  // initial values to suppress warning
197  real BX0 = 0, BY0 = 0, BZ0 = 0, BX1 = 0, BY1 = 0, BZ1 = 0;
198  real BXc = 0, BYc = 0, BZc = 0;
199  _harm[n](X, Y, Z, BX0, BY0, BZ0);
200  _harm[n + 1](X, Y, Z, BX1, BY1, BZ1);
201  if (_Nconstants)
202  _harm[_Nmodels + 1](X, Y, Z, BXc, BYc, BZc);
203  if (interpolate) {
204  // Convert to a time derivative
205  BX1 = (BX1 - BX0) / _dt0;
206  BY1 = (BY1 - BY0) / _dt0;
207  BZ1 = (BZ1 - BZ0) / _dt0;
208  }
209  BX0 += t * BX1 + BXc;
210  BY0 += t * BY1 + BYc;
211  BZ0 += t * BZ1 + BZc;
212  if (diffp) {
213  Geocentric::Unrotate(M, BX1, BY1, BZ1, Bxt, Byt, Bzt);
214  Bxt *= - _a;
215  Byt *= - _a;
216  Bzt *= - _a;
217  }
218  Geocentric::Unrotate(M, BX0, BY0, BZ0, Bx, By, Bz);
219  Bx *= - _a;
220  By *= - _a;
221  Bz *= - _a;
222  }
223 
224  MagneticCircle MagneticModel::Circle(real t, real lat, real h) const {
225  real t1 = t - _t0;
226  int n = max(min(int(floor(t1 / _dt0)), _Nmodels - 1), 0);
227  bool interpolate = n + 1 < _Nmodels;
228  t1 -= n * _dt0;
229  real X, Y, Z, M[Geocentric::dim2_];
230  _earth.IntForward(lat, 0, h, X, Y, Z, M);
231  // Y = 0, cphi = M[7], sphi = M[8];
232 
233  return (_Nconstants == 0 ?
234  MagneticCircle(_a, _earth._f, lat, h, t,
235  M[7], M[8], t1, _dt0, interpolate,
236  _harm[n].Circle(X, Z, true),
237  _harm[n + 1].Circle(X, Z, true)) :
238  MagneticCircle(_a, _earth._f, lat, h, t,
239  M[7], M[8], t1, _dt0, interpolate,
240  _harm[n].Circle(X, Z, true),
241  _harm[n + 1].Circle(X, Z, true),
242  _harm[_Nmodels + 1].Circle(X, Z, true)));
243  }
244 
245  void MagneticModel::FieldComponents(real Bx, real By, real Bz,
246  real Bxt, real Byt, real Bzt,
247  real& H, real& F, real& D, real& I,
248  real& Ht, real& Ft,
249  real& Dt, real& It) {
250  H = Math::hypot(Bx, By);
251  Ht = H != 0 ? (Bx * Bxt + By * Byt) / H : Math::hypot(Bxt, Byt);
252  D = H != 0 ? Math::atan2d(Bx, By) : Math::atan2d(Bxt, Byt);
253  Dt = (H != 0 ? (By * Bxt - Bx * Byt) / Math::sq(H) : 0) / Math::degree();
254  F = Math::hypot(H, Bz);
255  Ft = F != 0 ? (H * Ht + Bz * Bzt) / F : Math::hypot(Ht, Bzt);
256  I = F != 0 ? Math::atan2d(-Bz, H) : Math::atan2d(-Bzt, Ht);
257  It = (F != 0 ? (Bz * Ht - H * Bzt) / Math::sq(F) : 0) / Math::degree();
258  }
259 
261  string path;
262  char* magneticpath = getenv("GEOGRAPHICLIB_MAGNETIC_PATH");
263  if (magneticpath)
264  path = string(magneticpath);
265  if (!path.empty())
266  return path;
267  char* datapath = getenv("GEOGRAPHICLIB_DATA");
268  if (datapath)
269  path = string(datapath);
270  return (!path.empty() ? path : string(GEOGRAPHICLIB_DATA)) + "/magnetic";
271  }
272 
274  string name;
275  char* magneticname = getenv("GEOGRAPHICLIB_MAGNETIC_NAME");
276  if (magneticname)
277  name = string(magneticname);
278  return !name.empty() ? name : string(GEOGRAPHICLIB_MAGNETIC_DEFAULT_NAME);
279  }
280 
281 } // namespace GeographicLib
static void readcoeffs(std::istream &stream, int &N, int &M, std::vector< real > &C, std::vector< real > &S, bool truncate=false)
static bool isfinite(T x)
Definition: Math.cpp:372
Header for GeographicLib::Utility class.
Header for GeographicLib::MagneticCircle class.
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:98
#define GEOGRAPHICLIB_MAGNETIC_DEFAULT_NAME
Geomagnetic field on a circle of latitude.
static void FieldComponents(real Bx, real By, real Bz, real &H, real &F, real &D, real &I)
Geocentric coordinates
Definition: Geocentric.hpp:67
static T atan2d(T y, T x)
Definition: Math.cpp:310
static T sq(T x)
Definition: Math.hpp:201
#define GEOGRAPHICLIB_DATA
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
Header for GeographicLib::MagneticModel class.
static T degree()
Definition: Math.hpp:185
static std::string DefaultMagneticName()
Exception handling for GeographicLib.
Definition: Constants.hpp:390
static T hypot(T x, T y)
Definition: Math.cpp:58
static std::string DefaultMagneticPath()
Spherical harmonic series.
static bool ParseLine(const std::string &line, std::string &key, std::string &val)
Definition: Utility.cpp:22
Header for GeographicLib::SphericalEngine class.
MagneticCircle Circle(real t, real lat, real h) const