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/TransverseMercator.hpp"
00043
00044 #define GEOGRAPHICLIB_TRANSVERSEMERCATOR_CPP "$Id: TransverseMercator.cpp 6807 2010-02-01 11:26:34Z karney $"
00045
00046 RCSID_DECL(GEOGRAPHICLIB_TRANSVERSEMERCATOR_CPP)
00047 RCSID_DECL(GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP)
00048
00049 namespace GeographicLib {
00050
00051 using namespace std;
00052
00053 const Math::real TransverseMercator::tol =
00054 real(0.1)*sqrt(numeric_limits<real>::epsilon());
00055
00056 const Math::real TransverseMercator::overflow =
00057 1 / sq(numeric_limits<real>::epsilon());
00058
00059 TransverseMercator::TransverseMercator(real a, real r, real k0)
00060 : _a(a)
00061 , _r(r)
00062 , _f(_r != 0 ? 1 / _r : 0)
00063 , _k0(k0)
00064 , _e2(_f * (2 - _f))
00065 , _e(sqrt(abs(_e2)))
00066 , _e2m(1 - _e2)
00067
00068
00069 , _c( sqrt(_e2m) * exp(eatanhe(real(1))) )
00070 , _n(_f / (2 - _f))
00071 {
00072 if (!(_a > 0))
00073 throw GeographicErr("Major radius is not positive");
00074 if (!(_k0 > 0))
00075 throw GeographicErr("Scale is not positive");
00076
00077
00078 real nx = sq(_n);
00079 switch (maxpow) {
00080 case 4:
00081 _b1 = 1/(1+_n)*(nx*(nx+16)+64)/64;
00082 _alp[1] = _n*(_n*(_n*(164*_n+225)-480)+360)/720;
00083 _bet[1] = _n*(_n*((555-4*_n)*_n-960)+720)/1440;
00084 _alp[2] = nx*(_n*(557*_n-864)+390)/1440;
00085 _bet[2] = nx*((96-437*_n)*_n+30)/1440;
00086 nx *= _n;
00087 _alp[3] = (427-1236*_n)*nx/1680;
00088 _bet[3] = (119-148*_n)*nx/3360;
00089 nx *= _n;
00090 _alp[4] = 49561*nx/161280;
00091 _bet[4] = 4397*nx/161280;
00092 break;
00093 case 5:
00094 _b1 = 1/(1+_n)*(nx*(nx+16)+64)/64;
00095 _alp[1] = _n*(_n*(_n*((328-635*_n)*_n+450)-960)+720)/1440;
00096 _bet[1] = _n*(_n*(_n*((-3645*_n-64)*_n+8880)-15360)+11520)/23040;
00097 _alp[2] = nx*(_n*(_n*(4496*_n+3899)-6048)+2730)/10080;
00098 _bet[2] = nx*(_n*(_n*(4416*_n-3059)+672)+210)/10080;
00099 nx *= _n;
00100 _alp[3] = nx*(_n*(15061*_n-19776)+6832)/26880;
00101 _bet[3] = nx*((-627*_n-592)*_n+476)/13440;
00102 nx *= _n;
00103 _alp[4] = (49561-171840*_n)*nx/161280;
00104 _bet[4] = (4397-3520*_n)*nx/161280;
00105 nx *= _n;
00106 _alp[5] = 34729*nx/80640;
00107 _bet[5] = 4583*nx/161280;
00108 break;
00109 case 6:
00110 _b1 = 1/(1+_n)*(nx*(nx*(nx+4)+64)+256)/256;
00111 _alp[1] = _n*(_n*(_n*(_n*(_n*(31564*_n-66675)+34440)+47250)-100800)+
00112 75600)/151200;
00113 _bet[1] = _n*(_n*(_n*(_n*(_n*(384796*_n-382725)-6720)+932400)-1612800)+
00114 1209600)/2419200;
00115 _alp[2] = nx*(_n*(_n*((863232-1983433*_n)*_n+748608)-1161216)+524160)/
00116 1935360;
00117 _bet[2] = nx*(_n*(_n*((1695744-1118711*_n)*_n-1174656)+258048)+80640)/
00118 3870720;
00119 nx *= _n;
00120 _alp[3] = nx*(_n*(_n*(670412*_n+406647)-533952)+184464)/725760;
00121 _bet[3] = nx*(_n*(_n*(22276*_n-16929)-15984)+12852)/362880;
00122 nx *= _n;
00123 _alp[4] = nx*(_n*(6601661*_n-7732800)+2230245)/7257600;
00124 _bet[4] = nx*((-830251*_n-158400)*_n+197865)/7257600;
00125 nx *= _n;
00126 _alp[5] = (3438171-13675556*_n)*nx/7983360;
00127 _bet[5] = (453717-435388*_n)*nx/15966720;
00128 nx *= _n;
00129 _alp[6] = 212378941*nx/319334400;
00130 _bet[6] = 20648693*nx/638668800;
00131 break;
00132 case 7:
00133 _b1 = 1/(1+_n)*(nx*(nx*(nx+4)+64)+256)/256;
00134 _alp[1] = _n*(_n*(_n*(_n*(_n*(_n*(1804025*_n+2020096)-4267200)+2204160)+
00135 3024000)-6451200)+4838400)/9676800;
00136 _bet[1] = _n*(_n*(_n*(_n*(_n*((6156736-5406467*_n)*_n-6123600)-107520)+
00137 14918400)-25804800)+19353600)/38707200;
00138 _alp[2] = nx*(_n*(_n*(_n*(_n*(4626384*_n-9917165)+4316160)+3743040)-
00139 5806080)+2620800)/9676800;
00140 _bet[2] = nx*(_n*(_n*(_n*(_n*(829456*_n-5593555)+8478720)-5873280)+
00141 1290240)+403200)/19353600;
00142 nx *= _n;
00143 _alp[3] = nx*(_n*(_n*((26816480-67102379*_n)*_n+16265880)-21358080)+
00144 7378560)/29030400;
00145 _bet[3] = nx*(_n*(_n*(_n*(9261899*_n+3564160)-2708640)-2557440)+
00146 2056320)/58060800;
00147 nx *= _n;
00148 _alp[4] = nx*(_n*(_n*(155912000*_n+72618271)-85060800)+24532695)/
00149 79833600;
00150 _bet[4] = nx*(_n*(_n*(14928352*_n-9132761)-1742400)+2176515)/79833600;
00151 nx *= _n;
00152 _alp[5] = nx*(_n*(102508609*_n-109404448)+27505368)/63866880;
00153 _bet[5] = nx*((-8005831*_n-1741552)*_n+1814868)/63866880;
00154 nx *= _n;
00155 _alp[6] = (2760926233.0-12282192400.0*_n)*nx/4151347200.0;
00156 _bet[6] = (268433009-261810608*_n)*nx/8302694400.0;
00157 nx *= _n;
00158 _alp[7] = 1522256789.0*nx/1383782400.0;
00159 _bet[7] = 219941297*nx/5535129600.0;
00160 break;
00161 case 8:
00162 _b1 = 1/(1+_n)*(nx*(nx*(nx*(25*nx+64)+256)+4096)+16384)/16384;
00163 _alp[1] = _n*(_n*(_n*(_n*(_n*(_n*((37884525-75900428*_n)*_n+42422016)-
00164 89611200)+46287360)+63504000)-135475200)+
00165 101606400)/203212800;
00166 _bet[1] = _n*(_n*(_n*(_n*(_n*(_n*(_n*(31777436*_n-37845269)+43097152)-
00167 42865200)-752640)+104428800)-180633600)+
00168 135475200)/270950400;
00169 _alp[2] = nx*(_n*(_n*(_n*(_n*(_n*(148003883*_n+83274912)-178508970)+
00170 77690880)+67374720)-104509440)+47174400)/
00171 174182400;
00172 _bet[2] = nx*(_n*(_n*(_n*(_n*(_n*(24749483*_n+14930208)-100683990)+
00173 152616960)-105719040)+23224320)+7257600)/
00174 348364800;
00175 nx *= _n;
00176 _alp[3] = nx*(_n*(_n*(_n*(_n*(318729724*_n-738126169)+294981280)+
00177 178924680)-234938880)+81164160)/319334400;
00178 _bet[3] = nx*(_n*(_n*(_n*((101880889-232468668*_n)*_n+39205760)-
00179 29795040)-28131840)+22619520)/638668800;
00180 nx *= _n;
00181 _alp[4] = nx*(_n*(_n*((14967552000.0-40176129013.0*_n)*_n+6971354016.0)-
00182 8165836800.0)+2355138720.0)/7664025600.0;
00183 _bet[4] = nx*(_n*(_n*(_n*(324154477*_n+1433121792.0)-876745056)-
00184 167270400)+208945440)/7664025600.0;
00185 nx *= _n;
00186 _alp[5] = nx*(_n*(_n*(10421654396.0*_n+3997835751.0)-4266773472.0)+
00187 1072709352.0)/2490808320.0;
00188 _bet[5] = nx*(_n*(_n*(457888660*_n-312227409)-67920528)+70779852)/
00189 2490808320.0;
00190 nx *= _n;
00191 _alp[6] = nx*(_n*(175214326799.0*_n-171950693600.0)+38652967262.0)/
00192 58118860800.0;
00193 _bet[6] = nx*((-19841813847.0*_n-3665348512.0)*_n+3758062126.0)/
00194 116237721600.0;
00195 nx *= _n;
00196 _alp[7] = (13700311101.0-67039739596.0*_n)*nx/12454041600.0;
00197 _bet[7] = (1979471673.0-1989295244.0*_n)*nx/49816166400.0;
00198 nx *= _n;
00199 _alp[8] = 1424729850961.0*nx/743921418240.0;
00200 _bet[8] = 191773887257.0*nx/3719607091200.0;
00201 break;
00202 default:
00203 STATIC_ASSERT(maxpow >= 4 && maxpow <= 8, "Bad value of maxpow");
00204 }
00205
00206
00207 _a1 = _b1 * _a;
00208 }
00209
00210 const TransverseMercator
00211 TransverseMercator::UTM(Constants::WGS84_a(), Constants::WGS84_r(),
00212 Constants::UTM_k0());
00213
00214 void TransverseMercator::Forward(real lon0, real lat, real lon,
00215 real& x, real& y, real& gamma, real& k)
00216 const throw() {
00217
00218 if (lon - lon0 > 180)
00219 lon -= lon0 - 360;
00220 else if (lon - lon0 <= -180)
00221 lon -= lon0 + 360;
00222 else
00223 lon -= lon0;
00224
00225
00226 int
00227 latsign = lat < 0 ? -1 : 1,
00228 lonsign = lon < 0 ? -1 : 1;
00229 lon *= lonsign;
00230 lat *= latsign;
00231 bool backside = lon > 90;
00232 if (backside) {
00233 if (lat == 0)
00234 latsign = -1;
00235 lon = 180 - lon;
00236 }
00237 real
00238 phi = lat * Constants::degree(),
00239 lam = lon * Constants::degree();
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257 real etap, xip;
00258 if (lat < 90) {
00259 real
00260 c = max(real(0), cos(lam)),
00261 tau = tanx(phi),
00262 secphi = Math::hypot(real(1), tau),
00263 sig = sinh( eatanhe(sin(phi)) ),
00264 taup = (Math::hypot(real(1), sig) * tau - sig * secphi);
00265 xip = atan2(taup, c);
00266
00267
00268 etap = Math::asinh(sin(lam) / Math::hypot(taup, c));
00269
00270
00271
00272 gamma = atan(tanx(lam) *
00273 taup / Math::hypot(real(1), taup));
00274
00275
00276
00277
00278
00279
00280
00281 k = sqrt(_e2m + _e2 * sq(cos(phi))) * secphi / Math::hypot(taup, c);
00282 } else {
00283 xip = Constants::pi()/2;
00284 etap = 0;
00285 gamma = lam;
00286 k = _c;
00287 }
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341 real
00342 c0 = cos(2 * xip), ch0 = cosh(2 * etap),
00343 s0 = sin(2 * xip), sh0 = sinh(2 * etap),
00344 ar = 2 * c0 * ch0, ai = -2 * s0 * sh0;
00345 int n = maxpow;
00346 real
00347 xi0 = (n & 1 ? _alp[n] : 0), eta0 = 0,
00348 xi1 = 0, eta1 = 0;
00349 real
00350 yr0 = (n & 1 ? 2 * maxpow * _alp[n--] : 0), yi0 = 0,
00351 yr1 = 0, yi1 = 0;
00352 while (n) {
00353 xi1 = ar * xi0 - ai * eta0 - xi1 + _alp[n];
00354 eta1 = ai * xi0 + ar * eta0 - eta1;
00355 yr1 = ar * yr0 - ai * yi0 - yr1 + 2 * n * _alp[n];
00356 yi1 = ai * yr0 + ar * yi0 - yi1;
00357 --n;
00358 xi0 = ar * xi1 - ai * eta1 - xi0 + _alp[n];
00359 eta0 = ai * xi1 + ar * eta1 - eta0;
00360 yr0 = ar * yr1 - ai * yi1 - yr0 + 2 * n * _alp[n];
00361 yi0 = ai * yr1 + ar * yi1 - yi0;
00362 --n;
00363 }
00364 ar /= 2; ai /= 2;
00365 yr1 = 1 - yr1 + ar * yr0 - ai * yi0;
00366 yi1 = - yi1 + ai * yr0 + ar * yi0;
00367 ar = s0 * ch0; ai = c0 * sh0;
00368 real
00369 xi = xip + ar * xi0 - ai * eta0,
00370 eta = etap + ai * xi0 + ar * eta0;
00371
00372
00373 gamma -= atan2(yi1, yr1);
00374 k *= _b1 * Math::hypot(yr1, yi1);
00375 gamma /= Constants::degree();
00376 y = _a1 * _k0 * (backside ? Constants::pi() - xi : xi) * latsign;
00377 x = _a1 * _k0 * eta * lonsign;
00378 if (backside)
00379 gamma = 180 - gamma;
00380 gamma *= latsign * lonsign;
00381 k *= _k0;
00382 }
00383
00384 void TransverseMercator::Reverse(real lon0, real x, real y,
00385 real& lat, real& lon, real& gamma, real& k)
00386 const throw() {
00387
00388
00389
00390 real
00391 xi = y / (_a1 * _k0),
00392 eta = x / (_a1 * _k0);
00393
00394 int
00395 xisign = xi < 0 ? -1 : 1,
00396 etasign = eta < 0 ? -1 : 1;
00397 xi *= xisign;
00398 eta *= etasign;
00399 bool backside = xi > Constants::pi()/2;
00400 if (backside)
00401 xi = Constants::pi() - xi;
00402 real
00403 c0 = cos(2 * xi), ch0 = cosh(2 * eta),
00404 s0 = sin(2 * xi), sh0 = sinh(2 * eta),
00405 ar = 2 * c0 * ch0, ai = -2 * s0 * sh0;
00406 int n = maxpow;
00407 real
00408 xip0 = (n & 1 ? -_bet[n] : 0), etap0 = 0,
00409 xip1 = 0, etap1 = 0;
00410 real
00411 yr0 = (n & 1 ? - 2 * maxpow * _bet[n--] : 0), yi0 = 0,
00412 yr1 = 0, yi1 = 0;
00413 while (n) {
00414 xip1 = ar * xip0 - ai * etap0 - xip1 - _bet[n];
00415 etap1 = ai * xip0 + ar * etap0 - etap1;
00416 yr1 = ar * yr0 - ai * yi0 - yr1 - 2 * n * _bet[n];
00417 yi1 = ai * yr0 + ar * yi0 - yi1;
00418 --n;
00419 xip0 = ar * xip1 - ai * etap1 - xip0 - _bet[n];
00420 etap0 = ai * xip1 + ar * etap1 - etap0;
00421 yr0 = ar * yr1 - ai * yi1 - yr0 - 2 * n * _bet[n];
00422 yi0 = ai * yr1 + ar * yi1 - yi0;
00423 --n;
00424 }
00425 ar /= 2; ai /= 2;
00426 yr1 = 1 - yr1 + ar * yr0 - ai * yi0;
00427 yi1 = - yi1 + ai * yr0 + ar * yi0;
00428 ar = s0 * ch0; ai = c0 * sh0;
00429 real
00430 xip = xi + ar * xip0 - ai * etap0,
00431 etap = eta + ai * xip0 + ar * etap0;
00432
00433 gamma = atan2(yi1, yr1);
00434 k = _b1 / Math::hypot(yr1, yi1);
00435
00436
00437
00438
00439
00440 real lam, phi;
00441 real
00442 s = sinh(etap),
00443 c = max(real(0), cos(xip)),
00444 r = Math::hypot(s, c);
00445 if (r > 0) {
00446 lam = atan2(s, c);
00447
00448 real
00449 taup = sin(xip)/r,
00450 tau = taup;
00451
00452 for (int i = 0; i < numit; ++i) {
00453 real
00454 tau1 = Math::hypot(real(1), tau),
00455 sig = sinh( eatanhe( tau / tau1 ) ),
00456 sig1 = Math::hypot(real(1), sig),
00457 dtau = - (sig1 * tau - sig * tau1 - taup) * (1 + _e2m * sq(tau)) /
00458 ( (sig1 * tau1 - sig * tau) * _e2m * tau1 );
00459 tau += dtau;
00460 if (abs(dtau) < tol * max(real(1), tau))
00461 break;
00462 }
00463 phi = atan(tau);
00464 gamma += atan(tanx(xip) * tanh(etap));
00465
00466 k *= sqrt(_e2m + _e2 * sq(cos(phi))) * Math::hypot(real(1), tau) * r;
00467 } else {
00468 phi = Constants::pi()/2;
00469 lam = 0;
00470 k *= _c;
00471 }
00472 lat = phi / Constants::degree() * xisign;
00473 lon = lam / Constants::degree();
00474 if (backside)
00475 lon = 180 - lon;
00476 lon *= etasign;
00477
00478 if (lon + lon0 >= 180)
00479 lon += lon0 - 360;
00480 else if (lon + lon0 < -180)
00481 lon += lon0 + 360;
00482 else
00483 lon += lon0;
00484 gamma /= Constants::degree();
00485 if (backside)
00486 gamma = 180 - gamma;
00487 gamma *= xisign * etasign;
00488 k *= _k0;
00489 }
00490
00491 }