GeographicLib  1.43
CassiniSoldner.cpp
Go to the documentation of this file.
1 /**
2  * \file CassiniSoldner.cpp
3  * \brief Implementation for GeographicLib::CassiniSoldner class
4  *
7  * http://geographiclib.sourceforge.net/
8  **********************************************************************/
9
11
12 namespace GeographicLib {
13
14  using namespace std;
15
17  : eps1_(real(0.01) * sqrt(numeric_limits<real>::epsilon()))
18  , tiny_(sqrt(numeric_limits<real>::min()))
19  , _earth(earth) {}
20
21  CassiniSoldner::CassiniSoldner(real lat0, real lon0, const Geodesic& earth)
22  : eps1_(real(0.01) * sqrt(numeric_limits<real>::epsilon()))
23  , tiny_(sqrt(numeric_limits<real>::min()))
24  , _earth(earth)
25  { Reset(lat0, lon0); }
26
27  void CassiniSoldner::Reset(real lat0, real lon0) {
28  _meridian = _earth.Line(lat0, lon0, real(0),
32  real
33  phi = LatitudeOrigin() * Math::degree(),
34  f = _earth.Flattening();
35  _sbet0 = (1 - f) * sin(phi);
36  _cbet0 = abs(LatitudeOrigin()) == 90 ? 0 : cos(phi);
37  Math::norm(_sbet0, _cbet0);
38  }
39
40  void CassiniSoldner::Forward(real lat, real lon, real& x, real& y,
41  real& azi, real& rk) const {
42  if (!Init())
43  return;
45  real sig12, s12, azi1, azi2;
46  lat = Math::AngRound(lat);
47  sig12 = _earth.Inverse(lat, -abs(dlon), lat, abs(dlon), s12, azi1, azi2);
48  if (sig12 < 100 * tiny_)
49  sig12 = s12 = 0;
50  sig12 *= real(0.5);
51  s12 *= real(0.5);
52  if (s12 == 0) {
53  real da = (azi2 - azi1)/2;
54  if (abs(dlon) <= 90) {
55  azi1 = 90 - da;
56  azi2 = 90 + da;
57  } else {
58  azi1 = -90 - da;
59  azi2 = -90 + da;
60  }
61  }
62  if (dlon < 0) {
63  azi2 = azi1;
64  s12 = -s12;
65  sig12 = -sig12;
66  }
67  x = s12;
68  azi = Math::AngNormalize(azi2);
69  GeodesicLine perp(_earth.Line(lat, dlon, azi, Geodesic::GEODESICSCALE));
70  real t;
71  perp.GenPosition(true, -sig12,
73  t, t, t, t, t, t, rk, t);
74
75  real
76  alp0 = perp.EquatorialAzimuth() * Math::degree(),
77  calp0 = cos(alp0), salp0 = sin(alp0),
78  sbet1 = lat >=0 ? calp0 : -calp0,
79  cbet1 = abs(dlon) <= 90 ? abs(salp0) : -abs(salp0),
80  sbet01 = sbet1 * _cbet0 - cbet1 * _sbet0,
81  cbet01 = cbet1 * _cbet0 + sbet1 * _sbet0,
82  sig01 = atan2(sbet01, cbet01) / Math::degree();
83  _meridian.GenPosition(true, sig01,
85  t, t, t, y, t, t, t, t);
86  }
87
88  void CassiniSoldner::Reverse(real x, real y, real& lat, real& lon,
89  real& azi, real& rk) const {
90  if (!Init())
91  return;
92  real lat1, lon1;
93  real azi0, t;
94  _meridian.Position(y, lat1, lon1, azi0);
95  _earth.Direct(lat1, lon1, azi0 + 90, x, lat, lon, azi, rk, t);
96  }
97
98 } // namespace GeographicLib
static T AngNormalize(T x)
Definition: Math.hpp:445
CassiniSoldner(const Geodesic &earth=Geodesic::WGS84())
Math::real LatitudeOrigin() const
Math::real Direct(real lat1, real lon1, real azi1, real s12, real &lat2, real &lon2, real &azi2, real &m12, real &M12, real &M21, real &S12) const
Definition: Geodesic.hpp:379
GeodesicLine Line(real lat1, real lon1, real azi1, unsigned caps=ALL) const
Definition: Geodesic.cpp:118
void Forward(real lat, real lon, real &x, real &y, real &azi, real &rk) const
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, real &s12, real &azi1, real &azi2, real &m12, real &M12, real &M21, real &S12) const
Definition: Geodesic.hpp:679
Math::real Position(real s12, real &lat2, real &lon2, real &azi2, real &m12, real &M12, real &M21, real &S12) const
static void norm(T &x, T &y)
Definition: Math.hpp:392
void Reverse(real x, real y, real &lat, real &lon, real &azi, real &rk) const
Math::real GenPosition(bool arcmode, real s12_a12, unsigned outmask, real &lat2, real &lon2, real &azi2, real &s12, real &m12, real &M12, real &M21, real &S12) const
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
static T degree()
Definition: Math.hpp:228
static T AngDiff(T x, T y)
Definition: Math.hpp:475
Math::real Flattening() const
Definition: Geodesic.hpp:863
Geodesic calculations
Definition: Geodesic.hpp:171
void Reset(real lat0, real lon0)
static T AngRound(T x)
Definition: Math.hpp:498
Math::real LongitudeOrigin() const