MagneticField.cpp
Go to the documentation of this file.
1 
12 #include <iostream>
13 #include <string>
14 #include <sstream>
15 #include <fstream>
18 #include <GeographicLib/DMS.hpp>
20 
21 #if defined(_MSC_VER)
22 // Squelch warnings about constant conditional expressions
23 # pragma warning (disable: 4127)
24 #endif
25 
26 #include "MagneticField.usage"
27 
28 int main(int argc, const char* const argv[]) {
29  try {
30  using namespace GeographicLib;
31  typedef Math::real real;
33  bool verbose = false, longfirst = false;
34  std::string dir;
36  std::string istring, ifile, ofile, cdelim;
37  char lsep = ';';
38  real time = 0, lat = 0, h = 0;
39  bool timeset = false, circle = false, rate = false;
40  real hguard = 500000, tguard = 50;
41  int prec = 1;
42 
43  for (int m = 1; m < argc; ++m) {
44  std::string arg(argv[m]);
45  if (arg == "-n") {
46  if (++m == argc) return usage(1, true);
47  model = argv[m];
48  } else if (arg == "-d") {
49  if (++m == argc) return usage(1, true);
50  dir = argv[m];
51  } else if (arg == "-t") {
52  if (++m == argc) return usage(1, true);
53  try {
54  time = Utility::fractionalyear<real>(std::string(argv[m]));
55  timeset = true;
56  circle = false;
57  }
58  catch (const std::exception& e) {
59  std::cerr << "Error decoding argument of " << arg << ": "
60  << e.what() << "\n";
61  return 1;
62  }
63  } else if (arg == "-c") {
64  if (m + 3 >= argc) return usage(1, true);
65  try {
66  using std::abs;
67  time = Utility::fractionalyear<real>(std::string(argv[++m]));
68  DMS::flag ind;
69  lat = DMS::Decode(std::string(argv[++m]), ind);
70  if (ind == DMS::LONGITUDE)
71  throw GeographicErr("Bad hemisphere letter on latitude");
72  if (!(abs(lat) <= 90))
73  throw GeographicErr("Latitude not in [-90d, 90d]");
74  h = Utility::val<real>(std::string(argv[++m]));
75  timeset = false;
76  circle = true;
77  }
78  catch (const std::exception& e) {
79  std::cerr << "Error decoding argument of " << arg << ": "
80  << e.what() << "\n";
81  return 1;
82  }
83  } else if (arg == "-r")
84  rate = !rate;
85  else if (arg == "-w")
86  longfirst = !longfirst;
87  else if (arg == "-p") {
88  if (++m == argc) return usage(1, true);
89  try {
90  prec = Utility::val<int>(std::string(argv[m]));
91  }
92  catch (const std::exception&) {
93  std::cerr << "Precision " << argv[m] << " is not a number\n";
94  return 1;
95  }
96  } else if (arg == "-T") {
97  if (++m == argc) return usage(1, true);
98  try {
99  tguard = Utility::val<real>(std::string(argv[m]));
100  }
101  catch (const std::exception& e) {
102  std::cerr << "Error decoding argument of " << arg << ": "
103  << e.what() << "\n";
104  return 1;
105  }
106  } else if (arg == "-H") {
107  if (++m == argc) return usage(1, true);
108  try {
109  hguard = Utility::val<real>(std::string(argv[m]));
110  }
111  catch (const std::exception& e) {
112  std::cerr << "Error decoding argument of " << arg << ": "
113  << e.what() << "\n";
114  return 1;
115  }
116  } else if (arg == "-v")
117  verbose = true;
118  else if (arg == "--input-string") {
119  if (++m == argc) return usage(1, true);
120  istring = argv[m];
121  } else if (arg == "--input-file") {
122  if (++m == argc) return usage(1, true);
123  ifile = argv[m];
124  } else if (arg == "--output-file") {
125  if (++m == argc) return usage(1, true);
126  ofile = argv[m];
127  } else if (arg == "--line-separator") {
128  if (++m == argc) return usage(1, true);
129  if (std::string(argv[m]).size() != 1) {
130  std::cerr << "Line separator must be a single character\n";
131  return 1;
132  }
133  lsep = argv[m][0];
134  } else if (arg == "--comment-delimiter") {
135  if (++m == argc) return usage(1, true);
136  cdelim = argv[m];
137  } else if (arg == "--version") {
138  std::cout << argv[0] << ": GeographicLib version "
139  << GEOGRAPHICLIB_VERSION_STRING << "\n";
140  return 0;
141  } else {
142  int retval = usage(!(arg == "-h" || arg == "--help"), arg != "--help");
143  if (arg == "-h")
144  std::cout<< "\nDefault magnetic path = \""
146  << "\"\nDefault magnetic name = \""
148  << "\"\n";
149  return retval;
150  }
151  }
152 
153  if (!ifile.empty() && !istring.empty()) {
154  std::cerr << "Cannot specify --input-string and --input-file together\n";
155  return 1;
156  }
157  if (ifile == "-") ifile.clear();
158  std::ifstream infile;
159  std::istringstream instring;
160  if (!ifile.empty()) {
161  infile.open(ifile.c_str());
162  if (!infile.is_open()) {
163  std::cerr << "Cannot open " << ifile << " for reading\n";
164  return 1;
165  }
166  } else if (!istring.empty()) {
167  std::string::size_type m = 0;
168  while (true) {
169  m = istring.find(lsep, m);
170  if (m == std::string::npos)
171  break;
172  istring[m] = '\n';
173  }
174  instring.str(istring);
175  }
176  std::istream* input = !ifile.empty() ? &infile :
177  (!istring.empty() ? &instring : &std::cin);
178 
179  std::ofstream outfile;
180  if (ofile == "-") ofile.clear();
181  if (!ofile.empty()) {
182  outfile.open(ofile.c_str());
183  if (!outfile.is_open()) {
184  std::cerr << "Cannot open " << ofile << " for writing\n";
185  return 1;
186  }
187  }
188  std::ostream* output = !ofile.empty() ? &outfile : &std::cout;
189 
190  tguard = std::max(real(0), tguard);
191  hguard = std::max(real(0), hguard);
192  prec = std::min(10 + Math::extra_digits(), std::max(0, prec));
193  int retval = 0;
194  try {
195  const MagneticModel m(model, dir);
196  if ((timeset || circle)
197  && (!Math::isfinite(time) ||
198  time < m.MinTime() - tguard ||
199  time > m.MaxTime() + tguard))
200  throw GeographicErr("Time " + Utility::str(time) +
201  " too far outside allowed range [" +
202  Utility::str(m.MinTime()) + "," +
203  Utility::str(m.MaxTime()) + "]");
204  if (circle
205  && (!Math::isfinite(h) ||
206  h < m.MinHeight() - hguard ||
207  h > m.MaxHeight() + hguard))
208  throw GeographicErr("Height " + Utility::str(h/1000) +
209  "km too far outside allowed range [" +
210  Utility::str(m.MinHeight()/1000) + "km," +
211  Utility::str(m.MaxHeight()/1000) + "km]");
212  if (verbose) {
213  std::cerr << "Magnetic file: " << m.MagneticFile() << "\n"
214  << "Name: " << m.MagneticModelName() << "\n"
215  << "Description: " << m.Description() << "\n"
216  << "Date & Time: " << m.DateTime() << "\n"
217  << "Time range: ["
218  << m.MinTime() << ","
219  << m.MaxTime() << "]\n"
220  << "Height range: ["
221  << m.MinHeight()/1000 << "km,"
222  << m.MaxHeight()/1000 << "km]\n";
223  }
224  if ((timeset || circle) && (time < m.MinTime() || time > m.MaxTime()))
225  std::cerr << "WARNING: Time " << time
226  << " outside allowed range ["
227  << m.MinTime() << "," << m.MaxTime() << "]\n";
228  if (circle && (h < m.MinHeight() || h > m.MaxHeight()))
229  std::cerr << "WARNING: Height " << h/1000
230  << "km outside allowed range ["
231  << m.MinHeight()/1000 << "km,"
232  << m.MaxHeight()/1000 << "km]\n";
233  const MagneticCircle c(circle ? m.Circle(time, lat, h) :
234  MagneticCircle());
235  std::string s, eol, stra, strb;
236  std::istringstream str;
237  while (std::getline(*input, s)) {
238  try {
239  eol = "\n";
240  if (!cdelim.empty()) {
241  std::string::size_type n = s.find(cdelim);
242  if (n != std::string::npos) {
243  eol = " " + s.substr(n) + "\n";
244  s = s.substr(0, n);
245  }
246  }
247  str.clear(); str.str(s);
248  if (!(timeset || circle)) {
249  if (!(str >> stra))
250  throw GeographicErr("Incomplete input: " + s);
251  time = Utility::fractionalyear<real>(stra);
252  if (time < m.MinTime() - tguard || time > m.MaxTime() + tguard)
253  throw GeographicErr("Time " + Utility::str(time) +
254  " too far outside allowed range [" +
255  Utility::str(m.MinTime()) + "," +
256  Utility::str(m.MaxTime()) +
257  "]");
258  if (time < m.MinTime() || time > m.MaxTime())
259  std::cerr << "WARNING: Time " << time
260  << " outside allowed range ["
261  << m.MinTime() << "," << m.MaxTime() << "]\n";
262  }
263  real lon;
264  if (circle) {
265  if (!(str >> strb))
266  throw GeographicErr("Incomplete input: " + s);
267  DMS::flag ind;
268  lon = DMS::Decode(strb, ind);
269  if (ind == DMS::LATITUDE)
270  throw GeographicErr("Bad hemisphere letter on " + strb);
271  } else {
272  if (!(str >> stra >> strb))
273  throw GeographicErr("Incomplete input: " + s);
274  DMS::DecodeLatLon(stra, strb, lat, lon, longfirst);
275  h = 0; // h is optional
276  if (str >> h) {
277  if (h < m.MinHeight() - hguard || h > m.MaxHeight() + hguard)
278  throw GeographicErr("Height " + Utility::str(h/1000) +
279  "km too far outside allowed range [" +
280  Utility::str(m.MinHeight()/1000) + "km," +
281  Utility::str(m.MaxHeight()/1000) + "km]");
282  if (h < m.MinHeight() || h > m.MaxHeight())
283  std::cerr << "WARNING: Height " << h/1000
284  << "km outside allowed range ["
285  << m.MinHeight()/1000 << "km,"
286  << m.MaxHeight()/1000 << "km]\n";
287  }
288  else
289  str.clear();
290  }
291  if (str >> stra)
292  throw GeographicErr("Extra junk in input: " + s);
293  real bx, by, bz, bxt, byt, bzt;
294  if (circle)
295  c(lon, bx, by, bz, bxt, byt, bzt);
296  else
297  m(time, lat, lon, h, bx, by, bz, bxt, byt, bzt);
298  real H, F, D, I, Ht, Ft, Dt, It;
299  MagneticModel::FieldComponents(bx, by, bz, bxt, byt, bzt,
300  H, F, D, I, Ht, Ft, Dt, It);
301 
302  *output << DMS::Encode(D, prec + 1, DMS::NUMBER) << " "
303  << DMS::Encode(I, prec + 1, DMS::NUMBER) << " "
304  << Utility::str(H, prec) << " "
305  << Utility::str(by, prec) << " "
306  << Utility::str(bx, prec) << " "
307  << Utility::str(-bz, prec) << " "
308  << Utility::str(F, prec) << eol;
309  if (rate)
310  *output << DMS::Encode(Dt, prec + 1, DMS::NUMBER) << " "
311  << DMS::Encode(It, prec + 1, DMS::NUMBER) << " "
312  << Utility::str(Ht, prec) << " "
313  << Utility::str(byt, prec) << " "
314  << Utility::str(bxt, prec) << " "
315  << Utility::str(-bzt, prec) << " "
316  << Utility::str(Ft, prec) << eol;
317  }
318  catch (const std::exception& e) {
319  *output << "ERROR: " << e.what() << "\n";
320  retval = 1;
321  }
322  }
323  }
324  catch (const std::exception& e) {
325  std::cerr << "Error reading " << model << ": " << e.what() << "\n";
326  retval = 1;
327  }
328  return retval;
329  }
330  catch (const std::exception& e) {
331  std::cerr << "Caught exception: " << e.what() << "\n";
332  return 1;
333  }
334  catch (...) {
335  std::cerr << "Caught unknown exception\n";
336  return 1;
337  }
338 }
Matrix3f m
Math::real MinHeight() const
#define max(a, b)
Definition: datatypes.h:20
#define I
Definition: main.h:112
float real
Definition: datatypes.h:10
const std::string & Description() const
Key F(std::uint64_t j)
static const double lat
noiseModel::Diagonal::shared_ptr model
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
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
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
Model of the earth&#39;s magnetic field.
Geomagnetic field on a circle of latitude.
static int extra_digits()
Definition: Math.hpp:187
static void FieldComponents(real Bx, real By, real Bz, real &H, real &F, real &D, real &I)
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
static std::string Encode(real angle, component trailing, unsigned prec, flag ind=NONE, char dmssep=char(0))
Definition: src/DMS.cpp:299
int main(int argc, const char *const argv[])
std::vector< int > ind
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const ArgReturnType arg() const
#define GEOGRAPHICLIB_VERSION_STRING
Definition: Config.h:3
static Math::real Decode(const std::string &dms, flag &ind)
Definition: src/DMS.cpp:28
Namespace for GeographicLib.
#define time
Array< double, 1, 3 > e(1./3., 0.5, 2.)
RealScalar s
Header for GeographicLib::MagneticModel class.
static std::string str(T x, int p=-1)
Definition: Utility.hpp:276
static std::string DefaultMagneticName()
const std::string & MagneticModelName() const
static void DecodeLatLon(const std::string &dmsa, const std::string &dmsb, real &lat, real &lon, bool longfirst=false)
Definition: src/DMS.cpp:252
static int set_digits(int ndigits=0)
Definition: Utility.cpp:48
const double h
Math::real MaxHeight() const
Exception handling for GeographicLib.
Definition: Constants.hpp:389
const std::string & MagneticFile() const
static std::string DefaultMagneticPath()
static const double lon
const std::string & DateTime() const
#define abs(x)
Definition: datatypes.h:17
MagneticCircle Circle(real t, real lat, real h) const
Header for GeographicLib::DMS class.


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