LambertConformalConic.cpp

Go to the documentation of this file.
00001 /**
00002  * \file LambertConformalConic.cpp
00003  * \brief Implementation for GeographicLib::LambertConformalConic class
00004  *
00005  * Copyright (c) Charles Karney (2009, 2010) <charles@karney.com>
00006  * and licensed under the LGPL.  For more information, see
00007  * http://geographiclib.sourceforge.net/
00008  **********************************************************************/
00009 
00010 #include "GeographicLib/LambertConformalConic.hpp"
00011 
00012 #define GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_CPP "$Id: LambertConformalConic.cpp 6816 2010-02-05 21:03:10Z karney $"
00013 
00014 RCSID_DECL(GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_CPP)
00015 RCSID_DECL(GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_HPP)
00016 
00017 namespace GeographicLib {
00018 
00019   using namespace std;
00020 
00021   const Math::real LambertConformalConic::eps = numeric_limits<real>::epsilon();
00022   const Math::real LambertConformalConic::eps2 = sqrt(eps);
00023   const Math::real LambertConformalConic::epsx = sq(eps);
00024   const Math::real LambertConformalConic::tol = real(0.1) * eps2;
00025   // Large enough value so that atan(sinh(ahypover)) = pi/2
00026   const Math::real LambertConformalConic::ahypover =
00027     real(numeric_limits<real>::digits) * log(real(numeric_limits<real>::radix))
00028     + 2;
00029 
00030   LambertConformalConic::LambertConformalConic(real a, real r,
00031                                                real stdlat, real k0)
00032     : _a(a)
00033     , _r(r)
00034     , _f(_r != 0 ? 1 / _r : 0)
00035     , _e2(_f * (2 - _f))
00036     , _e(sqrt(abs(_e2)))
00037     , _e2m(1 - _e2)
00038   {
00039     if (!(_a > 0))
00040       throw GeographicErr("Major radius is not positive");
00041     if (!(k0 > 0))
00042       throw GeographicErr("Scale is not positive");
00043     if (!(abs(stdlat) <= 90))
00044       throw GeographicErr("Standard latitude not in [-90, 90]");
00045     real
00046       phi = stdlat * Constants::degree(),
00047       sphi = sin(phi),
00048       cphi = abs(stdlat) != 90 ? cos(phi) : 0;
00049     Init(sphi, cphi, sphi, cphi, k0);
00050   }
00051 
00052   LambertConformalConic::LambertConformalConic(real a, real r,
00053                                                real stdlat1, real stdlat2,
00054                                                real k1)
00055     : _a(a)
00056     , _r(r)
00057     , _f(_r != 0 ? 1 / _r : 0)
00058     , _e2(_f * (2 - _f))
00059     , _e(sqrt(abs(_e2)))
00060     , _e2m(1 - _e2)
00061   {
00062     if (!(_a > 0))
00063       throw GeographicErr("Major radius is not positive");
00064     if (!(k1 > 0))
00065       throw GeographicErr("Scale is not positive");
00066     if (!(abs(stdlat1) <= 90))
00067       throw GeographicErr("Standard latitude 1 not in [-90, 90]");
00068     if (!(abs(stdlat2) <= 90))
00069       throw GeographicErr("Standard latitude 2 not in [-90, 90]");
00070     if (abs(stdlat1) == 90 || abs(stdlat2) == 90)
00071       if (!(stdlat1 == stdlat2))
00072         throw GeographicErr
00073           ("Standard latitudes must be equal is either is a pole");
00074     real
00075       phi1 = stdlat1 * Constants::degree(),
00076       phi2 = stdlat2 * Constants::degree();
00077     Init(sin(phi1), abs(stdlat1) != 90 ? cos(phi1) : 0,
00078          sin(phi2), abs(stdlat2) != 90 ? cos(phi2) : 0, k1);
00079   }
00080 
00081   LambertConformalConic::LambertConformalConic(real a, real r,
00082                                                real sinlat1, real coslat1,
00083                                                real sinlat2, real coslat2,
00084                                                real k1)
00085     : _a(a)
00086     , _r(r)
00087     , _f(_r != 0 ? 1 / _r : 0)
00088     , _e2(_f * (2 - _f))
00089     , _e(sqrt(abs(_e2)))
00090     , _e2m(1 - _e2)
00091   {
00092     if (!(_a > 0))
00093       throw GeographicErr("Major radius is not positive");
00094     if (!(k1 > 0))
00095       throw GeographicErr("Scale is not positive");
00096     if (coslat1 == 0 || coslat2 == 0)
00097       if (!(coslat1 == coslat2 && sinlat1 == sinlat2))
00098         throw GeographicErr
00099           ("Standard latitudes must be equal is either is a pole");
00100     Init(sinlat1, coslat1, sinlat2, coslat2, k1);
00101   }
00102 
00103   void LambertConformalConic::Init(real sphi1, real cphi1,
00104                                    real sphi2, real cphi2, real k1) throw() {
00105     // Determine hemisphere of tangent latitude
00106     _sign = sphi1 + sphi2 > 0 ? 1 : sphi1 + sphi2 < 0  ? -1 :
00107       atan2(sphi1, cphi1) + atan2(sphi2, cphi2) >= 0 ? 1 : -1;
00108     // Internally work with tangent latitude positive
00109     sphi1 *= _sign;
00110     sphi2 *= _sign;
00111     real
00112       m1 = mf(sphi1, cphi1), lt1 = logtf(sphi1, cphi1),
00113       m2 = mf(sphi2, cphi2), lt2 = logtf(sphi2, cphi2),
00114       sindiff = abs(sphi1 * cphi2 - cphi1 * sphi2);
00115     real sphi0, cphi0;          // Use phi0 = tangent latitude
00116     if (sindiff > eps2) {
00117       // sphi0 = Snyder's n, p 108, eq 15-8
00118       sphi0 = (log(m1) - log (m2)) / (lt1 - lt2);
00119       cphi0 = sindiff < real(0.7) ? sqrt(1 - sq(sphi0)) :
00120         // cos(phi0) = sqrt( - (sin(phi0) - 1) * (sin(phi0) + 1) )
00121         sqrt( -(logmtf(sphi1) - logmtf(sphi2)) / (lt1 - lt2)
00122               * (sphi0 + 1));
00123     } else {
00124       // Set phi0 = (phi1 + phi2)/2
00125       sphi0 = (sphi1 * sqrt((1 + cphi2)/(1 + cphi1)) +
00126                sphi2 * sqrt((1 + cphi1)/(1 + cphi2)))/2;
00127       cphi0 = (cphi1 * sqrt((1 + sphi2)/(1 + sphi1)) +
00128                cphi2 * sqrt((1 + sphi1)/(1 + sphi2)))/2;
00129     }
00130     _n = sphi0;                 // Snyder's n
00131     _nc = cphi0;                // sqrt(1 - sq(n))
00132     _lat0 = atan2(_sign * sphi0, cphi0) / Constants::degree();
00133     _lt0 = logtf(sphi0, cphi0); // Snyder's log(t0)
00134     _t0n = exp(_n * _lt0);      // Snyder's t0^n
00135     _t0nm1 = Math::expm1(_n * _lt0);      // Snyder's t0^n - 1
00136     // k1 * m1/t1^n = k1 * m2/t2^n = k1 * n * (Snyder's F)
00137     _scale = k1 *
00138       (_nc == 0 ?               // if phi0 = pi/2, n = t = m = 0, so take limit
00139        2/( sqrt(_e2m) * exp(eatanhe(real(1))) ) :
00140        (lt1 > lt2 ? m1 / exp(_n * lt1) : m2 / exp(_n * lt2)));
00141     // Scale at phi0
00142     _k0 = _scale * _t0n / mf(sphi0, cphi0);
00143   }
00144 
00145   const LambertConformalConic
00146   LambertConformalConic::Mercator(Constants::WGS84_a(), Constants::WGS84_r(),
00147                                   real(0), real(1));
00148 
00149   void LambertConformalConic::Forward(real lon0, real lat, real lon,
00150                                       real& x, real& y, real& gamma, real& k)
00151     const throw() {
00152     if (lon - lon0 > 180)
00153       lon -= lon0 - 360;
00154     else if (lon - lon0 <= -180)
00155       lon -= lon0 + 360;
00156     else
00157       lon -= lon0;
00158     lat *= _sign;
00159     real
00160       phi = lat * Constants::degree(),
00161       lam = lon * Constants::degree(),
00162       sphi = sin(phi), cphi = abs(lat) != 90 ? cos(phi) : 0,
00163       m = mf(sphi, cphi),
00164       lt = logtf(sphi, cphi),
00165       tn = exp(_n * lt),
00166       theta = _n * lam, stheta = sin(theta), ctheta = cos(theta);
00167     x = _a * _scale * tn * (_n > eps2 ? stheta/_n : lam);
00168     if (_n > real(0.5))
00169       y = (_t0n - tn * ctheta)/_n;
00170     else {
00171       // write as
00172       // _t0n * ( (1 - ctheta)/_n - (tn/_t0n - 1)/_n * ctheta )
00173       // _t0n * ( stheta^2/(1 + ctheta) / _n
00174       //          - ((t/_t0)^n - 1)/_n * ctheta )
00175       // _t0n * ( stheta^2/(1 + ctheta) / _n
00176       //          - (exp(_n * log(t/_t0)) - 1)/_n * ctheta )
00177       // _t0n * (ax - bx)
00178       real
00179         ax = (_n > eps2 ? sq(stheta) / _n : lam * theta) / (1 + ctheta),
00180         fx = lt - _lt0,
00181         bx = ctheta * (abs(_n * fx) > eps ? Math::expm1(_n * fx) / _n : fx);
00182       y = _t0n * (ax - bx);
00183     }
00184     y *= _a * _scale * _sign;
00185     k = _scale * (m == 0 && _nc == 0 && sphi > 0 ?
00186                   sqrt(_e2m) * exp(eatanhe(real(1))) / 2 :
00187                   tn/m);        // infinite if pole and _n < 1
00188     gamma = theta * _sign;
00189   }
00190 
00191   void LambertConformalConic::Reverse(real lon0, real x, real y,
00192                                       real& lat, real& lon,
00193                                       real& gamma, real& k)
00194     const throw() {
00195     y *= _sign;
00196     x /= _a * _scale;
00197     y /= _a * _scale;
00198     real
00199       /*
00200         rho = hypot(x, rho0 - y)
00201         rho0 = _a * _scale * _t0n / _n
00202         _n * rho = hypot(x * _n, _a * _scale * _t0n - y * _n)
00203         tn    = _n * rho / (_a * _scale)
00204               = hypot(_n * x/(_a * _scale), _t0n - _n * y/(_a * _scale))
00205         theta = atan2(_n * x/(_a * _scale), _t0n - _n * y/(_a * _scale))
00206         lam = theta/_n
00207       */
00208       x1 = _n * x, y1 = _n * y,
00209       tn = Math::hypot(x1, _t0n - y1), theta = atan2(x1, _t0n - y1),
00210       lam = _n != 0 ? theta/_n : x/_t0n;
00211     real q;                     // -log(t)
00212     if (_n != 0 && tn < real(0.5))
00213       q = -log(tn) / _n;
00214     else {
00215       real
00216         b1 = _t0nm1 - _n * y,
00217         tnm1 = (sq(x1) +  2 * b1 + sq(b1)) / (tn + 1); // tn - 1
00218       q = _n != 0 ? - Math::log1p(tnm1)/_n :
00219         -2 * (_lt0  - y) / (tn + 1); // _t0nm1 -> _n * _lt0
00220     }
00221     // Clip to [-ahypover, ahypover] to avoid overflow later
00222     q = q < ahypover ? (q > -ahypover ? q : -ahypover) : ahypover;
00223     // Recast Snyder's 15-9 as
00224     // q = q' - eatanh(tanh(q'))
00225     // where q = -log(t) and q' = asinh(tan(phi))
00226     // q is known and we solve for q' by Newton's method.
00227     // Write f(q') = q' - eatanh(tanh(q')) - q
00228     // f'(q') = (1 - e^2)/(1 - e^2 * tanh(q')^2)
00229     // Starting guess is q' = q.
00230     real qp = q;
00231     for (int i = 0; i < numit; ++i) {
00232       real
00233         t = tanh(qp),
00234         dqp = -(qp - eatanhe(t) - q) * (1 - _e2 * sq(t)) / _e2m;
00235       qp += dqp;
00236       if (abs(dqp) < tol)
00237         break;
00238     }
00239     double
00240       phi = _sign * atan(sinh(qp)),
00241       m = mf(tanh(qp), 1/cosh(qp));
00242     lat = phi / Constants::degree();
00243     lon = lam / Constants::degree();
00244     // Avoid losing a bit of accuracy in lon (assuming lon0 is an integer)
00245     if (lon + lon0 >= 180)
00246       lon += lon0 - 360;
00247     else if (lon + lon0 < -180)
00248       lon += lon0 + 360;
00249     else
00250       lon += lon0;
00251     gamma = _sign * theta / Constants::degree();
00252     k = _scale * (m == 0 && _nc == 0 && qp > 0 ?
00253                   sqrt(_e2m) * exp(eatanhe(real(1))) / 2 :
00254                   tn/m);        // infinite if pole and _n < 1
00255   }
00256 
00257   void LambertConformalConic::SetScale(real lat, real k) {
00258     if (!(k > 0))
00259       throw GeographicErr("Scale is not positive");
00260     real x, y, gamma, kold;
00261     Forward(0, lat, 0, x, y, gamma, kold);
00262     k /= kold;
00263     _scale *= k;
00264     _k0 *= k;
00265   }
00266 
00267 } // namespace GeographicLib

Generated on 21 May 2010 for GeographicLib by  doxygen 1.6.1