TransverseMercator.cpp

Go to the documentation of this file.
00001 /**
00002  * \file TransverseMercator.cpp
00003  * \brief Implementation for GeographicLib::TransverseMercator class
00004  *
00005  * Copyright (c) Charles Karney (2008, 2009, 2010) <charles@karney.com>
00006  * and licensed under the LGPL.  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&auml;rjestelm&auml;&auml;n liittyv&auml;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 "$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   // Overflow value s.t. atan(overflow) = pi/2
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       // _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 (!(_a > 0))
00073       throw GeographicErr("Major radius is not positive");
00074     if (!(_k0 > 0))
00075       throw GeographicErr("Scale is not positive");
00076     // If coefficents might overflow an int, convert them to double (and they
00077     // are all exactly representable as doubles).
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     // _a1 is the equivalent radius for computing the circumference of
00206     // ellipse.
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     // Avoid losing a bit of accuracy in lon (assuming lon0 is an integer)
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     // Now lon in (-180, 180]
00225     // Explicitly enforce the parity
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     // phi = latitude
00241     // phi' = conformal latitude
00242     // psi = isometric latitude
00243     // tau = tan(phi)
00244     // tau' = tan(phi')
00245     // [xi', eta'] = Gauss-Schreiber TM coordinates
00246     // [xi, eta] = Gauss-Krueger TM coordinates
00247     //
00248     // We use
00249     //   tan(phi') = sinh(psi)
00250     //   sin(phi') = tanh(psi)
00251     //   cos(phi') = sech(psi)
00252     //   denom^2    = 1-cos(phi')^2*sin(lam)^2 = 1-sech(psi)^2*sin(lam)^2
00253     //   sin(xip)   = sin(phi')/denom          = tanh(psi)/denom
00254     //   cos(xip)   = cos(phi')*cos(lam)/denom = sech(psi)*cos(lam)/denom
00255     //   cosh(etap) = 1/denom                  = 1/denom
00256     //   sinh(etap) = cos(phi')*sin(lam)/denom = sech(psi)*sin(lam)/denom
00257     real etap, xip;
00258     if (lat < 90) {
00259       real
00260         c = max(real(0), cos(lam)), // cos(pi/2) might be negative
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       // Used to be
00267       //   etap = Math::atanh(sin(lam) / cosh(psi));
00268       etap = Math::asinh(sin(lam) / Math::hypot(taup, c));
00269       // convergence and scale for Gauss-Schreiber TM (xip, etap) -- gamma0 =
00270       // atan(tan(xip) * tanh(etap)) = atan(tan(lam) * sin(phi'));
00271       // sin(phi') = tau'/sqrt(1 + tau'^2)
00272       gamma = atan(tanx(lam) *
00273                    taup / Math::hypot(real(1), taup)); // Krueger p 22 (44)
00274       // k0 = sqrt(1 - _e2 * sin(phi)^2) * (cos(phi') / cos(phi)) * cosh(etap)
00275       // Note 1/cos(phi) = cosh(psip);
00276       // and cos(phi') * cosh(etap) = 1/hypot(sinh(psi), cos(lam))
00277       //
00278       // This form has cancelling errors.  This property is lost if cosh(psip)
00279       // is replaced by 1/cos(phi), even though it's using "primary" data (phi
00280       // instead of psip).
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     // {xi',eta'} is {northing,easting} for Gauss-Schreiber transverse Mercator
00289     // (for eta' = 0, xi' = bet). {xi,eta} is {northing,easting} for transverse
00290     // Mercator with constant scale on the central meridian (for eta = 0, xip =
00291     // rectifying latitude).  Define
00292     //
00293     //   zeta = xi + i*eta
00294     //   zeta' = xi' + i*eta'
00295     //
00296     // The conversion from conformal to rectifying latitude can be expresses as
00297     // a series in _n:
00298     //
00299     //   zeta = zeta' + sum(h[j-1]' * sin(2 * j * zeta'), j = 1..maxpow)
00300     //
00301     // where h[j]' = O(_n^j).  The reversion of this series gives
00302     //
00303     //   zeta' = zeta - sum(h[j-1] * sin(2 * j * zeta), j = 1..maxpow)
00304     //
00305     // which is used in Reverse.
00306     //
00307     // Evaluate sums via Clenshaw method.  See
00308     //    http://mathworld.wolfram.com/ClenshawRecurrenceFormula.html
00309     //
00310     // Let
00311     //
00312     //    S = sum(c[k] * F[k](x), k = 0..N)
00313     //    F[n+1](x) = alpha(n,x) * F[n](x) + beta(n,x) * F[n-1](x)
00314     //
00315     // Evaluate S with
00316     //
00317     //    y[N+2] = y[N+1] = 0
00318     //    y[k] = alpha(k,x) * y[k+1] + beta(k+1,x) * y[k+2] + c[k]
00319     //    S = c[0] * F[0](x) + y[1] * F[1](x) + beta(1,x) * F[0](x) * y[2]
00320     //
00321     // Here we have
00322     //
00323     //    x = 2 * zeta'
00324     //    F[n](x) = sin(n * x)
00325     //    a(n, x) = 2 * cos(x)
00326     //    b(n, x) = -1
00327     //    [ sin(A+B) - 2*cos(B)*sin(A) + sin(A-B) = 0, A = n*x, B = x ]
00328     //    N = maxpow
00329     //    c[k] = _alp[k]
00330     //    S = y[1] * sin(x)
00331     //
00332     // For the derivative we have
00333     //
00334     //    x = 2 * zeta'
00335     //    F[n](x) = cos(n * x)
00336     //    a(n, x) = 2 * cos(x)
00337     //    b(n, x) = -1
00338     //    [ cos(A+B) - 2*cos(B)*cos(A) + cos(A-B) = 0, A = n*x, B = x ]
00339     //    c[0] = 1; c[k] = 2*k*_alp[k]
00340     //    S = (c[0] - y[2]) + y[1] * cos(x)
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; // 2 * cos(2*zeta')
00345     int n = maxpow;
00346     real
00347       xi0 = (n & 1 ? _alp[n] : 0), eta0 = 0,
00348       xi1 = 0, eta1 = 0;
00349     real                        // Accumulators for dzeta/dzeta'
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;           // cos(2*zeta')
00365     yr1 = 1 - yr1 + ar * yr0 - ai * yi0;
00366     yi1 =   - yi1 + ai * yr0 + ar * yi0;
00367     ar = s0 * ch0; ai = c0 * sh0; // sin(2*zeta')
00368     real
00369       xi  = xip  + ar * xi0 - ai * eta0,
00370       eta = etap + ai * xi0 + ar * eta0;
00371     // Fold in change in convergence and scale for Gauss-Schreiber TM to
00372     // Gauss-Krueger TM.
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     // This undoes the steps in Forward.  The wrinkles are: (1) Use of the
00388     // reverted series to express zeta' in terms of zeta. (2) Newton's method
00389     // to solve for phi in terms of psi.
00390     real
00391       xi = y / (_a1 * _k0),
00392       eta = x / (_a1 * _k0);
00393     // Explicitly enforce the parity
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; // 2 * cos(2*zeta)
00406     int n = maxpow;
00407     real                        // Accumulators for zeta'
00408       xip0 = (n & 1 ? -_bet[n] : 0), etap0 = 0,
00409       xip1 = 0, etap1 = 0;
00410     real                        // Accumulators for dzeta'/dzeta
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;           // cos(2*zeta')
00426     yr1 = 1 - yr1 + ar * yr0 - ai * yi0;
00427     yi1 =   - yi1 + ai * yr0 + ar * yi0;
00428     ar = s0 * ch0; ai = c0 * sh0; // sin(2*zeta)
00429     real
00430       xip  = xi  + ar * xip0 - ai * etap0,
00431       etap = eta + ai * xip0 + ar * etap0;
00432     // Convergence and scale for Gauss-Schreiber TM to Gauss-Krueger TM.
00433     gamma = atan2(yi1, yr1);
00434     k = _b1 / Math::hypot(yr1, yi1);
00435     // JHS 154 has
00436     //
00437     //   phi' = asin(sin(xi') / cosh(eta')) (Krueger p 17 (25))
00438     //   lam = asin(tanh(eta') / cos(phi')
00439     //   psi = asinh(tan(phi'))
00440     real lam, phi;
00441     real
00442       s = sinh(etap),
00443       c = max(real(0), cos(xip)), // cos(pi/2) might be negative
00444       r = Math::hypot(s, c);
00445     if (r > 0) {
00446       lam = atan2(s, c);        // Krueger p 17 (25)
00447       // Use Newton's< method to solve for tau
00448       real
00449         taup = sin(xip)/r,
00450         tau = taup;
00451       // min iterations = 1, max iterations = 2; mean = 1.99
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)); // Krueger p 19 (31)
00465       // Note cos(phi') * cosh(eta') = r
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     // Avoid losing a bit of accuracy in lon (assuming lon0 is an integer)
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 } // namespace GeographicLib

Generated on 21 May 2010 for GeographicLib by  doxygen 1.6.1