GeographicLib
1.21
|
00001 /** 00002 * \file CassiniSoldner.cpp 00003 * \brief Implementation for GeographicLib::CassiniSoldner class 00004 * 00005 * Copyright (c) Charles Karney (2009-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/CassiniSoldner.hpp> 00011 00012 #define GEOGRAPHICLIB_CASSINISOLDNER_CPP \ 00013 "$Id: 2823df38f3e31fd8b882e2f55ad72d6419b03246 $" 00014 00015 RCSID_DECL(GEOGRAPHICLIB_CASSINISOLDNER_CPP) 00016 RCSID_DECL(GEOGRAPHICLIB_CASSINISOLDNER_HPP) 00017 00018 namespace GeographicLib { 00019 00020 using namespace std; 00021 00022 const Math::real CassiniSoldner::eps1_ = 00023 real(0.01) * sqrt(numeric_limits<real>::epsilon()); 00024 const Math::real CassiniSoldner::tiny_ = sqrt(numeric_limits<real>::min()); 00025 00026 void CassiniSoldner::Reset(real lat0, real lon0) throw() { 00027 _meridian = _earth.Line(lat0, lon0, real(0), 00028 Geodesic::LATITUDE | Geodesic::LONGITUDE | 00029 Geodesic::DISTANCE | Geodesic::DISTANCE_IN | 00030 Geodesic::AZIMUTH); 00031 real 00032 phi = LatitudeOrigin() * Math::degree<real>(), 00033 f = _earth.Flattening(); 00034 _sbet0 = (1 - f) * sin(phi); 00035 _cbet0 = abs(LatitudeOrigin()) == 90 ? 0 : cos(phi); 00036 SinCosNorm(_sbet0, _cbet0); 00037 } 00038 00039 void CassiniSoldner::Forward(real lat, real lon, real& x, real& y, 00040 real& azi, real& rk) const throw() { 00041 if (!Init()) 00042 return; 00043 real dlon = AngNormalize(lon - LongitudeOrigin()); 00044 real sig12, s12, azi1, azi2; 00045 lat = AngRound(lat); 00046 sig12 = _earth.Inverse(lat, -abs(dlon), lat, abs(dlon), s12, azi1, azi2); 00047 if (sig12 < 100 * tiny_) 00048 sig12 = s12 = 0; 00049 sig12 *= real(0.5); 00050 s12 *= real(0.5); 00051 if (s12 == 0) { 00052 real da = (azi2 - azi1)/2; 00053 if (abs(dlon) <= 90) { 00054 azi1 = 90 - da; 00055 azi2 = 90 + da; 00056 } else { 00057 azi1 = -90 - da; 00058 azi2 = -90 + da; 00059 } 00060 } 00061 if (dlon < 0) { 00062 azi2 = azi1; 00063 s12 = -s12; 00064 sig12 = -sig12; 00065 } 00066 x = s12; 00067 azi = AngNormalize(azi2); 00068 GeodesicLine perp(_earth.Line(lat, dlon, azi2, Geodesic::GEODESICSCALE)); 00069 real t; 00070 perp.GenPosition(true, -sig12, 00071 Geodesic::GEODESICSCALE, 00072 t, t, t, t, t, t, rk, t); 00073 00074 real 00075 alp0 = perp.EquatorialAzimuth() * Math::degree<real>(), 00076 calp0 = cos(alp0), salp0 = sin(alp0), 00077 sbet1 = lat >=0 ? calp0 : -calp0, 00078 cbet1 = abs(dlon) <= 90 ? abs(salp0) : -abs(salp0), 00079 sbet01 = sbet1 * _cbet0 - cbet1 * _sbet0, 00080 cbet01 = cbet1 * _cbet0 + sbet1 * _sbet0, 00081 sig01 = atan2(sbet01, cbet01) / Math::degree<real>(); 00082 _meridian.GenPosition(true, sig01, 00083 Geodesic::DISTANCE, 00084 t, t, t, y, t, t, t, t); 00085 } 00086 00087 void CassiniSoldner::Reverse(real x, real y, real& lat, real& lon, 00088 real& azi, real& rk) const throw() { 00089 if (!Init()) 00090 return; 00091 real lat1, lon1; 00092 real azi0, t; 00093 _meridian.Position(y, lat1, lon1, azi0); 00094 _earth.Direct(lat1, lon1, azi0 + 90, x, lat, lon, azi, rk, t); 00095 } 00096 00097 } // namespace GeographicLib