src/MagneticModel.cpp
Go to the documentation of this file.
1 
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())
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, 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::val<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::val<real>(val);
120  else if (key == "DeltaEpoch")
121  _dt0 = Utility::val<real>(val);
122  else if (key == "NumModels")
123  _Nmodels = Utility::val<int>(val);
124  else if (key == "NumConstants")
125  _Nconstants = Utility::val<int>(val);
126  else if (key == "MinTime")
127  _tmin = Utility::val<real>(val);
128  else if (key == "MaxTime")
129  _tmax = Utility::val<real>(val);
130  else if (key == "MinHeight")
131  _hmin = Utility::val<real>(val);
132  else if (key == "MaxHeight")
133  _hmax = Utility::val<real>(val);
134  else if (key == "Normalization") {
135  if (val == "FULL" || val == "Full" || val == "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;
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 
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 
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 != 0 ? (Bx * Bxt + By * Byt) / H : Math::hypot(Bxt, Byt);
240  D = H != 0 ? Math::atan2d(Bx, By) : Math::atan2d(Bxt, Byt);
241  Dt = (H != 0 ? (By * Bxt - Bx * Byt) / Math::sq(H) : 0) / Math::degree();
242  F = Math::hypot(H, Bz);
243  Ft = F != 0 ? (H * Ht + Bz * Bzt) / F : Math::hypot(Ht, Bzt);
244  I = F != 0 ? Math::atan2d(-Bz, H) : Math::atan2d(-Bzt, Ht);
245  It = (F != 0 ? (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
const gtsam::Symbol key('X', 0)
const char Y
#define max(a, b)
Definition: datatypes.h:20
#define I
Definition: main.h:112
Matrix< RealScalar, Dynamic, Dynamic > M
Definition: bench_gemm.cpp:51
Key F(std::uint64_t j)
static const double lat
Header for GeographicLib::Utility class.
#define min(a, b)
Definition: datatypes.h:19
static bool isfinite(T x)
Definition: Math.hpp:806
Header for GeographicLib::MagneticCircle class.
int n
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:102
Definition: BFloat16.h:88
std::vector< std::vector< real > > _G
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy y set format x g set format y g set format x2 g set format y2 g set format z g set angles radians set nogrid set key title set key left top Right noreverse box linetype linewidth samplen spacing width set nolabel set noarrow set nologscale set logscale x set set pointsize set encoding default set nopolar set noparametric set set set set surface set nocontour set clabel set mapping cartesian set nohidden3d set cntrparam order set cntrparam linear set cntrparam levels auto set cntrparam points set size set set xzeroaxis lt lw set x2zeroaxis lt lw set yzeroaxis lt lw set y2zeroaxis lt lw set tics in set ticslevel set tics set mxtics default set mytics default set mx2tics default set my2tics default set xtics border mirror norotate autofreq set ytics border mirror norotate autofreq set ztics border nomirror norotate autofreq set nox2tics set noy2tics set timestamp bottom norotate set rrange [*:*] noreverse nowriteback set trange [*:*] noreverse nowriteback set urange [*:*] noreverse nowriteback set vrange [*:*] noreverse nowriteback set xlabel matrix size set x2label set timefmt d m y n H
#define N
Definition: gksort.c:12
Geomagnetic field on a circle of latitude.
MagneticModel(const MagneticModel &)
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
EIGEN_DEVICE_FUNC const FloorReturnType floor() const
#define Z
Definition: icosphere.cpp:21
static T hypot(T x, T y)
Definition: Math.hpp:243
static T sq(T x)
Definition: Math.hpp:232
static T atan2d(T y, T x)
Definition: Math.hpp:691
Definition: main.h:100
Namespace for GeographicLib.
void IntForward(real lat, real lon, real h, real &X, real &Y, real &Z, real M[dim2_]) const
void Field(real t, real lat, real lon, real h, bool diffp, real &Bx, real &By, real &Bz, real &Bxt, real &Byt, real &Bzt) const
void ReadMetadata(const std::string &name)
Header for GeographicLib::MagneticModel class.
static T degree()
Definition: Math.hpp:216
static std::string DefaultMagneticName()
std::vector< SphericalHarmonic > _harm
const double h
static const size_t dim2_
Definition: Geocentric.hpp:77
Exception handling for GeographicLib.
Definition: Constants.hpp:389
static void Unrotate(real M[dim2_], real X, real Y, real Z, real &x, real &y, real &z)
Definition: Geocentric.hpp:89
static EIGEN_DEPRECATED const end_t end
static std::string DefaultMagneticPath()
static const double lon
#define GEOGRAPHICLIB_MAGNETIC_DEFAULT_NAME
Annotation for function names.
Definition: attr.h:48
#define X
Definition: icosphere.cpp:20
Spherical harmonic series.
T interpolate(const T &X, const T &Y, double t, typename MakeOptionalJacobian< T, T >::type Hx={}, typename MakeOptionalJacobian< T, T >::type Hy={})
Definition: Lie.h:327
static bool ParseLine(const std::string &line, std::string &key, std::string &val)
Definition: Utility.cpp:22
Header for GeographicLib::SphericalEngine class.
std::vector< std::vector< real > > _H
SphericalHarmonic::normalization _norm
Point2 t(10, 10)
MagneticCircle Circle(real t, real lat, real h) const
#define GEOGRAPHICLIB_DATA


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:35