46 # pragma warning (disable: 4127)
49 namespace GeographicLib {
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)))
60 , _f(f <= 1 ? f : 1/f)
98 stol = tol_ * max(
real(1), abs(taup));
105 dtau = (taup - taupa) * (1 + _mv *
Math::sq(tau)) /
108 if (!(abs(dtau) >= stol))
124 t1 = (d1 ? snu * dnv / d1 : (snu < 0 ? -overflow() : overflow())),
125 t2 = (d2 ? sinh( _e * Math::asinh(_e * snu / d2) ) :
126 (snu < 0 ? -overflow() : overflow()));
130 lam = (d1 != 0 && d2 != 0) ?
131 atan2(dnu * snv, cnu * cnv) - _e * atan2(_e * cnu * snv, dnu * cnv) :
135 void TransverseMercatorExact::dwdzeta(
real ,
148 bool TransverseMercatorExact::zetainv0(
real psi,
real lam,
real& u,
real& v)
152 lam > (1 - 2 * _e) *
Math::pi()/2 &&
153 psi < lam - (1 - _e) *
Math::pi()/2) {
168 v = atan2(cos(lamx), sinh(psix)) * (1 + _mu/2);
171 }
else if (psi < _e *
Math::pi()/2 &&
172 lam > (1 - 2 * _e) *
Math::pi()/2) {
185 dlam = lam - (1 - _e) *
Math::pi()/2,
192 ang = atan2(dlam-psi, psi+dlam) -
real(0.75) *
Math::pi();
194 retval = rad < _e * taytol_;
198 v = rad * sin(ang) + _Ev.
K();
204 u = atan2(sinh(psi), cos(lam));
213 void TransverseMercatorExact::zetainv(
real taup,
real lam,
real& u,
real& v)
218 if (zetainv0(psi, lam, u, v))
223 real snu, cnu, dnu, snv, cnv, dnv;
224 _Eu.
sncndn(u, snu, cnu, dnu);
225 _Ev.
sncndn(v, snv, cnv, dnv);
226 real tau1, lam1, du1, dv1;
227 zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau1, lam1);
228 dwdzeta(u, snu, cnu, dnu, v, snv, cnv, dnv, du1, dv1);
233 delu = tau1 * du1 - lam1 * dv1,
234 delv = tau1 * dv1 + lam1 * du1;
240 if (!(delw2 >= stol2))
251 xi = _Eu.
E(snu, cnu, dnu) - _mu * snu * cnu * dnu / d;
252 eta = v - _Ev.
E(snv, cnv, dnv) + _mv * snv * cnv * dnv / d;
255 void TransverseMercatorExact::dwdsigma(
real ,
264 dnr = dnu * cnv * dnv,
265 dni = - _mu * snu * cnu * snv;
267 dv = 2 * dnr * dni / d;
271 bool TransverseMercatorExact::sigmainv0(
real xi,
real eta,
real& u,
real& v)
274 if (eta >
real(1.25) * _Ev.
KE() ||
275 (xi < -
real(0.25) * _Eu.
E() && xi < eta - _Ev.
KE())) {
286 }
else if ((eta >
real(0.75) * _Ev.
KE() && xi <
real(0.25) * _Eu.
E())
301 deta = eta - _Ev.
KE(),
305 ang = atan2(deta-xi, xi+deta) -
real(0.75) *
Math::pi();
307 retval = rad < 2 * taytol_;
311 v = rad * sin(ang) + _Ev.
K();
314 u = xi * _Eu.
K()/_Eu.
E();
315 v = eta * _Eu.
K()/_Eu.
E();
323 if (sigmainv0(xi, eta, u, v))
327 real snu, cnu, dnu, snv, cnv, dnv;
328 _Eu.
sncndn(u, snu, cnu, dnu);
329 _Ev.
sncndn(v, snv, cnv, dnv);
330 real xi1, eta1, du1, dv1;
331 sigma(u, snu, cnu, dnu, v, snv, cnv, dnv, xi1, eta1);
332 dwdsigma(u, snu, cnu, dnu, v, snv, cnv, dnv, du1, dv1);
336 delu = xi1 * du1 - eta1 * dv1,
337 delv = xi1 * dv1 + eta1 * du1;
343 if (!(delw2 >= tol2_))
348 void TransverseMercatorExact::Scale(
real tau,
real ,
355 gamma = atan2(_mv * snu * snv * cnv, cnu * dnu * dnv);
369 k = sqrt(_mv + _mu / sec2) * sqrt(sec2) *
375 real& x, real& y, real& gamma, real& k)
380 latsign = (!_extendp && lat < 0) ? -1 : 1,
381 lonsign = (!_extendp && lon < 0) ? -1 : 1;
384 bool backside = !_extendp && lon > 90;
400 }
else if (lat == 0 && lon == 90 * (1 - _e)) {
404 zetainv(taup(tau), lam, u, v);
406 real snu, cnu, dnu, snv, cnv, dnv;
407 _Eu.
sncndn(u, snu, cnu, dnu);
408 _Ev.
sncndn(v, snv, cnv, dnv);
411 sigma(u, snu, cnu, dnu, v, snv, cnv, dnv, xi, eta);
413 xi = 2 * _Eu.
E() - xi;
414 y = xi * _a * _k0 * latsign;
415 x = eta * _a * _k0 * lonsign;
422 zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau, lam);
424 Scale(tau, lam, snu, cnu, dnu, snv, cnv, dnv, gamma, k);
429 gamma *= latsign * lonsign;
434 real& lat, real& lon,
435 real& gamma, real& k)
440 eta = x / (_a * _k0);
443 latsign = !_extendp && y < 0 ? -1 : 1,
444 lonsign = !_extendp && x < 0 ? -1 : 1;
447 bool backside = !_extendp && xi > _Eu.
E();
449 xi = 2 * _Eu.
E()- xi;
453 if (xi == 0 && eta == _Ev.
KE()) {
457 sigmainv(xi, eta, u, v);
459 real snu, cnu, dnu, snv, cnv, dnv;
460 _Eu.
sncndn(u, snu, cnu, dnu);
461 _Ev.
sncndn(v, snv, cnv, dnv);
463 if (v != 0 || u != _Eu.
K()) {
464 zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau, lam);
469 Scale(tau, lam, snu, cnu, dnu, snv, cnv, dnv, gamma, k);
473 lon = lam = gamma = 0;
484 gamma *= latsign * lonsign;