GeographicLib  1.35
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 namespace GeographicLib {
13 
14  using namespace std;
15 
16  const Math::real PolarStereographic::tol_ =
17  real(0.1)*sqrt(numeric_limits<real>::epsilon());
18  // Overflow value s.t. atan(overflow_) = pi/2
19  const Math::real PolarStereographic::overflow_ =
20  1 / Math::sq(numeric_limits<real>::epsilon());
21 
22  PolarStereographic::PolarStereographic(real a, real f, real k0)
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 
40  const PolarStereographic
41  PolarStereographic::UPS(Constants::WGS84_a<real>(),
42  Constants::WGS84_f<real>(),
43  Constants::UPS_k0<real>());
44 
45  // This formulation converts to conformal coordinates by tau = tan(phi) and
46  // tau' = tan(phi') where phi' is the conformal latitude. The formulas are:
47  // tau = tan(phi)
48  // secphi = hypot(1, tau)
49  // sig = sinh(e * atanh(e * tau / secphi))
50  // taup = tan(phip) = tau * hypot(1, sig) - sig * hypot(1, tau)
51  // c = (1 - f) * exp(e * atanh(e))
52  //
53  // Forward:
54  // rho = (2*k0*a/c) / (hypot(1, taup) + taup) (taup >= 0)
55  // = (2*k0*a/c) * (hypot(1, taup) - taup) (taup < 0)
56  //
57  // Reverse:
58  // taup = ((2*k0*a/c) / rho - rho / (2*k0*a/c))/2
59  //
60  // Scale:
61  // k = (rho/a) * secphi * sqrt((1-e2) + e2 / secphi^2)
62  //
63  // In limit rho -> 0, tau -> inf, taup -> inf, secphi -> inf, secphip -> inf
64  // secphip = taup = exp(-e * atanh(e)) * tau = exp(-e * atanh(e)) * secphi
65 
66  void PolarStereographic::Forward(bool northp, real lat, real lon,
67  real& x, real& y, real& gamma, real& k)
68  const throw() {
69  lat *= northp ? 1 : -1;
70  real
71  phi = lat * Math::degree<real>(),
72  tau = lat != -90 ? tanx(phi) : -overflow_,
73  secphi = Math::hypot(real(1), tau),
74  sig = sinh( eatanhe(tau / secphi) ),
75  taup = Math::hypot(real(1), sig) * tau - sig * secphi,
76  rho = Math::hypot(real(1), taup) + abs(taup);
77  rho = taup >= 0 ? (lat != 90 ? 1/rho : 0) : rho;
78  rho *= 2 * _k0 * _a / _c;
79  k = lat != 90 ? (rho / _a) * secphi * sqrt(_e2m + _e2 / Math::sq(secphi)) :
80  _k0;
81  lon = Math::AngNormalize(lon);
82  real
83  lam = lon * Math::degree<real>();
84  x = rho * (lon == -180 ? 0 : sin(lam));
85  y = (northp ? -rho : rho) * (abs(lon) == 90 ? 0 : cos(lam));
86  gamma = northp ? lon : -lon;
87  }
88 
89  void PolarStereographic::Reverse(bool northp, real x, real y,
90  real& lat, real& lon, real& gamma, real& k)
91  const throw() {
92  real
93  rho = Math::hypot(x, y),
94  t = rho / (2 * _k0 * _a / _c),
95  taup = (1 / t - t) / 2,
96  tau = taup * _Cx,
97  stol = tol_ * max(real(1), abs(taup));
98  if (abs(tau) < overflow_) {
99  // min iterations = 1, max iterations = 2; mean = 1.99
100  for (int i = 0; i < numit_; ++i) {
101  real
102  tau1 = Math::hypot(real(1), tau),
103  sig = sinh( eatanhe( tau / tau1 ) ),
104  taupa = Math::hypot(real(1), sig) * tau - sig * tau1,
105  dtau = (taup - taupa) * (1 + _e2m * Math::sq(tau)) /
106  ( _e2m * tau1 * Math::hypot(real(1), taupa) );
107  tau += dtau;
108  if (!(abs(dtau) >= stol))
109  break;
110  }
111  }
112  real
113  phi = atan(tau),
114  secphi = Math::hypot(real(1), tau);
115  k = rho != 0 ?
116  (rho / _a) * secphi * sqrt(_e2m + _e2 / Math::sq(secphi)) : _k0;
117  lat = (northp ? 1 : -1) * (rho != 0 ? phi / Math::degree<real>() : 90);
118  lon = -atan2( -x, northp ? -y : y ) / Math::degree<real>();
119  gamma = northp ? lon : -lon;
120  }
121 
122  void PolarStereographic::SetScale(real lat, real k) {
123  if (!(Math::isfinite(k) && k > 0))
124  throw GeographicErr("Scale is not positive");
125  if (!(-90 < lat && lat <= 90))
126  throw GeographicErr("Latitude must be in (-90d, 90d]");
127  real x, y, gamma, kold;
128  _k0 = 1;
129  Forward(true, lat, 0, x, y, gamma, kold);
130  _k0 *= k/kold;
131  }
132 
133 } // namespace GeographicLib
static T AngNormalize(T x)
Definition: Math.hpp:388
GeographicLib::Math::real real
Definition: GeodSolve.cpp:40
static bool isfinite(T x)
Definition: Math.hpp:435
static const PolarStereographic UPS
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:165
static T sq(T x)
Definition: Math.hpp:153
void SetScale(real lat, real k=real(1))
Polar stereographic projection.
Exception handling for GeographicLib.
Definition: Constants.hpp:320
Header for GeographicLib::PolarStereographic class.