46 # pragma warning (disable: 4127)
55 : tol_(numeric_limits<real>::epsilon())
56 , tol2_(real(0.1) * tol_)
57 , taytol_(pow(tol_, real(0.6)))
68 if (!(isfinite(_a) && _a > 0))
74 if (!(isfinite(_k0) && _k0 > 0))
94 overflow = 1 /
Math::sq(numeric_limits<real>::epsilon());
98 t1 = (d1 != 0 ? snu * dnv / d1 : (signbit(snu) ? -overflow : overflow)),
99 t2 = (d2 != 0 ? sinh( _e * asinh(_e * snu / d2) ) :
100 (signbit(snu) ? -overflow : overflow));
103 taup = t1 * hypot(
real(1), t2) - t2 * hypot(
real(1), t1);
104 lam = (d1 != 0 && d2 != 0) ?
105 atan2(dnu * snv, cnu * cnv) - _e * atan2(_e * cnu * snv, dnu * cnv) :
109 void TransverseMercatorExact::dwdzeta(
real ,
122 bool TransverseMercatorExact::zetainv0(
real psi,
real lam,
126 lam > (1 - 2 * _e) *
Math::pi()/2 &&
127 psi < lam - (1 - _e) *
Math::pi()/2) {
140 u = asinh(sin(lamx) / hypot(cos(lamx), sinh(psix))) *
142 v = atan2(cos(lamx), sinh(psix)) * (1 + _mu/2);
145 }
else if (psi < _e *
Math::pi()/2 &&
146 lam > (1 - 2 * _e) *
Math::pi()/2) {
159 dlam = lam - (1 - _e) *
Math::pi()/2,
160 rad = hypot(psi, dlam),
166 ang = atan2(dlam-psi, psi+dlam) -
real(0.75) *
Math::pi();
168 retval = rad < _e * taytol_;
169 rad = cbrt(3 / (_mv * _e) * rad);
172 v = rad * sin(ang) + _eEv.
K();
177 v = asinh(sin(lam) / hypot(cos(lam), sinh(psi)));
178 u = atan2(sinh(psi), cos(lam));
187 void TransverseMercatorExact::zetainv(
real taup,
real lam,
191 scal = 1/hypot(
real(1), taup);
192 if (zetainv0(psi, lam, u, v))
197 real snu, cnu, dnu, snv, cnv, dnv;
198 _eEu.
sncndn(u, snu, cnu, dnu);
199 _eEv.
sncndn(v, snv, cnv, dnv);
200 real tau1, lam1, du1, dv1;
201 zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau1, lam1);
202 dwdzeta(u, snu, cnu, dnu, v, snv, cnv, dnv, du1, dv1);
207 delu = tau1 * du1 - lam1 * dv1,
208 delv = tau1 * dv1 + lam1 * du1;
214 if (!(delw2 >= stol2))
225 xi = _eEu.
E(snu, cnu, dnu) - _mu * snu * cnu * dnu / d;
226 eta = v - _eEv.
E(snv, cnv, dnv) + _mv * snv * cnv * dnv / d;
229 void TransverseMercatorExact::dwdsigma(
real ,
238 dnr = dnu * cnv * dnv,
239 dni = - _mu * snu * cnu * snv;
241 dv = 2 * dnr * dni / d;
245 bool TransverseMercatorExact::sigmainv0(
real xi,
real eta,
248 if (eta >
real(1.25) * _eEv.
KE() ||
249 (xi < -
real(0.25) * _eEu.
E() && xi < eta - _eEv.
KE())) {
260 }
else if ((eta >
real(0.75) * _eEv.
KE() && xi <
real(0.25) * _eEu.
E())
261 || eta > _eEv.
KE()) {
275 deta = eta - _eEv.
KE(),
276 rad = hypot(xi, deta),
279 ang = atan2(deta-xi, xi+deta) -
real(0.75) *
Math::pi();
281 retval = rad < 2 * taytol_;
282 rad = cbrt(3 / _mv * rad);
285 v = rad * sin(ang) + _eEv.
K();
288 u = xi * _eEu.
K()/_eEu.
E();
289 v = eta * _eEu.
K()/_eEu.
E();
295 void TransverseMercatorExact::sigmainv(
real xi,
real eta,
297 if (sigmainv0(xi, eta, u, v))
301 real snu, cnu, dnu, snv, cnv, dnv;
302 _eEu.
sncndn(u, snu, cnu, dnu);
303 _eEv.
sncndn(v, snv, cnv, dnv);
304 real xi1, eta1, du1, dv1;
305 sigma(u, snu, cnu, dnu, v, snv, cnv, dnv, xi1, eta1);
306 dwdsigma(u, snu, cnu, dnu, v, snv, cnv, dnv, du1, dv1);
310 delu = xi1 * du1 - eta1 * dv1,
311 delv = xi1 * dv1 + eta1 * du1;
317 if (!(delw2 >= tol2_))
322 void TransverseMercatorExact::Scale(
real tau,
real ,
329 gamma = atan2(_mv * snu * snv * cnv, cnu * dnu * dnv);
343 k = sqrt(_mv + _mu / sec2) * sqrt(sec2) *
350 real& gamma, real& k)
const {
355 latsign = (!_extendp && signbit(lat)) ? -1 : 1,
356 lonsign = (!_extendp && signbit(lon)) ? -1 : 1;
359 bool backside = !_extendp && lon >
Math::qd;
374 }
else if (lat == 0 && lon ==
Math::qd * (1 - _e)) {
381 real snu, cnu, dnu, snv, cnv, dnv;
382 _eEu.
sncndn(u, snu, cnu, dnu);
383 _eEv.
sncndn(v, snv, cnv, dnv);
386 sigma(u, snu, cnu, dnu, v, snv, cnv, dnv, xi, eta);
388 xi = 2 * _eEu.
E() - xi;
389 y = xi * _a * _k0 * latsign;
390 x = eta * _a * _k0 * lonsign;
397 zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau, lam);
399 Scale(tau, lam, snu, cnu, dnu, snv, cnv, dnv, gamma, k);
404 gamma *= latsign * lonsign;
409 real& lat, real& lon,
410 real& gamma, real& k)
const {
414 eta = x / (_a * _k0);
417 xisign = (!_extendp && signbit(xi)) ? -1 : 1,
418 etasign = (!_extendp && signbit(eta)) ? -1 : 1;
421 bool backside = !_extendp && xi > _eEu.
E();
423 xi = 2 * _eEu.
E()- xi;
427 if (xi == 0 && eta == _eEv.
KE()) {
431 sigmainv(xi, eta, u, v);
433 real snu, cnu, dnu, snv, cnv, dnv;
434 _eEu.
sncndn(u, snu, cnu, dnu);
435 _eEv.
sncndn(v, snv, cnv, dnv);
437 if (v != 0 || u != _eEu.
K()) {
438 zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau, lam);
443 Scale(tau, lam, snu, cnu, dnu, snv, cnv, dnv, gamma, k);
447 lon = lam = gamma = 0;
458 gamma *= xisign * etasign;
GeographicLib::Math::real real
#define GEOGRAPHICLIB_PANIC
Header for GeographicLib::TransverseMercatorExact class.
void sncndn(real x, real &sn, real &cn, real &dn) const
Exception handling for GeographicLib.
static T tauf(T taup, T es)
static T AngNormalize(T x)
static T taupf(T tau, T es)
static T AngDiff(T x, T y, T &e)
@ hd
degrees per half turn
@ qd
degrees per quarter turn
An exact implementation of the transverse Mercator projection.
void Forward(real lon0, real lat, real lon, real &x, real &y, real &gamma, real &k) const
TransverseMercatorExact(real a, real f, real k0, bool extendp=false)
static const TransverseMercatorExact & UTM()
void Reverse(real lon0, real x, real y, real &lat, real &lon, real &gamma, real &k) const
Namespace for GeographicLib.