00001
00002
00003
00004
00005
00006
00007
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
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
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
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
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 }