GeographicLib  1.21
PolarStereographic.hpp
Go to the documentation of this file.
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&ndash;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