00001 /** 00002 * \file PolarStereographic.hpp 00003 * \brief Header for GeographicLib::PolarStereographic class 00004 * 00005 * Copyright (c) Charles Karney (2008, 2009, 2010) <charles@karney.com> 00006 * and licensed under the LGPL. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #if !defined(GEOGRAPHICLIB_POLARSTEREOGRAPHIC_HPP) 00011 #define GEOGRAPHICLIB_POLARSTEREOGRAPHIC_HPP "$Id: PolarStereographic.hpp 6906 2010-12-02 22:10:56Z karney $" 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 class PolarStereographic { 00030 private: 00031 typedef Math::real real; 00032 // _Cx used to be _C but g++ 3.4 has a macro of that name 00033 const real _a, _r, _f, _e2, _e, _e2m, _Cx, _c; 00034 real _k0; 00035 static const real tol, overflow; 00036 static const int numit = 5; 00037 static inline real sq(real x) throw() { return x * x; } 00038 // Return e * atanh(e * x) for f >= 0, else return 00039 // - sqrt(-e2) * atan( sqrt(-e2) * x) for f < 0 00040 inline real eatanhe(real x) const throw() { 00041 return _f >= 0 ? _e * Math::atanh(_e * x) : - _e * std::atan(_e * x); 00042 } 00043 public: 00044 00045 /** 00046 * Constructor for a ellipsoid with 00047 * 00048 * @param[in] a equatorial radius (meters) 00049 * @param[in] r reciprocal flattening. Setting \e r = 0 implies \e r = inf 00050 * or flattening = 0 (i.e., a sphere). Negative \e r indicates a prolate 00051 * ellipsoid. 00052 * @param[in] k0 central scale factor. 00053 * 00054 * An exception is thrown if either of the axes of the ellipsoid is 00055 * not positive \e a or if \e k0 is not positive. 00056 **********************************************************************/ 00057 PolarStereographic(real a, real r, real k0); 00058 00059 /** 00060 * Set the scale for the projection. 00061 * 00062 * @param[in] lat (degrees) assuming \e northp = true. 00063 * @param[in] k scale at latitude \e lat (default 1). 00064 * 00065 * This allows a "latitude of true scale" to be specified. An exception is 00066 * thrown if \e k is not positive or if \e lat is not in the range (-90, 00067 * 90]. 00068 **********************************************************************/ 00069 void SetScale(real lat, real k = real(1)); 00070 00071 /** 00072 * Forward projection, from geographic to polar stereographic. 00073 * 00074 * @param[in] northp the pole which is the center of projection (true means 00075 * north, false means south). 00076 * @param[in] lat latitude of point (degrees). 00077 * @param[in] lon longitude of point (degrees). 00078 * @param[out] x easting of point (meters). 00079 * @param[out] y northing of point (meters). 00080 * @param[out] gamma meridian convergence at point (degrees). 00081 * @param[out] k scale of projection at point. 00082 * 00083 * No false easting or northing is added. \e lat should be in the range 00084 * (-90, 90] for \e northp = true and in the range [-90, 90) for \e northp 00085 * = false; \e lon should be in the range [-180, 360]. 00086 **********************************************************************/ 00087 void Forward(bool northp, real lat, real lon, 00088 real& x, real& y, real& gamma, real& k) const throw(); 00089 00090 /** 00091 * Reverse projection, from polar stereographic to geographic. 00092 * 00093 * @param[in] northp the pole which is the center of projection (true means 00094 * north, false means south). 00095 * @param[in] x easting of point (meters). 00096 * @param[in] y northing of point (meters). 00097 * @param[out] lat latitude of point (degrees). 00098 * @param[out] lon longitude of point (degrees). 00099 * @param[out] gamma meridian convergence at point (degrees). 00100 * @param[out] k scale of projection at point. 00101 * 00102 * No false easting or northing is added. The value of \e lon returned is 00103 * in the range [-180, 180). 00104 **********************************************************************/ 00105 void Reverse(bool northp, real x, real y, 00106 real& lat, real& lon, real& gamma, real& k) const throw(); 00107 00108 /** 00109 * PolarStereographic::Forward without returning the convergence and scale. 00110 **********************************************************************/ 00111 void Forward(bool northp, real lat, real lon, 00112 real& x, real& y) const throw() { 00113 real gamma, k; 00114 Forward(northp, lat, lon, x, y, gamma, k); 00115 } 00116 00117 /** 00118 * PolarStereographic::Reverse without returning the convergence and scale. 00119 **********************************************************************/ 00120 void Reverse(bool northp, real x, real y, 00121 real& lat, real& lon) const throw() { 00122 real gamma, k; 00123 Reverse(northp, x, y, lat, lon, gamma, k); 00124 } 00125 00126 /** \name Inspector functions 00127 **********************************************************************/ 00128 ///@{ 00129 /** 00130 * @return \e a the equatorial radius of the ellipsoid (meters). This is 00131 * the value used in the constructor. 00132 **********************************************************************/ 00133 Math::real MajorRadius() const throw() { return _a; } 00134 00135 /** 00136 * @return \e r the inverse flattening of the ellipsoid. This is the 00137 * value used in the constructor. A value of 0 is returned for a sphere 00138 * (infinite inverse flattening). 00139 **********************************************************************/ 00140 Math::real InverseFlattening() const throw() { return _r; } 00141 00142 /** 00143 * The central scale for the projection. This is the value of \e k0 used 00144 * in the constructor and is the scale at the pole unless overridden by 00145 * PolarStereographic::SetScale. 00146 **********************************************************************/ 00147 Math::real CentralScale() const throw() { return _k0; } 00148 ///@} 00149 00150 /** 00151 * A global instantiation of PolarStereographic with the WGS84 ellipsoid 00152 * and the UPS scale factor. However, unlike UPS, no false easting or 00153 * northing is added. 00154 **********************************************************************/ 00155 static const PolarStereographic UPS; 00156 }; 00157 00158 } // namespace GeographicLib 00159 00160 #endif