GeographicLib
1.21
|
00001 /** 00002 * \file PolarStereographic.cpp 00003 * \brief Implementation for GeographicLib::PolarStereographic 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 00010 #include <GeographicLib/PolarStereographic.hpp> 00011 00012 #define GEOGRAPHICLIB_POLARSTEREOGRAPHIC_CPP \ 00013 "$Id: 3a2dee07e6ef1c55ddcdc2178d818c8edd4d1cd4 $" 00014 00015 RCSID_DECL(GEOGRAPHICLIB_POLARSTEREOGRAPHIC_CPP) 00016 RCSID_DECL(GEOGRAPHICLIB_POLARSTEREOGRAPHIC_HPP) 00017 00018 namespace GeographicLib { 00019 00020 using namespace std; 00021 00022 const Math::real PolarStereographic::tol_ = 00023 real(0.1)*sqrt(numeric_limits<real>::epsilon()); 00024 // Overflow value s.t. atan(overflow_) = pi/2 00025 const Math::real PolarStereographic::overflow_ = 00026 1 / Math::sq(numeric_limits<real>::epsilon()); 00027 00028 PolarStereographic::PolarStereographic(real a, real f, real k0) 00029 : _a(a) 00030 , _f(f <= 1 ? f : 1/f) 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 (!(Math::isfinite(_a) && _a > 0)) 00039 throw GeographicErr("Major radius is not positive"); 00040 if (!(Math::isfinite(_f) && _f < 1)) 00041 throw GeographicErr("Minor radius is not positive"); 00042 if (!(Math::isfinite(_k0) && _k0 > 0)) 00043 throw GeographicErr("Scale is not positive"); 00044 } 00045 00046 const PolarStereographic 00047 PolarStereographic::UPS(Constants::WGS84_a<real>(), 00048 Constants::WGS84_f<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 ? tanx(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 / Math::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 * Math::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 ? 00120 (rho / _a) * secphi * sqrt(_e2m + _e2 / Math::sq(secphi)) : _k0; 00121 lat = (northp ? 1 : -1) * (rho != 0 ? phi / Math::degree<real>() : 90); 00122 lon = -atan2( -x, northp ? -y : y ) / Math::degree<real>(); 00123 gamma = northp ? lon : -lon; 00124 } 00125 00126 void PolarStereographic::SetScale(real lat, real k) { 00127 if (!(Math::isfinite(k) && k > 0)) 00128 throw GeographicErr("Scale is not positive"); 00129 if (!(-90 < lat && lat <= 90)) 00130 throw GeographicErr("Latitude must be in (-90d, 90d]"); 00131 real x, y, gamma, kold; 00132 _k0 = 1; 00133 Forward(true, lat, 0, x, y, gamma, kold); 00134 _k0 *= k/kold; 00135 } 00136 00137 } // namespace GeographicLib