GeographicLib  1.37
PolarStereographic.cpp
Go to the documentation of this file.
1 /**
2  * \file PolarStereographic.cpp
3  * \brief Implementation for GeographicLib::PolarStereographic class
4  *
5  * Copyright (c) Charles Karney (2008-2012) <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 constant conditional expressions
14 # pragma warning (disable: 4127)
15 #endif
16 
17 namespace GeographicLib {
18 
19  using namespace std;
20 
21  PolarStereographic::PolarStereographic(real a, real f, real k0)
22  : tol_(real(0.1)*sqrt(numeric_limits<real>::epsilon()))
23  , _a(a)
24  , _f(f <= 1 ? f : 1/f)
25  , _e2(_f * (2 - _f))
26  , _e(sqrt(abs(_e2)))
27  , _e2m(1 - _e2)
28  , _Cx(exp(eatanhe(real(1))))
29  , _c( (1 - _f) * _Cx )
30  , _k0(k0)
31  {
32  if (!(Math::isfinite(_a) && _a > 0))
33  throw GeographicErr("Major radius is not positive");
34  if (!(Math::isfinite(_f) && _f < 1))
35  throw GeographicErr("Minor radius is not positive");
36  if (!(Math::isfinite(_k0) && _k0 > 0))
37  throw GeographicErr("Scale is not positive");
38  }
39 
41  static const PolarStereographic ups(Constants::WGS84_a(),
44  return ups;
45  }
46 
47  // This formulation converts to conformal coordinates by tau = tan(phi) and
48  // tau' = tan(phi') where phi' is the conformal latitude. The formulas are:
49  // tau = tan(phi)
50  // secphi = hypot(1, tau)
51  // sig = sinh(e * atanh(e * tau / secphi))
52  // taup = tan(phip) = tau * hypot(1, sig) - sig * hypot(1, tau)
53  // c = (1 - f) * exp(e * atanh(e))
54  //
55  // Forward:
56  // rho = (2*k0*a/c) / (hypot(1, taup) + taup) (taup >= 0)
57  // = (2*k0*a/c) * (hypot(1, taup) - taup) (taup < 0)
58  //
59  // Reverse:
60  // taup = ((2*k0*a/c) / rho - rho / (2*k0*a/c))/2
61  //
62  // Scale:
63  // k = (rho/a) * secphi * sqrt((1-e2) + e2 / secphi^2)
64  //
65  // In limit rho -> 0, tau -> inf, taup -> inf, secphi -> inf, secphip -> inf
66  // secphip = taup = exp(-e * atanh(e)) * tau = exp(-e * atanh(e)) * secphi
67 
68  void PolarStereographic::Forward(bool northp, real lat, real lon,
69  real& x, real& y, real& gamma, real& k)
70  const {
71  lat *= northp ? 1 : -1;
72  real
73  phi = lat * Math::degree(),
74  tau = lat != -90 ? tanx(phi) : -overflow(),
75  secphi = Math::hypot(real(1), tau),
76  sig = sinh( eatanhe(tau / secphi) ),
77  taup = Math::hypot(real(1), sig) * tau - sig * secphi,
78  rho = Math::hypot(real(1), taup) + abs(taup);
79  rho = taup >= 0 ? (lat != 90 ? 1/rho : 0) : rho;
80  rho *= 2 * _k0 * _a / _c;
81  k = lat != 90 ? (rho / _a) * secphi * sqrt(_e2m + _e2 / Math::sq(secphi)) :
82  _k0;
83  lon = Math::AngNormalize(lon);
84  real
85  lam = lon * Math::degree();
86  x = rho * (lon == -180 ? 0 : sin(lam));
87  y = (northp ? -rho : rho) * (abs(lon) == 90 ? 0 : cos(lam));
88  gamma = northp ? lon : -lon;
89  }
90 
91  void PolarStereographic::Reverse(bool northp, real x, real y,
92  real& lat, real& lon, real& gamma, real& k)
93  const {
94  real
95  rho = Math::hypot(x, y),
96  t = rho / (2 * _k0 * _a / _c),
97  taup = (1 / t - t) / 2,
98  tau = taup * _Cx,
99  stol = tol_ * max(real(1), abs(taup));
100  if (abs(tau) < overflow()) {
101  // min iterations = 1, max iterations = 2; mean = 1.99
102  for (int i = 0; i < numit_ || GEOGRAPHICLIB_PANIC; ++i) {
103  real
104  tau1 = Math::hypot(real(1), tau),
105  sig = sinh( eatanhe( tau / tau1 ) ),
106  taupa = Math::hypot(real(1), sig) * tau - sig * tau1,
107  dtau = (taup - taupa) * (1 + _e2m * Math::sq(tau)) /
108  ( _e2m * tau1 * Math::hypot(real(1), taupa) );
109  tau += dtau;
110  if (!(abs(dtau) >= stol))
111  break;
112  }
113  }
114  real
115  phi = atan(tau),
116  secphi = Math::hypot(real(1), tau);
117  k = rho ? (rho / _a) * secphi * sqrt(_e2m + _e2 / Math::sq(secphi)) : _k0;
118  lat = (northp ? 1 : -1) * (rho ? phi / Math::degree() : 90);
119  lon = -atan2( -x, northp ? -y : y ) / Math::degree();
120  gamma = northp ? lon : -lon;
121  }
122 
123  void PolarStereographic::SetScale(real lat, real k) {
124  if (!(Math::isfinite(k) && k > 0))
125  throw GeographicErr("Scale is not positive");
126  if (!(-90 < lat && lat <= 90))
127  throw GeographicErr("Latitude must be in (-90d, 90d]");
128  real x, y, gamma, kold;
129  _k0 = 1;
130  Forward(true, lat, 0, x, y, gamma, kold);
131  _k0 *= k/kold;
132  }
133 
134 } // namespace GeographicLib
static T AngNormalize(T x)
Definition: Math.hpp:399
static bool isfinite(T x)
Definition: Math.hpp:445
void Forward(bool northp, real lat, real lon, real &x, real &y, real &gamma, real &k) const
PolarStereographic(real a, real f, real k0)
void Reverse(bool northp, real x, real y, real &lat, real &lon, real &gamma, real &k) const
static T hypot(T x, T y)
Definition: Math.hpp:254
static T sq(T x)
Definition: Math.hpp:243
static const PolarStereographic & UPS()
void SetScale(real lat, real k=real(1))
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
static T degree()
Definition: Math.hpp:227
Polar stereographic projection.
Exception handling for GeographicLib.
Definition: Constants.hpp:362
Header for GeographicLib::PolarStereographic class.
#define GEOGRAPHICLIB_PANIC
Definition: Math.hpp:86