GeographicLib
1.21
|
00001 /** 00002 * \file Gravity.cpp 00003 * \brief Command line utility for evaluating gravity 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 Gravity \ 00011 * Gravity.cpp \ 00012 * ../src/CircularEngine.cpp \ 00013 * ../src/DMS.cpp \ 00014 * ../src/Geocentric.cpp \ 00015 * ../src/GravityCircle.cpp \ 00016 * ../src/GravityModel.cpp \ 00017 * ../src/NormalGravity.cpp \ 00018 * ../src/SphericalEngine.cpp \ 00019 * ../src/Utility.cpp 00020 * 00021 * See the <a href="Gravity.1.html">man page</a> for usage 00022 * information. 00023 **********************************************************************/ 00024 00025 #include <iostream> 00026 #include <string> 00027 #include <sstream> 00028 #include <fstream> 00029 #include <GeographicLib/GravityModel.hpp> 00030 #include <GeographicLib/GravityCircle.hpp> 00031 #include <GeographicLib/DMS.hpp> 00032 #include <GeographicLib/Utility.hpp> 00033 00034 #include "Gravity.usage" 00035 00036 int main(int argc, char* argv[]) { 00037 try { 00038 using namespace GeographicLib; 00039 typedef Math::real real; 00040 bool verbose = false; 00041 std::string dir; 00042 std::string model = GravityModel::DefaultGravityName(); 00043 std::string istring, ifile, ofile, cdelim; 00044 char lsep = ';'; 00045 real lat = 0, h = 0; 00046 bool circle = false; 00047 int prec = -1; 00048 enum { 00049 GRAVITY = 0, 00050 DISTURBANCE = 1, 00051 ANOMALY = 2, 00052 UNDULATION = 3, 00053 }; 00054 unsigned mode = GRAVITY; 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 == "-G") 00064 mode = GRAVITY; 00065 else if (arg == "-D") 00066 mode = DISTURBANCE; 00067 else if (arg == "-A") 00068 mode = ANOMALY; 00069 else if (arg == "-H") 00070 mode = UNDULATION; 00071 else if (arg == "-c") { 00072 if (m + 2 >= argc) return usage(1, true); 00073 try { 00074 DMS::flag ind; 00075 lat = DMS::Decode(std::string(argv[++m]), ind); 00076 if (ind == DMS::LONGITUDE) 00077 throw GeographicErr("Bad hemisphere letter on latitude"); 00078 if (!(std::abs(lat) <= 90)) 00079 throw GeographicErr("Latitude not in [-90d, 90d]"); 00080 h = Utility::num<real>(std::string(argv[++m])); 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 == "-p") { 00089 if (++m == argc) return usage(1, true); 00090 try { 00091 prec = Utility::num<int>(std::string(argv[m])); 00092 } 00093 catch (const std::exception&) { 00094 std::cerr << "Precision " << argv[m] << " is not a number\n"; 00095 return 1; 00096 } 00097 } else if (arg == "-v") 00098 verbose = true; 00099 else if (arg == "--input-string") { 00100 if (++m == argc) return usage(1, true); 00101 istring = argv[m]; 00102 } else if (arg == "--input-file") { 00103 if (++m == argc) return usage(1, true); 00104 ifile = argv[m]; 00105 } else if (arg == "--output-file") { 00106 if (++m == argc) return usage(1, true); 00107 ofile = argv[m]; 00108 } else if (arg == "--line-separator") { 00109 if (++m == argc) return usage(1, true); 00110 if (std::string(argv[m]).size() != 1) { 00111 std::cerr << "Line separator must be a single character\n"; 00112 return 1; 00113 } 00114 lsep = argv[m][0]; 00115 } else if (arg == "--comment-delimiter") { 00116 if (++m == argc) return usage(1, true); 00117 cdelim = argv[m]; 00118 } else if (arg == "--version") { 00119 std::cout 00120 << argv[0] 00121 << ": $Id: c87c647c3e973929010cdb2fd5d1eaa6aa739eca $\n" 00122 << "GeographicLib version " << GEOGRAPHICLIB_VERSION_STRING << "\n"; 00123 return 0; 00124 } else { 00125 int retval = usage(!(arg == "-h" || arg == "--help"), arg != "--help"); 00126 if (arg == "-h") 00127 std::cout<< "\nDefault gravity path = \"" 00128 << GravityModel::DefaultGravityPath() 00129 << "\"\nDefault gravity name = \"" 00130 << GravityModel::DefaultGravityName() 00131 << "\"\n"; 00132 return retval; 00133 } 00134 } 00135 00136 if (!ifile.empty() && !istring.empty()) { 00137 std::cerr << "Cannot specify --input-string and --input-file together\n"; 00138 return 1; 00139 } 00140 if (ifile == "-") ifile.clear(); 00141 std::ifstream infile; 00142 std::istringstream instring; 00143 if (!ifile.empty()) { 00144 infile.open(ifile.c_str()); 00145 if (!infile.is_open()) { 00146 std::cerr << "Cannot open " << ifile << " for reading\n"; 00147 return 1; 00148 } 00149 } else if (!istring.empty()) { 00150 std::string::size_type m = 0; 00151 while (true) { 00152 m = istring.find(lsep, m); 00153 if (m == std::string::npos) 00154 break; 00155 istring[m] = '\n'; 00156 } 00157 instring.str(istring); 00158 } 00159 std::istream* input = !ifile.empty() ? &infile : 00160 (!istring.empty() ? &instring : &std::cin); 00161 00162 std::ofstream outfile; 00163 if (ofile == "-") ofile.clear(); 00164 if (!ofile.empty()) { 00165 outfile.open(ofile.c_str()); 00166 if (!outfile.is_open()) { 00167 std::cerr << "Cannot open " << ofile << " for writing\n"; 00168 return 1; 00169 } 00170 } 00171 std::ostream* output = !ofile.empty() ? &outfile : &std::cout; 00172 00173 switch (mode) { 00174 case GRAVITY: 00175 prec = std::min(16, prec < 0 ? 5 : prec); 00176 break; 00177 case DISTURBANCE: 00178 case ANOMALY: 00179 prec = std::min(14, prec < 0 ? 3 : prec); 00180 break; 00181 case UNDULATION: 00182 default: 00183 prec = std::min(12, prec < 0 ? 4 : prec); 00184 break; 00185 } 00186 int retval = 0; 00187 try { 00188 const GravityModel g(model, dir); 00189 if (circle) { 00190 if (!Math::isfinite<real>(h)) 00191 throw GeographicErr("Bad height"); 00192 else if (mode == UNDULATION && h != 0) 00193 throw GeographicErr("Height should be zero for geoid undulations"); 00194 } 00195 if (verbose) { 00196 std::cerr << "Gravity file: " << g.GravityFile() << "\n" 00197 << "Name: " << g.GravityModelName() << "\n" 00198 << "Description: " << g.Description() << "\n" 00199 << "Date & Time: " << g.DateTime() << "\n"; 00200 } 00201 unsigned mask = (mode == GRAVITY ? GravityModel::GRAVITY : 00202 (mode == DISTURBANCE ? GravityModel::DISTURBANCE : 00203 (mode == ANOMALY ? GravityModel::SPHERICAL_ANOMALY : 00204 GravityModel::GEOID_HEIGHT))); // mode == UNDULATION 00205 const GravityCircle c(circle ? g.Circle(lat, h, mask) : GravityCircle()); 00206 std::string s, stra, strb; 00207 while (std::getline(*input, s)) { 00208 try { 00209 std::string eol("\n"); 00210 if (!cdelim.empty()) { 00211 std::string::size_type m = s.find(cdelim); 00212 if (m != std::string::npos) { 00213 eol = " " + s.substr(m) + "\n"; 00214 s = s.substr(0, m); 00215 } 00216 } 00217 std::istringstream str(s); 00218 real lon; 00219 if (circle) { 00220 if (!(str >> strb)) 00221 throw GeographicErr("Incomplete input: " + s); 00222 DMS::flag ind; 00223 lon = DMS::Decode(strb, ind); 00224 if (ind == DMS::LATITUDE) 00225 throw GeographicErr("Bad hemisphere letter on " + strb); 00226 if (lon < -180 || lon > 360) 00227 throw GeographicErr("Longitude " + strb + "not in [-180d, 360d]"); 00228 } else { 00229 if (!(str >> stra >> strb)) 00230 throw GeographicErr("Incomplete input: " + s); 00231 DMS::DecodeLatLon(stra, strb, lat, lon); 00232 h = 0; 00233 if (!(str >> h)) // h is optional 00234 str.clear(); 00235 if (mode == UNDULATION && h != 0) 00236 throw GeographicErr("Height must be zero for geoid heights"); 00237 } 00238 if (str >> stra) 00239 throw GeographicErr("Extra junk in input: " + s); 00240 switch (mode) { 00241 case GRAVITY: 00242 { 00243 real gx, gy, gz; 00244 if (circle) { 00245 c.Gravity(lon, gx, gy, gz); 00246 } else { 00247 g.Gravity(lat, lon, h, gx, gy, gz); 00248 } 00249 *output << Utility::str<real>(gx, prec) << " " 00250 << Utility::str<real>(gy, prec) << " " 00251 << Utility::str<real>(gz, prec) << eol; 00252 } 00253 break; 00254 case DISTURBANCE: 00255 { 00256 real deltax, deltay, deltaz; 00257 if (circle) { 00258 c.Disturbance(lon, deltax, deltay, deltaz); 00259 } else { 00260 g.Disturbance(lat, lon, h, deltax, deltay, deltaz); 00261 } 00262 // Convert to mGals 00263 *output << Utility::str<real>(deltax * 1e5, prec) << " " 00264 << Utility::str<real>(deltay * 1e5, prec) << " " 00265 << Utility::str<real>(deltaz * 1e5, prec) 00266 << eol; 00267 } 00268 break; 00269 case ANOMALY: 00270 { 00271 real Dg01, xi, eta; 00272 if (circle) 00273 c.SphericalAnomaly(lon, Dg01, xi, eta); 00274 else 00275 g.SphericalAnomaly(lat, lon, h, Dg01, xi, eta); 00276 Dg01 *= 1e5; // Convert to mGals 00277 xi *= 3600; // Convert to arcsecs 00278 eta *= 3600; 00279 *output << Utility::str<real>(Dg01, prec) << " " 00280 << Utility::str<real>(xi, prec) << " " 00281 << Utility::str<real>(eta, prec) << eol; 00282 } 00283 break; 00284 case UNDULATION: 00285 default: 00286 { 00287 real N = circle ? c.GeoidHeight(lon) : g.GeoidHeight(lat, lon); 00288 *output << Utility::str<real>(N, prec) << eol; 00289 } 00290 break; 00291 } 00292 } 00293 catch (const std::exception& e) { 00294 *output << "ERROR: " << e.what() << "\n"; 00295 retval = 1; 00296 } 00297 } 00298 } 00299 catch (const std::exception& e) { 00300 std::cerr << "Error reading " << model << ": " << e.what() << "\n"; 00301 retval = 1; 00302 } 00303 return retval; 00304 } 00305 catch (const std::exception& e) { 00306 std::cerr << "Caught exception: " << e.what() << "\n"; 00307 return 1; 00308 } 00309 catch (...) { 00310 std::cerr << "Caught unknown exception\n"; 00311 return 1; 00312 } 00313 }