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