GeographicLib
1.21
|
00001 /** 00002 * \file TransverseMercator.cpp 00003 * \brief Implementation for GeographicLib::TransverseMercator class 00004 * 00005 * Copyright (c) Charles Karney (2008-2011) <charles@karney.com> and licensed 00006 * under the MIT/X11 License. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 * 00009 * This implementation follows closely 00010 * <a href="http://www.jhs-suositukset.fi/suomi/jhs154"> JHS 154, ETRS89 - 00011 * järjestelmään liittyvät karttaprojektiot, 00012 * tasokoordinaatistot ja karttalehtijako</a> (Map projections, plane 00013 * coordinates, and map sheet index for ETRS89), published by JUHTA, Finnish 00014 * Geodetic Institute, and the National Land Survey of Finland (2006). 00015 * 00016 * The relevant section is available as the 2008 PDF file 00017 * http://docs.jhs-suositukset.fi/jhs-suositukset/JHS154/JHS154_liite1.pdf 00018 * 00019 * This is a straight transcription of the formulas in this paper with the 00020 * following exceptions: 00021 * - use of 6th order series instead of 4th order series. This reduces the 00022 * error to about 5nm for the UTM range of coordinates (instead of 200nm), 00023 * with a speed penalty of only 1%; 00024 * - use Newton's method instead of plain iteration to solve for latitude in 00025 * terms of isometric latitude in the Reverse method; 00026 * - use of Horner's representation for evaluating polynomials and Clenshaw's 00027 * method for summing trigonometric series; 00028 * - several modifications of the formulas to improve the numerical accuracy; 00029 * - evaluating the convergence and scale using the expression for the 00030 * projection or its inverse. 00031 * 00032 * If the preprocessor variable TM_TX_MAXPOW is set to an integer between 4 and 00033 * 8, then this specifies the order of the series used for the forward and 00034 * reverse transformations. The default value is 6. (The series accurate to 00035 * 12th order is given in \ref tmseries.) 00036 * 00037 * Other equivalent implementations are given in 00038 * - http://www.ign.fr/DISPLAY/000/526/702/5267021/NTG_76.pdf 00039 * - http://www.lantmateriet.se/upload/filer/kartor/geodesi_gps_och_detaljmatning/geodesi/Formelsamling/Gauss_Conformal_Projection.pdf 00040 **********************************************************************/ 00041 00042 #include <GeographicLib/TransverseMercator.hpp> 00043 00044 #define GEOGRAPHICLIB_TRANSVERSEMERCATOR_CPP \ 00045 "$Id: 7b5a1854a015da061b8fdad0a4b35be7e06fcb9a $" 00046 00047 RCSID_DECL(GEOGRAPHICLIB_TRANSVERSEMERCATOR_CPP) 00048 RCSID_DECL(GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP) 00049 00050 namespace GeographicLib { 00051 00052 using namespace std; 00053 00054 const Math::real TransverseMercator::tol_ = 00055 real(0.1)*sqrt(numeric_limits<real>::epsilon()); 00056 // Overflow value s.t. atan(overflow_) = pi/2 00057 const Math::real TransverseMercator::overflow_ = 00058 1 / Math::sq(numeric_limits<real>::epsilon()); 00059 00060 TransverseMercator::TransverseMercator(real a, real f, real k0) 00061 : _a(a) 00062 , _f(f <= 1 ? f : 1/f) 00063 , _k0(k0) 00064 , _e2(_f * (2 - _f)) 00065 , _e(sqrt(abs(_e2))) 00066 , _e2m(1 - _e2) 00067 // _c = sqrt( pow(1 + _e, 1 + _e) * pow(1 - _e, 1 - _e) ) ) 00068 // See, for example, Lee (1976), p 100. 00069 , _c( sqrt(_e2m) * exp(eatanhe(real(1))) ) 00070 , _n(_f / (2 - _f)) 00071 { 00072 if (!(Math::isfinite(_a) && _a > 0)) 00073 throw GeographicErr("Major radius is not positive"); 00074 if (!(Math::isfinite(_f) && _f < 1)) 00075 throw GeographicErr("Minor radius is not positive"); 00076 if (!(Math::isfinite(_k0) && _k0 > 0)) 00077 throw GeographicErr("Scale is not positive"); 00078 // If coefficents might overflow_ an int, convert them to double (and they 00079 // are all exactly representable as doubles). 00080 real nx = Math::sq(_n); 00081 switch (maxpow_) { 00082 case 4: 00083 _b1 = 1/(1+_n)*(nx*(nx+16)+64)/64; 00084 _alp[1] = _n*(_n*(_n*(164*_n+225)-480)+360)/720; 00085 _bet[1] = _n*(_n*((555-4*_n)*_n-960)+720)/1440; 00086 _alp[2] = nx*(_n*(557*_n-864)+390)/1440; 00087 _bet[2] = nx*((96-437*_n)*_n+30)/1440; 00088 nx *= _n; 00089 _alp[3] = (427-1236*_n)*nx/1680; 00090 _bet[3] = (119-148*_n)*nx/3360; 00091 nx *= _n; 00092 _alp[4] = 49561*nx/161280; 00093 _bet[4] = 4397*nx/161280; 00094 break; 00095 case 5: 00096 _b1 = 1/(1+_n)*(nx*(nx+16)+64)/64; 00097 _alp[1] = _n*(_n*(_n*((328-635*_n)*_n+450)-960)+720)/1440; 00098 _bet[1] = _n*(_n*(_n*((-3645*_n-64)*_n+8880)-15360)+11520)/23040; 00099 _alp[2] = nx*(_n*(_n*(4496*_n+3899)-6048)+2730)/10080; 00100 _bet[2] = nx*(_n*(_n*(4416*_n-3059)+672)+210)/10080; 00101 nx *= _n; 00102 _alp[3] = nx*(_n*(15061*_n-19776)+6832)/26880; 00103 _bet[3] = nx*((-627*_n-592)*_n+476)/13440; 00104 nx *= _n; 00105 _alp[4] = (49561-171840*_n)*nx/161280; 00106 _bet[4] = (4397-3520*_n)*nx/161280; 00107 nx *= _n; 00108 _alp[5] = 34729*nx/80640; 00109 _bet[5] = 4583*nx/161280; 00110 break; 00111 case 6: 00112 _b1 = 1/(1+_n)*(nx*(nx*(nx+4)+64)+256)/256; 00113 _alp[1] = _n*(_n*(_n*(_n*(_n*(31564*_n-66675)+34440)+47250)-100800)+ 00114 75600)/151200; 00115 _bet[1] = _n*(_n*(_n*(_n*(_n*(384796*_n-382725)-6720)+932400)-1612800)+ 00116 1209600)/2419200; 00117 _alp[2] = nx*(_n*(_n*((863232-1983433*_n)*_n+748608)-1161216)+524160)/ 00118 1935360; 00119 _bet[2] = nx*(_n*(_n*((1695744-1118711*_n)*_n-1174656)+258048)+80640)/ 00120 3870720; 00121 nx *= _n; 00122 _alp[3] = nx*(_n*(_n*(670412*_n+406647)-533952)+184464)/725760; 00123 _bet[3] = nx*(_n*(_n*(22276*_n-16929)-15984)+12852)/362880; 00124 nx *= _n; 00125 _alp[4] = nx*(_n*(6601661*_n-7732800)+2230245)/7257600; 00126 _bet[4] = nx*((-830251*_n-158400)*_n+197865)/7257600; 00127 nx *= _n; 00128 _alp[5] = (3438171-13675556*_n)*nx/7983360; 00129 _bet[5] = (453717-435388*_n)*nx/15966720; 00130 nx *= _n; 00131 _alp[6] = 212378941*nx/319334400; 00132 _bet[6] = 20648693*nx/638668800; 00133 break; 00134 case 7: 00135 _b1 = 1/(1+_n)*(nx*(nx*(nx+4)+64)+256)/256; 00136 _alp[1] = _n*(_n*(_n*(_n*(_n*(_n*(1804025*_n+2020096)-4267200)+2204160)+ 00137 3024000)-6451200)+4838400)/9676800; 00138 _bet[1] = _n*(_n*(_n*(_n*(_n*((6156736-5406467*_n)*_n-6123600)-107520)+ 00139 14918400)-25804800)+19353600)/38707200; 00140 _alp[2] = nx*(_n*(_n*(_n*(_n*(4626384*_n-9917165)+4316160)+3743040)- 00141 5806080)+2620800)/9676800; 00142 _bet[2] = nx*(_n*(_n*(_n*(_n*(829456*_n-5593555)+8478720)-5873280)+ 00143 1290240)+403200)/19353600; 00144 nx *= _n; 00145 _alp[3] = nx*(_n*(_n*((26816480-67102379*_n)*_n+16265880)-21358080)+ 00146 7378560)/29030400; 00147 _bet[3] = nx*(_n*(_n*(_n*(9261899*_n+3564160)-2708640)-2557440)+ 00148 2056320)/58060800; 00149 nx *= _n; 00150 _alp[4] = nx*(_n*(_n*(155912000*_n+72618271)-85060800)+24532695)/ 00151 79833600; 00152 _bet[4] = nx*(_n*(_n*(14928352*_n-9132761)-1742400)+2176515)/79833600; 00153 nx *= _n; 00154 _alp[5] = nx*(_n*(102508609*_n-109404448)+27505368)/63866880; 00155 _bet[5] = nx*((-8005831*_n-1741552)*_n+1814868)/63866880; 00156 nx *= _n; 00157 _alp[6] = (2760926233.0-12282192400.0*_n)*nx/4151347200.0; 00158 _bet[6] = (268433009-261810608*_n)*nx/8302694400.0; 00159 nx *= _n; 00160 _alp[7] = 1522256789.0*nx/1383782400.0; 00161 _bet[7] = 219941297*nx/5535129600.0; 00162 break; 00163 case 8: 00164 _b1 = 1/(1+_n)*(nx*(nx*(nx*(25*nx+64)+256)+4096)+16384)/16384; 00165 _alp[1] = _n*(_n*(_n*(_n*(_n*(_n*((37884525-75900428*_n)*_n+42422016)- 00166 89611200)+46287360)+63504000)-135475200)+ 00167 101606400)/203212800; 00168 _bet[1] = _n*(_n*(_n*(_n*(_n*(_n*(_n*(31777436*_n-37845269)+43097152)- 00169 42865200)-752640)+104428800)-180633600)+ 00170 135475200)/270950400; 00171 _alp[2] = nx*(_n*(_n*(_n*(_n*(_n*(148003883*_n+83274912)-178508970)+ 00172 77690880)+67374720)-104509440)+47174400)/ 00173 174182400; 00174 _bet[2] = nx*(_n*(_n*(_n*(_n*(_n*(24749483*_n+14930208)-100683990)+ 00175 152616960)-105719040)+23224320)+7257600)/ 00176 348364800; 00177 nx *= _n; 00178 _alp[3] = nx*(_n*(_n*(_n*(_n*(318729724*_n-738126169)+294981280)+ 00179 178924680)-234938880)+81164160)/319334400; 00180 _bet[3] = nx*(_n*(_n*(_n*((101880889-232468668*_n)*_n+39205760)- 00181 29795040)-28131840)+22619520)/638668800; 00182 nx *= _n; 00183 _alp[4] = nx*(_n*(_n*((14967552000.0-40176129013.0*_n)*_n+6971354016.0)- 00184 8165836800.0)+2355138720.0)/7664025600.0; 00185 _bet[4] = nx*(_n*(_n*(_n*(324154477*_n+1433121792.0)-876745056)- 00186 167270400)+208945440)/7664025600.0; 00187 nx *= _n; 00188 _alp[5] = nx*(_n*(_n*(10421654396.0*_n+3997835751.0)-4266773472.0)+ 00189 1072709352.0)/2490808320.0; 00190 _bet[5] = nx*(_n*(_n*(457888660*_n-312227409)-67920528)+70779852)/ 00191 2490808320.0; 00192 nx *= _n; 00193 _alp[6] = nx*(_n*(175214326799.0*_n-171950693600.0)+38652967262.0)/ 00194 58118860800.0; 00195 _bet[6] = nx*((-19841813847.0*_n-3665348512.0)*_n+3758062126.0)/ 00196 116237721600.0; 00197 nx *= _n; 00198 _alp[7] = (13700311101.0-67039739596.0*_n)*nx/12454041600.0; 00199 _bet[7] = (1979471673.0-1989295244.0*_n)*nx/49816166400.0; 00200 nx *= _n; 00201 _alp[8] = 1424729850961.0*nx/743921418240.0; 00202 _bet[8] = 191773887257.0*nx/3719607091200.0; 00203 break; 00204 default: 00205 STATIC_ASSERT(maxpow_ >= 4 && maxpow_ <= 8, "Bad value of maxpow_"); 00206 } 00207 // _a1 is the equivalent radius for computing the circumference of 00208 // ellipse. 00209 _a1 = _b1 * _a; 00210 } 00211 00212 const TransverseMercator 00213 TransverseMercator::UTM(Constants::WGS84_a<real>(), 00214 Constants::WGS84_f<real>(), 00215 Constants::UTM_k0<real>()); 00216 00217 void TransverseMercator::Forward(real lon0, real lat, real lon, 00218 real& x, real& y, real& gamma, real& k) 00219 const throw() { 00220 // Avoid losing a bit of accuracy in lon (assuming lon0 is an integer) 00221 if (lon - lon0 > 180) 00222 lon -= lon0 + 360; 00223 else if (lon - lon0 <= -180) 00224 lon -= lon0 - 360; 00225 else 00226 lon -= lon0; 00227 // Now lon in (-180, 180] 00228 // Explicitly enforce the parity 00229 int 00230 latsign = lat < 0 ? -1 : 1, 00231 lonsign = lon < 0 ? -1 : 1; 00232 lon *= lonsign; 00233 lat *= latsign; 00234 bool backside = lon > 90; 00235 if (backside) { 00236 if (lat == 0) 00237 latsign = -1; 00238 lon = 180 - lon; 00239 } 00240 real 00241 phi = lat * Math::degree<real>(), 00242 lam = lon * Math::degree<real>(); 00243 // phi = latitude 00244 // phi' = conformal latitude 00245 // psi = isometric latitude 00246 // tau = tan(phi) 00247 // tau' = tan(phi') 00248 // [xi', eta'] = Gauss-Schreiber TM coordinates 00249 // [xi, eta] = Gauss-Krueger TM coordinates 00250 // 00251 // We use 00252 // tan(phi') = sinh(psi) 00253 // sin(phi') = tanh(psi) 00254 // cos(phi') = sech(psi) 00255 // denom^2 = 1-cos(phi')^2*sin(lam)^2 = 1-sech(psi)^2*sin(lam)^2 00256 // sin(xip) = sin(phi')/denom = tanh(psi)/denom 00257 // cos(xip) = cos(phi')*cos(lam)/denom = sech(psi)*cos(lam)/denom 00258 // cosh(etap) = 1/denom = 1/denom 00259 // sinh(etap) = cos(phi')*sin(lam)/denom = sech(psi)*sin(lam)/denom 00260 real etap, xip; 00261 if (lat != 90) { 00262 real 00263 c = max(real(0), cos(lam)), // cos(pi/2) might be negative 00264 tau = tan(phi), 00265 secphi = Math::hypot(real(1), tau), 00266 sig = sinh( eatanhe(tau / secphi) ), 00267 taup = Math::hypot(real(1), sig) * tau - sig * secphi; 00268 xip = atan2(taup, c); 00269 // Used to be 00270 // etap = Math::atanh(sin(lam) / cosh(psi)); 00271 etap = Math::asinh(sin(lam) / Math::hypot(taup, c)); 00272 // convergence and scale for Gauss-Schreiber TM (xip, etap) -- gamma0 = 00273 // atan(tan(xip) * tanh(etap)) = atan(tan(lam) * sin(phi')); 00274 // sin(phi') = tau'/sqrt(1 + tau'^2) 00275 gamma = atan(tanx(lam) * 00276 taup / Math::hypot(real(1), taup)); // Krueger p 22 (44) 00277 // k0 = sqrt(1 - _e2 * sin(phi)^2) * (cos(phi') / cos(phi)) * cosh(etap) 00278 // Note 1/cos(phi) = cosh(psip); 00279 // and cos(phi') * cosh(etap) = 1/hypot(sinh(psi), cos(lam)) 00280 // 00281 // This form has cancelling errors. This property is lost if cosh(psip) 00282 // is replaced by 1/cos(phi), even though it's using "primary" data (phi 00283 // instead of psip). 00284 k = sqrt(_e2m + _e2 * Math::sq(cos(phi))) * secphi / Math::hypot(taup, c); 00285 } else { 00286 xip = Math::pi<real>()/2; 00287 etap = 0; 00288 gamma = lam; 00289 k = _c; 00290 } 00291 // {xi',eta'} is {northing,easting} for Gauss-Schreiber transverse Mercator 00292 // (for eta' = 0, xi' = bet). {xi,eta} is {northing,easting} for transverse 00293 // Mercator with constant scale on the central meridian (for eta = 0, xip = 00294 // rectifying latitude). Define 00295 // 00296 // zeta = xi + i*eta 00297 // zeta' = xi' + i*eta' 00298 // 00299 // The conversion from conformal to rectifying latitude can be expresses as 00300 // a series in _n: 00301 // 00302 // zeta = zeta' + sum(h[j-1]' * sin(2 * j * zeta'), j = 1..maxpow_) 00303 // 00304 // where h[j]' = O(_n^j). The reversion of this series gives 00305 // 00306 // zeta' = zeta - sum(h[j-1] * sin(2 * j * zeta), j = 1..maxpow_) 00307 // 00308 // which is used in Reverse. 00309 // 00310 // Evaluate sums via Clenshaw method. See 00311 // http://mathworld.wolfram.com/ClenshawRecurrenceFormula.html 00312 // 00313 // Let 00314 // 00315 // S = sum(c[k] * F[k](x), k = 0..N) 00316 // F[n+1](x) = alpha(n,x) * F[n](x) + beta(n,x) * F[n-1](x) 00317 // 00318 // Evaluate S with 00319 // 00320 // y[N+2] = y[N+1] = 0 00321 // y[k] = alpha(k,x) * y[k+1] + beta(k+1,x) * y[k+2] + c[k] 00322 // S = c[0] * F[0](x) + y[1] * F[1](x) + beta(1,x) * F[0](x) * y[2] 00323 // 00324 // Here we have 00325 // 00326 // x = 2 * zeta' 00327 // F[n](x) = sin(n * x) 00328 // a(n, x) = 2 * cos(x) 00329 // b(n, x) = -1 00330 // [ sin(A+B) - 2*cos(B)*sin(A) + sin(A-B) = 0, A = n*x, B = x ] 00331 // N = maxpow_ 00332 // c[k] = _alp[k] 00333 // S = y[1] * sin(x) 00334 // 00335 // For the derivative we have 00336 // 00337 // x = 2 * zeta' 00338 // F[n](x) = cos(n * x) 00339 // a(n, x) = 2 * cos(x) 00340 // b(n, x) = -1 00341 // [ cos(A+B) - 2*cos(B)*cos(A) + cos(A-B) = 0, A = n*x, B = x ] 00342 // c[0] = 1; c[k] = 2*k*_alp[k] 00343 // S = (c[0] - y[2]) + y[1] * cos(x) 00344 real 00345 c0 = cos(2 * xip), ch0 = cosh(2 * etap), 00346 s0 = sin(2 * xip), sh0 = sinh(2 * etap), 00347 ar = 2 * c0 * ch0, ai = -2 * s0 * sh0; // 2 * cos(2*zeta') 00348 int n = maxpow_; 00349 real 00350 xi0 = (n & 1 ? _alp[n] : 0), eta0 = 0, 00351 xi1 = 0, eta1 = 0; 00352 real // Accumulators for dzeta/dzeta' 00353 yr0 = (n & 1 ? 2 * maxpow_ * _alp[n--] : 0), yi0 = 0, 00354 yr1 = 0, yi1 = 0; 00355 while (n) { 00356 xi1 = ar * xi0 - ai * eta0 - xi1 + _alp[n]; 00357 eta1 = ai * xi0 + ar * eta0 - eta1; 00358 yr1 = ar * yr0 - ai * yi0 - yr1 + 2 * n * _alp[n]; 00359 yi1 = ai * yr0 + ar * yi0 - yi1; 00360 --n; 00361 xi0 = ar * xi1 - ai * eta1 - xi0 + _alp[n]; 00362 eta0 = ai * xi1 + ar * eta1 - eta0; 00363 yr0 = ar * yr1 - ai * yi1 - yr0 + 2 * n * _alp[n]; 00364 yi0 = ai * yr1 + ar * yi1 - yi0; 00365 --n; 00366 } 00367 ar /= 2; ai /= 2; // cos(2*zeta') 00368 yr1 = 1 - yr1 + ar * yr0 - ai * yi0; 00369 yi1 = - yi1 + ai * yr0 + ar * yi0; 00370 ar = s0 * ch0; ai = c0 * sh0; // sin(2*zeta') 00371 real 00372 xi = xip + ar * xi0 - ai * eta0, 00373 eta = etap + ai * xi0 + ar * eta0; 00374 // Fold in change in convergence and scale for Gauss-Schreiber TM to 00375 // Gauss-Krueger TM. 00376 gamma -= atan2(yi1, yr1); 00377 k *= _b1 * Math::hypot(yr1, yi1); 00378 gamma /= Math::degree<real>(); 00379 y = _a1 * _k0 * (backside ? Math::pi<real>() - xi : xi) * latsign; 00380 x = _a1 * _k0 * eta * lonsign; 00381 if (backside) 00382 gamma = 180 - gamma; 00383 gamma *= latsign * lonsign; 00384 k *= _k0; 00385 } 00386 00387 void TransverseMercator::Reverse(real lon0, real x, real y, 00388 real& lat, real& lon, real& gamma, real& k) 00389 const throw() { 00390 // This undoes the steps in Forward. The wrinkles are: (1) Use of the 00391 // reverted series to express zeta' in terms of zeta. (2) Newton's method 00392 // to solve for phi in terms of tan(phi). 00393 real 00394 xi = y / (_a1 * _k0), 00395 eta = x / (_a1 * _k0); 00396 // Explicitly enforce the parity 00397 int 00398 xisign = xi < 0 ? -1 : 1, 00399 etasign = eta < 0 ? -1 : 1; 00400 xi *= xisign; 00401 eta *= etasign; 00402 bool backside = xi > Math::pi<real>()/2; 00403 if (backside) 00404 xi = Math::pi<real>() - xi; 00405 real 00406 c0 = cos(2 * xi), ch0 = cosh(2 * eta), 00407 s0 = sin(2 * xi), sh0 = sinh(2 * eta), 00408 ar = 2 * c0 * ch0, ai = -2 * s0 * sh0; // 2 * cos(2*zeta) 00409 int n = maxpow_; 00410 real // Accumulators for zeta' 00411 xip0 = (n & 1 ? -_bet[n] : 0), etap0 = 0, 00412 xip1 = 0, etap1 = 0; 00413 real // Accumulators for dzeta'/dzeta 00414 yr0 = (n & 1 ? - 2 * maxpow_ * _bet[n--] : 0), yi0 = 0, 00415 yr1 = 0, yi1 = 0; 00416 while (n) { 00417 xip1 = ar * xip0 - ai * etap0 - xip1 - _bet[n]; 00418 etap1 = ai * xip0 + ar * etap0 - etap1; 00419 yr1 = ar * yr0 - ai * yi0 - yr1 - 2 * n * _bet[n]; 00420 yi1 = ai * yr0 + ar * yi0 - yi1; 00421 --n; 00422 xip0 = ar * xip1 - ai * etap1 - xip0 - _bet[n]; 00423 etap0 = ai * xip1 + ar * etap1 - etap0; 00424 yr0 = ar * yr1 - ai * yi1 - yr0 - 2 * n * _bet[n]; 00425 yi0 = ai * yr1 + ar * yi1 - yi0; 00426 --n; 00427 } 00428 ar /= 2; ai /= 2; // cos(2*zeta') 00429 yr1 = 1 - yr1 + ar * yr0 - ai * yi0; 00430 yi1 = - yi1 + ai * yr0 + ar * yi0; 00431 ar = s0 * ch0; ai = c0 * sh0; // sin(2*zeta) 00432 real 00433 xip = xi + ar * xip0 - ai * etap0, 00434 etap = eta + ai * xip0 + ar * etap0; 00435 // Convergence and scale for Gauss-Schreiber TM to Gauss-Krueger TM. 00436 gamma = atan2(yi1, yr1); 00437 k = _b1 / Math::hypot(yr1, yi1); 00438 // JHS 154 has 00439 // 00440 // phi' = asin(sin(xi') / cosh(eta')) (Krueger p 17 (25)) 00441 // lam = asin(tanh(eta') / cos(phi') 00442 // psi = asinh(tan(phi')) 00443 real lam, phi; 00444 real 00445 s = sinh(etap), 00446 c = max(real(0), cos(xip)), // cos(pi/2) might be negative 00447 r = Math::hypot(s, c); 00448 if (r != 0) { 00449 lam = atan2(s, c); // Krueger p 17 (25) 00450 // Use Newton's method to solve for tau 00451 real 00452 taup = sin(xip)/r, 00453 // To lowest order in e^2, taup = (1 - e^2) * tau = _e2m * tau; so use 00454 // tau = taup/_e2m as a starting guess. Only 1 iteration is needed for 00455 // |lat| < 3.35 deg, otherwise 2 iterations are needed. If, instead, 00456 // tau = taup is used the mean number of iterations increases to 1.99 00457 // (2 iterations are needed except near tau = 0). 00458 tau = taup/_e2m, 00459 stol = tol_ * max(real(1), abs(taup)); 00460 // min iterations = 1, max iterations = 2; mean = 1.94 00461 for (int i = 0; i < numit_; ++i) { 00462 real 00463 tau1 = Math::hypot(real(1), tau), 00464 sig = sinh( eatanhe( tau / tau1 ) ), 00465 taupa = Math::hypot(real(1), sig) * tau - sig * tau1, 00466 dtau = (taup - taupa) * (1 + _e2m * Math::sq(tau)) / 00467 ( _e2m * tau1 * Math::hypot(real(1), taupa) ); 00468 tau += dtau; 00469 if (!(abs(dtau) >= stol)) 00470 break; 00471 } 00472 phi = atan(tau); 00473 gamma += atan(tanx(xip) * tanh(etap)); // Krueger p 19 (31) 00474 // Note cos(phi') * cosh(eta') = r 00475 k *= sqrt(_e2m + _e2 * Math::sq(cos(phi))) * 00476 Math::hypot(real(1), tau) * r; 00477 } else { 00478 phi = Math::pi<real>()/2; 00479 lam = 0; 00480 k *= _c; 00481 } 00482 lat = phi / Math::degree<real>() * xisign; 00483 lon = lam / Math::degree<real>(); 00484 if (backside) 00485 lon = 180 - lon; 00486 lon *= etasign; 00487 // Avoid losing a bit of accuracy in lon (assuming lon0 is an integer) 00488 if (lon + lon0 >= 180) 00489 lon += lon0 - 360; 00490 else if (lon + lon0 < -180) 00491 lon += lon0 + 360; 00492 else 00493 lon += lon0; 00494 gamma /= Math::degree<real>(); 00495 if (backside) 00496 gamma = 180 - gamma; 00497 gamma *= xisign * etasign; 00498 k *= _k0; 00499 } 00500 00501 } // namespace GeographicLib