Gnomonic.cpp

Go to the documentation of this file.
00001 /**
00002  * \file Gnomonic.cpp
00003  * \brief Implementation for GeographicLib::Gnomonic class
00004  *
00005  * Copyright (c) Charles Karney (2010) <charles@karney.com> and licensed under
00006  * the LGPL.  For more information, see http://geographiclib.sourceforge.net/
00007  **********************************************************************/
00008 
00009 #include "GeographicLib/Gnomonic.hpp"
00010 
00011 #define GEOGRAPHICLIB_GNOMONIC_CPP "$Id: Gnomonic.cpp 6921 2010-12-31 14:34:50Z karney $"
00012 
00013 RCSID_DECL(GEOGRAPHICLIB_GNOMONIC_CPP)
00014 RCSID_DECL(GEOGRAPHICLIB_GNOMONIC_HPP)
00015 
00016 namespace GeographicLib {
00017 
00018   using namespace std;
00019 
00020   const Math::real Gnomonic::eps0 = numeric_limits<real>::epsilon();
00021   const Math::real Gnomonic::eps = real(0.01) * sqrt(eps0);
00022 
00023   void Gnomonic::Forward(real lat0, real lon0, real lat, real lon,
00024                          real& x, real& y, real& azi, real& rk)
00025     const throw() {
00026     real azi0, m, M, t;
00027     _earth.GenInverse(lat0, lon0, lat, lon,
00028                       Geodesic::AZIMUTH | Geodesic::REDUCEDLENGTH |
00029                       Geodesic::GEODESICSCALE,
00030                       t, azi0, azi, m, M, t, t);
00031     rk = M;
00032     if (M <= 0)
00033       x = y = Math::NaN();
00034     else {
00035       real rho = m/M;
00036       azi0 *= Math::degree<real>();
00037       x = rho * sin(azi0);
00038       y = rho * cos(azi0);
00039     }
00040   }
00041 
00042   void Gnomonic::Reverse(real lat0, real lon0, real x, real y,
00043                          real& lat, real& lon, real& azi, real& rk)
00044     const throw() {
00045     real
00046       azi0 = atan2(x, y) / Math::degree<real>(),
00047       rho = Math::hypot(x, y),
00048       s = _a * atan(rho/_a);
00049     bool little = rho <= _a;
00050     if (!little)
00051       rho = 1/rho;
00052     GeodesicLine line(_earth.Line(lat0, lon0, azi0,
00053                                   Geodesic::LATITUDE | Geodesic::LONGITUDE |
00054                                   Geodesic::AZIMUTH | Geodesic::DISTANCE_IN |
00055                                   Geodesic::REDUCEDLENGTH |
00056                                   Geodesic::GEODESICSCALE));
00057     int count = numit, trip = 0;
00058     real lat1, lon1, azi1, M;
00059     while (count--) {
00060       real m, t;
00061       line.Position(s, lat1, lon1, azi1, m, M, t);
00062       if (trip)
00063         break;
00064       // If little, solve rho(s) = rho with drho(s)/ds = 1/M^2
00065       // else solve 1/rho(s) = 1/rho with d(1/rho(s))/ds = -1/m^2
00066       real ds = little ? (m/M - rho) * M * M : (rho - M/m) * m * m;
00067       s -= ds;
00068       if (!(abs(ds) >= eps * _a))
00069         ++trip;
00070     }
00071     if (trip) {
00072       lat = lat1; lon = lon1; azi = azi1; rk = M;
00073     } else
00074       lat = lon = azi = rk = Math::NaN();
00075     return;
00076   }
00077 } // namespace GeographicLib