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