GeographicLib  1.43
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-2015) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * http://geographiclib.sourceforge.net/
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 "wmm2015"
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)
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  , _norm(SphericalHarmonic::SCHMIDT)
53  , _earth(earth)
54  {
55  if (_dir.empty())
56  _dir = DefaultMagneticPath();
57  ReadMetadata(_name);
58  _G.resize(_Nmodels + 1 + _Nconstants);
59  _H.resize(_Nmodels + 1 + _Nconstants);
60  {
61  string coeff = _filename + ".cof";
62  ifstream coeffstr(coeff.c_str(), ios::binary);
63  if (!coeffstr.good())
64  throw GeographicErr("Error opening " + coeff);
65  char id[idlength_ + 1];
66  coeffstr.read(id, idlength_);
67  if (!coeffstr.good())
68  throw GeographicErr("No header in " + coeff);
69  id[idlength_] = '\0';
70  if (_id != string(id))
71  throw GeographicErr("ID mismatch: " + _id + " vs " + id);
72  for (int i = 0; i < _Nmodels + 1 + _Nconstants; ++i) {
73  int N, M;
74  SphericalEngine::coeff::readcoeffs(coeffstr, N, M, _G[i], _H[i]);
75  if (!(M < 0 || _G[i][0] == 0))
76  throw GeographicErr("A degree 0 term is not permitted");
77  _harm.push_back(SphericalHarmonic(_G[i], _H[i], N, N, M, _a, _norm));
78  }
79  int pos = int(coeffstr.tellg());
80  coeffstr.seekg(0, ios::end);
81  if (pos != coeffstr.tellg())
82  throw GeographicErr("Extra data in " + coeff);
83  }
84  }
85 
86  void MagneticModel::ReadMetadata(const std::string& name) {
87  const char* spaces = " \t\n\v\f\r";
88  _filename = _dir + "/" + name + ".wmm";
89  ifstream metastr(_filename.c_str());
90  if (!metastr.good())
91  throw GeographicErr("Cannot open " + _filename);
92  string line;
93  getline(metastr, line);
94  if (!(line.size() >= 6 && line.substr(0,5) == "WMMF-"))
95  throw GeographicErr(_filename + " does not contain WMMF-n signature");
96  string::size_type n = line.find_first_of(spaces, 5);
97  if (n != string::npos)
98  n -= 5;
99  string version = line.substr(5, n);
100  if (!(version == "1" || version == "2"))
101  throw GeographicErr("Unknown version in " + _filename + ": " + version);
102  string key, val;
103  while (getline(metastr, line)) {
104  if (!Utility::ParseLine(line, key, val))
105  continue;
106  // Process key words
107  if (key == "Name")
108  _name = val;
109  else if (key == "Description")
110  _description = val;
111  else if (key == "ReleaseDate")
112  _date = val;
113  else if (key == "Radius")
114  _a = Utility::num<real>(val);
115  else if (key == "Type") {
116  if (!(val == "Linear" || val == "linear"))
117  throw GeographicErr("Only linear models are supported");
118  } else if (key == "Epoch")
119  _t0 = Utility::num<real>(val);
120  else if (key == "DeltaEpoch")
121  _dt0 = Utility::num<real>(val);
122  else if (key == "NumModels")
123  _Nmodels = Utility::num<int>(val);
124  else if (key == "NumConstants")
125  _Nconstants = Utility::num<int>(val);
126  else if (key == "MinTime")
127  _tmin = Utility::num<real>(val);
128  else if (key == "MaxTime")
129  _tmax = Utility::num<real>(val);
130  else if (key == "MinHeight")
131  _hmin = Utility::num<real>(val);
132  else if (key == "MaxHeight")
133  _hmax = Utility::num<real>(val);
134  else if (key == "Normalization") {
135  if (val == "FULL" || val == "Full" || val == "full")
136  _norm = SphericalHarmonic::FULL;
137  else if (val == "SCHMIDT" || val == "Schmidt" || val == "schmidt")
139  else
140  throw GeographicErr("Unknown normalization " + val);
141  } else if (key == "ByteOrder") {
142  if (val == "Big" || val == "big")
143  throw GeographicErr("Only little-endian ordering is supported");
144  else if (!(val == "Little" || val == "little"))
145  throw GeographicErr("Unknown byte ordering " + val);
146  } else if (key == "ID")
147  _id = val;
148  // else unrecognized keywords are skipped
149  }
150  // Check values
151  if (!(Math::isfinite(_a) && _a > 0))
152  throw GeographicErr("Reference radius must be positive");
153  if (!(_t0 > 0))
154  throw GeographicErr("Epoch time not defined");
155  if (_tmin >= _tmax)
156  throw GeographicErr("Min time exceeds max time");
157  if (_hmin >= _hmax)
158  throw GeographicErr("Min height exceeds max height");
159  if (int(_id.size()) != idlength_)
160  throw GeographicErr("Invalid ID");
161  if (_Nmodels < 1)
162  throw GeographicErr("NumModels must be positive");
163  if (!(_Nconstants == 0 || _Nconstants == 1))
164  throw GeographicErr("NumConstants must be 0 or 1");
165  if (!(_dt0 > 0)) {
166  if (_Nmodels > 1)
167  throw GeographicErr("DeltaEpoch must be positive");
168  else
169  _dt0 = 1;
170  }
171  }
172 
173  void MagneticModel::Field(real t, real lat, real lon, real h, bool diffp,
174  real& Bx, real& By, real& Bz,
175  real& Bxt, real& Byt, real& Bzt) const {
176  t -= _t0;
177  int n = max(min(int(floor(t / _dt0)), _Nmodels - 1), 0);
178  bool interpolate = n + 1 < _Nmodels;
179  t -= n * _dt0;
180  real X, Y, Z;
181  real M[Geocentric::dim2_];
182  _earth.IntForward(lat, lon, h, X, Y, Z, M);
183  // Components in geocentric basis
184  // initial values to suppress warning
185  real BX0 = 0, BY0 = 0, BZ0 = 0, BX1 = 0, BY1 = 0, BZ1 = 0;
186  real BXc = 0, BYc = 0, BZc = 0;
187  _harm[n](X, Y, Z, BX0, BY0, BZ0);
188  _harm[n + 1](X, Y, Z, BX1, BY1, BZ1);
189  if (_Nconstants)
190  _harm[_Nmodels + 1](X, Y, Z, BXc, BYc, BZc);
191  if (interpolate) {
192  // Convert to a time derivative
193  BX1 = (BX1 - BX0) / _dt0;
194  BY1 = (BY1 - BY0) / _dt0;
195  BZ1 = (BZ1 - BZ0) / _dt0;
196  }
197  BX0 += t * BX1 + BXc;
198  BY0 += t * BY1 + BYc;
199  BZ0 += t * BZ1 + BZc;
200  if (diffp) {
201  Geocentric::Unrotate(M, BX1, BY1, BZ1, Bxt, Byt, Bzt);
202  Bxt *= - _a;
203  Byt *= - _a;
204  Bzt *= - _a;
205  }
206  Geocentric::Unrotate(M, BX0, BY0, BZ0, Bx, By, Bz);
207  Bx *= - _a;
208  By *= - _a;
209  Bz *= - _a;
210  }
211 
212  MagneticCircle MagneticModel::Circle(real t, real lat, real h) const {
213  real t1 = t - _t0;
214  int n = max(min(int(floor(t1 / _dt0)), _Nmodels - 1), 0);
215  bool interpolate = n + 1 < _Nmodels;
216  t1 -= n * _dt0;
217  real X, Y, Z, M[Geocentric::dim2_];
218  _earth.IntForward(lat, 0, h, X, Y, Z, M);
219  // Y = 0, cphi = M[7], sphi = M[8];
220 
221  return (_Nconstants == 0 ?
222  MagneticCircle(_a, _earth._f, lat, h, t,
223  M[7], M[8], t1, _dt0, interpolate,
224  _harm[n].Circle(X, Z, true),
225  _harm[n + 1].Circle(X, Z, true)) :
226  MagneticCircle(_a, _earth._f, lat, h, t,
227  M[7], M[8], t1, _dt0, interpolate,
228  _harm[n].Circle(X, Z, true),
229  _harm[n + 1].Circle(X, Z, true),
230  _harm[_Nmodels + 1].Circle(X, Z, true)));
231  }
232 
233  void MagneticModel::FieldComponents(real Bx, real By, real Bz,
234  real Bxt, real Byt, real Bzt,
235  real& H, real& F, real& D, real& I,
236  real& Ht, real& Ft,
237  real& Dt, real& It) {
238  H = Math::hypot(Bx, By);
239  Ht = H ? (Bx * Bxt + By * Byt) / H : Math::hypot(Bxt, Byt);
240  D = H ? Math::atan2d(Bx, By) : Math::atan2d(Bxt, Byt);
241  Dt = (H ? (By * Bxt - Bx * Byt) / Math::sq(H) : 0) / Math::degree();
242  F = Math::hypot(H, Bz);
243  Ft = F ? (H * Ht + Bz * Bzt) / F : Math::hypot(Ht, Bzt);
244  I = F ? Math::atan2d(-Bz, H) : Math::atan2d(-Bzt, Ht);
245  It = (F ? (Bz * Ht - H * Bzt) / Math::sq(F) : 0) / Math::degree();
246  }
247 
249  string path;
250  char* magneticpath = getenv("GEOGRAPHICLIB_MAGNETIC_PATH");
251  if (magneticpath)
252  path = string(magneticpath);
253  if (!path.empty())
254  return path;
255  char* datapath = getenv("GEOGRAPHICLIB_DATA");
256  if (datapath)
257  path = string(datapath);
258  return (!path.empty() ? path : string(GEOGRAPHICLIB_DATA)) + "/magnetic";
259  }
260 
262  string name;
263  char* magneticname = getenv("GEOGRAPHICLIB_MAGNETIC_NAME");
264  if (magneticname)
265  name = string(magneticname);
266  return !name.empty() ? name : string(GEOGRAPHICLIB_MAGNETIC_DEFAULT_NAME);
267  }
268 
269 } // namespace GeographicLib
Header for GeographicLib::Utility class.
static bool isfinite(T x)
Definition: Math.hpp:614
Header for GeographicLib::MagneticCircle class.
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:102
#define GEOGRAPHICLIB_MAGNETIC_DEFAULT_NAME
Geomagnetic field on a circle of latitude.
MagneticCircle Circle(real t, real lat, real h) const
static void FieldComponents(real Bx, real By, real Bz, real &H, real &F, real &D, real &I)
static void readcoeffs(std::istream &stream, int &N, int &M, std::vector< real > &C, std::vector< real > &S)
Geocentric coordinates
Definition: Geocentric.hpp:67
static T hypot(T x, T y)
Definition: Math.hpp:255
static T sq(T x)
Definition: Math.hpp:244
static T atan2d(T y, T x)
Definition: Math.hpp:551
#define GEOGRAPHICLIB_DATA
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
Header for GeographicLib::MagneticModel class.
static T degree()
Definition: Math.hpp:228
static std::string DefaultMagneticName()
Exception handling for GeographicLib.
Definition: Constants.hpp:382
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.