MagneticField.cpp
Go to the documentation of this file.
1 
12 #include <iostream>
13 #include <string>
14 #include <sstream>
15 #include <fstream>
16 #include <GeographicLib/MagneticModel.hpp>
17 #include <GeographicLib/MagneticCircle.hpp>
18 #include <GeographicLib/DMS.hpp>
19 #include <GeographicLib/Utility.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;
32  Utility::set_digits();
33  bool verbose = false, longfirst = false;
34  std::string dir;
35  std::string model = MagneticModel::DefaultMagneticName();
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 = \""
145  << MagneticModel::DefaultMagneticPath()
146  << "\"\nDefault magnetic name = \""
147  << MagneticModel::DefaultMagneticName()
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 }
ind
std::vector< int > ind
Definition: Slicing_stdvector_cxx11.cpp:1
H
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
Definition: gnuplot_common_settings.hh:74
GeographicLib::Math::extra_digits
static int extra_digits()
Definition: Math.hpp:187
MagneticModel.hpp
Header for GeographicLib::MagneticModel class.
D
MatrixXcd D
Definition: EigenSolver_EigenSolver_MatrixType.cpp:14
s
RealScalar s
Definition: level1_cplx_impl.h:126
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
GeographicLib::Utility::str
static std::string str(T x, int p=-1)
Definition: Utility.hpp:276
GeographicLib
Namespace for GeographicLib.
Definition: JacobiConformal.hpp:15
GeographicLib::DMS::DecodeLatLon
static void DecodeLatLon(const std::string &dmsa, const std::string &dmsb, real &lat, real &lon, bool longfirst=false)
Definition: src/DMS.cpp:252
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
GeographicLib::Math::isfinite
static bool isfinite(T x)
Definition: Math.hpp:806
GeographicLib::MagneticModel
Model of the earth's magnetic field.
Definition: MagneticModel.hpp:69
real
float real
Definition: datatypes.h:10
h
const double h
Definition: testSimpleHelicopter.cpp:19
GeographicLib::DMS::Decode
static Math::real Decode(const std::string &dms, flag &ind)
Definition: src/DMS.cpp:28
GeographicLib::GeographicErr
Exception handling for GeographicLib.
Definition: Constants.hpp:389
size
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
n
int n
Definition: BiCGSTAB_simple.cpp:1
GeographicLib::Math::real
double real
Definition: Math.hpp:129
GeographicLib::DMS::NUMBER
@ NUMBER
Definition: DMS.hpp:67
GeographicLib::MagneticModel::DefaultMagneticName
static std::string DefaultMagneticName()
Definition: src/MagneticModel.cpp:261
GeographicLib::DMS::Encode
static std::string Encode(real angle, component trailing, unsigned prec, flag ind=NONE, char dmssep=char(0))
Definition: src/DMS.cpp:299
I
#define I
Definition: main.h:112
arg
Definition: cast.h:1419
Utility.hpp
Header for GeographicLib::Utility class.
main
int main(int argc, const char *const argv[])
Definition: MagneticField.cpp:28
GeographicLib::DMS::flag
flag
Definition: DMS.hpp:41
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
GeographicLib::DMS::LONGITUDE
@ LONGITUDE
Definition: DMS.hpp:56
arg
EIGEN_DEVICE_FUNC const EIGEN_STRONG_INLINE ArgReturnType arg() const
Definition: ArrayCwiseUnaryOps.h:66
GeographicLib::MagneticModel::FieldComponents
static void FieldComponents(real Bx, real By, real Bz, real &H, real &F, real &D, real &I)
Definition: MagneticModel.hpp:202
str::str
str(const char *c, const SzType &n)
Definition: pytypes.h:1565
gtsam::symbol_shorthand::F
Key F(std::uint64_t j)
Definition: inference/Symbol.h:153
time
#define time
Definition: timeAdaptAutoDiff.cpp:31
model
noiseModel::Diagonal::shared_ptr model
Definition: doc/Code/Pose2SLAMExample.cpp:7
GeographicLib::DMS::LATITUDE
@ LATITUDE
Definition: DMS.hpp:51
str
Definition: pytypes.h:1560
MagneticCircle.hpp
Header for GeographicLib::MagneticCircle class.
GEOGRAPHICLIB_VERSION_STRING
#define GEOGRAPHICLIB_VERSION_STRING
Definition: Config.h:3
min
#define min(a, b)
Definition: datatypes.h:19
abs
#define abs(x)
Definition: datatypes.h:17
lon
static const double lon
Definition: testGeographicLib.cpp:34
max
#define max(a, b)
Definition: datatypes.h:20
gtsam.examples.ShonanAveragingCLI.str
str
Definition: ShonanAveragingCLI.py:115
real
Definition: main.h:100
GeographicLib::Utility::set_digits
static int set_digits(int ndigits=0)
Definition: Utility.cpp:48
GeographicLib::MagneticCircle
Geomagnetic field on a circle of latitude.
Definition: MagneticCircle.hpp:37
DMS.hpp
Header for GeographicLib::DMS class.
GeographicLib::MagneticModel::DefaultMagneticPath
static std::string DefaultMagneticPath()
Definition: src/MagneticModel.cpp:248
lat
static const double lat
Definition: testGeographicLib.cpp:34

gtsam
Author(s):
autogenerated on Wed May 28 2025 03:01:44