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