GeographicLib
1.21
|
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