GeographicLib  2.1.1
CassiniSoldner.cpp
Go to the documentation of this file.
1 /**
2  * \file CassiniSoldner.cpp
3  * \brief Implementation for GeographicLib::CassiniSoldner class
4  *
5  * Copyright (c) Charles Karney (2009-2022) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * https://geographiclib.sourceforge.io/
8  **********************************************************************/
9 
11 
12 #if defined(_MSC_VER)
13 // Squelch warnings about enum-float expressions
14 # pragma warning (disable: 5055)
15 #endif
16 
17 namespace GeographicLib {
18 
19  using namespace std;
20 
22  : _earth(earth) {}
23 
24  CassiniSoldner::CassiniSoldner(real lat0, real lon0, const Geodesic& earth)
25  : _earth(earth)
26  { Reset(lat0, lon0); }
27 
28  void CassiniSoldner::Reset(real lat0, real lon0) {
29  _meridian = _earth.Line(lat0, lon0, real(0),
33  real f = _earth.Flattening();
34  Math::sincosd(LatitudeOrigin(), _sbet0, _cbet0);
35  _sbet0 *= (1 - f);
36  Math::norm(_sbet0, _cbet0);
37  }
38 
39  void CassiniSoldner::Forward(real lat, real lon, real& x, real& y,
40  real& azi, real& rk) const {
41  if (!Init())
42  return;
43  real dlon = Math::AngDiff(LongitudeOrigin(), lon);
44  real sig12, s12, azi1, azi2;
45  sig12 = _earth.Inverse(lat, -fabs(dlon), lat, fabs(dlon), s12, azi1, azi2);
46  sig12 *= real(0.5);
47  s12 *= real(0.5);
48  if (s12 == 0) {
49  real da = Math::AngDiff(azi1, azi2)/2;
50  if (fabs(dlon) <= Math::qd) {
51  azi1 = Math::qd - da;
52  azi2 = Math::qd + da;
53  } else {
54  azi1 = -Math::qd - da;
55  azi2 = -Math::qd + da;
56  }
57  }
58  if (signbit(dlon)) {
59  azi2 = azi1;
60  s12 = -s12;
61  sig12 = -sig12;
62  }
63  x = s12;
64  azi = Math::AngNormalize(azi2);
65  GeodesicLine perp(_earth.Line(lat, dlon, azi, Geodesic::GEODESICSCALE));
66  real t;
67  perp.GenPosition(true, -sig12,
69  t, t, t, t, t, t, rk, t);
70 
71  real salp0, calp0;
72  Math::sincosd(perp.EquatorialAzimuth(), salp0, calp0);
73  real
74  sbet1 = lat >=0 ? calp0 : -calp0,
75  cbet1 = fabs(dlon) <= Math::qd ? fabs(salp0) : -fabs(salp0),
76  sbet01 = sbet1 * _cbet0 - cbet1 * _sbet0,
77  cbet01 = cbet1 * _cbet0 + sbet1 * _sbet0,
78  sig01 = atan2(sbet01, cbet01) / Math::degree();
79  _meridian.GenPosition(true, sig01,
81  t, t, t, y, t, t, t, t);
82  }
83 
84  void CassiniSoldner::Reverse(real x, real y, real& lat, real& lon,
85  real& azi, real& rk) const {
86  if (!Init())
87  return;
88  real lat1, lon1;
89  real azi0, t;
90  _meridian.Position(y, lat1, lon1, azi0);
91  _earth.Direct(lat1, lon1, azi0 + Math::qd, x, lat, lon, azi, rk, t);
92  }
93 
94 } // namespace GeographicLib
Header for GeographicLib::CassiniSoldner class.
Math::real LatitudeOrigin() const
CassiniSoldner(const Geodesic &earth=Geodesic::WGS84())
Math::real LongitudeOrigin() const
void Reverse(real x, real y, real &lat, real &lon, real &azi, real &rk) const
void Forward(real lat, real lon, real &x, real &y, real &azi, real &rk) const
void Reset(real lat0, real lon0)
Math::real EquatorialAzimuth() const
Math::real Position(real s12, real &lat2, real &lon2, real &azi2, real &m12, real &M12, real &M21, real &S12) 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
Geodesic calculations
Definition: Geodesic.hpp:172
Math::real Flattening() const
Definition: Geodesic.hpp:955
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:385
GeodesicLine Line(real lat1, real lon1, real azi1, unsigned caps=ALL) const
Definition: Geodesic.cpp:118
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:680
static T degree()
Definition: Math.hpp:200
static void sincosd(T x, T &sinx, T &cosx)
Definition: Math.cpp:106
static void norm(T &x, T &y)
Definition: Math.hpp:222
static T AngNormalize(T x)
Definition: Math.cpp:71
static T AngDiff(T x, T y, T &e)
Definition: Math.cpp:82
@ qd
degrees per quarter turn
Definition: Math.hpp:141
Namespace for GeographicLib.
Definition: Accumulator.cpp:12