PolarStereographic.cpp

Go to the documentation of this file.
00001 /**
00002  * \file PolarStereographic.cpp
00003  * \brief Implementation for GeographicLib::PolarStereographic class
00004  *
00005  * Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.com>
00006  * and licensed under the LGPL.  For more information, see
00007  * http://geographiclib.sourceforge.net/
00008  **********************************************************************/
00009 
00010 #include "GeographicLib/PolarStereographic.hpp"
00011 
00012 #define GEOGRAPHICLIB_POLARSTEREOGRAPHIC_CPP "$Id: PolarStereographic.cpp 6937 2011-02-01 20:17:13Z karney $"
00013 
00014 RCSID_DECL(GEOGRAPHICLIB_POLARSTEREOGRAPHIC_CPP)
00015 RCSID_DECL(GEOGRAPHICLIB_POLARSTEREOGRAPHIC_HPP)
00016 
00017 namespace GeographicLib {
00018 
00019   using namespace std;
00020 
00021   const Math::real PolarStereographic::tol =
00022     real(0.1)*sqrt(numeric_limits<real>::epsilon());
00023   // Overflow value s.t. atan(overflow) = pi/2
00024   const Math::real PolarStereographic::overflow =
00025     1 / sq(numeric_limits<real>::epsilon());
00026 
00027   PolarStereographic::PolarStereographic(real a, real r, real k0)
00028     : _a(a)
00029     , _r(r)
00030     , _f(_r != 0 ? 1 / _r : 0)
00031     , _e2(_f * (2 - _f))
00032     , _e(sqrt(abs(_e2)))
00033     , _e2m(1 - _e2)
00034     , _Cx(exp(eatanhe(real(1))))
00035     , _c( (1 - _f) * _Cx )
00036     , _k0(k0)
00037   {
00038     if (!(_a > 0))
00039       throw GeographicErr("Major radius is not positive");
00040     if (!(_f < 1))
00041       throw GeographicErr("Minor radius is not positive");
00042     if (!(_k0 > 0))
00043       throw GeographicErr("Scale is not positive");
00044   }
00045 
00046   const PolarStereographic
00047   PolarStereographic::UPS(Constants::WGS84_a<real>(),
00048                           Constants::WGS84_r<real>(),
00049                           Constants::UPS_k0<real>());
00050 
00051   // This formulation converts to conformal coordinates by tau = tan(phi) and
00052   // tau' = tan(phi') where phi' is the conformal latitude.  The formulas are:
00053   //    tau = tan(phi)
00054   //    secphi = hypot(1, tau)
00055   //    sig = sinh(e * atanh(e * tau / secphi))
00056   //    taup = tan(phip) = tau * hypot(1, sig) - sig * hypot(1, tau)
00057   //    c = (1 - f) * exp(e * atanh(e))
00058   //
00059   // Forward:
00060   //   rho = (2*k0*a/c) / (hypot(1, taup) + taup)  (taup >= 0)
00061   //       = (2*k0*a/c) * (hypot(1, taup) - taup)  (taup <  0)
00062   //
00063   // Reverse:
00064   //   taup = ((2*k0*a/c) / rho - rho / (2*k0*a/c))/2
00065   //
00066   // Scale:
00067   //   k = (rho/a) * secphi * sqrt((1-e2) + e2 / secphi^2)
00068   //
00069   // In limit rho -> 0, tau -> inf, taup -> inf, secphi -> inf, secphip -> inf
00070   //   secphip = taup = exp(-e * atanh(e)) * tau = exp(-e * atanh(e)) * secphi
00071 
00072   void PolarStereographic::Forward(bool northp, real lat, real lon,
00073                                    real& x, real& y, real& gamma, real& k)
00074     const throw() {
00075     lat *= northp ? 1 : -1;
00076     real
00077       phi = lat * Math::degree<real>(),
00078       tau = lat != -90 ? tan(phi) : -overflow,
00079       secphi = Math::hypot(real(1), tau),
00080       sig = sinh( eatanhe(tau / secphi) ),
00081       taup = Math::hypot(real(1), sig) * tau - sig * secphi,
00082       rho =  Math::hypot(real(1), taup) + abs(taup);
00083     rho = taup >= 0 ? (lat != 90 ? 1/rho : 0) : rho;
00084     rho *= 2 * _k0 * _a / _c;
00085     k = lat != 90 ? (rho / _a) * secphi * sqrt(_e2m + _e2 / sq(secphi)) :
00086       _k0;
00087     lon = lon >= 180 ? lon - 360 : lon < -180 ? lon + 360 : lon;
00088     real
00089       lam = lon * Math::degree<real>();
00090     x = rho * (lon == -180 ? 0 : sin(lam));
00091     y = (northp ? -rho : rho) * (abs(lon) == 90 ? 0 : cos(lam));
00092     gamma = northp ? lon : -lon;
00093   }
00094 
00095   void PolarStereographic::Reverse(bool northp, real x, real y,
00096                                    real& lat, real& lon, real& gamma, real& k)
00097     const throw() {
00098     real
00099       rho = Math::hypot(x, y),
00100       t = rho / (2 * _k0 * _a / _c),
00101       taup = (1 / t - t) / 2,
00102       tau = taup * _Cx,
00103       stol = tol * max(real(1), abs(taup));
00104     // min iterations = 1, max iterations = 2; mean = 1.99
00105     for (int i = 0; i < numit; ++i) {
00106       real
00107         tau1 = Math::hypot(real(1), tau),
00108         sig = sinh( eatanhe( tau / tau1 ) ),
00109         taupa = Math::hypot(real(1), sig) * tau - sig * tau1,
00110         dtau = (taup - taupa) * (1 + _e2m * sq(tau)) /
00111         ( _e2m * tau1 * Math::hypot(real(1), taupa) );
00112       tau += dtau;
00113       if (!(abs(dtau) >= stol))
00114         break;
00115     }
00116     real
00117       phi = atan(tau),
00118       secphi = Math::hypot(real(1), tau);
00119     k = rho != 0 ? (rho / _a) * secphi * sqrt(_e2m + _e2 / sq(secphi)) : _k0;
00120     lat = (northp ? 1 : -1) * (rho != 0 ? phi / Math::degree<real>() : 90);
00121     lon = -atan2( -x, northp ? -y : y ) / Math::degree<real>();
00122     gamma = northp ? lon : -lon;
00123   }
00124 
00125   void PolarStereographic::SetScale(real lat, real k) {
00126     if (!(k > 0))
00127       throw GeographicErr("Scale is not positive");
00128     if (!(-90 < lat && lat <= 90))
00129       throw GeographicErr("Latitude must be in (-90d, 90d]");
00130     real x, y, gamma, kold;
00131     _k0 = 1;
00132     Forward(true, lat, 0, x, y, gamma, kold);
00133     _k0 *= k/kold;
00134   }
00135 
00136 } // namespace GeographicLib