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