00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "GeographicLib/PolarStereographic.hpp"
00011
00012 #define GEOGRAPHICLIB_POLARSTEREOGRAPHIC_CPP "$Id: PolarStereographic.cpp 6785 2010-01-05 22:15:42Z 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 PolarStereographic::PolarStereographic(real a, real r, real k0)
00025 : _a(a)
00026 , _r(r)
00027 , _f(_r != 0 ? 1 / _r : 0)
00028 , _e2(_f * (2 - _f))
00029 , _e(sqrt(abs(_e2)))
00030 , _e2m(1 - _e2)
00031
00032 , _c( sqrt(_e2m) * exp(eatanhe(real(1))) )
00033 , _k0(k0)
00034 {
00035 if (!(_a > 0))
00036 throw GeographicErr("Major radius is not positive");
00037 if (!(_k0 > 0))
00038 throw GeographicErr("Scale is not positive");
00039 }
00040
00041 const PolarStereographic
00042 PolarStereographic::UPS(Constants::WGS84_a(), Constants::WGS84_r(),
00043 Constants::UPS_k0());
00044
00045 void PolarStereographic::Forward(bool northp, real lat, real lon,
00046 real& x, real& y, real& gamma, real& k)
00047 const throw() {
00048 real theta = 90 - (northp ? lat : -lat);
00049 real rho;
00050 theta *= Constants::degree();
00051 real
00052
00053 ctheta = cos(theta),
00054 f = exp(eatanhe(ctheta)),
00055 t2 = 2 * tan(theta/2) * f,
00056 m = sin(theta) / sqrt(1 - _e2 * sq(ctheta));
00057 rho = _a * _k0 * t2 / _c;
00058 k = m < numeric_limits<real>::epsilon() ? _k0 : rho / (_a * m);
00059 lon = lon >= 180 ? lon - 360 : lon < -180 ? lon + 360 : lon;
00060 real
00061 lam = lon * Constants::degree();
00062 x = rho * (lon == -180 ? 0 : sin(lam));
00063 y = (northp ? -rho : rho) * (abs(lon) == 90 ? 0 : cos(lam));
00064 gamma = northp ? lon : -lon;
00065 }
00066
00067 void PolarStereographic::Reverse(bool northp, real x, real y,
00068 real& lat, real& lon, real& gamma, real& k)
00069 const throw() {
00070 real
00071 rho = Math::hypot(x, y),
00072 t2 = rho * _c / (_a * _k0),
00073 theta = 2 * atan(t2/2);
00074
00075
00076
00077
00078 for (int i = 0; i < numit; ++i) {
00079 real
00080 ctheta = cos(theta),
00081 f = exp(eatanhe(ctheta)),
00082 v = 2 * tan(theta/2) * f - t2,
00083
00084 dv = _e2m * f / ((1 - _e2 * sq(ctheta)) * sq(cos(theta/2))),
00085 dtheta = -v/dv;
00086 theta += dtheta;
00087 if (abs(dtheta) < tol)
00088 break;
00089 }
00090 lat = (northp ? 1 : -1) * (90 - theta / Constants::degree());
00091
00092 lon = -atan2( -x, northp ? -y : y ) / Constants::degree();
00093 real m = sin(theta) / sqrt(1 - _e2 * sq(cos(theta)));
00094 k = m == 0 ? _k0 : rho / (_a * m);
00095 gamma = northp ? lon : -lon;
00096 }
00097
00098 void PolarStereographic::SetScale(real lat, real k) {
00099 if (!(k > 0))
00100 throw GeographicErr("Scale is not positive");
00101 real x, y, gamma, kold;
00102 Forward(true, lat, 0, x, y, gamma, kold);
00103 _k0 *= k/kold;
00104 }
00105
00106 }