44 namespace GeographicLib {
48 const Math::real TransverseMercatorExact::tol_ =
49 numeric_limits<real>::epsilon();
50 const Math::real TransverseMercatorExact::tol1_ =
real(0.1) * sqrt(tol_);
51 const Math::real TransverseMercatorExact::tol2_ =
real(0.1) * tol_;
52 const Math::real TransverseMercatorExact::taytol_ = pow(tol_,
real(0.6));
59 , _f(f <= 1 ? f : 1/f)
80 Constants::WGS84_f<real>(),
81 Constants::UTM_k0<real>());
91 Math::real TransverseMercatorExact::taupinv(
real taup)
const throw() {
95 stol = tol_ * max(
real(1), abs(taup));
97 for (
int i = 0; i < numit_; ++i) {
102 dtau = (taup - taupa) * (1 + _mv *
Math::sq(tau)) /
105 if (!(abs(dtau) >= stol))
113 real& taup,
real& lam)
const throw() {
121 t1 = (d1 ? snu * dnv / d1 : (snu < 0 ? -overflow_ : overflow_)),
122 t2 = (d2 ? sinh( _e *
Math::asinh(_e * snu / d2) ) :
123 (snu < 0 ? -overflow_ : overflow_));
127 lam = (d1 != 0 && d2 != 0) ?
128 atan2(dnu * snv, cnu * cnv) - _e * atan2(_e * cnu * snv, dnu * cnv) :
132 void TransverseMercatorExact::dwdzeta(
real ,
145 bool TransverseMercatorExact::zetainv0(
real psi,
real lam,
real& u,
real& v)
148 if (psi < -_e * Math::pi<real>()/4 &&
149 lam > (1 - 2 * _e) * Math::pi<real>()/2 &&
150 psi < lam - (1 - _e) * Math::pi<real>()/2) {
162 lamx = (Math::pi<real>()/2 - lam) / _e;
165 v = atan2(cos(lamx), sinh(psix)) * (1 + _mu/2);
168 }
else if (psi < _e * Math::pi<real>()/2 &&
169 lam > (1 - 2 * _e) * Math::pi<real>()/2) {
182 dlam = lam - (1 - _e) * Math::pi<real>()/2,
189 ang = atan2(dlam-psi, psi+dlam) -
real(0.75) * Math::pi<real>();
191 retval = rad < _e * taytol_;
195 v = rad * sin(ang) + _Ev.K();
201 u = atan2(sinh(psi), cos(lam));
203 u *= _Eu.K() / (Math::pi<real>()/2);
204 v *= _Eu.K() / (Math::pi<real>()/2);
210 void TransverseMercatorExact::zetainv(
real taup,
real lam,
real& u,
real& v)
215 if (zetainv0(psi, lam, u, v))
219 for (
int i = 0, trip = 0; i < numit_; ++i) {
220 real snu, cnu, dnu, snv, cnv, dnv;
221 _Eu.sncndn(u, snu, cnu, dnu);
222 _Ev.sncndn(v, snv, cnv, dnv);
223 real tau1, lam1, du1, dv1;
224 zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau1, lam1);
225 dwdzeta(u, snu, cnu, dnu, v, snv, cnv, dnv, du1, dv1);
230 delu = tau1 * du1 - lam1 * dv1,
231 delv = tau1 * dv1 + lam1 * du1;
237 if (!(delw2 >= stol2))
244 real& xi,
real& eta)
const throw() {
248 xi = _Eu.E(snu, cnu, dnu) - _mu * snu * cnu * dnu / d;
249 eta = v - _Ev.E(snv, cnv, dnv) + _mv * snv * cnv * dnv / d;
252 void TransverseMercatorExact::dwdsigma(
real ,
261 dnr = dnu * cnv * dnv,
262 dni = - _mu * snu * cnu * snv;
264 dv = 2 * dnr * dni / d;
268 bool TransverseMercatorExact::sigmainv0(
real xi,
real eta,
real& u,
real& v)
271 if (eta >
real(1.25) * _Ev.KE() ||
272 (xi < -
real(0.25) * _Eu.E() && xi < eta - _Ev.KE())) {
283 }
else if ((eta >
real(0.75) * _Ev.KE() && xi <
real(0.25) * _Eu.E())
298 deta = eta - _Ev.KE(),
302 ang = atan2(deta-xi, xi+deta) -
real(0.75) * Math::pi<real>();
304 retval = rad < 2 * taytol_;
308 v = rad * sin(ang) + _Ev.K();
311 u = xi * _Eu.K()/_Eu.E();
312 v = eta * _Eu.K()/_Eu.E();
320 if (sigmainv0(xi, eta, u, v))
323 for (
int i = 0, trip = 0; i < numit_; ++i) {
324 real snu, cnu, dnu, snv, cnv, dnv;
325 _Eu.sncndn(u, snu, cnu, dnu);
326 _Ev.sncndn(v, snv, cnv, dnv);
327 real xi1, eta1, du1, dv1;
328 sigma(u, snu, cnu, dnu, v, snv, cnv, dnv, xi1, eta1);
329 dwdsigma(u, snu, cnu, dnu, v, snv, cnv, dnv, du1, dv1);
333 delu = xi1 * du1 - eta1 * dv1,
334 delv = xi1 * dv1 + eta1 * du1;
340 if (!(delw2 >= tol2_))
345 void TransverseMercatorExact::Scale(
real tau,
real ,
348 real& gamma,
real& k)
const throw() {
352 gamma = atan2(_mv * snu * snv * cnv, cnu * dnu * dnv);
366 k = sqrt(_mv + _mu / sec2) * sqrt(sec2) *
372 real& x, real& y, real& gamma, real& k)
377 latsign = (!_extendp && lat < 0) ? -1 : 1,
378 lonsign = (!_extendp && lon < 0) ? -1 : 1;
381 bool backside = !_extendp && lon > 90;
388 phi = lat * Math::degree<real>(),
389 lam = lon * Math::degree<real>(),
397 }
else if (lat == 0 && lon == 90 * (1 - _e)) {
401 zetainv(taup(tau), lam, u, v);
403 real snu, cnu, dnu, snv, cnv, dnv;
404 _Eu.sncndn(u, snu, cnu, dnu);
405 _Ev.sncndn(v, snv, cnv, dnv);
408 sigma(u, snu, cnu, dnu, v, snv, cnv, dnv, xi, eta);
410 xi = 2 * _Eu.E() - xi;
411 y = xi * _a * _k0 * latsign;
412 x = eta * _a * _k0 * lonsign;
419 zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau, lam);
421 Scale(tau, lam, snu, cnu, dnu, snv, cnv, dnv, gamma, k);
422 gamma /= Math::degree<real>();
426 gamma *= latsign * lonsign;
431 real& lat, real& lon,
432 real& gamma, real& k)
437 eta = x / (_a * _k0);
440 latsign = !_extendp && y < 0 ? -1 : 1,
441 lonsign = !_extendp && x < 0 ? -1 : 1;
444 bool backside = !_extendp && xi > _Eu.E();
446 xi = 2 * _Eu.E()- xi;
450 if (xi == 0 && eta == _Ev.KE()) {
454 sigmainv(xi, eta, u, v);
456 real snu, cnu, dnu, snv, cnv, dnv;
457 _Eu.sncndn(u, snu, cnu, dnu);
458 _Ev.sncndn(v, snv, cnv, dnv);
460 if (v != 0 || u != _Eu.K()) {
461 zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau, lam);
464 lat = phi / Math::degree<real>();
465 lon = lam / Math::degree<real>();
466 Scale(tau, lam, snu, cnu, dnu, snv, cnv, dnv, gamma, k);
467 gamma /= Math::degree<real>();
470 phi = Math::pi<real>()/2;
472 lon = lam = gamma = 0;
483 y *= _a * _k0 * latsign;
484 x *= _a * _k0 * lonsign;
487 gamma *= latsign * lonsign;
static T AngNormalize(T x)
An exact implementation of the transverse Mercator projection.
GeographicLib::Math::real real
static bool isfinite(T x)
Header for GeographicLib::TransverseMercatorExact class.
TransverseMercatorExact(real a, real f, real k0, bool extendp=false)
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
static T AngDiff(T x, T y)
Exception handling for GeographicLib.
static const TransverseMercatorExact UTM