00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "GeographicLib/CassiniSoldner.hpp"
00012
00013 #define GEOGRAPHICLIB_CASSINISOLDNER_CPP "$Id: CassiniSoldner.cpp 6921 2010-12-31 14:34:50Z 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 Geodesic::LATITUDE | Geodesic::LONGITUDE |
00029 Geodesic::DISTANCE | Geodesic::DISTANCE_IN |
00030 Geodesic::AZIMUTH);
00031 real
00032 phi = LatitudeOrigin() * Math::degree<real>(),
00033 f = _earth.InverseFlattening() != 0 ? 1 / _earth.InverseFlattening() : 0;
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 * eps2)
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 }