GeographicLib  1.43
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-2015) <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  PolarStereographic::PolarStereographic(real a, real f, real k0)
17  : _a(a)
18  , _f(f <= 1 ? f : 1/f)
19  , _e2(_f * (2 - _f))
20  , _es((_f < 0 ? -1 : 1) * sqrt(abs(_e2)))
21  , _e2m(1 - _e2)
22  , _c( (1 - _f) * exp(Math::eatanhe(real(1), _es)) )
23  , _k0(k0)
24  {
25  if (!(Math::isfinite(_a) && _a > 0))
26  throw GeographicErr("Major radius is not positive");
27  if (!(Math::isfinite(_f) && _f < 1))
28  throw GeographicErr("Minor radius is not positive");
29  if (!(Math::isfinite(_k0) && _k0 > 0))
30  throw GeographicErr("Scale is not positive");
31  }
32 
34  static const PolarStereographic ups(Constants::WGS84_a(),
37  return ups;
38  }
39 
40  // This formulation converts to conformal coordinates by tau = tan(phi) and
41  // tau' = tan(phi') where phi' is the conformal latitude. The formulas are:
42  // tau = tan(phi)
43  // secphi = hypot(1, tau)
44  // sig = sinh(e * atanh(e * tau / secphi))
45  // taup = tan(phip) = tau * hypot(1, sig) - sig * hypot(1, tau)
46  // c = (1 - f) * exp(e * atanh(e))
47  //
48  // Forward:
49  // rho = (2*k0*a/c) / (hypot(1, taup) + taup) (taup >= 0)
50  // = (2*k0*a/c) * (hypot(1, taup) - taup) (taup < 0)
51  //
52  // Reverse:
53  // taup = ((2*k0*a/c) / rho - rho / (2*k0*a/c))/2
54  //
55  // Scale:
56  // k = (rho/a) * secphi * sqrt((1-e2) + e2 / secphi^2)
57  //
58  // In limit rho -> 0, tau -> inf, taup -> inf, secphi -> inf, secphip -> inf
59  // secphip = taup = exp(-e * atanh(e)) * tau = exp(-e * atanh(e)) * secphi
60 
61  void PolarStereographic::Forward(bool northp, real lat, real lon,
62  real& x, real& y, real& gamma, real& k)
63  const {
64  lat *= northp ? 1 : -1;
65  real
66  tau = Math::tand(lat),
67  secphi = Math::hypot(real(1), tau),
68  taup = Math::taupf(tau, _es),
69  rho = Math::hypot(real(1), taup) + abs(taup);
70  rho = taup >= 0 ? (lat != 90 ? 1/rho : 0) : rho;
71  rho *= 2 * _k0 * _a / _c;
72  k = lat != 90 ? (rho / _a) * secphi * sqrt(_e2m + _e2 / Math::sq(secphi)) :
73  _k0;
74  lon = Math::AngNormalize(lon);
75  real
76  lam = lon * Math::degree();
77  x = rho * (lon == -180 ? 0 : sin(lam));
78  y = (northp ? -rho : rho) * (abs(lon) == 90 ? 0 : cos(lam));
79  gamma = northp ? lon : -lon;
80  }
81 
82  void PolarStereographic::Reverse(bool northp, real x, real y,
83  real& lat, real& lon, real& gamma, real& k)
84  const {
85  real
86  rho = Math::hypot(x, y),
87  t = rho / (2 * _k0 * _a / _c),
88  taup = (1 / t - t) / 2,
89  tau = Math::tauf(taup, _es),
90  phi = atan(tau),
91  secphi = Math::hypot(real(1), tau);
92  k = rho ? (rho / _a) * secphi * sqrt(_e2m + _e2 / Math::sq(secphi)) : _k0;
93  lat = (northp ? 1 : -1) * (rho ? phi / Math::degree() : 90);
94  lon = 0 - atan2( -x, northp ? -y : y ) / Math::degree();
95  gamma = northp ? lon : -lon;
96  }
97 
98  void PolarStereographic::SetScale(real lat, real k) {
99  if (!(Math::isfinite(k) && k > 0))
100  throw GeographicErr("Scale is not positive");
101  if (!(-90 < lat && lat <= 90))
102  throw GeographicErr("Latitude must be in (-90d, 90d]");
103  real x, y, gamma, kold;
104  _k0 = 1;
105  Forward(true, lat, 0, x, y, gamma, kold);
106  _k0 *= k/kold;
107  }
108 
109 } // namespace GeographicLib
static T AngNormalize(T x)
Definition: Math.hpp:445
static bool isfinite(T x)
Definition: Math.hpp:614
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:102
void Reverse(bool northp, real x, real y, real &lat, real &lon, real &gamma, real &k) const
PolarStereographic(real a, real f, real k0)
static T hypot(T x, T y)
Definition: Math.hpp:255
static T sq(T x)
Definition: Math.hpp:244
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:228
Polar stereographic projection.
static T tand(T x)
Definition: Math.hpp:517
static T tauf(T taup, T es)
Exception handling for GeographicLib.
Definition: Constants.hpp:382
Header for GeographicLib::PolarStereographic class.
static T taupf(T tau, T es)
void Forward(bool northp, real lat, real lon, real &x, real &y, real &gamma, real &k) const