1 /** 2 * \file CassiniSoldner.cpp 3 * \brief Implementation for GeographicLib::CassiniSoldner class 4 * 5 * Copyright (c) Charles Karney (2009-2015) <charles@karney.com> and licensed 6 * under the MIT/X11 License. For more information, see 7 * https://geographiclib.sourceforge.io/ 8 **********************************************************************/ 9 10 #include <GeographicLib/CassiniSoldner.hpp> 11 12 namespace GeographicLib { 13 14 using namespace std; 15 CassiniSoldner(const Geodesic & earth)16 CassiniSoldner::CassiniSoldner(const Geodesic& earth) 17 : _earth(earth) {} 18 CassiniSoldner(real lat0,real lon0,const Geodesic & earth)19 CassiniSoldner::CassiniSoldner(real lat0, real lon0, const Geodesic& earth) 20 : _earth(earth) 21 { Reset(lat0, lon0); } 22 Reset(real lat0,real lon0)23 void CassiniSoldner::Reset(real lat0, real lon0) { 24 _meridian = _earth.Line(lat0, lon0, real(0), 25 Geodesic::LATITUDE | Geodesic::LONGITUDE | 26 Geodesic::DISTANCE | Geodesic::DISTANCE_IN | 27 Geodesic::AZIMUTH); 28 real f = _earth.Flattening(); 29 Math::sincosd(LatitudeOrigin(), _sbet0, _cbet0); 30 _sbet0 *= (1 - f); 31 Math::norm(_sbet0, _cbet0); 32 } 33 Forward(real lat,real lon,real & x,real & y,real & azi,real & rk) const34 void CassiniSoldner::Forward(real lat, real lon, real& x, real& y, 35 real& azi, real& rk) const { 36 if (!Init()) 37 return; 38 real dlon = Math::AngDiff(LongitudeOrigin(), lon); 39 real sig12, s12, azi1, azi2; 40 sig12 = _earth.Inverse(lat, -abs(dlon), lat, abs(dlon), s12, azi1, azi2); 41 sig12 *= real(0.5); 42 s12 *= real(0.5); 43 if (s12 == 0) { 44 real da = Math::AngDiff(azi1, azi2)/2; 45 if (abs(dlon) <= 90) { 46 azi1 = 90 - da; 47 azi2 = 90 + da; 48 } else { 49 azi1 = -90 - da; 50 azi2 = -90 + da; 51 } 52 } 53 if (dlon < 0) { 54 azi2 = azi1; 55 s12 = -s12; 56 sig12 = -sig12; 57 } 58 x = s12; 59 azi = Math::AngNormalize(azi2); 60 GeodesicLine perp(_earth.Line(lat, dlon, azi, Geodesic::GEODESICSCALE)); 61 real t; 62 perp.GenPosition(true, -sig12, 63 Geodesic::GEODESICSCALE, 64 t, t, t, t, t, t, rk, t); 65 66 real salp0, calp0; 67 Math::sincosd(perp.EquatorialAzimuth(), salp0, calp0); 68 real 69 sbet1 = lat >=0 ? calp0 : -calp0, 70 cbet1 = abs(dlon) <= 90 ? abs(salp0) : -abs(salp0), 71 sbet01 = sbet1 * _cbet0 - cbet1 * _sbet0, 72 cbet01 = cbet1 * _cbet0 + sbet1 * _sbet0, 73 sig01 = atan2(sbet01, cbet01) / Math::degree(); 74 _meridian.GenPosition(true, sig01, 75 Geodesic::DISTANCE, 76 t, t, t, y, t, t, t, t); 77 } 78 Reverse(real x,real y,real & lat,real & lon,real & azi,real & rk) const79 void CassiniSoldner::Reverse(real x, real y, real& lat, real& lon, 80 real& azi, real& rk) const { 81 if (!Init()) 82 return; 83 real lat1, lon1; 84 real azi0, t; 85 _meridian.Position(y, lat1, lon1, azi0); 86 _earth.Direct(lat1, lon1, azi0 + 90, x, lat, lon, azi, rk, t); 87 } 88 89 } // namespace GeographicLib 90