TransverseMercatorExact.hpp

Go to the documentation of this file.
00001 /**
00002  * \file TransverseMercatorExact.hpp
00003  * \brief Header for GeographicLib::TransverseMercatorExact 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_TRANSVERSEMERCATOREXACT_HPP)
00011 #define GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_HPP "$Id: TransverseMercatorExact.hpp 6824 2010-04-19 14:25:10Z karney $"
00012 
00013 #include "GeographicLib/Constants.hpp"
00014 #include "GeographicLib/EllipticFunction.hpp"
00015 
00016 namespace GeographicLib {
00017 
00018   /**
00019    * \brief An exact implementation of the Transverse Mercator Projection
00020    *
00021    * Implementation of the Transverse Mercator Projection given in
00022    *  - L. P. Lee,
00023    *    <a href="http://dx.doi.org/10.3138/X687-1574-4325-WM62"> Conformal
00024    *    Projections Based On Jacobian Elliptic Functions</a>, Part V of
00025    *    Conformal Projections Based on Elliptic Functions,
00026    *    (B. V. Gutsell, Toronto, 1976), 128pp.,
00027    *    ISBN: 0919870163
00028    *    (also appeared as:
00029    *    Monograph 16, Suppl. No. 1 to Canadian Cartographer, Vol 13).
00030    *
00031    * This method gives the correct results for forward and reverse
00032    * transformations subject to the branch cut rules (see the description of
00033    * the \e extendp argument to the constructor).  The maximum error is about 8
00034    * nm (ground distance) for the forward and reverse transformations.  The
00035    * error in the convergence is 2e-15&quot;, the relative error in the scale
00036    * is 7e-12%%.  (See \ref tmerrors for the weasel words.)  The method is
00037    * "exact" in the sense that the errors are close to the round-off limit and
00038    * that no changes are needed in the algorithms for them to be used with
00039    * reals of a higher precision.  Thus the errors using long double (with a
00040    * 64-bit fraction) are about 2000 times smaller than using double (with a
00041    * 53-bit fraction).
00042    *
00043    * This algorithm is about 4.5 times slower than the 6th-order Kr&uuml;ger
00044    * method, GeographicLib::TransverseMercator, taking about 11 us for a
00045    * combined forward and reverse projection on a 2.6 GHz Intel machine (g++,
00046    * version 4.3.0, -O3).
00047    *
00048    * The ellipsoid parameters and the central scale are set in the constructor.
00049    * The central meridian (which is a trivial shift of the longitude) is
00050    * specified as the \e lon0 argument of the TransverseMercatorExact::Forward
00051    * and TransverseMercatorExact::Reverse functions.  The latitude of origin is
00052    * taken to be the equator.  See the documentation on
00053    * GeographicLib::TransverseMercator for how to include a false easting,
00054    * false northing, or a latitude of origin.
00055    *
00056    * See TransverseMercatorExact.cpp for more information on the
00057    * implementation.
00058    *
00059    * See \ref transversemercator for a discussion of this projection.
00060    **********************************************************************/
00061 
00062   class TransverseMercatorExact {
00063   private:
00064     typedef Math::real real;
00065     static const real tol, tol1, tol2, taytol, overflow;
00066     static const int numit = 10;
00067     const real _a, _r, _f, _k0, _mu, _mv, _e, _ep2;
00068     const bool _extendp;
00069     const EllipticFunction _Eu, _Ev;
00070     static inline real sq(real x) throw() { return x * x; }
00071     // tan(x) for x in [-pi/2, pi/2] ensuring that the sign is right
00072     static inline real tanx(real x) throw() {
00073       real t = std::tan(x);
00074       return x >= 0 ? (t >= 0 ? t : overflow) : (t < 0 ? t : -overflow);
00075     }
00076 
00077     real taup(real tau) const throw();
00078     real taupinv(real taup) const throw();
00079 
00080     void zeta(real u, real snu, real cnu, real dnu,
00081               real v, real snv, real cnv, real dnv,
00082               real& taup, real& lam) const throw();
00083 
00084     void dwdzeta(real u, real snu, real cnu, real dnu,
00085                  real v, real snv, real cnv, real dnv,
00086                  real& du, real& dv) const throw();
00087 
00088     bool zetainv0(real psi, real lam, real& u, real& v) const throw();
00089     void zetainv(real taup, real lam, real& u, real& v) const throw();
00090 
00091     void sigma(real u, real snu, real cnu, real dnu,
00092                real v, real snv, real cnv, real dnv,
00093                real& xi, real& eta) const throw();
00094 
00095     void dwdsigma(real u, real snu, real cnu, real dnu,
00096                   real v, real snv, real cnv, real dnv,
00097                   real& du, real& dv) const throw();
00098 
00099     bool sigmainv0(real xi, real eta, real& u, real& v) const throw();
00100     void sigmainv(real xi, real eta, real& u, real& v) const throw();
00101 
00102     void Scale(real tau, real lam,
00103                real snu, real cnu, real dnu,
00104                real snv, real cnv, real dnv,
00105                real& gamma, real& k) const throw();
00106 
00107   public:
00108 
00109     /**
00110      * Constructor for a ellipsoid radius \e a (meters), reciprocal flattening
00111      * \e r, and central scale factor \e k0.  The transverse Mercator
00112      * projection has a branch point singularity at \e lat = 0 and \e lon - \e
00113      * lon0 = 90 (1 - \e e) or (for TransverseMercatorExact::UTM) x = 18381 km,
00114      * y = 0m.  The \e extendp argument governs where the branch cut is placed.
00115      * With \e extendp = false, the "standard" convention is followed, namely
00116      * the cut is placed along x > 18381 km, y = 0m.  Forward can be called
00117      * with any \e lat and \e lon then produces the transformation shown in
00118      * Lee, Fig 46.  Reverse analytically continues this in the +/- \e x
00119      * direction.  As a consequence, Reverse may map multiple points to the
00120      * same geographic location; for example, for TransverseMercatorExact::UTM,
00121      * \e x = 22051449.037349 m, \e y = -7131237.022729 m and \e x =
00122      * 29735142.378357 m, \e y = 4235043.607933 m both map to \e lat = -2 deg,
00123      * \e lon = 88 deg.
00124      *
00125      * With \e extendp = true, the branch cut is moved to the lower left
00126      * quadrant.  The various symmetries of the transverse Mercator projection
00127      * can be used to explore the projection on any sheet.  In this mode the
00128      * domains of \e lat, \e lon, \e x, and \e y are restricted to
00129      * - the union of
00130      *   - \e lat in [0, 90] and \e lon - \e lon0 in [0, 90]
00131      *   - \e lat in (-90, 0] and \e lon - \e lon0 in [90 (1 - \e e), 90]
00132      * - the union of
00133      *   - \e x/(\e k0 \e a) in [0, inf) and
00134      *     \e y/(\e k0 \e a) in [0, E(\e e^2)]
00135      *   - \e x/(\e k0 \e a) in [K(1 - \e e^2) - E(1 - \e e^2), inf) and
00136      *     \e y/(\e k0 \e a) in (-inf, 0]
00137      * .
00138      * See \ref extend for a full discussion of the treatment of the branch
00139      * cut.
00140      *
00141      * The method will work for all ellipsoids used in terrestial geodesy.  The
00142      * method cannot be applied directly to the case of a sphere (\e r = inf)
00143      * because some the constants characterizing this method diverge in that
00144      * limit, and in practise, \e r should be smaller than about
00145      * 1/std::numeric_limits<real>::epsilon().  However,
00146      * GeographicLib::TransverseMercator treats the sphere exactly.  An
00147      * exception is thrown if \e a, \e r, or \e k0 is not positive.
00148      **********************************************************************/
00149     TransverseMercatorExact(real a, real r, real k0, bool extendp = false);
00150 
00151     /**
00152      * Convert from latitude \e lat (degrees) and longitude \e lon (degrees) to
00153      * transverse Mercator easting \e x (meters) and northing \e y (meters).
00154      * The central meridian of the transformation is \e lon0 (degrees).  Also
00155      * return the meridian convergence \e gamma (degrees) and the scale \e k.
00156      * No false easting or northing is added.  \e lat should be in the range
00157      * [-90, 90]; \e lon and \e lon0 should be in the range [-180, 360].
00158      **********************************************************************/
00159     void Forward(real lon0, real lat, real lon,
00160                  real& x, real& y, real& gamma, real& k) const throw();
00161 
00162     /**
00163      * Convert from transverse Mercator easting \e x (meters) and northing \e y
00164      * (meters) to latitude \e lat (degrees) and longitude \e lon (degrees) .
00165      * The central meridian of the transformation is \e lon0 (degrees).  Also
00166      * return the meridian convergence \e gamma (degrees) and the scale \e k.
00167      * No false easting or northing is added.  \e lon0 should be in the range
00168      * [-180, 360].  The value of \e lon returned is in the range [-180, 180).
00169      **********************************************************************/
00170     void Reverse(real lon0, real x, real y,
00171                  real& lat, real& lon, real& gamma, real& k) const throw();
00172 
00173     /**
00174      * The major radius of the ellipsoid (meters).  This is that value of \e a
00175      * used in the constructor.
00176      **********************************************************************/
00177     Math::real MajorRadius() const throw() { return _a; }
00178 
00179     /**
00180      * The inverse flattening of the ellipsoid.  This is that value of \e r
00181      * used in the constructor.
00182      **********************************************************************/
00183     Math::real InverseFlattening() const throw() { return _r; }
00184 
00185     /**
00186      * The central scale for the projection.  This is that value of \e k0 used
00187      * in the constructor and is the scale on the central meridian.
00188      **********************************************************************/
00189     Math::real CentralScale() const throw() { return _k0; }
00190 
00191     /**
00192      * A global instantiation of TransverseMercatorExact with the WGS84
00193      * ellipsoid and the UTM scale factor.  However, unlike UTM, no false
00194      * easting or northing is added.
00195      **********************************************************************/
00196     const static TransverseMercatorExact UTM;
00197   };
00198 
00199 } // namespace GeographicLib
00200 
00201 #endif

Generated on 21 May 2010 for GeographicLib by  doxygen 1.6.1