GeographicLib  1.35
Gnomonic.cpp
Go to the documentation of this file.
1 /**
2  * \file Gnomonic.cpp
3  * \brief Implementation for GeographicLib::Gnomonic class
4  *
5  * Copyright (c) Charles Karney (2010-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 #if defined(_MSC_VER)
13 // Squelch warnings about potentially uninitialized local variables
14 # pragma warning (disable: 4701)
15 #endif
16 
17 namespace GeographicLib {
18 
19  using namespace std;
20 
21  const Math::real Gnomonic::eps0_ = numeric_limits<real>::epsilon();
22  const Math::real Gnomonic::eps_ = real(0.01) * sqrt(eps0_);
23 
24  void Gnomonic::Forward(real lat0, real lon0, real lat, real lon,
25  real& x, real& y, real& azi, real& rk)
26  const throw() {
27  real azi0, m, M, t;
28  _earth.GenInverse(lat0, lon0, lat, lon,
31  t, azi0, azi, m, M, t, t);
32  rk = M;
33  if (M <= 0)
34  x = y = Math::NaN<real>();
35  else {
36  real rho = m/M;
37  azi0 *= Math::degree<real>();
38  x = rho * sin(azi0);
39  y = rho * cos(azi0);
40  }
41  }
42 
43  void Gnomonic::Reverse(real lat0, real lon0, real x, real y,
44  real& lat, real& lon, real& azi, real& rk)
45  const throw() {
46  real
47  azi0 = atan2(x, y) / Math::degree<real>(),
48  rho = Math::hypot(x, y),
49  s = _a * atan(rho/_a);
50  bool little = rho <= _a;
51  if (!little)
52  rho = 1/rho;
53  GeodesicLine line(_earth.Line(lat0, lon0, azi0,
58  int count = numit_, trip = 0;
59  real lat1, lon1, azi1, M;
60  while (count--) {
61  real m, t;
62  line.Position(s, lat1, lon1, azi1, m, M, t);
63  if (trip)
64  break;
65  // If little, solve rho(s) = rho with drho(s)/ds = 1/M^2
66  // else solve 1/rho(s) = 1/rho with d(1/rho(s))/ds = -1/m^2
67  real ds = little ? (m/M - rho) * M * M : (rho - M/m) * m * m;
68  s -= ds;
69  // Reversed test to allow escape with NaNs
70  if (!(abs(ds) >= eps_ * _a))
71  ++trip;
72  }
73  if (trip) {
74  lat = lat1; lon = lon1; azi = azi1; rk = M;
75  } else
76  lat = lon = azi = rk = Math::NaN<real>();
77  return;
78  }
79 
80 } // namespace GeographicLib
GeographicLib::Math::real real
Definition: GeodSolve.cpp:40
void Reverse(real lat0, real lon0, real x, real y, real &lat, real &lon, real &azi, real &rk) const
Definition: Gnomonic.cpp:43
Math::real Position(real s12, real &lat2, real &lon2, real &azi2, real &m12, real &M12, real &M21, real &S12) const
Header for GeographicLib::Gnomonic class.
static T hypot(T x, T y)
Definition: Math.hpp:165
void Forward(real lat0, real lon0, real lat, real lon, real &x, real &y, real &azi, real &rk) const
Definition: Gnomonic.cpp:24