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