00001 /** 00002 * \file PolarStereographic.hpp 00003 * \brief Header for GeographicLib::PolarStereographic class 00004 * 00005 * Copyright (c) Charles Karney (2008-2014) <charles@karney.com> and licensed 00006 * under the MIT/X11 License. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #if !defined(GEOGRAPHICLIB_POLARSTEREOGRAPHIC_HPP) 00011 #define GEOGRAPHICLIB_POLARSTEREOGRAPHIC_HPP 1 00012 00013 #include <GeographicLib/Constants.hpp> 00014 00015 namespace GeographicLib { 00016 00017 /** 00018 * \brief Polar stereographic projection 00019 * 00020 * Implementation taken from the report, 00021 * - J. P. Snyder, 00022 * <a href="http://pubs.er.usgs.gov/usgspubs/pp/pp1395"> Map Projections: A 00023 * Working Manual</a>, USGS Professional Paper 1395 (1987), 00024 * pp. 160--163. 00025 * 00026 * This is a straightforward implementation of the equations in Snyder except 00027 * that Newton's method is used to invert the projection. 00028 * 00029 * Example of use: 00030 * \include example-PolarStereographic.cpp 00031 **********************************************************************/ 00032 class GEOGRAPHICLIB_EXPORT PolarStereographic { 00033 private: 00034 typedef Math::real real; 00035 real tol_; 00036 // _Cx used to be _C but g++ 3.4 has a macro of that name 00037 real _a, _f, _e2, _e, _e2m, _Cx, _c; 00038 real _k0; 00039 static const int numit_ = 5; 00040 static inline real overflow() { 00041 // Overflow value s.t. atan(overflow_) = pi/2 00042 static const real 00043 overflow = 1 / Math::sq(std::numeric_limits<real>::epsilon()); 00044 return overflow; 00045 } 00046 // tan(x) for x in [-pi/2, pi/2] ensuring that the sign is right 00047 static inline real tanx(real x) { 00048 using std::tan; 00049 real t = tan(x); 00050 // Write the tests this way to ensure that tanx(NaN()) is NaN() 00051 return x >= 0 ? 00052 (!(t < 0) ? t : overflow()) : 00053 (!(t >= 0) ? t : -overflow()); 00054 } 00055 // Return e * atanh(e * x) for f >= 0, else return 00056 // - sqrt(-e2) * atan( sqrt(-e2) * x) for f < 0 00057 inline real eatanhe(real x) const { 00058 using std::atan; 00059 return _f >= 0 ? _e * Math::atanh(_e * x) : - _e * atan(_e * x); 00060 } 00061 public: 00062 00063 /** 00064 * Constructor for a ellipsoid with 00065 * 00066 * @param[in] a equatorial radius (meters). 00067 * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere. 00068 * Negative \e f gives a prolate ellipsoid. If \e f > 1, set 00069 * flattening to 1/\e f. 00070 * @param[in] k0 central scale factor. 00071 * @exception GeographicErr if \e a, (1 − \e f) \e a, or \e k0 is 00072 * not positive. 00073 **********************************************************************/ 00074 PolarStereographic(real a, real f, real k0); 00075 00076 /** 00077 * Set the scale for the projection. 00078 * 00079 * @param[in] lat (degrees) assuming \e northp = true. 00080 * @param[in] k scale at latitude \e lat (default 1). 00081 * @exception GeographicErr \e k is not positive. 00082 * @exception GeographicErr if \e lat is not in (−90°, 00083 * 90°]. 00084 **********************************************************************/ 00085 void SetScale(real lat, real k = real(1)); 00086 00087 /** 00088 * Forward projection, from geographic to polar stereographic. 00089 * 00090 * @param[in] northp the pole which is the center of projection (true means 00091 * north, false means south). 00092 * @param[in] lat latitude of point (degrees). 00093 * @param[in] lon longitude of point (degrees). 00094 * @param[out] x easting of point (meters). 00095 * @param[out] y northing of point (meters). 00096 * @param[out] gamma meridian convergence at point (degrees). 00097 * @param[out] k scale of projection at point. 00098 * 00099 * No false easting or northing is added. \e lat should be in the range 00100 * (−90°, 90°] for \e northp = true and in the range 00101 * [−90°, 90°) for \e northp = false; \e lon should 00102 * be in the range [−540°, 540°). 00103 **********************************************************************/ 00104 void Forward(bool northp, real lat, real lon, 00105 real& x, real& y, real& gamma, real& k) const; 00106 00107 /** 00108 * Reverse projection, from polar stereographic to geographic. 00109 * 00110 * @param[in] northp the pole which is the center of projection (true means 00111 * north, false means south). 00112 * @param[in] x easting of point (meters). 00113 * @param[in] y northing of point (meters). 00114 * @param[out] lat latitude of point (degrees). 00115 * @param[out] lon longitude of point (degrees). 00116 * @param[out] gamma meridian convergence at point (degrees). 00117 * @param[out] k scale of projection at point. 00118 * 00119 * No false easting or northing is added. The value of \e lon returned is 00120 * in the range [−180°, 180°). 00121 **********************************************************************/ 00122 void Reverse(bool northp, real x, real y, 00123 real& lat, real& lon, real& gamma, real& k) const; 00124 00125 /** 00126 * PolarStereographic::Forward without returning the convergence and scale. 00127 **********************************************************************/ 00128 void Forward(bool northp, real lat, real lon, 00129 real& x, real& y) const { 00130 real gamma, k; 00131 Forward(northp, lat, lon, x, y, gamma, k); 00132 } 00133 00134 /** 00135 * PolarStereographic::Reverse without returning the convergence and scale. 00136 **********************************************************************/ 00137 void Reverse(bool northp, real x, real y, 00138 real& lat, real& lon) const { 00139 real gamma, k; 00140 Reverse(northp, x, y, lat, lon, gamma, k); 00141 } 00142 00143 /** \name Inspector functions 00144 **********************************************************************/ 00145 ///@{ 00146 /** 00147 * @return \e a the equatorial radius of the ellipsoid (meters). This is 00148 * the value used in the constructor. 00149 **********************************************************************/ 00150 Math::real MajorRadius() const { return _a; } 00151 00152 /** 00153 * @return \e f the flattening of the ellipsoid. This is the value used in 00154 * the constructor. 00155 **********************************************************************/ 00156 Math::real Flattening() const { return _f; } 00157 00158 /// \cond SKIP 00159 /** 00160 * <b>DEPRECATED</b> 00161 * @return \e r the inverse flattening of the ellipsoid. 00162 **********************************************************************/ 00163 Math::real InverseFlattening() const { return 1/_f; } 00164 /// \endcond 00165 00166 /** 00167 * The central scale for the projection. This is the value of \e k0 used 00168 * in the constructor and is the scale at the pole unless overridden by 00169 * PolarStereographic::SetScale. 00170 **********************************************************************/ 00171 Math::real CentralScale() const { return _k0; } 00172 ///@} 00173 00174 /** 00175 * A global instantiation of PolarStereographic with the WGS84 ellipsoid 00176 * and the UPS scale factor. However, unlike UPS, no false easting or 00177 * northing is added. 00178 **********************************************************************/ 00179 static const PolarStereographic& UPS(); 00180 }; 00181 00182 } // namespace GeographicLib 00183 00184 #endif // GEOGRAPHICLIB_POLARSTEREOGRAPHIC_HPP