GeographicLib: MagneticModel.cpp Source File (original) (raw)

11#include

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

24#if !defined(GEOGRAPHICLIB_MAGNETIC_DEFAULT_NAME)

25# define GEOGRAPHICLIB_MAGNETIC_DEFAULT_NAME "wmm2025"

26#endif

28#if defined(_MSC_VER)

30# pragma warning (disable: 4996)

31#endif

35 using namespace std;

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)

55 , _earth(earth)

56 {

57 if (_dir.empty())

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::max();

63 if (Mmax < 0) Mmax = numeric_limits::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())

73 char id[idlength_ + 1];

74 coeffstr.read(id, idlength_);

75 if (!coeffstr.good())

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; }

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())

95 }

96 }

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())

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

113

114

115

116

117 if (!(version == "1" || version == "2"))

118 throw GeographicErr("Unknown version in " + _filename + ": " + version);

119 string key, val;

120 while (getline(metastr, line)) {

122 continue;

123

124 if (key == "Name")

125 _name = val;

126 else if (key == "Description")

127 _description = val;

128 else if (key == "ReleaseDate")

129 _date = val;

130 else if (key == "Radius")

131 _a = Utility::val(val);

132 else if (key == "Type") {

133 if (!(val == "Linear" || val == "linear"))

134 throw GeographicErr("Only linear models are supported");

135 } else if (key == "Epoch")

136 _t0 = Utility::val(val);

137 else if (key == "DeltaEpoch")

138 _dt0 = Utility::val(val);

139 else if (key == "NumModels")

140 _nNmodels = Utility::val(val);

141 else if (key == "NumConstants")

142 _nNconstants = Utility::val(val);

143 else if (key == "MinTime")

144 _tmin = Utility::val(val);

145 else if (key == "MaxTime")

146 _tmax = Utility::val(val);

147 else if (key == "MinHeight")

148 _hmin = Utility::val(val);

149 else if (key == "MaxHeight")

150 _hmax = Utility::val(val);

151 else if (key == "Normalization") {

152 if (val == "FULL" || val == "Full" || val == "full")

154 else if (val == "SCHMIDT" || val == "Schmidt" || val == "schmidt")

156 else

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")

164 _id = val;

165

166 }

167

168 if (!(isfinite(_a) && _a > 0))

169 throw GeographicErr("Reference radius must be positive");

170 if (!(_t0 > 0))

171 throw GeographicErr("Epoch time not defined");

172 if (_tmin >= _tmax)

173 throw GeographicErr("Min time exceeds max time");

174 if (_hmin >= _hmax)

175 throw GeographicErr("Min height exceeds max height");

176 if (int(_id.size()) != idlength_)

177 throw GeographicErr("Invalid ID");

178 if (_nNmodels < 1)

179 throw GeographicErr("NumModels must be positive");

180 if (!(_nNconstants == 0 || _nNconstants == 1))

181 throw GeographicErr("NumConstants must be 0 or 1");

182 if (!(_dt0 > 0)) {

183 if (_nNmodels > 1)

184 throw GeographicErr("DeltaEpoch must be positive");

185 else

186 _dt0 = 1;

187 }

188 }

189

191 real& BX, real& BY, real& BZ,

192 real& BXt, real& BYt, real& BZt) const {

193 t -= _t0;

194 int n = max(min(int(floor(t / _dt0)), _nNmodels - 1), 0);

195 bool interpolate = n + 1 < _nNmodels;

196 t -= n * _dt0;

197

198

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);

202 if (_nNconstants)

203 _harm[_nNmodels + 1](X, Y, Z, BXc, BYc, BZc);

204 if (interpolate) {

205

206 BXt = (BXt - BX) / _dt0;

207 BYt = (BYt - BY) / _dt0;

208 BZt = (BZt - BZ) / _dt0;

209 }

210 BX += t * BXt + BXc;

211 BY += t * BYt + BYc;

212 BZ += t * BZt + BZc;

213

214 BXt = BXt * - _a;

215 BYt = BYt * - _a;

216 BZt = BZt * - _a;

217

218 BX *= - _a;

219 BY *= - _a;

220 BZ *= - _a;

221 }

222

223 void MagneticModel::Field(real t, real lat, real lon, real h, bool diffp,

226 real X, Y, Z;

227 real M[Geocentric::dim2_];

228 _earth.IntForward(lat, lon, h, X, Y, Z, M);

229

230

231 real BX = 0, BY = 0, BZ = 0, BXt = 0, BYt = 0, BZt = 0;

233 if (diffp)

234 Geocentric::Unrotate(M, BXt, BYt, BZt, Bxt, Byt, Bzt);

235 Geocentric::Unrotate(M, BX, BY, BZ, Bx, By, Bz);

236 }

237

239 real t1 = t - _t0;

240 int n = max(min(int(floor(t1 / _dt0)), _nNmodels - 1), 0);

241 bool interpolate = n + 1 < _nNmodels;

242 t1 -= n * _dt0;

243 real X, Y, Z, M[Geocentric::dim2_];

244 _earth.IntForward(lat, 0, h, X, Y, Z, M);

245

246

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)));

257 }

258

260 real Bxt, real Byt, real Bzt,

261 real& H, real& F, real& D, real& I,

262 real& Ht, real& Ft,

263 real& Dt, real& It) {

264 H = hypot(Bx, By);

265 Ht = H != 0 ? (Bx * Bxt + By * Byt) / H : hypot(Bxt, Byt);

268 F = hypot(H, Bz);

269 Ft = F != 0 ? (H * Ht + Bz * Bzt) / F : hypot(Ht, Bzt);

272 }

273

275 string path;

276 char* magneticpath = getenv("GEOGRAPHICLIB_MAGNETIC_PATH");

277 if (magneticpath)

278 path = string(magneticpath);

279 if (!path.empty())

280 return path;

281 char* datapath = getenv("GEOGRAPHICLIB_DATA");

282 if (datapath)

283 path = string(datapath);

284 return (!path.empty() ? path : string(GEOGRAPHICLIB_DATA)) + "/magnetic";

285 }

286

288 string name;

289 char* magneticname = getenv("GEOGRAPHICLIB_MAGNETIC_NAME");

290 if (magneticname)

291 name = string(magneticname);

293 }

294

295}

GeographicLib::Math::real real

#define GEOGRAPHICLIB_DATA

Header for GeographicLib::MagneticCircle class.

#define GEOGRAPHICLIB_MAGNETIC_DEFAULT_NAME

Definition MagneticModel.cpp:25

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

Definition MagneticModel.cpp:190

MagneticCircle Circle(real t, real lat, real h) const

Definition MagneticModel.cpp:238

static std::string DefaultMagneticPath()

Definition MagneticModel.cpp:274

static std::string DefaultMagneticName()

Definition MagneticModel.cpp:287

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.