CassiniSoldner.cpp

Go to the documentation of this file.
00001 
00002 /**
00003  * \file CassiniSoldner.cpp
00004  * \brief Implementation for GeographicLib::CassiniSoldner class
00005  *
00006  * Copyright (c) Charles Karney (2009, 2010) <charles@karney.com>
00007  * and licensed under the LGPL.  For more information, see
00008  * http://geographiclib.sourceforge.net/
00009  **********************************************************************/
00010 
00011 #include "GeographicLib/CassiniSoldner.hpp"
00012 
00013 #define GEOGRAPHICLIB_CASSINISOLDNER_CPP "$Id: CassiniSoldner.cpp 6827 2010-05-20 19:56:18Z karney $"
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::eps2 = sqrt(numeric_limits<real>::min());
00025 
00026   void CassiniSoldner::Reset(real lat0, real lon0) throw() {
00027     _meridian = _earth.Line(lat0, lon0, real(0));
00028     real
00029       phi = LatitudeOrigin() * Constants::degree(),
00030       f = _earth.InverseFlattening() != 0 ? 1 / _earth.InverseFlattening() : 0;
00031     _sbet0 = (1 - f) * sin(phi);
00032     _cbet0 = abs(LatitudeOrigin()) == 90 ? 0 : cos(phi);
00033     SinCosNorm(_sbet0, _cbet0);
00034   }
00035 
00036   void CassiniSoldner::Forward(real lat, real lon, real& x, real& y,
00037                                real& azi, real& m) const throw() {
00038     if (!Init())
00039       return;
00040     real dlon = AngNormalize(lon - LongitudeOrigin());
00041     real sig12, s12, azi1, azi2, m12;
00042     lat = AngRound(lat);
00043     sig12 = _earth.Inverse(lat, -abs(dlon), lat, abs(dlon),
00044                            s12, azi1, azi2, m12);
00045     if (sig12 < 100 * eps2)
00046       sig12 = s12 = 0;
00047     sig12 *= real(0.5);
00048     s12 *= real(0.5);
00049     if (s12 == 0) {
00050       real da = (azi2 - azi1)/2;
00051       if (abs(dlon) <= 90) {
00052         azi1 = 90 - da;
00053         azi2 = 90 + da;
00054       } else {
00055         azi1 = -90 - da;
00056         azi2 = -90 + da;
00057       }
00058     }
00059     if (dlon < 0) {
00060       azi2 = azi1;
00061       s12 = -s12;
00062       sig12 = -sig12;
00063     }
00064     x = s12;
00065     azi = AngNormalize(azi2);
00066     GeodesicLine perp = _earth.Line(lat, dlon, azi2);
00067     real M12;
00068     perp.Scale(-sig12, M12, m);
00069 
00070     real
00071       alp0 = perp.EquatorialAzimuth() * Constants::degree(),
00072       calp0 = cos(alp0), salp0 = sin(alp0),
00073       sbet1 = lat >=0 ? calp0 : -calp0,
00074       cbet1 = abs(dlon) <= 90 ? abs(salp0) : -abs(salp0),
00075       sbet01 = sbet1 * _cbet0 - cbet1 * _sbet0,
00076       cbet01 = cbet1 * _cbet0 + sbet1 * _sbet0,
00077       sig01 = atan2(sbet01, cbet01) / Constants::degree();
00078     real latx, lonx, azix, m12x;
00079     y = _meridian.Position(sig01, latx, lonx, azix, m12x, true);
00080   }
00081 
00082   void CassiniSoldner::Reverse(real x, real y, real& lat, real& lon,
00083                                real& azi, real& m) const throw() {
00084     if (!Init())
00085       return;
00086     real lat1, lon1;
00087     real azi0, m0;
00088     _meridian.Position(y, lat1, lon1, azi0, m0);
00089     GeodesicLine perp = _earth.Line(lat1, lon1, azi0 + 90);
00090     real sig12 = perp.Position(x, lat, lon, azi, m0);
00091     real M21;
00092     perp.Scale(sig12, m, M21);
00093   }
00094 
00095 } // namespace GeographicLib

Generated on 21 May 2010 for GeographicLib by  doxygen 1.6.1