Constants.hpp

Go to the documentation of this file.
00001 /**
00002  * \file Constants.hpp
00003  * \brief Header for GeographicLib::Constants 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_CONSTANTS_HPP)
00011 #define GEOGRAPHICLIB_CONSTANTS_HPP "$Id: Constants.hpp 6827 2010-05-20 19:56:18Z karney $"
00012 
00013 /**
00014  * A simple compile-time assert.  This is designed to be compatible with the
00015  * C++0X static_assert.
00016  **********************************************************************/
00017 #if !defined(STATIC_ASSERT)
00018 #define STATIC_ASSERT(cond,reason) { enum{ STATIC_ASSERT_ENUM=1/int(cond) }; }
00019 #endif
00020 
00021 #if defined(__GNUC__)
00022 // Suppress "defined but not used" warnings
00023 #define RCSID_DECL(x) namespace \
00024 { char VAR_ ## x [] __attribute__((unused)) = x; }
00025 #else
00026 /**
00027  * Insertion of RCS Id strings into the object file.
00028  **********************************************************************/
00029 #define RCSID_DECL(x) namespace { char VAR_ ## x [] = x; }
00030 #endif
00031 
00032 RCSID_DECL(GEOGRAPHICLIB_CONSTANTS_HPP)
00033 
00034 #if !defined(GEOGRAPHICLIB_PREC)
00035 /**
00036  * The precision of floating point numbers used in %GeographicLib.  0 means
00037  * float; 1 (default) means double; 2 means long double.  Nearly all the
00038  * testing has been carried out with doubles and that's the recommended
00039  * configuration.  Note that with Microsoft Visual Studio, long double is the
00040  * same as double.
00041  **********************************************************************/
00042 #define GEOGRAPHICLIB_PREC 1
00043 #endif
00044 
00045 #include <cmath>
00046 #include <limits>
00047 #include <algorithm>
00048 #include <stdexcept>
00049 
00050 /**
00051  * \brief Namespace for %GeographicLib
00052  *
00053  * All of %GeographicLib is defined within the GeographicLib namespace.  In
00054  * addtion all the header files are included via %GeographicLib/filename.  This
00055  * minimizes the likelihood of conflicts with other packages.
00056  **********************************************************************/
00057 namespace GeographicLib {
00058 
00059   /**
00060    * \brief Mathematical functions needed by %GeographicLib
00061    *
00062    * Define mathematical functions in order to localize system dependencies and
00063    * to provide generic versions of the functions.  In addition define a real
00064    * type to be used by %GeographicLib.
00065    **********************************************************************/
00066   class Math {
00067   private:
00068     void dummy() {
00069       STATIC_ASSERT((GEOGRAPHICLIB_PREC) >= 0 && (GEOGRAPHICLIB_PREC) <= 2,
00070                     "Bad value of precision");
00071     }
00072     Math();                     // Disable constructor
00073   public:
00074 
00075 #if GEOGRAPHICLIB_PREC == 1
00076     /**
00077      * The real type for %GeographicLib. Nearly all the testing has been done
00078      * with \e real = double.  However, the algorithms should also work with
00079      * float and long double (where available).
00080      **********************************************************************/
00081     typedef double real;
00082 #elif GEOGRAPHICLIB_PREC == 0
00083     typedef float real;
00084 #elif GEOGRAPHICLIB_PREC == 2
00085     typedef long double real;
00086 #else
00087     typedef double real;
00088 #endif
00089 
00090 #if defined(DOXYGEN)
00091     /**
00092      * Return sqrt(\e x<sup>2</sup> + \e y<sup>2</sup>) avoiding underflow and
00093      * overflow.
00094      **********************************************************************/
00095     static inline real hypot(real x, real y) throw() {
00096       x = std::abs(x);
00097       y = std::abs(y);
00098       real
00099         a = std::max(x, y),
00100         b = std::min(x, y) / a;
00101       return a * std::sqrt(1 + b * b);
00102     }
00103 #elif defined(_MSC_VER)
00104     static inline double hypot(double x, double y) throw()
00105     { return _hypot(x, y); }
00106     static inline float hypot(float x, float y) throw()
00107     { return _hypotf(x, y); }
00108 #else
00109     // Use overloading to define generic versions
00110     static inline double hypot(double x, double y) throw()
00111     { return ::hypot(x, y); }
00112     static inline float hypot(float x, float y) throw()
00113     { return ::hypotf(x, y); }
00114 #if !defined(__NO_LONG_DOUBLE_MATH)
00115     static inline long double hypot(long double x, long double y) throw()
00116     { return ::hypotl(x, y); }
00117 #endif
00118 #endif
00119 
00120 #if defined(DOXYGEN) || defined(_MSC_VER)
00121     /**
00122      * Return exp(\e x) - 1 accurate near \e x = 0.  This is taken from
00123      * N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd
00124      * Edition (SIAM, 2002), Sec 1.14.1, p 19.
00125      **********************************************************************/
00126     static inline real expm1(real x) throw() {
00127       volatile real
00128         y = std::exp(x),
00129         z = y - 1;
00130       // The reasoning here is similar to that for log1p.  The expression
00131       // mathematically reduces to exp(x) - 1, and the factor z/log(y) = (y -
00132       // 1)/log(y) is a slowly varying quantity near y = 1 and is accurately
00133       // computed.
00134       return std::abs(x) > 1 ? z : z == 0 ?  x : x * z / std::log(y);
00135     }
00136 #else
00137     static inline double expm1(double x) throw() { return ::expm1(x); }
00138     static inline float expm1(float x) throw() { return ::expm1f(x); }
00139 #if !defined(__NO_LONG_DOUBLE_MATH)
00140     static inline long double expm1(long double x) throw()
00141     { return ::expm1l(x); }
00142 #endif
00143 #endif
00144 
00145 #if defined(DOXYGEN) || defined(_MSC_VER)
00146     /**
00147      * Return log(\e x + 1) accurate near \e x = 0.  This is taken See
00148      * D. Goldberg,
00149      * <a href="http://docs.sun.com/source/806-3568/ncg_goldberg.html"> What
00150      * every computer scientist should know about floating-point arithmetic</a>
00151      * (1991), Theorem 4.  See also, Higham (op. cit.), Answer to Problem 1.5,
00152      * p 528.
00153      **********************************************************************/
00154     static inline real log1p(real x) throw() {
00155       volatile real
00156         y = 1 + x,
00157         z = y - 1;
00158       // Here's the explanation for this magic: y = 1 + z, exactly, and z
00159       // approx x, thus log(y)/z (which is nearly constant near z = 0) returns
00160       // a good approximation to the true log(1 + x)/x.  The multiplication x *
00161       // (log(y)/z) introduces little additional error.
00162       return z == 0 ? x : x * std::log(y) / z;
00163     }
00164 #else
00165     static inline double log1p(double x) throw() { return ::log1p(x); }
00166     static inline float log1p(float x) throw() { return ::log1pf(x); }
00167 #if !defined(__NO_LONG_DOUBLE_MATH)
00168     static inline long double log1p(long double x) throw()
00169     { return ::log1pl(x); }
00170 #endif
00171 #endif
00172 
00173 #if defined(DOXYGEN) || defined(_MSC_VER)
00174     /**
00175      * Return asinh(\e x).  This is defined in terms of Math::log1p(\e x) in
00176      * order to maintain accuracy near \e x = 0.  In addition, the odd parity
00177      * of the function is enforced.
00178      **********************************************************************/
00179     static inline real asinh(real x) throw() {
00180       real y = std::abs(x);     // Enforce odd parity
00181       y = log1p(y * (1 + y/(hypot(real(1), y) + 1)));
00182       return x < 0 ? -y : y;
00183     }
00184 #else
00185     static inline double asinh(double x) throw() { return ::asinh(x); }
00186     static inline float asinh(float x) throw() { return ::asinhf(x); }
00187 #if !defined(__NO_LONG_DOUBLE_MATH)
00188     static inline long double asinh(long double x) throw()
00189     { return ::asinhl(x); }
00190 #endif
00191 #endif
00192 
00193 #if defined(DOXYGEN) || defined(_MSC_VER)
00194     /**
00195      * Return atanh(\e x).  This is defined in terms of Math::log1p(\e x) in
00196      * order to maintain accuracy near \e x = 0.  In addition, the odd parity
00197      * of the function is enforced.
00198      **********************************************************************/
00199     static inline real atanh(real x) throw() {
00200       real y = std::abs(x);     // Enforce odd parity
00201       y = log1p(2 * y/(1 - y))/2;
00202       return x < 0 ? -y : y;
00203     }
00204 #else
00205     static inline double atanh(double x) throw() { return ::atanh(x); }
00206     static inline float atanh(float x) throw() { return ::atanhf(x); }
00207 #if !defined(__NO_LONG_DOUBLE_MATH)
00208     static inline long double atanh(long double x) throw()
00209     { return ::atanhl(x); }
00210 #endif
00211 #endif
00212 
00213 #if defined(DOXYGEN) || defined(_MSC_VER)
00214     /**
00215      * Return the real cube root of \e x.
00216      **********************************************************************/
00217     static inline real cbrt(real x) throw() {
00218       real y = std::pow(std::abs(x), 1/real(3)); // Return the real cube root
00219       return x < 0 ? -y : y;
00220     }
00221 #else
00222     static inline double cbrt(double x) throw() { return ::cbrt(x); }
00223     static inline float cbrt(float x) throw() { return ::cbrtf(x); }
00224 #if !defined(__NO_LONG_DOUBLE_MATH)
00225     static inline long double cbrt(long double x) throw() { return ::cbrtl(x); }
00226 #endif
00227 #endif
00228 
00229     /**
00230      * Return true if number is finite, false if NaN or infinite.
00231      **********************************************************************/
00232     static inline bool isfinite(real x) throw() {
00233 #if defined(DOXYGEN)
00234       return std::abs(x) <= std::numeric_limits<real>::max();
00235 #elif defined(_MSC_VER)
00236       return _finite(x) != 0;
00237 #else
00238       return std::isfinite(x) != 0;
00239 #endif
00240     }
00241   };
00242 
00243   /**
00244    * \brief %Constants needed by %GeographicLib
00245    *
00246    * Define constants specifying the WGS84 ellipsoid, the UTM and UPS
00247    * projections, and various unit conversions.
00248    **********************************************************************/
00249   class Constants {
00250   private:
00251     typedef Math::real real;
00252     Constants();                // Disable constructor
00253 
00254   public:
00255     /**
00256      * pi
00257      **********************************************************************/
00258     static inline Math::real pi() throw()
00259     // good for about 123-bit accuracy
00260     { return real(3.141592653589793238462643383279502884L); }
00261     /**
00262      * Factor to convert from degrees to radians
00263      **********************************************************************/
00264     static inline Math::real degree() throw() { return pi() / 180; }
00265     /**
00266      * Factor to convert from minutes to radians
00267      **********************************************************************/
00268     static inline Math::real arcminute() throw() { return degree() / 60; }
00269     /**
00270      * Factor to convert from seconds to radians
00271      **********************************************************************/
00272     static inline Math::real arcsecond() throw() { return arcminute() / 60; }
00273 
00274     /** \name Ellipsoid parameters
00275      **********************************************************************/
00276     ///@{
00277     /**
00278      * Major radius of WGS84 ellipsoid
00279      **********************************************************************/
00280     static inline Math::real WGS84_a() throw() { return 6378137 * meter(); }
00281     /**
00282      * Reciprocal flattening of WGS84 ellipsoid
00283      **********************************************************************/
00284     static inline Math::real WGS84_r() throw() { return real(298.257223563L); }
00285     /**
00286      * Central scale factor for UTM
00287      **********************************************************************/
00288     static inline Math::real UTM_k0() throw() {return real(0.9996L); }
00289     /**
00290      * Central scale factor for UPS
00291      **********************************************************************/
00292     static inline Math::real UPS_k0() throw() { return real(0.994L); }
00293     ///@}
00294 
00295     /** \name SI units
00296      **********************************************************************/
00297     ///@{
00298     /**
00299      * Factor to convert from meters to meters (i.e., 1, but this lets the
00300      * internal system of units be changed if necessary).
00301      **********************************************************************/
00302     static inline Math::real meter() throw() { return real(1); }
00303     /**
00304      * Factor to convert from kilometers to meters.
00305      **********************************************************************/
00306     static inline Math::real kilometer() throw() { return 1000 * meter(); }
00307     ///@}
00308 
00309     /**
00310      * Factor to convert from nautical miles (approximately 1 arc minute) to
00311      * meters.
00312      **********************************************************************/
00313     static inline Math::real nauticalmile() throw() { return 1852 * meter(); }
00314 
00315     /** \name Anachronistic British units
00316      **********************************************************************/
00317     ///@{
00318     /**
00319      * Factor to convert from international feet to meters.
00320      **********************************************************************/
00321     static inline Math::real foot() throw()
00322     { return real(0.0254L) * 12 * meter(); }
00323     /**
00324      * Factor to convert from yards to meters.
00325      **********************************************************************/
00326     static inline Math::real yard() throw() { return 3 * foot(); }
00327     /**
00328      * Factor to convert from fathoms to meters.
00329      **********************************************************************/
00330     static inline Math::real fathom() throw() { return 2 * yard(); }
00331     /**
00332      * Factor to convert from chains to meters.
00333      **********************************************************************/
00334     static inline Math::real chain() throw() { return 22 * yard(); }
00335     /**
00336      * Factor to convert from furlongs to meters.
00337      **********************************************************************/
00338     static inline Math::real furlong() throw() { return 10 * chain(); }
00339     /**
00340      * Factor to convert from statute miles to meters.
00341      **********************************************************************/
00342     static inline Math::real mile() throw() { return 8 * furlong(); }
00343     ///@}
00344 
00345     /** \name Anachronistic US units
00346      **********************************************************************/
00347     ///@{
00348     /**
00349      * Factor to convert from US survery feet to meters.
00350      **********************************************************************/
00351     static inline Math::real surveyfoot() throw()
00352     { return real(1200) / real(3937) * meter(); }
00353     ///@}
00354   };
00355 
00356   /**
00357    * \brief %Exception handling for %GeographicLib
00358    *
00359    * A class to handle exceptions.  It's derived off std::runtime_error so it
00360    * can be caught by the usual catch clauses.
00361    **********************************************************************/
00362   class GeographicErr : public std::runtime_error {
00363   public:
00364 
00365     /**
00366      * Constructor takes a string message, \e msg, which is accessible in the
00367      * catch clause, via what().
00368      **********************************************************************/
00369     GeographicErr(const std::string& msg) : std::runtime_error(msg) {}
00370   };
00371 
00372 } // namespace GeographicLib
00373 
00374 #endif

Generated on 21 May 2010 for GeographicLib by  doxygen 1.6.1