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