GeographicLib
1.21
|
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 }