GeographicLib  1.21
CassiniSoldner.cpp
Go to the documentation of this file.
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