src/PolarStereographic.cpp
Go to the documentation of this file.
1 
11 
12 namespace GeographicLib {
13 
14  using namespace std;
15 
17  : _a(a)
18  , _f(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("Equatorial radius is not positive");
27  if (!(Math::isfinite(_f) && _f < 1))
28  throw GeographicErr("Polar semi-axis 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 
62  real& x, real& y,
63  real& gamma, real& k) const {
64  lat = Math::LatFix(lat);
65  lat *= northp ? 1 : -1;
66  real
67  tau = Math::tand(lat),
68  secphi = Math::hypot(real(1), tau),
69  taup = Math::taupf(tau, _es),
70  rho = Math::hypot(real(1), taup) + abs(taup);
71  rho = taup >= 0 ? (lat != 90 ? 1/rho : 0) : rho;
72  rho *= 2 * _k0 * _a / _c;
73  k = lat != 90 ? (rho / _a) * secphi * sqrt(_e2m + _e2 / Math::sq(secphi)) :
74  _k0;
75  Math::sincosd(lon, x, y);
76  x *= rho;
77  y *= (northp ? -rho : rho);
78  gamma = Math::AngNormalize(northp ? lon : -lon);
79  }
80 
81  void PolarStereographic::Reverse(bool northp, real x, real y,
82  real& lat, real& lon,
83  real& gamma, real& k) const {
84  real
85  rho = Math::hypot(x, y),
86  t = rho != 0 ? rho / (2 * _k0 * _a / _c) :
88  taup = (1 / t - t) / 2,
89  tau = Math::tauf(taup, _es),
90  secphi = Math::hypot(real(1), tau);
91  k = rho != 0 ? (rho / _a) * secphi * sqrt(_e2m + _e2 / Math::sq(secphi)) :
92  _k0;
93  lat = (northp ? 1 : -1) * Math::atand(tau);
94  lon = Math::atan2d(x, northp ? -y : y );
95  gamma = Math::AngNormalize(northp ? lon : -lon);
96  }
97 
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:440
Scalar * y
static T atand(T x)
Definition: Math.hpp:723
static const double lat
static bool isfinite(T x)
Definition: Math.hpp:806
static T LatFix(T x)
Definition: Math.hpp:467
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:102
Definition: BFloat16.h:88
static void sincosd(T x, T &sinx, T &cosx)
Definition: Math.hpp:558
void Reverse(bool northp, real x, real y, real &lat, real &lon, real &gamma, real &k) const
static double epsilon
Definition: testRot3.cpp:37
PolarStereographic(real a, real f, real k0)
EIGEN_DEVICE_FUNC const ExpReturnType exp() const
static T hypot(T x, T y)
Definition: Math.hpp:243
static T sq(T x)
Definition: Math.hpp:232
static T atan2d(T y, T x)
Definition: Math.hpp:691
static const PolarStereographic & UPS()
void SetScale(real lat, real k=real(1))
Definition: main.h:100
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
Namespace for GeographicLib.
Polar stereographic projection.
static T tand(T x)
Definition: Math.hpp:671
Exception handling for GeographicLib.
Definition: Constants.hpp:389
static const double lon
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
static T tauf(T taup, T es)
Definition: Math.cpp:31
Header for GeographicLib::PolarStereographic class.
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
#define abs(x)
Definition: datatypes.h:17
Point2 t(10, 10)
static T taupf(T tau, T es)
Definition: Math.cpp:25
void Forward(bool northp, real lat, real lon, real &x, real &y, real &gamma, real &k) const


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:35:14