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
10#include <GeographicLib/MagneticModel.hpp>
11#include <fstream>
12#include <GeographicLib/SphericalEngine.hpp>
13#include <GeographicLib/MagneticCircle.hpp>
14#include <GeographicLib/Utility.hpp>
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
33namespace 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")
150 _norm = SphericalHarmonic::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
269 string MagneticModel::DefaultMagneticPath() {
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
282 string MagneticModel::DefaultMagneticName() {
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
real
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
MagneticCircle.hpp
Header for GeographicLib::MagneticCircle class.
GEOGRAPHICLIB_DATA
#define GEOGRAPHICLIB_DATA
Definition: MagneticModel.cpp:20
GEOGRAPHICLIB_MAGNETIC_DEFAULT_NAME
#define GEOGRAPHICLIB_MAGNETIC_DEFAULT_NAME
Definition: MagneticModel.cpp:25
MagneticModel.hpp
Header for GeographicLib::MagneticModel class.
SphericalEngine.hpp
Header for GeographicLib::SphericalEngine class.
Utility.hpp
Header for GeographicLib::Utility class.
GeographicLib::Geocentric
Geocentric coordinates
Definition: Geocentric.hpp:67
GeographicLib::GeographicErr
Exception handling for GeographicLib.
Definition: Constants.hpp:316
GeographicLib::MagneticCircle
Geomagnetic field on a circle of latitude.
Definition: MagneticCircle.hpp:37
GeographicLib::MagneticModel::FieldGeocentric
void FieldGeocentric(real t, real X, real Y, real Z, real &BX, real &BY, real &BZ, real &BXt, real &BYt, real &BZt) const
Definition: MagneticModel.cpp:185
GeographicLib::MagneticModel::Circle
MagneticCircle Circle(real t, real lat, real h) const
Definition: MagneticModel.cpp:233
GeographicLib::MagneticModel::DefaultMagneticPath
static std::string DefaultMagneticPath()
Definition: MagneticModel.cpp:269
GeographicLib::MagneticModel::DefaultMagneticName
static std::string DefaultMagneticName()
Definition: MagneticModel.cpp:282
GeographicLib::MagneticModel::FieldComponents
static void FieldComponents(real Bx, real By, real Bz, real &H, real &F, real &D, real &I)
Definition: MagneticModel.hpp:249
GeographicLib::Math
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:76
GeographicLib::Math::degree
static T degree()
Definition: Math.hpp:200
GeographicLib::Math::atan2d
static T atan2d(T y, T x)
Definition: Math.cpp:183
GeographicLib::Math::sq
static T sq(T x)
Definition: Math.hpp:212
GeographicLib::SphericalEngine::coeff::readcoeffs
static void readcoeffs(std::istream &stream, int &N, int &M, std::vector< real > &C, std::vector< real > &S, bool truncate=false)
Definition: SphericalEngine.cpp:385
GeographicLib::SphericalHarmonic
Spherical harmonic series.
Definition: SphericalHarmonic.hpp:69
GeographicLib::Utility::ParseLine
static bool ParseLine(const std::string &line, std::string &key, std::string &value, char equals='\0', char comment='#')
Definition: Utility.cpp:170
GeographicLib
Namespace for GeographicLib.
Definition: Accumulator.cpp:12