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