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