00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #include <iostream>
00025 #include <string>
00026 #include <sstream>
00027 #include <fstream>
00028 #include <GeographicLib/MagneticModel.hpp>
00029 #include <GeographicLib/MagneticCircle.hpp>
00030 #include <GeographicLib/DMS.hpp>
00031 #include <GeographicLib/Utility.hpp>
00032
00033 #if defined(_MSC_VER)
00034
00035 # pragma warning (disable: 4127)
00036 #endif
00037
00038 #include "MagneticField.usage"
00039
00040 int main(int argc, char* argv[]) {
00041 try {
00042 using namespace GeographicLib;
00043 typedef Math::real real;
00044 Utility::set_digits();
00045 bool verbose = false;
00046 std::string dir;
00047 std::string model = MagneticModel::DefaultMagneticName();
00048 std::string istring, ifile, ofile, cdelim;
00049 char lsep = ';';
00050 real time = 0, lat = 0, h = 0;
00051 bool timeset = false, circle = false, rate = false;
00052 real hguard = 500000, tguard = 50;
00053 int prec = 1;
00054
00055 for (int m = 1; m < argc; ++m) {
00056 std::string arg(argv[m]);
00057 if (arg == "-n") {
00058 if (++m == argc) return usage(1, true);
00059 model = argv[m];
00060 } else if (arg == "-d") {
00061 if (++m == argc) return usage(1, true);
00062 dir = argv[m];
00063 } else if (arg == "-t") {
00064 if (++m == argc) return usage(1, true);
00065 try {
00066 time = Utility::fractionalyear<real>(std::string(argv[m]));
00067 timeset = true;
00068 circle = false;
00069 }
00070 catch (const std::exception& e) {
00071 std::cerr << "Error decoding argument of " << arg << ": "
00072 << e.what() << "\n";
00073 return 1;
00074 }
00075 } else if (arg == "-c") {
00076 if (m + 3 >= argc) return usage(1, true);
00077 try {
00078 using std::abs;
00079 time = Utility::fractionalyear<real>(std::string(argv[++m]));
00080 DMS::flag ind;
00081 lat = DMS::Decode(std::string(argv[++m]), ind);
00082 if (ind == DMS::LONGITUDE)
00083 throw GeographicErr("Bad hemisphere letter on latitude");
00084 if (!(abs(lat) <= 90))
00085 throw GeographicErr("Latitude not in [-90d, 90d]");
00086 h = Utility::num<real>(std::string(argv[++m]));
00087 timeset = false;
00088 circle = true;
00089 }
00090 catch (const std::exception& e) {
00091 std::cerr << "Error decoding argument of " << arg << ": "
00092 << e.what() << "\n";
00093 return 1;
00094 }
00095 } else if (arg == "-r")
00096 rate = !rate;
00097 else if (arg == "-p") {
00098 if (++m == argc) return usage(1, true);
00099 try {
00100 prec = Utility::num<int>(std::string(argv[m]));
00101 }
00102 catch (const std::exception&) {
00103 std::cerr << "Precision " << argv[m] << " is not a number\n";
00104 return 1;
00105 }
00106 } else if (arg == "-T") {
00107 if (++m == argc) return usage(1, true);
00108 try {
00109 tguard = Utility::num<real>(std::string(argv[m]));
00110 }
00111 catch (const std::exception& e) {
00112 std::cerr << "Error decoding argument of " << arg << ": "
00113 << e.what() << "\n";
00114 return 1;
00115 }
00116 } else if (arg == "-H") {
00117 if (++m == argc) return usage(1, true);
00118 try {
00119 hguard = Utility::num<real>(std::string(argv[m]));
00120 }
00121 catch (const std::exception& e) {
00122 std::cerr << "Error decoding argument of " << arg << ": "
00123 << e.what() << "\n";
00124 return 1;
00125 }
00126 } else if (arg == "-v")
00127 verbose = true;
00128 else if (arg == "--input-string") {
00129 if (++m == argc) return usage(1, true);
00130 istring = argv[m];
00131 } else if (arg == "--input-file") {
00132 if (++m == argc) return usage(1, true);
00133 ifile = argv[m];
00134 } else if (arg == "--output-file") {
00135 if (++m == argc) return usage(1, true);
00136 ofile = argv[m];
00137 } else if (arg == "--line-separator") {
00138 if (++m == argc) return usage(1, true);
00139 if (std::string(argv[m]).size() != 1) {
00140 std::cerr << "Line separator must be a single character\n";
00141 return 1;
00142 }
00143 lsep = argv[m][0];
00144 } else if (arg == "--comment-delimiter") {
00145 if (++m == argc) return usage(1, true);
00146 cdelim = argv[m];
00147 } else if (arg == "--version") {
00148 std::cout
00149 << argv[0] << ": GeographicLib version "
00150 << GEOGRAPHICLIB_VERSION_STRING << "\n";
00151 return 0;
00152 } else {
00153 int retval = usage(!(arg == "-h" || arg == "--help"), arg != "--help");
00154 if (arg == "-h")
00155 std::cout<< "\nDefault magnetic path = \""
00156 << MagneticModel::DefaultMagneticPath()
00157 << "\"\nDefault magnetic name = \""
00158 << MagneticModel::DefaultMagneticName()
00159 << "\"\n";
00160 return retval;
00161 }
00162 }
00163
00164 if (!ifile.empty() && !istring.empty()) {
00165 std::cerr << "Cannot specify --input-string and --input-file together\n";
00166 return 1;
00167 }
00168 if (ifile == "-") ifile.clear();
00169 std::ifstream infile;
00170 std::istringstream instring;
00171 if (!ifile.empty()) {
00172 infile.open(ifile.c_str());
00173 if (!infile.is_open()) {
00174 std::cerr << "Cannot open " << ifile << " for reading\n";
00175 return 1;
00176 }
00177 } else if (!istring.empty()) {
00178 std::string::size_type m = 0;
00179 while (true) {
00180 m = istring.find(lsep, m);
00181 if (m == std::string::npos)
00182 break;
00183 istring[m] = '\n';
00184 }
00185 instring.str(istring);
00186 }
00187 std::istream* input = !ifile.empty() ? &infile :
00188 (!istring.empty() ? &instring : &std::cin);
00189
00190 std::ofstream outfile;
00191 if (ofile == "-") ofile.clear();
00192 if (!ofile.empty()) {
00193 outfile.open(ofile.c_str());
00194 if (!outfile.is_open()) {
00195 std::cerr << "Cannot open " << ofile << " for writing\n";
00196 return 1;
00197 }
00198 }
00199 std::ostream* output = !ofile.empty() ? &outfile : &std::cout;
00200
00201 tguard = std::max(real(0), tguard);
00202 hguard = std::max(real(0), hguard);
00203 prec = std::min(10 + Math::extra_digits(), std::max(0, prec));
00204 int retval = 0;
00205 try {
00206 const MagneticModel m(model, dir);
00207 if ((timeset || circle)
00208 && (!Math::isfinite(time) ||
00209 time < m.MinTime() - tguard ||
00210 time > m.MaxTime() + tguard))
00211 throw GeographicErr("Time " + Utility::str(time) +
00212 " too far outside allowed range [" +
00213 Utility::str(m.MinTime()) + "," +
00214 Utility::str(m.MaxTime()) + "]");
00215 if (circle
00216 && (!Math::isfinite(h) ||
00217 h < m.MinHeight() - hguard ||
00218 h > m.MaxHeight() + hguard))
00219 throw GeographicErr("Height " + Utility::str(h/1000) +
00220 "km too far outside allowed range [" +
00221 Utility::str(m.MinHeight()/1000) + "km," +
00222 Utility::str(m.MaxHeight()/1000) + "km]");
00223 if (verbose) {
00224 std::cerr << "Magnetic file: " << m.MagneticFile() << "\n"
00225 << "Name: " << m.MagneticModelName() << "\n"
00226 << "Description: " << m.Description() << "\n"
00227 << "Date & Time: " << m.DateTime() << "\n"
00228 << "Time range: ["
00229 << m.MinTime() << ","
00230 << m.MaxTime() << "]\n"
00231 << "Height range: ["
00232 << m.MinHeight()/1000 << "km,"
00233 << m.MaxHeight()/1000 << "km]\n";
00234 }
00235 if ((timeset || circle) && (time < m.MinTime() || time > m.MaxTime()))
00236 std::cerr << "WARNING: Time " << time
00237 << " outside allowed range ["
00238 << m.MinTime() << "," << m.MaxTime() << "]\n";
00239 if (circle && (h < m.MinHeight() || h > m.MaxHeight()))
00240 std::cerr << "WARNING: Height " << h/1000
00241 << "km outside allowed range ["
00242 << m.MinHeight()/1000 << "km,"
00243 << m.MaxHeight()/1000 << "km]\n";
00244 const MagneticCircle c(circle ? m.Circle(time, lat, h) :
00245 MagneticCircle());
00246 std::string s, stra, strb;
00247 while (std::getline(*input, s)) {
00248 try {
00249 std::string eol("\n");
00250 if (!cdelim.empty()) {
00251 std::string::size_type m = s.find(cdelim);
00252 if (m != std::string::npos) {
00253 eol = " " + s.substr(m) + "\n";
00254 s = s.substr(0, m);
00255 }
00256 }
00257 std::istringstream str(s);
00258 if (!(timeset || circle)) {
00259 if (!(str >> stra))
00260 throw GeographicErr("Incomplete input: " + s);
00261 time = Utility::fractionalyear<real>(stra);
00262 if (time < m.MinTime() - tguard || time > m.MaxTime() + tguard)
00263 throw GeographicErr("Time " + Utility::str(time) +
00264 " too far outside allowed range [" +
00265 Utility::str(m.MinTime()) + "," +
00266 Utility::str(m.MaxTime()) +
00267 "]");
00268 if (time < m.MinTime() || time > m.MaxTime())
00269 std::cerr << "WARNING: Time " << time
00270 << " outside allowed range ["
00271 << m.MinTime() << "," << m.MaxTime() << "]\n";
00272 }
00273 real lon;
00274 if (circle) {
00275 if (!(str >> strb))
00276 throw GeographicErr("Incomplete input: " + s);
00277 DMS::flag ind;
00278 lon = DMS::Decode(strb, ind);
00279 if (ind == DMS::LATITUDE)
00280 throw GeographicErr("Bad hemisphere letter on " + strb);
00281 if (lon < -540 || lon >= 540)
00282 throw GeographicErr("Longitude " + strb +
00283 " not in [-540d, 540d)");
00284 } else {
00285 if (!(str >> stra >> strb))
00286 throw GeographicErr("Incomplete input: " + s);
00287 DMS::DecodeLatLon(stra, strb, lat, lon);
00288 h = 0;
00289 if (str >> h) {
00290 if (h < m.MinHeight() - hguard || h > m.MaxHeight() + hguard)
00291 throw GeographicErr("Height " + Utility::str(h/1000) +
00292 "km too far outside allowed range [" +
00293 Utility::str(m.MinHeight()/1000) + "km," +
00294 Utility::str(m.MaxHeight()/1000) + "km]");
00295 if (h < m.MinHeight() || h > m.MaxHeight())
00296 std::cerr << "WARNING: Height " << h/1000
00297 << "km outside allowed range ["
00298 << m.MinHeight()/1000 << "km,"
00299 << m.MaxHeight()/1000 << "km]\n";
00300 }
00301 else
00302 str.clear();
00303 }
00304 if (str >> stra)
00305 throw GeographicErr("Extra junk in input: " + s);
00306 real bx, by, bz, bxt, byt, bzt;
00307 if (circle)
00308 c(lon, bx, by, bz, bxt, byt, bzt);
00309 else
00310 m(time, lat, lon, h, bx, by, bz, bxt, byt, bzt);
00311 real H, F, D, I, Ht, Ft, Dt, It;
00312 MagneticModel::FieldComponents(bx, by, bz, bxt, byt, bzt,
00313 H, F, D, I, Ht, Ft, Dt, It);
00314
00315 *output << DMS::Encode(D, prec + 1, DMS::NUMBER) << " "
00316 << DMS::Encode(I, prec + 1, DMS::NUMBER) << " "
00317 << Utility::str(H, prec) << " "
00318 << Utility::str(by, prec) << " "
00319 << Utility::str(bx, prec) << " "
00320 << Utility::str(-bz, prec) << " "
00321 << Utility::str(F, prec) << eol;
00322 if (rate)
00323 *output << DMS::Encode(Dt, prec + 1, DMS::NUMBER) << " "
00324 << DMS::Encode(It, prec + 1, DMS::NUMBER) << " "
00325 << Utility::str(Ht, prec) << " "
00326 << Utility::str(byt, prec) << " "
00327 << Utility::str(bxt, prec) << " "
00328 << Utility::str(-bzt, prec) << " "
00329 << Utility::str(Ft, prec) << eol;
00330 }
00331 catch (const std::exception& e) {
00332 *output << "ERROR: " << e.what() << "\n";
00333 retval = 1;
00334 }
00335 }
00336 }
00337 catch (const std::exception& e) {
00338 std::cerr << "Error reading " << model << ": " << e.what() << "\n";
00339 retval = 1;
00340 }
00341 return retval;
00342 }
00343 catch (const std::exception& e) {
00344 std::cerr << "Caught exception: " << e.what() << "\n";
00345 return 1;
00346 }
00347 catch (...) {
00348 std::cerr << "Caught unknown exception\n";
00349 return 1;
00350 }
00351 }