00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include <GeographicLib/PolarStereographic.hpp>
00011
00012 #if defined(_MSC_VER)
00013
00014 # pragma warning (disable: 4127)
00015 #endif
00016
00017 namespace GeographicLib {
00018
00019 using namespace std;
00020
00021 PolarStereographic::PolarStereographic(real a, real f, real k0)
00022 : tol_(real(0.1)*sqrt(numeric_limits<real>::epsilon()))
00023 , _a(a)
00024 , _f(f <= 1 ? f : 1/f)
00025 , _e2(_f * (2 - _f))
00026 , _e(sqrt(abs(_e2)))
00027 , _e2m(1 - _e2)
00028 , _Cx(exp(eatanhe(real(1))))
00029 , _c( (1 - _f) * _Cx )
00030 , _k0(k0)
00031 {
00032 if (!(Math::isfinite(_a) && _a > 0))
00033 throw GeographicErr("Major radius is not positive");
00034 if (!(Math::isfinite(_f) && _f < 1))
00035 throw GeographicErr("Minor radius is not positive");
00036 if (!(Math::isfinite(_k0) && _k0 > 0))
00037 throw GeographicErr("Scale is not positive");
00038 }
00039
00040 const PolarStereographic& PolarStereographic::UPS() {
00041 static const PolarStereographic ups(Constants::WGS84_a(),
00042 Constants::WGS84_f(),
00043 Constants::UPS_k0());
00044 return ups;
00045 }
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068 void PolarStereographic::Forward(bool northp, real lat, real lon,
00069 real& x, real& y, real& gamma, real& k)
00070 const {
00071 lat *= northp ? 1 : -1;
00072 real
00073 phi = lat * Math::degree(),
00074 tau = lat != -90 ? tanx(phi) : -overflow(),
00075 secphi = Math::hypot(real(1), tau),
00076 sig = sinh( eatanhe(tau / secphi) ),
00077 taup = Math::hypot(real(1), sig) * tau - sig * secphi,
00078 rho = Math::hypot(real(1), taup) + abs(taup);
00079 rho = taup >= 0 ? (lat != 90 ? 1/rho : 0) : rho;
00080 rho *= 2 * _k0 * _a / _c;
00081 k = lat != 90 ? (rho / _a) * secphi * sqrt(_e2m + _e2 / Math::sq(secphi)) :
00082 _k0;
00083 lon = Math::AngNormalize(lon);
00084 real
00085 lam = lon * Math::degree();
00086 x = rho * (lon == -180 ? 0 : sin(lam));
00087 y = (northp ? -rho : rho) * (abs(lon) == 90 ? 0 : cos(lam));
00088 gamma = northp ? lon : -lon;
00089 }
00090
00091 void PolarStereographic::Reverse(bool northp, real x, real y,
00092 real& lat, real& lon, real& gamma, real& k)
00093 const {
00094 real
00095 rho = Math::hypot(x, y),
00096 t = rho / (2 * _k0 * _a / _c),
00097 taup = (1 / t - t) / 2,
00098 tau = taup * _Cx,
00099 stol = tol_ * max(real(1), abs(taup));
00100 if (abs(tau) < overflow()) {
00101
00102 for (int i = 0; i < numit_ || GEOGRAPHICLIB_PANIC; ++i) {
00103 real
00104 tau1 = Math::hypot(real(1), tau),
00105 sig = sinh( eatanhe( tau / tau1 ) ),
00106 taupa = Math::hypot(real(1), sig) * tau - sig * tau1,
00107 dtau = (taup - taupa) * (1 + _e2m * Math::sq(tau)) /
00108 ( _e2m * tau1 * Math::hypot(real(1), taupa) );
00109 tau += dtau;
00110 if (!(abs(dtau) >= stol))
00111 break;
00112 }
00113 }
00114 real
00115 phi = atan(tau),
00116 secphi = Math::hypot(real(1), tau);
00117 k = rho ? (rho / _a) * secphi * sqrt(_e2m + _e2 / Math::sq(secphi)) : _k0;
00118 lat = (northp ? 1 : -1) * (rho ? phi / Math::degree() : 90);
00119 lon = -atan2( -x, northp ? -y : y ) / Math::degree();
00120 gamma = northp ? lon : -lon;
00121 }
00122
00123 void PolarStereographic::SetScale(real lat, real k) {
00124 if (!(Math::isfinite(k) && k > 0))
00125 throw GeographicErr("Scale is not positive");
00126 if (!(-90 < lat && lat <= 90))
00127 throw GeographicErr("Latitude must be in (-90d, 90d]");
00128 real x, y, gamma, kold;
00129 _k0 = 1;
00130 Forward(true, lat, 0, x, y, gamma, kold);
00131 _k0 *= k/kold;
00132 }
00133
00134 }