GeographicLib  1.35
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-2011) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * http://geographiclib.sourceforge.net/
8  **********************************************************************/
9 
11 
12 namespace GeographicLib {
13 
14  using namespace std;
15 
16  const Math::real CassiniSoldner::eps1_ =
17  real(0.01) * sqrt(numeric_limits<real>::epsilon());
18  const Math::real CassiniSoldner::tiny_ = sqrt(numeric_limits<real>::min());
19 
20  void CassiniSoldner::Reset(real lat0, real lon0) throw() {
21  _meridian = _earth.Line(lat0, lon0, real(0),
25  real
26  phi = LatitudeOrigin() * Math::degree<real>(),
27  f = _earth.Flattening();
28  _sbet0 = (1 - f) * sin(phi);
29  _cbet0 = abs(LatitudeOrigin()) == 90 ? 0 : cos(phi);
30  SinCosNorm(_sbet0, _cbet0);
31  }
32 
33  void CassiniSoldner::Forward(real lat, real lon, real& x, real& y,
34  real& azi, real& rk) const throw() {
35  if (!Init())
36  return;
37  real dlon = Math::AngDiff(LongitudeOrigin(), Math::AngNormalize(lon));
38  real sig12, s12, azi1, azi2;
39  lat = AngRound(lat);
40  sig12 = _earth.Inverse(lat, -abs(dlon), lat, abs(dlon), s12, azi1, azi2);
41  if (sig12 < 100 * tiny_)
42  sig12 = s12 = 0;
43  sig12 *= real(0.5);
44  s12 *= real(0.5);
45  if (s12 == 0) {
46  real da = (azi2 - azi1)/2;
47  if (abs(dlon) <= 90) {
48  azi1 = 90 - da;
49  azi2 = 90 + da;
50  } else {
51  azi1 = -90 - da;
52  azi2 = -90 + da;
53  }
54  }
55  if (dlon < 0) {
56  azi2 = azi1;
57  s12 = -s12;
58  sig12 = -sig12;
59  }
60  x = s12;
61  azi = Math::AngNormalize(azi2);
62  GeodesicLine perp(_earth.Line(lat, dlon, azi, Geodesic::GEODESICSCALE));
63  real t;
64  perp.GenPosition(true, -sig12,
66  t, t, t, t, t, t, rk, t);
67 
68  real
69  alp0 = perp.EquatorialAzimuth() * Math::degree<real>(),
70  calp0 = cos(alp0), salp0 = sin(alp0),
71  sbet1 = lat >=0 ? calp0 : -calp0,
72  cbet1 = abs(dlon) <= 90 ? abs(salp0) : -abs(salp0),
73  sbet01 = sbet1 * _cbet0 - cbet1 * _sbet0,
74  cbet01 = cbet1 * _cbet0 + sbet1 * _sbet0,
75  sig01 = atan2(sbet01, cbet01) / Math::degree<real>();
76  _meridian.GenPosition(true, sig01,
78  t, t, t, y, t, t, t, t);
79  }
80 
81  void CassiniSoldner::Reverse(real x, real y, real& lat, real& lon,
82  real& azi, real& rk) const throw() {
83  if (!Init())
84  return;
85  real lat1, lon1;
86  real azi0, t;
87  _meridian.Position(y, lat1, lon1, azi0);
88  _earth.Direct(lat1, lon1, azi0 + 90, x, lat, lon, azi, rk, t);
89  }
90 
91 } // namespace GeographicLib