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