GeographicLib  2.1
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-2022) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * https://geographiclib.sourceforge.io/
8  **********************************************************************/
9 
11 
12 #if defined(_MSC_VER)
13 // Squelch warnings about enum-float expressions
14 # pragma warning (disable: 5055)
15 #endif
16 
17 namespace GeographicLib {
18 
19  using namespace std;
20 
21  PolarStereographic::PolarStereographic(real a, real f, real k0)
22  : _a(a)
23  , _f(f)
24  , _e2(_f * (2 - _f))
25  , _es((_f < 0 ? -1 : 1) * sqrt(fabs(_e2)))
26  , _e2m(1 - _e2)
27  , _c( (1 - _f) * exp(Math::eatanhe(real(1), _es)) )
28  , _k0(k0)
29  {
30  if (!(isfinite(_a) && _a > 0))
31  throw GeographicErr("Equatorial radius is not positive");
32  if (!(isfinite(_f) && _f < 1))
33  throw GeographicErr("Polar semi-axis is not positive");
34  if (!(isfinite(_k0) && _k0 > 0))
35  throw GeographicErr("Scale is not positive");
36  }
37 
39  static const PolarStereographic ups(Constants::WGS84_a(),
42  return ups;
43  }
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,
68  real& gamma, real& k) const {
69  lat = Math::LatFix(lat);
70  lat *= northp ? 1 : -1;
71  real
72  tau = Math::tand(lat),
73  secphi = hypot(real(1), tau),
74  taup = Math::taupf(tau, _es),
75  rho = hypot(real(1), taup) + fabs(taup);
76  rho = taup >= 0 ? (lat != Math::qd ? 1/rho : 0) : rho;
77  rho *= 2 * _k0 * _a / _c;
78  k = lat != Math::qd ?
79  (rho / _a) * secphi * sqrt(_e2m + _e2 / Math::sq(secphi)) : _k0;
80  Math::sincosd(lon, x, y);
81  x *= rho;
82  y *= (northp ? -rho : rho);
83  gamma = Math::AngNormalize(northp ? lon : -lon);
84  }
85 
86  void PolarStereographic::Reverse(bool northp, real x, real y,
87  real& lat, real& lon,
88  real& gamma, real& k) const {
89  real
90  rho = hypot(x, y),
91  t = rho != 0 ? rho / (2 * _k0 * _a / _c) :
92  Math::sq(numeric_limits<real>::epsilon()),
93  taup = (1 / t - t) / 2,
94  tau = Math::tauf(taup, _es),
95  secphi = hypot(real(1), tau);
96  k = rho != 0 ? (rho / _a) * secphi * sqrt(_e2m + _e2 / Math::sq(secphi)) :
97  _k0;
98  lat = (northp ? 1 : -1) * Math::atand(tau);
99  lon = Math::atan2d(x, northp ? -y : y );
100  gamma = Math::AngNormalize(northp ? lon : -lon);
101  }
102 
103  void PolarStereographic::SetScale(real lat, real k) {
104  if (!(isfinite(k) && k > 0))
105  throw GeographicErr("Scale is not positive");
106  if (!(-Math::qd < lat && lat <= Math::qd))
107  throw GeographicErr("Latitude must be in (-" + to_string(Math::qd)
108  + "d, " + to_string(Math::qd) + "d]");
109  real x, y, gamma, kold;
110  _k0 = 1;
111  Forward(true, lat, 0, x, y, gamma, kold);
112  _k0 *= k/kold;
113  }
114 
115 } // namespace GeographicLib
Header for GeographicLib::PolarStereographic class.
Exception handling for GeographicLib.
Definition: Constants.hpp:316
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:76
static T tand(T x)
Definition: Math.cpp:172
static T LatFix(T x)
Definition: Math.hpp:300
static void sincosd(T x, T &sinx, T &cosx)
Definition: Math.cpp:106
static T atan2d(T y, T x)
Definition: Math.cpp:183
static T sq(T x)
Definition: Math.hpp:212
static T tauf(T taup, T es)
Definition: Math.cpp:219
static T AngNormalize(T x)
Definition: Math.cpp:71
static T atand(T x)
Definition: Math.cpp:202
static T taupf(T tau, T es)
Definition: Math.cpp:209
@ qd
degrees per quarter turn
Definition: Math.hpp:141
Polar stereographic projection.
void Reverse(bool northp, real x, real y, real &lat, real &lon, real &gamma, real &k) const
PolarStereographic(real a, real f, real k0)
void SetScale(real lat, real k=real(1))
static const PolarStereographic & UPS()
void Forward(bool northp, real lat, real lon, real &x, real &y, real &gamma, real &k) const
Namespace for GeographicLib.
Definition: Accumulator.cpp:12