GeographicLib  1.21
Gravity.cpp
Go to the documentation of this file.
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 }