00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040 #include "GeographicLib/TransverseMercatorExact.hpp"
00041
00042 #define GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_CPP "$Id: TransverseMercatorExact.cpp 6813 2010-02-03 13:57:19Z karney $"
00043
00044 RCSID_DECL(GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_CPP)
00045 RCSID_DECL(GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_HPP)
00046
00047 namespace GeographicLib {
00048
00049 using namespace std;
00050
00051 const Math::real TransverseMercatorExact::tol =
00052 numeric_limits<real>::epsilon();
00053 const Math::real TransverseMercatorExact::tol1 = real(0.1) * sqrt(tol);
00054 const Math::real TransverseMercatorExact::tol2 = real(0.1) * tol;
00055 const Math::real TransverseMercatorExact::taytol = pow(tol, real(0.6));
00056
00057 const Math::real TransverseMercatorExact::overflow = 1 / sq(tol);
00058
00059 TransverseMercatorExact::TransverseMercatorExact(real a, real r, real k0,
00060 bool extendp)
00061 : _a(a)
00062 , _r(r)
00063 , _f(1 / _r)
00064 , _k0(k0)
00065 , _mu(_f * (2 - _f))
00066 , _mv(1 - _mu)
00067 , _e(sqrt(_mu))
00068 , _ep2(_mu / _mv)
00069 , _extendp(extendp)
00070 , _Eu(_mu)
00071 , _Ev(_mv)
00072 {
00073 if (!(_a > 0))
00074 throw GeographicErr("Major radius is not positive");
00075 if (!(_r > 0))
00076 throw GeographicErr("Inverse flattening is not positive");
00077 if (!(_k0 > 0))
00078 throw GeographicErr("Scale is not positive");
00079 }
00080
00081 const TransverseMercatorExact
00082 TransverseMercatorExact::UTM(Constants::WGS84_a(), Constants::WGS84_r(),
00083 Constants::UTM_k0());
00084
00085
00086 Math::real TransverseMercatorExact::taup(real tau) const throw() {
00087 real
00088 tau1 = Math::hypot(real(1), tau),
00089 sig = sinh( _e * Math::atanh(_e * tau / tau1) ),
00090 sig1 = Math::hypot(real(1), sig);
00091 return tau * sig1 - sig * tau1;
00092 }
00093
00094 Math::real TransverseMercatorExact::taupinv(real taup) const throw() {
00095 real tau = taup;
00096 for (int i = 0; i < numit; ++i) {
00097 real
00098 tau1 = Math::hypot(real(1), tau),
00099 sig = sinh( _e * Math::atanh(_e * tau / tau1 ) ),
00100 sig1 = Math::hypot(real(1), sig),
00101 dtau = - (sig1 * tau - sig * tau1 - taup) * (1 + _mv * sq(tau)) /
00102 ( (sig1 * tau1 - sig * tau) * _mv * tau1 );
00103 tau += dtau;
00104 if (abs(dtau) < tol * max(real(1), tau))
00105 break;
00106 }
00107 return tau;
00108 }
00109
00110 void TransverseMercatorExact::zeta(real u, real snu, real cnu, real dnu,
00111 real v, real snv, real cnv, real dnv,
00112 real& taup, real& lam) const throw() {
00113
00114
00115
00116
00117 real
00118 d1 = sqrt(sq(cnu) + _mv * sq(snu * snv)),
00119 d2 = sqrt(_mu * sq(cnu) + _mv * sq(cnv)),
00120 t1 = (d1 ? snu * dnv / d1 : snu < 0 ? -overflow : overflow),
00121 t2 = (d2 ? sinh( _e * Math::asinh(_e * snu / d2) ) :
00122 snu < 0 ? -overflow : overflow);
00123
00124
00125 taup = t1 * Math::hypot(real(1), t2) - t2 * Math::hypot(real(1), t1);
00126 lam = (d1 != 0 && d2 != 0) ?
00127 atan2(dnu * snv, cnu * cnv) - _e * atan2(_e * cnu * snv, dnu * cnv) :
00128 0;
00129 }
00130
00131 void TransverseMercatorExact::dwdzeta(real u, real snu, real cnu, real dnu,
00132 real v, real snv, real cnv, real dnv,
00133 real& du, real& dv) const throw() {
00134
00135
00136 real d = _mv * sq(sq(cnv) + _mu * sq(snu * snv));
00137 du = cnu * dnu * dnv * (sq(cnv) - _mu * sq(snu * snv)) / d;
00138 dv = -snu * snv * cnv * (sq(dnu * dnv) + _mu * sq(cnu)) / d;
00139 }
00140
00141
00142 bool TransverseMercatorExact::zetainv0(real psi, real lam, real& u, real& v)
00143 const throw() {
00144 bool retval = false;
00145 if (psi < -_e * Constants::pi()/4 &&
00146 lam > (1 - 2 * _e) * Constants::pi()/2 &&
00147 psi < lam - (1 - _e) * Constants::pi()/2) {
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157 real
00158 psix = 1 - psi / _e,
00159 lamx = (Constants::pi()/2 - lam) / _e;
00160 u = Math::asinh(sin(lamx) / Math::hypot(cos(lamx), sinh(psix))) *
00161 (1 + _mu/2);
00162 v = atan2(cos(lamx), sinh(psix)) * (1 + _mu/2);
00163 u = _Eu.K() - u;
00164 v = _Ev.K() - v;
00165 } else if (psi < _e * Constants::pi()/2 &&
00166 lam > (1 - 2 * _e) * Constants::pi()/2) {
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178 real
00179 dlam = lam - (1 - _e) * Constants::pi()/2,
00180 rad = Math::hypot(psi, dlam),
00181
00182
00183
00184
00185
00186 ang = atan2(dlam-psi, psi+dlam) - real(0.75) * Constants::pi();
00187
00188 retval = rad < _e * taytol;
00189 rad = Math::cbrt(3 / (_mv * _e) * rad);
00190 ang /= 3;
00191 u = rad * cos(ang);
00192 v = rad * sin(ang) + _Ev.K();
00193 } else {
00194
00195
00196
00197 v = Math::asinh(sin(lam) / Math::hypot(cos(lam), sinh(psi)));
00198 u = atan2(sinh(psi), cos(lam));
00199
00200 u *= _Eu.K() / (Constants::pi()/2);
00201 v *= _Eu.K() / (Constants::pi()/2);
00202 }
00203 return retval;
00204 }
00205
00206
00207 void TransverseMercatorExact::zetainv(real taup, real lam, real& u, real& v)
00208 const throw() {
00209 real
00210 psi = Math::asinh(taup),
00211 scal = 1/Math::hypot(real(1), taup);
00212 if (zetainv0(psi, lam, u, v))
00213 return;
00214 real stol2 = tol2 / sq(max(psi, real(1)));
00215
00216 for (int i = 0, trip = 0; i < numit; ++i) {
00217 real snu, cnu, dnu, snv, cnv, dnv;
00218 _Eu.sncndn(u, snu, cnu, dnu);
00219 _Ev.sncndn(v, snv, cnv, dnv);
00220 real tau1, lam1, du1, dv1;
00221 zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau1, lam1);
00222 dwdzeta(u, snu, cnu, dnu, v, snv, cnv, dnv, du1, dv1);
00223 tau1 -= taup;
00224 lam1 -= lam;
00225 tau1 *= scal;
00226 real
00227 delu = tau1 * du1 - lam1 * dv1,
00228 delv = tau1 * dv1 + lam1 * du1;
00229 u -= delu;
00230 v -= delv;
00231 if (trip)
00232 break;
00233 real delw2 = sq(delu) + sq(delv);
00234 if (delw2 < stol2)
00235 ++trip;
00236 }
00237 }
00238
00239 void TransverseMercatorExact::sigma(real u, real snu, real cnu, real dnu,
00240 real v, real snv, real cnv, real dnv,
00241 real& xi, real& eta) const throw() {
00242
00243
00244 real d = _mu * sq(cnu) + _mv * sq(cnv);
00245 xi = _Eu.E(snu, cnu, dnu) - _mu * snu * cnu * dnu / d;
00246 eta = v - _Ev.E(snv, cnv, dnv) + _mv * snv * cnv * dnv / d;
00247 }
00248
00249 void TransverseMercatorExact::dwdsigma(real u, real snu, real cnu, real dnu,
00250 real v, real snv, real cnv, real dnv,
00251 real& du, real& dv) const throw() {
00252
00253
00254 real d = _mv * sq(sq(cnv) + _mu * sq(snu * snv));
00255 real
00256 dnr = dnu * cnv * dnv,
00257 dni = - _mu * snu * cnu * snv;
00258 du = (sq(dnr) - sq(dni)) / d;
00259 dv = 2 * dnr * dni / d;
00260 }
00261
00262
00263 bool TransverseMercatorExact::sigmainv0(real xi, real eta, real& u, real& v)
00264 const throw() {
00265 bool retval = false;
00266 if (eta > real(1.25) * _Ev.KE() ||
00267 (xi < -real(0.25) * _Eu.E() && xi < eta - _Ev.KE())) {
00268
00269
00270
00271
00272 real
00273 x = xi - _Eu.E(),
00274 y = eta - _Ev.KE(),
00275 r2 = sq(x) + sq(y);
00276 u = _Eu.K() + x/r2;
00277 v = _Ev.K() - y/r2;
00278 } else if ((eta > real(0.75) * _Ev.KE() && xi < real(0.25) * _Eu.E())
00279 || eta > _Ev.KE()) {
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292 real
00293 deta = eta - _Ev.KE(),
00294 rad = Math::hypot(xi, deta),
00295
00296
00297 ang = atan2(deta-xi, xi+deta) - real(0.75) * Constants::pi();
00298
00299 retval = rad < 2 * taytol;
00300 rad = Math::cbrt(3 / _mv * rad);
00301 ang /= 3;
00302 u = rad * cos(ang);
00303 v = rad * sin(ang) + _Ev.K();
00304 } else {
00305
00306 u = xi * _Eu.K()/_Eu.E();
00307 v = eta * _Eu.K()/_Eu.E();
00308 }
00309 return retval;
00310 }
00311
00312
00313 void TransverseMercatorExact::sigmainv(real xi, real eta, real& u, real& v)
00314 const throw() {
00315 if (sigmainv0(xi, eta, u, v))
00316 return;
00317
00318 for (int i = 0, trip = 0; i < numit; ++i) {
00319 real snu, cnu, dnu, snv, cnv, dnv;
00320 _Eu.sncndn(u, snu, cnu, dnu);
00321 _Ev.sncndn(v, snv, cnv, dnv);
00322 real xi1, eta1, du1, dv1;
00323 sigma(u, snu, cnu, dnu, v, snv, cnv, dnv, xi1, eta1);
00324 dwdsigma(u, snu, cnu, dnu, v, snv, cnv, dnv, du1, dv1);
00325 xi1 -= xi;
00326 eta1 -= eta;
00327 real
00328 delu = xi1 * du1 - eta1 * dv1,
00329 delv = xi1 * dv1 + eta1 * du1;
00330 u -= delu;
00331 v -= delv;
00332 if (trip)
00333 break;
00334 real delw2 = sq(delu) + sq(delv);
00335 if (delw2 < tol2)
00336 ++trip;
00337 }
00338 }
00339
00340 void TransverseMercatorExact::Scale(real tau, real lam,
00341 real snu, real cnu, real dnu,
00342 real snv, real cnv, real dnv,
00343 real& gamma, real& k) const throw() {
00344 real sec2 = 1 + sq(tau);
00345
00346
00347 gamma = atan2(_mv * snu * snv * cnv, cnu * dnu * dnv);
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361 k = sqrt(_mv + _mu / sec2) * sqrt(sec2) *
00362 sqrt( (_mv * sq(snv) + sq(cnu * dnv)) /
00363 (_mu * sq(cnu) + _mv * sq(cnv)) );
00364 }
00365
00366 void TransverseMercatorExact::Forward(real lon0, real lat, real lon,
00367 real& x, real& y, real& gamma, real& k)
00368 const throw() {
00369
00370 if (lon - lon0 > 180)
00371 lon -= lon0 - 360;
00372 else if (lon - lon0 <= -180)
00373 lon -= lon0 + 360;
00374 else
00375 lon -= lon0;
00376
00377
00378 int
00379 latsign = !_extendp && lat < 0 ? -1 : 1,
00380 lonsign = !_extendp && lon < 0 ? -1 : 1;
00381 lon *= lonsign;
00382 lat *= latsign;
00383 bool backside = !_extendp && lon > 90;
00384 if (backside) {
00385 if (lat == 0)
00386 latsign = -1;
00387 lon = 180 - lon;
00388 }
00389 real
00390 phi = lat * Constants::degree(),
00391 lam = lon * Constants::degree(),
00392 tau = tanx(phi);
00393
00394
00395 real u, v;
00396 if (lat == 90) {
00397 u = _Eu.K();
00398 v = 0;
00399 } else if (lat == 0 && lon == 90 * (1 - _e)) {
00400 u = 0;
00401 v = _Ev.K();
00402 } else
00403 zetainv(taup(tau), lam, u, v);
00404
00405 real snu, cnu, dnu, snv, cnv, dnv;
00406 _Eu.sncndn(u, snu, cnu, dnu);
00407 _Ev.sncndn(v, snv, cnv, dnv);
00408
00409 real xi, eta;
00410 sigma(u, snu, cnu, dnu, v, snv, cnv, dnv, xi, eta);
00411 if (backside)
00412 xi = 2 * _Eu.E() - xi;
00413 y = xi * _a * _k0 * latsign;
00414 x = eta * _a * _k0 * lonsign;
00415
00416
00417 zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau, lam);
00418 tau=taupinv(tau);
00419 Scale(tau, lam, snu, cnu, dnu, snv, cnv, dnv, gamma, k);
00420 gamma /= Constants::degree();
00421 if (backside)
00422 gamma = 180 - gamma;
00423 gamma *= latsign * lonsign;
00424 k *= _k0;
00425 }
00426
00427 void TransverseMercatorExact::Reverse(real lon0, real x, real y,
00428 real& lat, real& lon,
00429 real& gamma, real& k)
00430 const throw() {
00431
00432 real
00433 xi = y / (_a * _k0),
00434 eta = x / (_a * _k0);
00435
00436 int
00437 latsign = !_extendp && y < 0 ? -1 : 1,
00438 lonsign = !_extendp && x < 0 ? -1 : 1;
00439 xi *= latsign;
00440 eta *= lonsign;
00441 bool backside = !_extendp && xi > _Eu.E();
00442 if (backside)
00443 xi = 2 * _Eu.E()- xi;
00444
00445
00446 real u, v;
00447 if (xi == 0 && eta == _Ev.KE()) {
00448 u = 0;
00449 v = _Ev.K();
00450 } else
00451 sigmainv(xi, eta, u, v);
00452
00453 real snu, cnu, dnu, snv, cnv, dnv;
00454 _Eu.sncndn(u, snu, cnu, dnu);
00455 _Ev.sncndn(v, snv, cnv, dnv);
00456 real phi, lam, tau;
00457 if (v != 0 || u != _Eu.K()) {
00458 zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau, lam);
00459 tau = taupinv(tau);
00460 phi = atan(tau);
00461 lat = phi / Constants::degree();
00462 lon = lam / Constants::degree();
00463 } else {
00464 tau = overflow;
00465 phi = Constants::pi()/2;
00466 lat = 90;
00467 lon = lam = 0;
00468 }
00469 Scale(tau, lam, snu, cnu, dnu, snv, cnv, dnv, gamma, k);
00470 gamma /= Constants::degree();
00471 if (backside)
00472 lon = 180 - lon;
00473 lon *= lonsign;
00474
00475 if (lon + lon0 >= 180)
00476 lon += lon0 - 360;
00477 else if (lon + lon0 < -180)
00478 lon += lon0 + 360;
00479 else
00480 lon += lon0;
00481 lat *= latsign;
00482 if (backside)
00483 y = 2 * _Eu.E() - y;
00484 y *= _a * _k0 * latsign;
00485 x *= _a * _k0 * lonsign;
00486 if (backside)
00487 gamma = 180 - gamma;
00488 gamma *= latsign * lonsign;
00489 k *= _k0;
00490 }
00491
00492 }