16#if !defined(GEOGRAPHICLIB_DATA)
18# define GEOGRAPHICLIB_DATA "C:/ProgramData/GeographicLib"
20# define GEOGRAPHICLIB_DATA "/usr/local/share/GeographicLib"
24#if !defined(GEOGRAPHICLIB_MAGNETIC_DEFAULT_NAME)
25# define GEOGRAPHICLIB_MAGNETIC_DEFAULT_NAME "wmm2025"
30# pragma warning (disable: 4996)
37 MagneticModel::MagneticModel(
const std::string& name,
const std::string& path,
41 , _description(
"NONE")
59 bool truncate = Nmax >= 0 || Mmax >= 0;
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();
66 _gG.resize(_nNmodels + 1 + _nNconstants);
67 _hH.resize(_nNmodels + 1 + _nNconstants);
69 string coeff = _filename +
".cof";
70 ifstream coeffstr(coeff.c_str(), ios::binary);
73 char id[idlength_ + 1];
74 coeffstr.read(
id, idlength_);
78 if (_id !=
string(
id))
80 for (
int i = 0; i < _nNmodels + 1 + _nNconstants; ++i) {
82 if (truncate) { N = Nmax; M = Mmax; }
85 if (!(M < 0 || _gG[i][0] == 0))
88 _nmx = max(_nmx, _harm.back().Coefficients().nmx());
89 _mmx = max(_mmx, _harm.back().Coefficients().mmx());
91 int pos = int(coeffstr.tellg());
92 coeffstr.seekg(0, ios::end);
93 if (pos != coeffstr.tellg())
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());
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)
111 string version(line, 5, n);
117 if (!(version ==
"1" || version ==
"2"))
118 throw GeographicErr(
"Unknown version in " + _filename +
": " + version);
120 while (getline(metastr, line)) {
126 else if (key ==
"Description")
128 else if (key ==
"ReleaseDate")
130 else if (key ==
"Radius")
131 _a = Utility::val<real>(val);
132 else if (key ==
"Type") {
133 if (!(val ==
"Linear" || val ==
"linear"))
135 }
else if (key ==
"Epoch")
136 _t0 = Utility::val<real>(val);
137 else if (key ==
"DeltaEpoch")
138 _dt0 = Utility::val<real>(val);
139 else if (key ==
"NumModels")
140 _nNmodels = Utility::val<int>(val);
141 else if (key ==
"NumConstants")
142 _nNconstants = Utility::val<int>(val);
143 else if (key ==
"MinTime")
144 _tmin = Utility::val<real>(val);
145 else if (key ==
"MaxTime")
146 _tmax = Utility::val<real>(val);
147 else if (key ==
"MinHeight")
148 _hmin = Utility::val<real>(val);
149 else if (key ==
"MaxHeight")
150 _hmax = Utility::val<real>(val);
151 else if (key ==
"Normalization") {
152 if (val ==
"FULL" || val ==
"Full" || val ==
"full")
154 else if (val ==
"SCHMIDT" || val ==
"Schmidt" || val ==
"schmidt")
157 throw GeographicErr(
"Unknown normalization " + val);
158 }
else if (key ==
"ByteOrder") {
159 if (val ==
"Big" || val ==
"big")
160 throw GeographicErr(
"Only little-endian ordering is supported");
161 else if (!(val ==
"Little" || val ==
"little"))
162 throw GeographicErr(
"Unknown byte ordering " + val);
163 }
else if (key ==
"ID")
168 if (!(isfinite(_a) && _a > 0))
169 throw GeographicErr(
"Reference radius must be positive");
171 throw GeographicErr(
"Epoch time not defined");
173 throw GeographicErr(
"Min time exceeds max time");
175 throw GeographicErr(
"Min height exceeds max height");
176 if (
int(_id.size()) != idlength_)
177 throw GeographicErr(
"Invalid ID");
179 throw GeographicErr(
"NumModels must be positive");
180 if (!(_nNconstants == 0 || _nNconstants == 1))
181 throw GeographicErr(
"NumConstants must be 0 or 1");
184 throw GeographicErr(
"DeltaEpoch must be positive");
191 real& BX, real& BY, real& BZ,
192 real& BXt, real& BYt, real& BZt)
const {
194 int n = max(min(
int(floor(t / _dt0)), _nNmodels - 1), 0);
195 bool interpolate = n + 1 < _nNmodels;
199 real BXc = 0, BYc = 0, BZc = 0;
200 _harm[n](X, Y, Z, BX, BY, BZ);
201 _harm[n + 1](X, Y, Z, BXt, BYt, BZt);
203 _harm[_nNmodels + 1](X, Y, Z, BXc, BYc, BZc);
206 BXt = (BXt - BX) / _dt0;
207 BYt = (BYt - BY) / _dt0;
208 BZt = (BZt - BZ) / _dt0;
227 real M[Geocentric::dim2_];
228 _earth.IntForward(lat, lon, h, X, Y, Z, M);
231 real BX = 0, BY = 0, BZ = 0, BXt = 0, BYt = 0, BZt = 0;
234 Geocentric::Unrotate(M, BXt, BYt, BZt, Bxt, Byt, Bzt);
235 Geocentric::Unrotate(M, BX, BY, BZ, Bx, By, Bz);
240 int n = max(min(
int(floor(t1 / _dt0)), _nNmodels - 1), 0);
241 bool interpolate = n + 1 < _nNmodels;
243 real X, Y, Z, M[Geocentric::dim2_];
244 _earth.IntForward(lat, 0, h, X, Y, Z, M);
247 return (_nNconstants == 0 ?
249 M[7], M[8], t1, _dt0, interpolate,
250 _harm[n].Circle(X, Z,
true),
251 _harm[n + 1].Circle(X, Z,
true)) :
253 M[7], M[8], t1, _dt0, interpolate,
254 _harm[n].Circle(X, Z,
true),
255 _harm[n + 1].Circle(X, Z,
true),
256 _harm[_nNmodels + 1].Circle(X, Z,
true)));
260 real Bxt, real Byt, real Bzt,
261 real& H, real& F, real& D, real& I,
263 real& Dt, real& It) {
265 Ht = H != 0 ? (Bx * Bxt + By * Byt) / H : hypot(Bxt, Byt);
269 Ft = F != 0 ? (H * Ht + Bz * Bzt) / F : hypot(Ht, Bzt);
276 char* magneticpath = getenv(
"GEOGRAPHICLIB_MAGNETIC_PATH");
278 path = string(magneticpath);
281 char* datapath = getenv(
"GEOGRAPHICLIB_DATA");
283 path = string(datapath);
289 char* magneticname = getenv(
"GEOGRAPHICLIB_MAGNETIC_NAME");
291 name = string(magneticname);
GeographicLib::Math::real real
#define GEOGRAPHICLIB_DATA
Header for GeographicLib::MagneticCircle class.
#define GEOGRAPHICLIB_MAGNETIC_DEFAULT_NAME
Header for GeographicLib::MagneticModel class.
Header for GeographicLib::SphericalEngine class.
Header for GeographicLib::Utility class.
Exception handling for GeographicLib.
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.
static T atan2d(T y, T x)
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='#')
Namespace for GeographicLib.