00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include <GeographicLib/Gnomonic.hpp>
00011
00012 #if defined(_MSC_VER)
00013
00014 # pragma warning (disable: 4701)
00015 #endif
00016
00017 namespace GeographicLib {
00018
00019 using namespace std;
00020
00021 Gnomonic::Gnomonic(const Geodesic& earth)
00022 : eps0_(numeric_limits<real>::epsilon())
00023 , eps_(real(0.01) * sqrt(eps0_))
00024 , _earth(earth)
00025 , _a(_earth.MajorRadius())
00026 , _f(_earth.Flattening())
00027 {}
00028
00029 void Gnomonic::Forward(real lat0, real lon0, real lat, real lon,
00030 real& x, real& y, real& azi, real& rk)
00031 const {
00032 real azi0, m, M, t;
00033 _earth.GenInverse(lat0, lon0, lat, lon,
00034 Geodesic::AZIMUTH | Geodesic::REDUCEDLENGTH |
00035 Geodesic::GEODESICSCALE,
00036 t, azi0, azi, m, M, t, t);
00037 rk = M;
00038 if (M <= 0)
00039 x = y = Math::NaN();
00040 else {
00041 real rho = m/M;
00042 azi0 *= Math::degree();
00043 x = rho * sin(azi0);
00044 y = rho * cos(azi0);
00045 }
00046 }
00047
00048 void Gnomonic::Reverse(real lat0, real lon0, real x, real y,
00049 real& lat, real& lon, real& azi, real& rk)
00050 const {
00051 real
00052 azi0 = atan2(x, y) / Math::degree(),
00053 rho = Math::hypot(x, y),
00054 s = _a * atan(rho/_a);
00055 bool little = rho <= _a;
00056 if (!little)
00057 rho = 1/rho;
00058 GeodesicLine line(_earth.Line(lat0, lon0, azi0,
00059 Geodesic::LATITUDE | Geodesic::LONGITUDE |
00060 Geodesic::AZIMUTH | Geodesic::DISTANCE_IN |
00061 Geodesic::REDUCEDLENGTH |
00062 Geodesic::GEODESICSCALE));
00063 int count = numit_, trip = 0;
00064 real lat1, lon1, azi1, M;
00065 while (count--) {
00066 real m, t;
00067 line.Position(s, lat1, lon1, azi1, m, M, t);
00068 if (trip)
00069 break;
00070
00071
00072 real ds = little ? (m/M - rho) * M * M : (rho - M/m) * m * m;
00073 s -= ds;
00074
00075 if (!(abs(ds) >= eps_ * _a))
00076 ++trip;
00077 }
00078 if (trip) {
00079 lat = lat1; lon = lon1; azi = azi1; rk = M;
00080 } else
00081 lat = lon = azi = rk = Math::NaN();
00082 return;
00083 }
00084
00085 }