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, 2011) <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 6967 2011-02-19 15:53:41Z karney $"
00012 
00013 /**
00014  * Are C++0X math functions available?
00015  **********************************************************************/
00016 #if !defined(GEOGRAPHICLIB_CPLUSPLUS0X_MATH)
00017 #if defined(__GXX_EXPERIMENTAL_CXX0X__)
00018 #define GEOGRAPHICLIB_CPLUSPLUS0X_MATH 1
00019 #else
00020 #define GEOGRAPHICLIB_CPLUSPLUS0X_MATH 0
00021 #endif
00022 #endif
00023 
00024 /**
00025  * A compile-time assert.  Use C++0X static_assert, if available.
00026  **********************************************************************/
00027 #if !defined(STATIC_ASSERT)
00028 #if defined(__GXX_EXPERIMENTAL_CXX0X__)
00029 #define STATIC_ASSERT static_assert
00030 #elif defined(_MSC_VER) && _MSC_VER >= 1600
00031 #define STATIC_ASSERT static_assert
00032 #else
00033 #define STATIC_ASSERT(cond,reason) { enum{ STATIC_ASSERT_ENUM=1/int(cond) }; }
00034 #endif
00035 #endif
00036 
00037 #if defined(__GNUC__)
00038 // Suppress "defined but not used" warnings
00039 #define RCSID_DECL(x) namespace \
00040 { char VAR_ ## x [] __attribute__((used)) = x; }
00041 #else
00042 /**
00043  * Insertion of RCS Id strings into the object file.
00044  **********************************************************************/
00045 #define RCSID_DECL(x) namespace { char VAR_ ## x [] = x; }
00046 #endif
00047 
00048 RCSID_DECL(GEOGRAPHICLIB_CONSTANTS_HPP)
00049 
00050 #if !defined(GEOGRAPHICLIB_PREC)
00051 /**
00052  * The precision of floating point numbers used in %GeographicLib.  0 means
00053  * float; 1 (default) means double; 2 means long double.  Nearly all the
00054  * testing has been carried out with doubles and that's the recommended
00055  * configuration.  Note that with Microsoft Visual Studio, long double is the
00056  * same as double.
00057  **********************************************************************/
00058 #define GEOGRAPHICLIB_PREC 1
00059 #endif
00060 
00061 #if defined(__CYGWIN__) && defined(__GNUC__) && __GNUC__ < 4
00062 // g++ 3.x under cygwin doesn't have long double
00063 #define __NO_LONG_DOUBLE_MATH 1
00064 #endif
00065 
00066 #include <cmath>
00067 #include <limits>
00068 #include <algorithm>
00069 #include <stdexcept>
00070 
00071 /**
00072  * \brief Namespace for %GeographicLib
00073  *
00074  * All of %GeographicLib is defined within the GeographicLib namespace.  In
00075  * addtion all the header files are included via %GeographicLib/filename.  This
00076  * minimizes the likelihood of conflicts with other packages.
00077  **********************************************************************/
00078 namespace GeographicLib {
00079 
00080   /**
00081    * \brief Mathematical functions needed by %GeographicLib
00082    *
00083    * Define mathematical functions in order to localize system dependencies and
00084    * to provide generic versions of the functions.  In addition define a real
00085    * type to be used by %GeographicLib.
00086    **********************************************************************/
00087   class Math {
00088   private:
00089     void dummy() {
00090       STATIC_ASSERT((GEOGRAPHICLIB_PREC) >= 0 && (GEOGRAPHICLIB_PREC) <= 2,
00091                     "Bad value of precision");
00092     }
00093     Math();                     // Disable constructor
00094   public:
00095 
00096 #if !defined(__NO_LONG_DOUBLE_MATH)
00097     /**
00098      * The extended precision type for real numbers, used for some testing.
00099      * This is long double on computers with this type; otherwise it is double.
00100      **********************************************************************/
00101     typedef long double extended;
00102 #else
00103     typedef double extended;
00104 #endif
00105 
00106 #if GEOGRAPHICLIB_PREC == 1
00107     /**
00108      * The real type for %GeographicLib. Nearly all the testing has been done
00109      * with \e real = double.  However, the algorithms should also work with
00110      * float and long double (where available).
00111      **********************************************************************/
00112     typedef double real;
00113 #elif GEOGRAPHICLIB_PREC == 0
00114     typedef float real;
00115 #elif GEOGRAPHICLIB_PREC == 2
00116     typedef extended real;
00117 #else
00118     typedef double real;
00119 #endif
00120 
00121     /**
00122      * @return \e pi
00123      **********************************************************************/
00124     template<typename T>
00125     static inline T pi() throw()
00126     // good for about 168-bit accuracy
00127     { return T(3.1415926535897932384626433832795028841971693993751L); }
00128     /**
00129      * A synonym for pi<real>().
00130      **********************************************************************/
00131     static inline real pi() throw() { return pi<real>(); }
00132     /**
00133      * <b>DEPRECATED</b> A synonym for pi<extened>().
00134      **********************************************************************/
00135     static inline extended epi() throw() { return pi<extended>(); }
00136 
00137     /**
00138      * @return the number of radians in a degree.
00139      **********************************************************************/
00140     template<typename T>
00141     static inline T degree() throw() { return pi<T>() / T(180); }
00142     /**
00143      * A synonym for degree<real>().
00144      **********************************************************************/
00145     static inline real degree() throw() { return degree<real>(); }
00146     /**
00147      * <b>DEPRECATED</b> A synonym for degree<extened>().
00148      **********************************************************************/
00149     static inline extended edegree() throw() { return degree<extended>(); }
00150 
00151 #if defined(DOXYGEN)
00152     /**
00153      * The hypotenuse function avoiding underflow and overflow.
00154      *
00155      * @param[in] x
00156      * @param[in] y
00157      * @return sqrt(\e x<sup>2</sup> + \e y<sup>2</sup>).
00158      **********************************************************************/
00159     template<typename T>
00160     static inline T hypot(T x, T y) throw() {
00161       x = std::abs(x);
00162       y = std::abs(y);
00163       T a = std::max(x, y),
00164         b = std::min(x, y) / a;
00165       return a * std::sqrt(1 + b * b);
00166     }
00167 #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH
00168     template<typename T>
00169     static inline T hypot(T x, T y) throw() { return std::hypot(x, y); }
00170 #elif defined(_MSC_VER)
00171     static inline double hypot(double x, double y) throw()
00172     { return _hypot(x, y); }
00173     static inline float hypot(float x, float y) throw()
00174     { return _hypotf(x, y); }
00175 #if !defined(__NO_LONG_DOUBLE_MATH)
00176     static inline long double hypot(long double x, long double y) throw()
00177     { return _hypot(x, y); }
00178 #endif
00179 #else
00180     // Use overloading to define generic versions
00181     static inline double hypot(double x, double y) throw()
00182     { return ::hypot(x, y); }
00183     static inline float hypot(float x, float y) throw()
00184     { return ::hypotf(x, y); }
00185 #if !defined(__NO_LONG_DOUBLE_MATH)
00186     static inline long double hypot(long double x, long double y) throw()
00187     { return ::hypotl(x, y); }
00188 #endif
00189 #endif
00190 
00191 #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MATH)
00192     /**
00193      * exp(\e x) - 1 accurate near \e x = 0.  This is taken from
00194      * N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd
00195      * Edition (SIAM, 2002), Sec 1.14.1, p 19.
00196      *
00197      * @param[in] x
00198      * @return exp(\e x) - 1.
00199      **********************************************************************/
00200     template<typename T>
00201     static inline T expm1(T x) throw() {
00202       volatile T
00203         y = std::exp(x),
00204         z = y - 1;
00205       // The reasoning here is similar to that for log1p.  The expression
00206       // mathematically reduces to exp(x) - 1, and the factor z/log(y) = (y -
00207       // 1)/log(y) is a slowly varying quantity near y = 1 and is accurately
00208       // computed.
00209       return std::abs(x) > 1 ? z : z == 0 ?  x : x * z / std::log(y);
00210     }
00211 #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH
00212     template<typename T>
00213     static inline T expm1(T x) throw() { return std::expm1(x); }
00214 #else
00215     static inline double expm1(double x) throw() { return ::expm1(x); }
00216     static inline float expm1(float x) throw() { return ::expm1f(x); }
00217 #if !defined(__NO_LONG_DOUBLE_MATH)
00218     static inline long double expm1(long double x) throw()
00219     { return ::expm1l(x); }
00220 #endif
00221 #endif
00222 
00223 #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MATH)
00224     /**
00225      * log(\e x + 1) accurate near \e x = 0.  This is taken See
00226      * D. Goldberg,
00227      * <a href="http://docs.sun.com/source/806-3568/ncg_goldberg.html"> What
00228      * every computer scientist should know about floating-point arithmetic</a>
00229      * (1991), Theorem 4.  See also, Higham (op. cit.), Answer to Problem 1.5,
00230      * p 528.
00231      *
00232      * @param[in] x
00233      * @return log(\e x + 1).
00234      **********************************************************************/
00235     template<typename T>
00236     static inline T log1p(T x) throw() {
00237       volatile T
00238         y = 1 + x,
00239         z = y - 1;
00240       // Here's the explanation for this magic: y = 1 + z, exactly, and z
00241       // approx x, thus log(y)/z (which is nearly constant near z = 0) returns
00242       // a good approximation to the true log(1 + x)/x.  The multiplication x *
00243       // (log(y)/z) introduces little additional error.
00244       return z == 0 ? x : x * std::log(y) / z;
00245     }
00246 #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH
00247     template<typename T>
00248     static inline T log1p(T x) throw() { return std::log1p(x); }
00249 #else
00250     static inline double log1p(double x) throw() { return ::log1p(x); }
00251     static inline float log1p(float x) throw() { return ::log1pf(x); }
00252 #if !defined(__NO_LONG_DOUBLE_MATH)
00253     static inline long double log1p(long double x) throw()
00254     { return ::log1pl(x); }
00255 #endif
00256 #endif
00257 
00258 #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MATH)
00259     /**
00260      * The inverse hyperbolic sine function.  This is defined in terms of
00261      * Math::log1p(\e x) in order to maintain accuracy near \e x = 0.  In
00262      * addition, the odd parity of the function is enforced.
00263      *
00264      * @param[in] x
00265      * @return asinh(\e x).
00266      **********************************************************************/
00267     template<typename T>
00268     static inline T asinh(T x) throw() {
00269       T y = std::abs(x);     // Enforce odd parity
00270       y = log1p(y * (1 + y/(hypot(T(1), y) + 1)));
00271       return x < 0 ? -y : y;
00272     }
00273 #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH
00274     template<typename T>
00275     static inline T asinh(T x) throw() { return std::asinh(x); }
00276 #else
00277     static inline double asinh(double x) throw() { return ::asinh(x); }
00278     static inline float asinh(float x) throw() { return ::asinhf(x); }
00279 #if !defined(__NO_LONG_DOUBLE_MATH)
00280     static inline long double asinh(long double x) throw()
00281     { return ::asinhl(x); }
00282 #endif
00283 #endif
00284 
00285 #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MATH)
00286     /**
00287      * The inverse hyperbolic tangent function.  This is defined in terms of
00288      * Math::log1p(\e x) in order to maintain accuracy near \e x = 0.  In
00289      * addition, the odd parity of the function is enforced.
00290      *
00291      * @param[in] x
00292      * @return atanh(\e x).
00293      **********************************************************************/
00294     template<typename T>
00295     static inline T atanh(T x) throw() {
00296       T y = std::abs(x);     // Enforce odd parity
00297       y = log1p(2 * y/(1 - y))/2;
00298       return x < 0 ? -y : y;
00299     }
00300 #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH
00301     template<typename T>
00302     static inline T atanh(T x) throw() { return std::atanh(x); }
00303 #else
00304     static inline double atanh(double x) throw() { return ::atanh(x); }
00305     static inline float atanh(float x) throw() { return ::atanhf(x); }
00306 #if !defined(__NO_LONG_DOUBLE_MATH)
00307     static inline long double atanh(long double x) throw()
00308     { return ::atanhl(x); }
00309 #endif
00310 #endif
00311 
00312 #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MATH)
00313     /**
00314      * The cube root function.
00315      *
00316      * @param[in] x
00317      * @return the real cube root of \e x.
00318      **********************************************************************/
00319     template<typename T>
00320     static inline T cbrt(T x) throw() {
00321       T y = std::pow(std::abs(x), 1/T(3)); // Return the real cube root
00322       return x < 0 ? -y : y;
00323     }
00324 #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH
00325     template<typename T>
00326     static inline T cbrt(T x) throw() { return std::cbrt(x); }
00327 #else
00328     static inline double cbrt(double x) throw() { return ::cbrt(x); }
00329     static inline float cbrt(float x) throw() { return ::cbrtf(x); }
00330 #if !defined(__NO_LONG_DOUBLE_MATH)
00331     static inline long double cbrt(long double x) throw() { return ::cbrtl(x); }
00332 #endif
00333 #endif
00334 
00335     /**
00336      * Test for finiteness.
00337      *
00338      * @param[in] x
00339      * @return true if number is finite, false if NaN or infinite.
00340      **********************************************************************/
00341     template<typename T>
00342     static inline bool isfinite(T x) throw() {
00343 #if defined(DOXYGEN)
00344       return std::abs(x) <= std::numeric_limits<T>::max();
00345 #elif (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MATH)
00346       return _finite(x) != 0;
00347 #else
00348       return std::isfinite(x);
00349 #endif
00350     }
00351 
00352     /**
00353      * The NaN (not a number)
00354      *
00355      * @return NaN if available, otherwise return the max real.
00356      **********************************************************************/
00357     template<typename T>
00358     static inline T NaN() throw() {
00359       return std::numeric_limits<T>::has_quiet_NaN ?
00360         std::numeric_limits<T>::quiet_NaN() :
00361         std::numeric_limits<T>::max();
00362     }
00363     /**
00364      * A synonym for NaN<real>().
00365      **********************************************************************/
00366     static inline real NaN() throw() { return NaN<real>(); }
00367 
00368     /**
00369      * Test for NaN.
00370      *
00371      * @param[in] x
00372      * @return true if argument is a NaN.
00373      **********************************************************************/
00374     template<typename T>
00375     static inline bool isnan(T x) throw() {
00376 #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MATH)
00377       return x != x;
00378 #else
00379       return std::isnan(x);
00380 #endif
00381     }
00382 
00383     /**
00384      * Infinity
00385      *
00386      * @return infinity if available, otherwise return the max real.
00387      **********************************************************************/
00388     template<typename T>
00389     static inline T infinity() throw() {
00390       return std::numeric_limits<T>::has_infinity ?
00391         std::numeric_limits<T>::infinity() :
00392         std::numeric_limits<T>::max();
00393     }
00394     /**
00395      * A synonym for infinity<real>().
00396      **********************************************************************/
00397     static inline real infinity() throw() { return infinity<real>(); }
00398   };
00399 
00400   /**
00401    * \brief %Constants needed by %GeographicLib
00402    *
00403    * Define constants specifying the WGS84 ellipsoid, the UTM and UPS
00404    * projections, and various unit conversions.
00405    **********************************************************************/
00406   class Constants {
00407   private:
00408     typedef Math::real real;
00409     Constants();                // Disable constructor
00410 
00411   public:
00412     /**
00413      * <b>DEPRECATED</b> A synonym for Math::pi<real>().
00414      **********************************************************************/
00415     static inline Math::real pi() throw() { return Math::pi<real>(); }
00416     /**
00417      * A synonym for Math::degree<real>().
00418      **********************************************************************/
00419     static inline Math::real degree() throw() { return Math::degree<real>(); }
00420     /**
00421      * @return the number of radians in an arcminute.
00422      **********************************************************************/
00423     static inline Math::real arcminute() throw()
00424     { return Math::degree<real>() / 60; }
00425     /**
00426      * @return the number of radians in an arcsecond.
00427      **********************************************************************/
00428     static inline Math::real arcsecond() throw()
00429     { return Math::degree<real>() / 3600; }
00430 
00431     /** \name Ellipsoid parameters
00432      **********************************************************************/
00433     ///@{
00434     /**
00435      * @return the equatorial radius of WGS84 ellipsoid
00436      **********************************************************************/
00437     template<typename T>
00438     static inline T WGS84_a() throw() { return T(6378137) * meter<T>(); }
00439     /**
00440      * A synonym for WGS84_a<real>().
00441      **********************************************************************/
00442     static inline Math::real WGS84_a() throw() { return WGS84_a<real>(); }
00443     /**
00444      * @return the reciprocal flattening of WGS84 ellipsoid
00445      **********************************************************************/
00446     template<typename T>
00447     static inline T WGS84_r() throw() {
00448       // 298.257223563 = 3393 * 87903691 / 1000000000
00449       return (T(3393) * T(87903691)) / T(1000000000);
00450     }
00451     /**
00452      * A synonym for WGS84_r<real>().
00453      **********************************************************************/
00454     static inline Math::real WGS84_r() throw() { return WGS84_r<real>(); }
00455     /**
00456      * @return the central scale factor for UTM
00457      **********************************************************************/
00458     template<typename T>
00459     static inline T UTM_k0() throw() {return T(9996) / T(10000); } // 0.9996
00460     /**
00461      * A synonym for UTM_k0<real>().
00462      **********************************************************************/
00463     static inline Math::real UTM_k0() throw() { return UTM_k0<real>(); }
00464     /**
00465      * @return the central scale factor for UPS
00466      **********************************************************************/
00467     template<typename T>
00468     static inline T UPS_k0() throw() { return T(994) / T(1000); } // 0.994
00469     /**
00470      * A synonym for UPS_k0<real>().
00471      **********************************************************************/
00472     static inline Math::real UPS_k0() throw() { return UPS_k0<real>(); }
00473     ///@}
00474 
00475     /** \name SI units
00476      **********************************************************************/
00477     ///@{
00478     /**
00479      * @return the number of meters in a meter.
00480      *
00481      * This is unity, but this lets the internal system of units be changed if
00482      * necessary.
00483      **********************************************************************/
00484     template<typename T>
00485     static inline T meter() throw() { return T(1); }
00486     /**
00487      * A synonym for meter<real>().
00488      **********************************************************************/
00489     static inline Math::real meter() throw() { return meter<real>(); }
00490     /**
00491      * @return the number of meters in a kilometer.
00492      **********************************************************************/
00493     static inline Math::real kilometer() throw()
00494     { return 1000 * meter<real>(); }
00495     /**
00496      * @return the number of meters in a nautical mile (approximately 1 arc
00497      *   minute)
00498      **********************************************************************/
00499     static inline Math::real nauticalmile() throw()
00500     { return 1852 * meter<real>(); }
00501     ///@}
00502 
00503     /** \name Anachronistic British units
00504      **********************************************************************/
00505     ///@{
00506     /**
00507      * @return the number of meters in an international foot.
00508      **********************************************************************/
00509     static inline Math::real foot() throw()
00510     { return real(0.0254L) * 12 * meter<real>(); }
00511     /**
00512      * @return the number of meters in a yard.
00513      **********************************************************************/
00514     static inline Math::real yard() throw() { return 3 * foot(); }
00515     /**
00516      * @return the number of meters in a fathom.
00517      **********************************************************************/
00518     static inline Math::real fathom() throw() { return 2 * yard(); }
00519     /**
00520      * @return the number of meters in a chain.
00521      **********************************************************************/
00522     static inline Math::real chain() throw() { return 22 * yard(); }
00523     /**
00524      * @return the number of meters in a furlong.
00525      **********************************************************************/
00526     static inline Math::real furlong() throw() { return 10 * chain(); }
00527     /**
00528      * @return the number of meters in a statute mile.
00529      **********************************************************************/
00530     static inline Math::real mile() throw() { return 8 * furlong(); }
00531     ///@}
00532 
00533     /** \name Anachronistic US units
00534      **********************************************************************/
00535     ///@{
00536     /**
00537      * @return the number of meters in a US survey foot.
00538      **********************************************************************/
00539     static inline Math::real surveyfoot() throw()
00540     { return real(1200) / real(3937) * meter<real>(); }
00541     ///@}
00542   };
00543 
00544   /**
00545    * \brief %Exception handling for %GeographicLib
00546    *
00547    * A class to handle exceptions.  It's derived off std::runtime_error so it
00548    * can be caught by the usual catch clauses.
00549    **********************************************************************/
00550   class GeographicErr : public std::runtime_error {
00551   public:
00552 
00553     /**
00554      * Constructor
00555      *
00556      * @param[in] msg a string message, which is accessible in the catch
00557      *   clause, via what().
00558      **********************************************************************/
00559     GeographicErr(const std::string& msg) : std::runtime_error(msg) {}
00560   };
00561 
00562 } // namespace GeographicLib
00563 
00564 #endif