46 # pragma warning (disable: 4127) 55 : tol_(numeric_limits<real>::epsilon())
56 , tol1_(real(0.1) * sqrt(tol_))
57 , tol2_(real(0.1) * tol_)
58 , taytol_(pow(tol_, real(0.6)))
86 void TransverseMercatorExact::zeta(real , real snu, real cnu, real dnu,
87 real , real snv, real cnv, real dnv,
88 real& taup, real& lam)
const {
95 overflow = 1 /
Math::sq(std::numeric_limits<real>::epsilon());
99 t1 = (d1 ? snu * dnv / d1 : (snu < 0 ? -overflow : overflow)),
100 t2 = (d2 ? sinh( _e *
Math::asinh(_e * snu / d2) ) :
101 (snu < 0 ? -overflow : overflow));
105 lam = (d1 != 0 && d2 != 0) ?
106 atan2(dnu * snv, cnu * cnv) - _e * atan2(_e * cnu * snv, dnu * cnv) :
110 void TransverseMercatorExact::dwdzeta(real ,
111 real snu, real cnu, real dnu,
113 real snv, real cnv, real dnv,
114 real& du, real& dv)
const {
123 bool TransverseMercatorExact::zetainv0(real psi, real lam, real& u, real& v)
127 lam > (1 - 2 * _e) *
Math::pi()/2 &&
128 psi < lam - (1 - _e) *
Math::pi()/2) {
143 v = atan2(cos(lamx), sinh(psix)) * (1 + _mu/2);
146 }
else if (psi < _e *
Math::pi()/2 &&
147 lam > (1 - 2 * _e) *
Math::pi()/2) {
160 dlam = lam - (1 - _e) *
Math::pi()/2,
167 ang = atan2(dlam-psi, psi+dlam) - real(0.75) *
Math::pi();
169 retval = rad < _e * taytol_;
173 v = rad * sin(ang) + _Ev.
K();
179 u = atan2(sinh(psi), cos(lam));
188 void TransverseMercatorExact::zetainv(real taup, real lam, real& u, real& v)
193 if (zetainv0(psi, lam, u, v))
195 real stol2 = tol2_ /
Math::sq(max(psi, real(1)));
198 real snu, cnu, dnu, snv, cnv, dnv;
199 _Eu.
sncndn(u, snu, cnu, dnu);
200 _Ev.
sncndn(v, snv, cnv, dnv);
201 real tau1, lam1, du1, dv1;
202 zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau1, lam1);
203 dwdzeta(u, snu, cnu, dnu, v, snv, cnv, dnv, du1, dv1);
208 delu = tau1 * du1 - lam1 * dv1,
209 delv = tau1 * dv1 + lam1 * du1;
215 if (!(delw2 >= stol2))
220 void TransverseMercatorExact::sigma(real , real snu, real cnu, real dnu,
221 real v, real snv, real cnv, real dnv,
222 real& xi, real& eta)
const {
226 xi = _Eu.
E(snu, cnu, dnu) - _mu * snu * cnu * dnu / d;
227 eta = v - _Ev.
E(snv, cnv, dnv) + _mv * snv * cnv * dnv / d;
230 void TransverseMercatorExact::dwdsigma(real ,
231 real snu, real cnu, real dnu,
233 real snv, real cnv, real dnv,
234 real& du, real& dv)
const {
239 dnr = dnu * cnv * dnv,
240 dni = - _mu * snu * cnu * snv;
242 dv = 2 * dnr * dni / d;
246 bool TransverseMercatorExact::sigmainv0(real xi, real eta, real& u, real& v)
249 if (eta > real(1.25) * _Ev.
KE() ||
250 (xi < -real(0.25) * _Eu.
E() && xi < eta - _Ev.
KE())) {
261 }
else if ((eta > real(0.75) * _Ev.
KE() && xi < real(0.25) * _Eu.
E())
276 deta = eta - _Ev.
KE(),
280 ang = atan2(deta-xi, xi+deta) - real(0.75) *
Math::pi();
282 retval = rad < 2 * taytol_;
286 v = rad * sin(ang) + _Ev.
K();
289 u = xi * _Eu.
K()/_Eu.
E();
290 v = eta * _Eu.
K()/_Eu.
E();
296 void TransverseMercatorExact::sigmainv(real xi, real eta, real& u, real& v)
298 if (sigmainv0(xi, eta, u, v))
302 real snu, cnu, dnu, snv, cnv, dnv;
303 _Eu.
sncndn(u, snu, cnu, dnu);
304 _Ev.
sncndn(v, snv, cnv, dnv);
305 real xi1, eta1, du1, dv1;
306 sigma(u, snu, cnu, dnu, v, snv, cnv, dnv, xi1, eta1);
307 dwdsigma(u, snu, cnu, dnu, v, snv, cnv, dnv, du1, dv1);
311 delu = xi1 * du1 - eta1 * dv1,
312 delv = xi1 * dv1 + eta1 * du1;
318 if (!(delw2 >= tol2_))
323 void TransverseMercatorExact::Scale(real tau, real ,
324 real snu, real cnu, real dnu,
325 real snv, real cnv, real dnv,
326 real& gamma, real& k)
const {
330 gamma = atan2(_mv * snu * snv * cnv, cnu * dnu * dnv);
344 k = sqrt(_mv + _mu / sec2) * sqrt(sec2) *
350 real& x, real& y, real& gamma, real& k)
356 latsign = (!_extendp && lat < 0) ? -1 : 1,
357 lonsign = (!_extendp && lon < 0) ? -1 : 1;
360 bool backside = !_extendp && lon > 90;
375 }
else if (lat == 0 && lon == 90 * (1 - _e)) {
382 real snu, cnu, dnu, snv, cnv, dnv;
383 _Eu.
sncndn(u, snu, cnu, dnu);
384 _Ev.
sncndn(v, snv, cnv, dnv);
387 sigma(u, snu, cnu, dnu, v, snv, cnv, dnv, xi, eta);
389 xi = 2 * _Eu.
E() - xi;
390 y = xi * _a * _k0 * latsign;
391 x = eta * _a * _k0 * lonsign;
398 zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau, lam);
400 Scale(tau, lam, snu, cnu, dnu, snv, cnv, dnv, gamma, k);
405 gamma *= latsign * lonsign;
410 real& lat, real& lon,
411 real& gamma, real& k)
416 eta = x / (_a * _k0);
419 latsign = !_extendp && y < 0 ? -1 : 1,
420 lonsign = !_extendp && x < 0 ? -1 : 1;
423 bool backside = !_extendp && xi > _Eu.
E();
425 xi = 2 * _Eu.
E()- xi;
429 if (xi == 0 && eta == _Ev.
KE()) {
433 sigmainv(xi, eta, u, v);
435 real snu, cnu, dnu, snv, cnv, dnv;
436 _Eu.
sncndn(u, snu, cnu, dnu);
437 _Ev.
sncndn(v, snv, cnv, dnv);
439 if (v != 0 || u != _Eu.
K()) {
440 zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau, lam);
445 Scale(tau, lam, snu, cnu, dnu, snv, cnv, dnv, gamma, k);
449 lon = lam = gamma = 0;
460 gamma *= latsign * lonsign;
static T AngNormalize(T x)
static const TransverseMercatorExact & UTM()
An exact implementation of the transverse Mercator projection.
static bool isfinite(T x)
Header for GeographicLib::TransverseMercatorExact class.
static T AngDiff(T x, T y, T &e)
TransverseMercatorExact(real a, real f, real k0, bool extendp=false)
Namespace for GeographicLib.
void sncndn(real x, real &sn, real &cn, real &dn) const
void Reverse(real lon0, real x, real y, real &lat, real &lon, real &gamma, real &k) const
void Forward(real lon0, real lat, real lon, real &x, real &y, real &gamma, real &k) const
Exception handling for GeographicLib.
static T tauf(T taup, T es)
static T taupf(T tau, T es)
#define GEOGRAPHICLIB_PANIC