| AlbersEqualArea.hpp | | AlbersEqualArea.hpp | |
| /** | | /** | |
| * \file AlbersEqualArea.hpp | | * \file AlbersEqualArea.hpp | |
| * \brief Header for GeographicLib::AlbersEqualArea class | | * \brief Header for GeographicLib::AlbersEqualArea class | |
| * | | * | |
|
| * Copyright (c) Charles Karney (2010) <charles@karney.com> and licensed un | | * Copyright (c) Charles Karney (2010, 2011) <charles@karney.com> and licen | |
| der | | sed | |
| * the LGPL. For more information, see http://geographiclib.sourceforge.ne | | * under the LGPL. For more information, see | |
| t/ | | * http://geographiclib.sourceforge.net/ | |
| **********************************************************************/ | | **********************************************************************/ | |
| | | | |
| #if !defined(GEOGRAPHICLIB_ALBERSEQUALAREA_HPP) | | #if !defined(GEOGRAPHICLIB_ALBERSEQUALAREA_HPP) | |
|
| #define GEOGRAPHICLIB_ALBERSEQUALAREA_HPP "$Id: AlbersEqualArea.hpp 6919 20
10-12-21 13:23:47Z karney $" | | #define GEOGRAPHICLIB_ALBERSEQUALAREA_HPP "$Id: 76631a35cadd11d3dcbfd2b2951
bd1528b8ae63b $" | |
| | | | |
|
| #include "GeographicLib/Constants.hpp" | | | |
| #include <algorithm> | | #include <algorithm> | |
|
| | | #include <GeographicLib/Constants.hpp> | |
| | | | |
| namespace GeographicLib { | | namespace GeographicLib { | |
| | | | |
| /** | | /** | |
| * \brief Albers Equal Area Conic Projection | | * \brief Albers Equal Area Conic Projection | |
| * | | * | |
| * Implementation taken from the report, | | * Implementation taken from the report, | |
| * - J. P. Snyder, | | * - J. P. Snyder, | |
| * <a href="http://pubs.er.usgs.gov/usgspubs/pp/pp1395"> Map Projection
s: A | | * <a href="http://pubs.er.usgs.gov/usgspubs/pp/pp1395"> Map Projection
s: A | |
| * Working Manual</a>, USGS Professional Paper 1395 (1987), | | * Working Manual</a>, USGS Professional Paper 1395 (1987), | |
| | | | |
| skipping to change at line 53 | | skipping to change at line 54 | |
| * AlbersEqualArea::Forward and AlbersEqualArea::Reverse functions. | | * AlbersEqualArea::Forward and AlbersEqualArea::Reverse functions. | |
| * AlbersEqualArea::Forward and AlbersEqualArea::Reverse also return the | | * AlbersEqualArea::Forward and AlbersEqualArea::Reverse also return the | |
| * meridian convergence, \e gamma, and azimuthal scale, \e k. A small sq
uare | | * meridian convergence, \e gamma, and azimuthal scale, \e k. A small sq
uare | |
| * aligned with the cardinal directions is projected to a rectangle with | | * aligned with the cardinal directions is projected to a rectangle with | |
| * dimensions \e k (in the E-W direction) and 1/\e k (in the N-S directio
n). | | * dimensions \e k (in the E-W direction) and 1/\e k (in the N-S directio
n). | |
| * The E-W sides of the rectangle are oriented \e gamma degrees | | * The E-W sides of the rectangle are oriented \e gamma degrees | |
| * counter-clockwise from the \e x axis. There is no provision in this c
lass | | * counter-clockwise from the \e x axis. There is no provision in this c
lass | |
| * for specifying a false easting or false northing or a different latitu
de | | * for specifying a false easting or false northing or a different latitu
de | |
| * of origin. | | * of origin. | |
| **********************************************************************/ | | **********************************************************************/ | |
|
| class AlbersEqualArea { | | class GEOGRAPHIC_EXPORT AlbersEqualArea { | |
| private: | | private: | |
| typedef Math::real real; | | typedef Math::real real; | |
| const real _a, _r, _f, _fm, _e2, _e, _e2m, _qZ, _qx; | | const real _a, _r, _f, _fm, _e2, _e, _e2m, _qZ, _qx; | |
| real _sign, _lat0, _k0; | | real _sign, _lat0, _k0; | |
| real _n0, _m02, _nrho0, _k2, _txi0, _scxi0, _sxi0; | | real _n0, _m02, _nrho0, _k2, _txi0, _scxi0, _sxi0; | |
|
| static const real eps, epsx, epsx2, tol, tol0, ahypover; | | static const real eps_; | |
| static const int numit = 5; // Newton iterations in Reverse | | static const real epsx_; | |
| static const int numit0 = 20; // Newton iterations in Init | | static const real epsx2_; | |
| static inline real sq(real x) throw() { return x * x; } | | static const real tol_; | |
| | | static const real tol0_; | |
| | | static const real ahypover_; | |
| | | static const int numit_ = 5; // Newton iterations in Reverse | |
| | | static const int numit0_ = 20; // Newton iterations in Init | |
| static inline real hyp(real x) throw() { return Math::hypot(real(1), x)
; } | | static inline real hyp(real x) throw() { return Math::hypot(real(1), x)
; } | |
| // atanh( e * x)/ e if f > 0 | | // atanh( e * x)/ e if f > 0 | |
| // atan (sqrt(-e2) * x)/sqrt(-e2) if f < 0 | | // atan (sqrt(-e2) * x)/sqrt(-e2) if f < 0 | |
| // x if f = 0 | | // x if f = 0 | |
| inline real atanhee(real x) const throw() { | | inline real atanhee(real x) const throw() { | |
| return _f > 0 ? Math::atanh(_e * x)/_e : | | return _f > 0 ? Math::atanh(_e * x)/_e : | |
| (_f < 0 ? std::atan(_e * x)/_e : x); | | (_f < 0 ? std::atan(_e * x)/_e : x); | |
| } | | } | |
| // return atanh(sqrt(x))/sqrt(x) - 1, accurate for small x | | // return atanh(sqrt(x))/sqrt(x) - 1, accurate for small x | |
| static real atanhxm1(real x) throw(); | | static real atanhxm1(real x) throw(); | |
| | | | |
| skipping to change at line 91 | | skipping to change at line 96 | |
| // | | // | |
| // General rules | | // General rules | |
| // h(x) = f(g(x)): Dh(x,y) = Df(g(x),g(y))*Dg(x,y) | | // h(x) = f(g(x)): Dh(x,y) = Df(g(x),g(y))*Dg(x,y) | |
| // h(x) = f(x)*g(x): | | // h(x) = f(x)*g(x): | |
| // Dh(x,y) = Df(x,y)*(g(x)+g(y))/2 + Dg(x,y)*(f(x)+f(y))/2 | | // Dh(x,y) = Df(x,y)*(g(x)+g(y))/2 + Dg(x,y)*(f(x)+f(y))/2 | |
| // | | // | |
| // sn(x) = x/sqrt(1+x^2): Dsn(x,y) = (x+y)/((sn(x)+sn(y))*(1+x^2)*(1+y^
2)) | | // sn(x) = x/sqrt(1+x^2): Dsn(x,y) = (x+y)/((sn(x)+sn(y))*(1+x^2)*(1+y^
2)) | |
| static inline real Dsn(real x, real y, real sx, real sy) throw() { | | static inline real Dsn(real x, real y, real sx, real sy) throw() { | |
| // sx = x/hyp(x) | | // sx = x/hyp(x) | |
| real t = x * y; | | real t = x * y; | |
|
| return t > 0 ? (x + y) * sq( (sx * sy)/t ) / (sx + sy) : | | return t > 0 ? (x + y) * Math::sq( (sx * sy)/t ) / (sx + sy) : | |
| (x - y != 0 ? (sx - sy) / (x - y) : 1); | | (x - y != 0 ? (sx - sy) / (x - y) : 1); | |
| } | | } | |
| // Datanhee(x,y) = atanhee((x-y)/(1-e^2*x*y))/(x-y) | | // Datanhee(x,y) = atanhee((x-y)/(1-e^2*x*y))/(x-y) | |
| inline real Datanhee(real x, real y) const throw() { | | inline real Datanhee(real x, real y) const throw() { | |
| real t = x - y, d = 1 - _e2 * x * y; | | real t = x - y, d = 1 - _e2 * x * y; | |
| return t != 0 ? atanhee(t / d) / t : 1 / d; | | return t != 0 ? atanhee(t / d) / t : 1 / d; | |
| } | | } | |
| // DDatanhee(x,y) = (Datanhee(1,y) - Datanhee(1,x))/(y-x) | | // DDatanhee(x,y) = (Datanhee(1,y) - Datanhee(1,x))/(y-x) | |
| real DDatanhee(real x, real y) const throw(); | | real DDatanhee(real x, real y) const throw(); | |
| void Init(real sphi1, real cphi1, real sphi2, real cphi2, real k1) thro
w(); | | void Init(real sphi1, real cphi1, real sphi2, real cphi2, real k1) thro
w(); | |
| | | | |
| skipping to change at line 292 | | skipping to change at line 297 | |
| /** | | /** | |
| * A global instantiation of AlbersEqualArea with the WGS84 ellipsoid,
\e | | * A global instantiation of AlbersEqualArea with the WGS84 ellipsoid,
\e | |
| * stdlat = -90<sup>o</sup>, and \e k0 = 1. This degenerates to the | | * stdlat = -90<sup>o</sup>, and \e k0 = 1. This degenerates to the | |
| * Lambert azimuthal equal area projection. | | * Lambert azimuthal equal area projection. | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static const AlbersEqualArea AzimuthalEqualAreaSouth; | | static const AlbersEqualArea AzimuthalEqualAreaSouth; | |
| }; | | }; | |
| | | | |
| } // namespace GeographicLib | | } // namespace GeographicLib | |
| | | | |
|
| #endif | | #endif // GEOGRAPHICLIB_ALBERSEQUALAREA_HPP | |
| | | | |
End of changes. 8 change blocks. |
| 12 lines changed or deleted | | 16 lines changed or added | |
|
| Constants.hpp | | Constants.hpp | |
| /** | | /** | |
| * \file Constants.hpp | | * \file Constants.hpp | |
| * \brief Header for GeographicLib::Constants class | | * \brief Header for GeographicLib::Constants class | |
| * | | * | |
| * Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co
m> | | * Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co
m> | |
| * and licensed under the LGPL. For more information, see | | * and licensed under the LGPL. For more information, see | |
| * http://geographiclib.sourceforge.net/ | | * http://geographiclib.sourceforge.net/ | |
| **********************************************************************/ | | **********************************************************************/ | |
| | | | |
| #if !defined(GEOGRAPHICLIB_CONSTANTS_HPP) | | #if !defined(GEOGRAPHICLIB_CONSTANTS_HPP) | |
|
| #define GEOGRAPHICLIB_CONSTANTS_HPP "$Id: Constants.hpp 6967 2011-02-19 15: | | #define GEOGRAPHICLIB_CONSTANTS_HPP "$Id: 484b008f7c73d568472abe2cc8fcbdfbc | |
| 53:41Z karney $" | | 81a63af $" | |
| | | | |
| | | #include <GeographicLib/Config.h> | |
| | | | |
| /** | | /** | |
| * Are C++0X math functions available? | | * Are C++0X math functions available? | |
| **********************************************************************/ | | **********************************************************************/ | |
| #if !defined(GEOGRAPHICLIB_CPLUSPLUS0X_MATH) | | #if !defined(GEOGRAPHICLIB_CPLUSPLUS0X_MATH) | |
|
| #if defined(__GXX_EXPERIMENTAL_CXX0X__) | | # if defined(__GXX_EXPERIMENTAL_CXX0X__) | |
| #define GEOGRAPHICLIB_CPLUSPLUS0X_MATH 1 | | # define GEOGRAPHICLIB_CPLUSPLUS0X_MATH 1 | |
| #else | | # else | |
| #define GEOGRAPHICLIB_CPLUSPLUS0X_MATH 0 | | # define GEOGRAPHICLIB_CPLUSPLUS0X_MATH 0 | |
| #endif | | # endif | |
| #endif | | #endif | |
| | | | |
| /** | | /** | |
| * A compile-time assert. Use C++0X static_assert, if available. | | * A compile-time assert. Use C++0X static_assert, if available. | |
| **********************************************************************/ | | **********************************************************************/ | |
| #if !defined(STATIC_ASSERT) | | #if !defined(STATIC_ASSERT) | |
|
| #if defined(__GXX_EXPERIMENTAL_CXX0X__) | | # if defined(__GXX_EXPERIMENTAL_CXX0X__) | |
| #define STATIC_ASSERT static_assert | | # define STATIC_ASSERT static_assert | |
| #elif defined(_MSC_VER) && _MSC_VER >= 1600 | | # elif defined(_MSC_VER) && _MSC_VER >= 1600 | |
| #define STATIC_ASSERT static_assert | | # define STATIC_ASSERT static_assert | |
| #else | | # else | |
| #define STATIC_ASSERT(cond,reason) { enum{ STATIC_ASSERT_ENUM=1/int(cond) } | | # define STATIC_ASSERT(cond,reason) \ | |
| ; } | | { enum{ STATIC_ASSERT_ENUM = 1/int(cond) }; } | |
| #endif | | # endif | |
| #endif | | #endif | |
| | | | |
| #if defined(__GNUC__) | | #if defined(__GNUC__) | |
| // Suppress "defined but not used" warnings | | // Suppress "defined but not used" warnings | |
|
| #define RCSID_DECL(x) namespace \ | | # define RCSID_DECL(x) namespace \ | |
| { char VAR_ ## x [] __attribute__((used)) = x; } | | { char VAR_ ## x [] __attribute__((used)) = x; } | |
| #else | | #else | |
| /** | | /** | |
| * Insertion of RCS Id strings into the object file. | | * Insertion of RCS Id strings into the object file. | |
| **********************************************************************/ | | **********************************************************************/ | |
|
| #define RCSID_DECL(x) namespace { char VAR_ ## x [] = x; } | | # define RCSID_DECL(x) namespace { char VAR_ ## x [] = x; } | |
| | | #endif | |
| | | | |
| | | #if defined(_WIN32) && defined(GEOGRAPHIC_SHARED_LIB) | |
| | | # if defined(Geographic_EXPORTS) | |
| | | # define GEOGRAPHIC_EXPORT __declspec(dllexport) | |
| | | # else | |
| | | # define GEOGRAPHIC_EXPORT __declspec(dllimport) | |
| | | # endif | |
| | | #else | |
| | | # define GEOGRAPHIC_EXPORT | |
| #endif | | #endif | |
| | | | |
| RCSID_DECL(GEOGRAPHICLIB_CONSTANTS_HPP) | | RCSID_DECL(GEOGRAPHICLIB_CONSTANTS_HPP) | |
| | | | |
| #if !defined(GEOGRAPHICLIB_PREC) | | #if !defined(GEOGRAPHICLIB_PREC) | |
| /** | | /** | |
| * The precision of floating point numbers used in %GeographicLib. 0 means | | * The precision of floating point numbers used in %GeographicLib. 0 means | |
| * float; 1 (default) means double; 2 means long double. Nearly all the | | * float; 1 (default) means double; 2 means long double. Nearly all the | |
| * testing has been carried out with doubles and that's the recommended | | * testing has been carried out with doubles and that's the recommended | |
| * configuration. Note that with Microsoft Visual Studio, long double is t
he | | * configuration. Note that with Microsoft Visual Studio, long double is t
he | |
| * same as double. | | * same as double. | |
| **********************************************************************/ | | **********************************************************************/ | |
| #define GEOGRAPHICLIB_PREC 1 | | #define GEOGRAPHICLIB_PREC 1 | |
| #endif | | #endif | |
| | | | |
|
| #if defined(__CYGWIN__) && defined(__GNUC__) && __GNUC__ < 4 | | | |
| // g++ 3.x under cygwin doesn't have long double | | | |
| #define __NO_LONG_DOUBLE_MATH 1 | | | |
| #endif | | | |
| | | | |
| #include <cmath> | | #include <cmath> | |
| #include <limits> | | #include <limits> | |
| #include <algorithm> | | #include <algorithm> | |
| #include <stdexcept> | | #include <stdexcept> | |
| | | | |
| /** | | /** | |
| * \brief Namespace for %GeographicLib | | * \brief Namespace for %GeographicLib | |
| * | | * | |
| * All of %GeographicLib is defined within the GeographicLib namespace. In | | * All of %GeographicLib is defined within the GeographicLib namespace. In | |
| * addtion all the header files are included via %GeographicLib/filename.
This | | * addtion all the header files are included via %GeographicLib/filename.
This | |
| | | | |
| skipping to change at line 87 | | skipping to change at line 95 | |
| **********************************************************************/ | | **********************************************************************/ | |
| namespace GeographicLib { | | namespace GeographicLib { | |
| | | | |
| /** | | /** | |
| * \brief Mathematical functions needed by %GeographicLib | | * \brief Mathematical functions needed by %GeographicLib | |
| * | | * | |
| * Define mathematical functions in order to localize system dependencies
and | | * Define mathematical functions in order to localize system dependencies
and | |
| * to provide generic versions of the functions. In addition define a re
al | | * to provide generic versions of the functions. In addition define a re
al | |
| * type to be used by %GeographicLib. | | * type to be used by %GeographicLib. | |
| **********************************************************************/ | | **********************************************************************/ | |
|
| class Math { | | class GEOGRAPHIC_EXPORT Math { | |
| private: | | private: | |
| void dummy() { | | void dummy() { | |
| STATIC_ASSERT((GEOGRAPHICLIB_PREC) >= 0 && (GEOGRAPHICLIB_PREC) <= 2, | | STATIC_ASSERT((GEOGRAPHICLIB_PREC) >= 0 && (GEOGRAPHICLIB_PREC) <= 2, | |
| "Bad value of precision"); | | "Bad value of precision"); | |
| } | | } | |
| Math(); // Disable constructor | | Math(); // Disable constructor | |
| public: | | public: | |
| | | | |
|
| #if !defined(__NO_LONG_DOUBLE_MATH) | | #if defined(HAVE_LONG_DOUBLE) | |
| /** | | /** | |
| * The extended precision type for real numbers, used for some testing. | | * The extended precision type for real numbers, used for some testing. | |
| * This is long double on computers with this type; otherwise it is dou
ble. | | * This is long double on computers with this type; otherwise it is dou
ble. | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| typedef long double extended; | | typedef long double extended; | |
| #else | | #else | |
| typedef double extended; | | typedef double extended; | |
| #endif | | #endif | |
| | | | |
| #if GEOGRAPHICLIB_PREC == 1 | | #if GEOGRAPHICLIB_PREC == 1 | |
| | | | |
| skipping to change at line 125 | | skipping to change at line 133 | |
| #elif GEOGRAPHICLIB_PREC == 2 | | #elif GEOGRAPHICLIB_PREC == 2 | |
| typedef extended real; | | typedef extended real; | |
| #else | | #else | |
| typedef double real; | | typedef double real; | |
| #endif | | #endif | |
| | | | |
| /** | | /** | |
| * @return \e pi | | * @return \e pi | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| template<typename T> | | template<typename T> | |
|
| static inline T pi() throw() | | static inline T pi() throw() { return std::atan2(T(0), -T(1)); } | |
| // good for about 168-bit accuracy | | | |
| { return T(3.1415926535897932384626433832795028841971693993751L); } | | | |
| /** | | /** | |
| * A synonym for pi<real>(). | | * A synonym for pi<real>(). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static inline real pi() throw() { return pi<real>(); } | | static inline real pi() throw() { return pi<real>(); } | |
| /** | | /** | |
| * <b>DEPRECATED</b> A synonym for pi<extened>(). | | * <b>DEPRECATED</b> A synonym for pi<extened>(). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static inline extended epi() throw() { return pi<extended>(); } | | static inline extended epi() throw() { return pi<extended>(); } | |
| | | | |
| /** | | /** | |
| | | | |
| skipping to change at line 151 | | skipping to change at line 157 | |
| static inline T degree() throw() { return pi<T>() / T(180); } | | static inline T degree() throw() { return pi<T>() / T(180); } | |
| /** | | /** | |
| * A synonym for degree<real>(). | | * A synonym for degree<real>(). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static inline real degree() throw() { return degree<real>(); } | | static inline real degree() throw() { return degree<real>(); } | |
| /** | | /** | |
| * <b>DEPRECATED</b> A synonym for degree<extened>(). | | * <b>DEPRECATED</b> A synonym for degree<extened>(). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static inline extended edegree() throw() { return degree<extended>(); } | | static inline extended edegree() throw() { return degree<extended>(); } | |
| | | | |
|
| | | /** | |
| | | * Square a number. | |
| | | * @param[in] x | |
| | | * @return \e x<sup>2</sup>. | |
| | | ********************************************************************** | |
| | | / | |
| | | template<typename T> | |
| | | static inline T sq(T x) throw() { return x * x; } | |
| | | | |
| #if defined(DOXYGEN) | | #if defined(DOXYGEN) | |
| /** | | /** | |
| * The hypotenuse function avoiding underflow and overflow. | | * The hypotenuse function avoiding underflow and overflow. | |
| * | | * | |
| * @param[in] x | | * @param[in] x | |
| * @param[in] y | | * @param[in] y | |
| * @return sqrt(\e x<sup>2</sup> + \e y<sup>2</sup>). | | * @return sqrt(\e x<sup>2</sup> + \e y<sup>2</sup>). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| template<typename T> | | template<typename T> | |
| static inline T hypot(T x, T y) throw() { | | static inline T hypot(T x, T y) throw() { | |
| x = std::abs(x); | | x = std::abs(x); | |
| y = std::abs(y); | | y = std::abs(y); | |
|
| T a = std::max(x, y), | | T a = (std::max)(x, y), | |
| b = std::min(x, y) / a; | | b = (std::min)(x, y) / (a ? a : 1); | |
| return a * std::sqrt(1 + b * b); | | return a * std::sqrt(1 + b * b); | |
| } | | } | |
| #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH | | #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH | |
| template<typename T> | | template<typename T> | |
| static inline T hypot(T x, T y) throw() { return std::hypot(x, y); } | | static inline T hypot(T x, T y) throw() { return std::hypot(x, y); } | |
| #elif defined(_MSC_VER) | | #elif defined(_MSC_VER) | |
| static inline double hypot(double x, double y) throw() | | static inline double hypot(double x, double y) throw() | |
| { return _hypot(x, y); } | | { return _hypot(x, y); } | |
| static inline float hypot(float x, float y) throw() | | static inline float hypot(float x, float y) throw() | |
| { return _hypotf(x, y); } | | { return _hypotf(x, y); } | |
|
| #if !defined(__NO_LONG_DOUBLE_MATH) | | #if defined(HAVE_LONG_DOUBLE) | |
| static inline long double hypot(long double x, long double y) throw() | | static inline long double hypot(long double x, long double y) throw() | |
| { return _hypot(x, y); } | | { return _hypot(x, y); } | |
| #endif | | #endif | |
| #else | | #else | |
| // Use overloading to define generic versions | | // Use overloading to define generic versions | |
| static inline double hypot(double x, double y) throw() | | static inline double hypot(double x, double y) throw() | |
| { return ::hypot(x, y); } | | { return ::hypot(x, y); } | |
| static inline float hypot(float x, float y) throw() | | static inline float hypot(float x, float y) throw() | |
| { return ::hypotf(x, y); } | | { return ::hypotf(x, y); } | |
|
| #if !defined(__NO_LONG_DOUBLE_MATH) | | #if defined(HAVE_LONG_DOUBLE) | |
| static inline long double hypot(long double x, long double y) throw() | | static inline long double hypot(long double x, long double y) throw() | |
| { return ::hypotl(x, y); } | | { return ::hypotl(x, y); } | |
| #endif | | #endif | |
| #endif | | #endif | |
| | | | |
| #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MA
TH) | | #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MA
TH) | |
| /** | | /** | |
| * exp(\e x) - 1 accurate near \e x = 0. This is taken from | | * exp(\e x) - 1 accurate near \e x = 0. This is taken from | |
| * N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd | | * N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd | |
| * Edition (SIAM, 2002), Sec 1.14.1, p 19. | | * Edition (SIAM, 2002), Sec 1.14.1, p 19. | |
| | | | |
| skipping to change at line 217 | | skipping to change at line 231 | |
| // 1)/log(y) is a slowly varying quantity near y = 1 and is accuratel
y | | // 1)/log(y) is a slowly varying quantity near y = 1 and is accuratel
y | |
| // computed. | | // computed. | |
| return std::abs(x) > 1 ? z : z == 0 ? x : x * z / std::log(y); | | return std::abs(x) > 1 ? z : z == 0 ? x : x * z / std::log(y); | |
| } | | } | |
| #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH | | #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH | |
| template<typename T> | | template<typename T> | |
| static inline T expm1(T x) throw() { return std::expm1(x); } | | static inline T expm1(T x) throw() { return std::expm1(x); } | |
| #else | | #else | |
| static inline double expm1(double x) throw() { return ::expm1(x); } | | static inline double expm1(double x) throw() { return ::expm1(x); } | |
| static inline float expm1(float x) throw() { return ::expm1f(x); } | | static inline float expm1(float x) throw() { return ::expm1f(x); } | |
|
| #if !defined(__NO_LONG_DOUBLE_MATH) | | #if defined(HAVE_LONG_DOUBLE) | |
| static inline long double expm1(long double x) throw() | | static inline long double expm1(long double x) throw() | |
| { return ::expm1l(x); } | | { return ::expm1l(x); } | |
| #endif | | #endif | |
| #endif | | #endif | |
| | | | |
| #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MA
TH) | | #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MA
TH) | |
| /** | | /** | |
|
| * log(\e x + 1) accurate near \e x = 0. This is taken See | | * log(\e x + 1) accurate near \e x = 0. | |
| * D. Goldberg, | | * | |
| * <a href="http://docs.sun.com/source/806-3568/ncg_goldberg.html"> Wha | | * This is taken from D. Goldberg, | |
| t | | * <a href="http://dx.doi.org/10.1145/103162.103163">What every compute | |
| * every computer scientist should know about floating-point arithmetic | | r | |
| </a> | | * scientist should know about floating-point arithmetic</a> (1991), | |
| * (1991), Theorem 4. See also, Higham (op. cit.), Answer to Problem 1 | | * Theorem 4. See also, Higham (op. cit.), Answer to Problem 1.5, p 52 | |
| .5, | | 8. | |
| * p 528. | | | |
| * | | * | |
| * @param[in] x | | * @param[in] x | |
| * @return log(\e x + 1). | | * @return log(\e x + 1). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| template<typename T> | | template<typename T> | |
| static inline T log1p(T x) throw() { | | static inline T log1p(T x) throw() { | |
| volatile T | | volatile T | |
| y = 1 + x, | | y = 1 + x, | |
| z = y - 1; | | z = y - 1; | |
| // Here's the explanation for this magic: y = 1 + z, exactly, and z | | // Here's the explanation for this magic: y = 1 + z, exactly, and z | |
| | | | |
| skipping to change at line 252 | | skipping to change at line 266 | |
| // a good approximation to the true log(1 + x)/x. The multiplication
x * | | // a good approximation to the true log(1 + x)/x. The multiplication
x * | |
| // (log(y)/z) introduces little additional error. | | // (log(y)/z) introduces little additional error. | |
| return z == 0 ? x : x * std::log(y) / z; | | return z == 0 ? x : x * std::log(y) / z; | |
| } | | } | |
| #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH | | #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH | |
| template<typename T> | | template<typename T> | |
| static inline T log1p(T x) throw() { return std::log1p(x); } | | static inline T log1p(T x) throw() { return std::log1p(x); } | |
| #else | | #else | |
| static inline double log1p(double x) throw() { return ::log1p(x); } | | static inline double log1p(double x) throw() { return ::log1p(x); } | |
| static inline float log1p(float x) throw() { return ::log1pf(x); } | | static inline float log1p(float x) throw() { return ::log1pf(x); } | |
|
| #if !defined(__NO_LONG_DOUBLE_MATH) | | #if defined(HAVE_LONG_DOUBLE) | |
| static inline long double log1p(long double x) throw() | | static inline long double log1p(long double x) throw() | |
| { return ::log1pl(x); } | | { return ::log1pl(x); } | |
| #endif | | #endif | |
| #endif | | #endif | |
| | | | |
| #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MA
TH) | | #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MA
TH) | |
| /** | | /** | |
| * The inverse hyperbolic sine function. This is defined in terms of | | * The inverse hyperbolic sine function. This is defined in terms of | |
| * Math::log1p(\e x) in order to maintain accuracy near \e x = 0. In | | * Math::log1p(\e x) in order to maintain accuracy near \e x = 0. In | |
| * addition, the odd parity of the function is enforced. | | * addition, the odd parity of the function is enforced. | |
| | | | |
| skipping to change at line 279 | | skipping to change at line 293 | |
| T y = std::abs(x); // Enforce odd parity | | T y = std::abs(x); // Enforce odd parity | |
| y = log1p(y * (1 + y/(hypot(T(1), y) + 1))); | | y = log1p(y * (1 + y/(hypot(T(1), y) + 1))); | |
| return x < 0 ? -y : y; | | return x < 0 ? -y : y; | |
| } | | } | |
| #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH | | #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH | |
| template<typename T> | | template<typename T> | |
| static inline T asinh(T x) throw() { return std::asinh(x); } | | static inline T asinh(T x) throw() { return std::asinh(x); } | |
| #else | | #else | |
| static inline double asinh(double x) throw() { return ::asinh(x); } | | static inline double asinh(double x) throw() { return ::asinh(x); } | |
| static inline float asinh(float x) throw() { return ::asinhf(x); } | | static inline float asinh(float x) throw() { return ::asinhf(x); } | |
|
| #if !defined(__NO_LONG_DOUBLE_MATH) | | #if defined(HAVE_LONG_DOUBLE) | |
| static inline long double asinh(long double x) throw() | | static inline long double asinh(long double x) throw() | |
| { return ::asinhl(x); } | | { return ::asinhl(x); } | |
| #endif | | #endif | |
| #endif | | #endif | |
| | | | |
| #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MA
TH) | | #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MA
TH) | |
| /** | | /** | |
| * The inverse hyperbolic tangent function. This is defined in terms o
f | | * The inverse hyperbolic tangent function. This is defined in terms o
f | |
| * Math::log1p(\e x) in order to maintain accuracy near \e x = 0. In | | * Math::log1p(\e x) in order to maintain accuracy near \e x = 0. In | |
| * addition, the odd parity of the function is enforced. | | * addition, the odd parity of the function is enforced. | |
| | | | |
| skipping to change at line 306 | | skipping to change at line 320 | |
| T y = std::abs(x); // Enforce odd parity | | T y = std::abs(x); // Enforce odd parity | |
| y = log1p(2 * y/(1 - y))/2; | | y = log1p(2 * y/(1 - y))/2; | |
| return x < 0 ? -y : y; | | return x < 0 ? -y : y; | |
| } | | } | |
| #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH | | #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH | |
| template<typename T> | | template<typename T> | |
| static inline T atanh(T x) throw() { return std::atanh(x); } | | static inline T atanh(T x) throw() { return std::atanh(x); } | |
| #else | | #else | |
| static inline double atanh(double x) throw() { return ::atanh(x); } | | static inline double atanh(double x) throw() { return ::atanh(x); } | |
| static inline float atanh(float x) throw() { return ::atanhf(x); } | | static inline float atanh(float x) throw() { return ::atanhf(x); } | |
|
| #if !defined(__NO_LONG_DOUBLE_MATH) | | #if defined(HAVE_LONG_DOUBLE) | |
| static inline long double atanh(long double x) throw() | | static inline long double atanh(long double x) throw() | |
| { return ::atanhl(x); } | | { return ::atanhl(x); } | |
| #endif | | #endif | |
| #endif | | #endif | |
| | | | |
| #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MA
TH) | | #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MA
TH) | |
| /** | | /** | |
| * The cube root function. | | * The cube root function. | |
| * | | * | |
| * @param[in] x | | * @param[in] x | |
| | | | |
| skipping to change at line 330 | | skipping to change at line 344 | |
| static inline T cbrt(T x) throw() { | | static inline T cbrt(T x) throw() { | |
| T y = std::pow(std::abs(x), 1/T(3)); // Return the real cube root | | T y = std::pow(std::abs(x), 1/T(3)); // Return the real cube root | |
| return x < 0 ? -y : y; | | return x < 0 ? -y : y; | |
| } | | } | |
| #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH | | #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH | |
| template<typename T> | | template<typename T> | |
| static inline T cbrt(T x) throw() { return std::cbrt(x); } | | static inline T cbrt(T x) throw() { return std::cbrt(x); } | |
| #else | | #else | |
| static inline double cbrt(double x) throw() { return ::cbrt(x); } | | static inline double cbrt(double x) throw() { return ::cbrt(x); } | |
| static inline float cbrt(float x) throw() { return ::cbrtf(x); } | | static inline float cbrt(float x) throw() { return ::cbrtf(x); } | |
|
| #if !defined(__NO_LONG_DOUBLE_MATH) | | #if defined(HAVE_LONG_DOUBLE) | |
| static inline long double cbrt(long double x) throw() { return ::cbrtl(
x); } | | static inline long double cbrt(long double x) throw() { return ::cbrtl(
x); } | |
| #endif | | #endif | |
| #endif | | #endif | |
| | | | |
| /** | | /** | |
| * Test for finiteness. | | * Test for finiteness. | |
| * | | * | |
| * @param[in] x | | * @param[in] x | |
| * @return true if number is finite, false if NaN or infinite. | | * @return true if number is finite, false if NaN or infinite. | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| template<typename T> | | template<typename T> | |
| static inline bool isfinite(T x) throw() { | | static inline bool isfinite(T x) throw() { | |
| #if defined(DOXYGEN) | | #if defined(DOXYGEN) | |
|
| return std::abs(x) <= std::numeric_limits<T>::max(); | | return std::abs(x) <= (std::numeric_limits<T>::max)(); | |
| #elif (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MATH) | | #elif (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MATH) | |
| return _finite(x) != 0; | | return _finite(x) != 0; | |
| #else | | #else | |
| return std::isfinite(x); | | return std::isfinite(x); | |
| #endif | | #endif | |
| } | | } | |
| | | | |
| /** | | /** | |
| * The NaN (not a number) | | * The NaN (not a number) | |
| * | | * | |
| * @return NaN if available, otherwise return the max real. | | * @return NaN if available, otherwise return the max real. | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| template<typename T> | | template<typename T> | |
| static inline T NaN() throw() { | | static inline T NaN() throw() { | |
| return std::numeric_limits<T>::has_quiet_NaN ? | | return std::numeric_limits<T>::has_quiet_NaN ? | |
| std::numeric_limits<T>::quiet_NaN() : | | std::numeric_limits<T>::quiet_NaN() : | |
|
| std::numeric_limits<T>::max(); | | (std::numeric_limits<T>::max)(); | |
| } | | } | |
| /** | | /** | |
| * A synonym for NaN<real>(). | | * A synonym for NaN<real>(). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static inline real NaN() throw() { return NaN<real>(); } | | static inline real NaN() throw() { return NaN<real>(); } | |
| | | | |
| /** | | /** | |
| * Test for NaN. | | * Test for NaN. | |
| * | | * | |
| * @param[in] x | | * @param[in] x | |
| | | | |
| skipping to change at line 392 | | skipping to change at line 406 | |
| | | | |
| /** | | /** | |
| * Infinity | | * Infinity | |
| * | | * | |
| * @return infinity if available, otherwise return the max real. | | * @return infinity if available, otherwise return the max real. | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| template<typename T> | | template<typename T> | |
| static inline T infinity() throw() { | | static inline T infinity() throw() { | |
| return std::numeric_limits<T>::has_infinity ? | | return std::numeric_limits<T>::has_infinity ? | |
| std::numeric_limits<T>::infinity() : | | std::numeric_limits<T>::infinity() : | |
|
| std::numeric_limits<T>::max(); | | (std::numeric_limits<T>::max)(); | |
| } | | } | |
| /** | | /** | |
| * A synonym for infinity<real>(). | | * A synonym for infinity<real>(). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static inline real infinity() throw() { return infinity<real>(); } | | static inline real infinity() throw() { return infinity<real>(); } | |
| }; | | }; | |
| | | | |
| /** | | /** | |
| * \brief %Constants needed by %GeographicLib | | * \brief %Constants needed by %GeographicLib | |
| * | | * | |
| * Define constants specifying the WGS84 ellipsoid, the UTM and UPS | | * Define constants specifying the WGS84 ellipsoid, the UTM and UPS | |
| * projections, and various unit conversions. | | * projections, and various unit conversions. | |
| **********************************************************************/ | | **********************************************************************/ | |
|
| class Constants { | | class GEOGRAPHIC_EXPORT Constants { | |
| private: | | private: | |
| typedef Math::real real; | | typedef Math::real real; | |
| Constants(); // Disable constructor | | Constants(); // Disable constructor | |
| | | | |
| public: | | public: | |
| /** | | /** | |
| * <b>DEPRECATED</b> A synonym for Math::pi<real>(). | | * <b>DEPRECATED</b> A synonym for Math::pi<real>(). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static inline Math::real pi() throw() { return Math::pi<real>(); } | | static inline Math::real pi() throw() { return Math::pi<real>(); } | |
| /** | | /** | |
| | | | |
| skipping to change at line 448 | | skipping to change at line 462 | |
| static inline T WGS84_a() throw() { return T(6378137) * meter<T>(); } | | static inline T WGS84_a() throw() { return T(6378137) * meter<T>(); } | |
| /** | | /** | |
| * A synonym for WGS84_a<real>(). | | * A synonym for WGS84_a<real>(). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static inline Math::real WGS84_a() throw() { return WGS84_a<real>(); } | | static inline Math::real WGS84_a() throw() { return WGS84_a<real>(); } | |
| /** | | /** | |
| * @return the reciprocal flattening of WGS84 ellipsoid | | * @return the reciprocal flattening of WGS84 ellipsoid | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| template<typename T> | | template<typename T> | |
| static inline T WGS84_r() throw() { | | static inline T WGS84_r() throw() { | |
|
| // 298.257223563 = 3393 * 87903691 / 1000000000 | | // 298.257223563 | |
| return (T(3393) * T(87903691)) / T(1000000000); | | return T(298) + T(257223563) / T(1000000000); | |
| } | | } | |
| /** | | /** | |
| * A synonym for WGS84_r<real>(). | | * A synonym for WGS84_r<real>(). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static inline Math::real WGS84_r() throw() { return WGS84_r<real>(); } | | static inline Math::real WGS84_r() throw() { return WGS84_r<real>(); } | |
| /** | | /** | |
| * @return the central scale factor for UTM | | * @return the central scale factor for UTM | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| template<typename T> | | template<typename T> | |
| static inline T UTM_k0() throw() {return T(9996) / T(10000); } // 0.999
6 | | static inline T UTM_k0() throw() {return T(9996) / T(10000); } // 0.999
6 | |
| | | | |
| skipping to change at line 545 | | skipping to change at line 559 | |
| ///@{ | | ///@{ | |
| /** | | /** | |
| * @return the number of meters in a US survey foot. | | * @return the number of meters in a US survey foot. | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static inline Math::real surveyfoot() throw() | | static inline Math::real surveyfoot() throw() | |
| { return real(1200) / real(3937) * meter<real>(); } | | { return real(1200) / real(3937) * meter<real>(); } | |
| ///@} | | ///@} | |
| }; | | }; | |
| | | | |
| /** | | /** | |
|
| * \brief %Exception handling for %GeographicLib | | * \brief An accumulator for sums. | |
| | | * | |
| | | * This allow many numbers of floating point type \e T to be added togeth | |
| | | er | |
| | | * with twice the normal precision. Thus if \e T is double, the effectiv | |
| | | e | |
| | | * precision of the sum is 106 bits or about 32 decimal places. The core | |
| | | * idea is the error free transformation of a sum, D. E. Knuth, TAOCP, Vo | |
| | | l 2, | |
| | | * 4.2.2, Theorem B. | |
| | | * | |
| | | * Two implementations are provided. The "fast" one is an implementation | |
| | | of | |
| | | * Algorithm 4.1, of T. Ogita, S. M. Rump, S. Oishi, | |
| | | * <a href="http://dx.doi.org/10.1137/030601818"> Accurate sum and dot | |
| | | * product</a>, SIAM J. Sci. Comp., 26(6) 1955-1988 (2005). The accumula | |
| | | tor | |
| | | * is represented by a two numbers _s + _t where _s is the result of the | |
| | | * normal sum and _t accumulates (approximately) the exact roundoff error | |
| | | s in | |
| | | * the summation of _s. | |
| | | * | |
| | | * The slow "non-fast" implementation follows J. R. Shewchuk, | |
| | | * <a href="http://dx.doi.org/10.1007/PL00009321"> Adaptive Precision | |
| | | * Floating-Point Arithmetic and Fast Robust Geometric Predicates</a>, | |
| | | * Discrete & Computational Geometry 18(3) 305-363 (1997). In this case, | |
| | | * with each addition of a number to the accumulator, _s and _t are adjus | |
| | | ted | |
| | | * so that _s represents the sum to an accuracy of 1 ulp. This may resul | |
| | | t in | |
| | | * considerably greater accuracy than the fast implementation (depending | |
| | | on | |
| | | * the input data). | |
| | | * | |
| | | * Approximate timings (summing vector<double> per addition) | |
| | | * - double: 2ns = 1 | |
| | | * - Accumulator<double, true>: 7ns = 3.5 | |
| | | * - Accumulator<double, true>: 21ns = 10 -- Optimize() on each addition | |
| | | * - Accumulator<double, false>: 23ns = 11 | |
| | | * . | |
| | | * Thus the slow method is about 3 times slower than the fast method. | |
| | | * However all times are negligible with the typical timings for geodesic | |
| | | * calculations which are roughly 2us. Thus the default mode is taken to | |
| | | be | |
| | | * slow, \e fast = false. | |
| | | * | |
| | | * In the documentation of the member functions, \e sum stands for the va | |
| | | lue | |
| | | * currently held in the accumulator. | |
| | | **********************************************************************/ | |
| | | template<typename T = Math::real, bool fast = false> | |
| | | class Accumulator { | |
| | | private: | |
| | | // _s accumulates for the straight sum | |
| | | // _t accumulates the errors. | |
| | | T _s, _t; | |
| | | // Error free transformation of a sum. Note that t can be the same as | |
| | | one | |
| | | // of the first two arguments. | |
| | | static inline T sum(T u, T v, T& t) { | |
| | | volatile T s = u + v; | |
| | | volatile T up = s - v; | |
| | | volatile T vpp = s - up; | |
| | | up -= u; | |
| | | vpp -= v; | |
| | | t = -(up + vpp); | |
| | | // u + v = s + t | |
| | | // = round(u + v) + t | |
| | | return s; | |
| | | } | |
| | | // Same as sum, but requires abs(u) >= abs(v). This isn't currently us | |
| | | ed. | |
| | | static inline T fastsum(T u, T v, T& t) { | |
| | | volatile T s = u + v; | |
| | | volatile T vp = s - u; | |
| | | t = v - vp; | |
| | | return s; | |
| | | } | |
| | | void Add(T y) throw() { | |
| | | if (fast) { | |
| | | _s = sum(_s, y, y); // Accumulate sum to _s with error in y; | |
| | | _t += y; // accumulate y in _t | |
| | | } else { // Here's Shewchuk's solution... | |
| | | T u; // hold exact sum as [s, t, u] | |
| | | y = sum(y, _t, u); // Accumulate starting at least significant | |
| | | end | |
| | | _s = sum(y, _s, _t); | |
| | | // Start is _s, _t decreasing and non-adjacent. Sum is now (s + t | |
| | | + u) | |
| | | // exactly with s, t, u non-adjacent and in decreasing order (excep | |
| | | t | |
| | | // for possible zeros). The following code tries to normalize the | |
| | | // result. Ideally, we want _s = round(s+t+u) and _u = round(s+t+u | |
| | | - | |
| | | // _s). The follow does an approximate job (and maintains the | |
| | | // decreasing non-adjacent property). Here are two "failures" usin | |
| | | g | |
| | | // 3-bit floats: | |
| | | // | |
| | | // Case 1: _s is not equal to round(s+t+u) -- off by 1 ulp | |
| | | // [12, -1] - 8 -> [4, 0, -1] -> [4, -1] = 3 should be [3, 0] = 3 | |
| | | // | |
| | | // Case 2: _s+_t is not as close to s+t+u as it shold be | |
| | | // [64, 5] + 4 -> [64, 8, 1] -> [64, 8] = 72 (off by 1) | |
| | | // should be [80, -7] = 73 (exact) | |
| | | // | |
| | | // "Fixing" these problems is probably not worth the expense. The | |
| | | // representation inevitably leads to small errors in the accumulat | |
| | | ed | |
| | | // values. The additional errors illustrated here amount to 1 ulp | |
| | | of | |
| | | // the less significant word during each addition to the Accumulato | |
| | | r | |
| | | // and an additional possible error of 1 ulp in the reported sum. | |
| | | // | |
| | | // Incidentally, the "ideal" representation described above is not | |
| | | // canonical, because _s = round(_s + _t) may not be true. For | |
| | | // example, with 3-bit floats: | |
| | | // | |
| | | // [128, 16] + 1 -> [160, -16] -- 160 = round(145). | |
| | | // But [160, 0] - 16 -> [128, 16] -- 128 = round(144). | |
| | | // | |
| | | if (_s == 0) // This implies t == 0, | |
| | | _s = u; // so result is u | |
| | | else | |
| | | _t += u; // otherwise just accumulate u to t. | |
| | | } | |
| | | } | |
| | | // Perhaps these should both be _s + _t? | |
| | | T Sum() const throw() { return fast ? _s + _t : _s; } | |
| | | T Sum(T y) const throw() { | |
| | | Accumulator a(*this); | |
| | | a.Add(y); | |
| | | return a.Sum(); | |
| | | } | |
| | | public: | |
| | | /** | |
| | | * Construct from a \e T. This is not declared explicit, so that you c | |
| | | an | |
| | | * write <code>Accumulator<double> a = 5;</code>. | |
| | | * | |
| | | * @param[in] y set \e sum = \e y. | |
| | | ********************************************************************** | |
| | | / | |
| | | Accumulator(T y = T(0)) throw() : _s(y), _t(0) { | |
| | | STATIC_ASSERT(!std::numeric_limits<T>::is_integer, | |
| | | "Accumulator type is not floating point"); | |
| | | }; | |
| | | /** | |
| | | * Set the accumulator to a number. | |
| | | * | |
| | | * @param[in] y set \e sum = \e y. | |
| | | ********************************************************************** | |
| | | / | |
| | | Accumulator& operator=(T y) throw() { _s = y; _t = 0; return *this; } | |
| | | /** | |
| | | * Return the value held in the accumulator. | |
| | | * | |
| | | * @return \e sum. | |
| | | ********************************************************************** | |
| | | / | |
| | | T operator()() const throw() { return Sum(); } | |
| | | /** | |
| | | * Return the result of adding a number to \e sum (but don't change \e | |
| | | sum). | |
| | | * | |
| | | * @param[in] y the number to be added to the sum. | |
| | | * @return \e sum + \e y. | |
| | | ********************************************************************** | |
| | | / | |
| | | T operator()(T y) const throw() { return Sum(y); } | |
| | | /** | |
| | | * Add a number to the accumulator. | |
| | | * | |
| | | * @param[in] y set \e sum += \e y. | |
| | | ********************************************************************** | |
| | | / | |
| | | Accumulator& operator+=(T y) throw() { Add(y); return *this; } | |
| | | /** | |
| | | * Subtract a number from the accumulator. | |
| | | * | |
| | | * @param[in] y set \e sum -= \e y. | |
| | | ********************************************************************** | |
| | | / | |
| | | Accumulator& operator-=(T y) throw() { Add(-y); return *this; } | |
| | | /** | |
| | | * Multiply accumulator by an integer. To avoid loss of accuracy, use | |
| | | only | |
| | | * integers such that \e n * \e T is exactly representable as a \e T (i | |
| | | .e., | |
| | | * +/- powers of two). Use \e n = -1 to negate \e sum. | |
| | | * | |
| | | * @param[in] n set \e sum *= \e n. | |
| | | ********************************************************************** | |
| | | / | |
| | | Accumulator& operator*=(int n) throw() { _s *= n; _t *= n; return *this | |
| | | ; } | |
| | | /** | |
| | | * Optimize how \e sum is stored in the accumulator. This is a no-op f | |
| | | or | |
| | | * the non-fast implementation. It is rarely necessary to do this; how | |
| | | ever | |
| | | * the accuracy might be improved if this function is called every time | |
| | | a | |
| | | * million numbers (for example) have been added to the accumulator. | |
| | | ********************************************************************** | |
| | | / | |
| | | void Optimize() throw() { if (fast) _s = sum(_s, _t, _t); } | |
| | | /** | |
| | | * Test equality of an Accumulator with a number. | |
| | | ********************************************************************** | |
| | | / | |
| | | bool operator==(T y) const throw() { return Sum() == y; } | |
| | | /** | |
| | | * Test inequality of an Accumulator with a number. | |
| | | ********************************************************************** | |
| | | / | |
| | | bool operator!=(T y) const throw() { return Sum() != y; } | |
| | | /** | |
| | | * Less operator on an Accumulator and a number. | |
| | | ********************************************************************** | |
| | | / | |
| | | bool operator<(T y) const throw() { return Sum() < y; } | |
| | | /** | |
| | | * Less or equal operator on an Accumulator and a number. | |
| | | ********************************************************************** | |
| | | / | |
| | | bool operator<=(T y) const throw() { return Sum() <= y; } | |
| | | /** | |
| | | * Greater operator on an Accumulator and a number. | |
| | | ********************************************************************** | |
| | | / | |
| | | bool operator>(T y) const throw() { return Sum() > y; } | |
| | | /** | |
| | | * Greater or equal operator on an Accumulator and a number. | |
| | | ********************************************************************** | |
| | | / | |
| | | bool operator>=(T y) const throw() { return Sum() >= y; } | |
| | | }; | |
| | | | |
| | | /** | |
| | | * \brief Exception handling for %GeographicLib | |
| * | | * | |
|
| * A class to handle exceptions. It's derived off std::runtime_error so
it | | * A class to handle exceptions. It's derived from std::runtime_error so
it | |
| * can be caught by the usual catch clauses. | | * can be caught by the usual catch clauses. | |
| **********************************************************************/ | | **********************************************************************/ | |
| class GeographicErr : public std::runtime_error { | | class GeographicErr : public std::runtime_error { | |
| public: | | public: | |
| | | | |
| /** | | /** | |
| * Constructor | | * Constructor | |
| * | | * | |
| * @param[in] msg a string message, which is accessible in the catch | | * @param[in] msg a string message, which is accessible in the catch | |
| * clause, via what(). | | * clause, via what(). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| GeographicErr(const std::string& msg) : std::runtime_error(msg) {} | | GeographicErr(const std::string& msg) : std::runtime_error(msg) {} | |
| }; | | }; | |
| | | | |
| } // namespace GeographicLib | | } // namespace GeographicLib | |
| | | | |
|
| #endif | | #endif // GEOGRAPHICLIB_CONSTANTS_HPP | |
| | | | |
End of changes. 27 change blocks. |
| 53 lines changed or deleted | | 307 lines changed or added | |
|
| DMS.hpp | | DMS.hpp | |
| /** | | /** | |
| * \file DMS.hpp | | * \file DMS.hpp | |
| * \brief Header for GeographicLib::DMS class | | * \brief Header for GeographicLib::DMS class | |
| * | | * | |
| * Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co
m> | | * Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co
m> | |
| * and licensed under the LGPL. For more information, see | | * and licensed under the LGPL. For more information, see | |
| * http://geographiclib.sourceforge.net/ | | * http://geographiclib.sourceforge.net/ | |
| **********************************************************************/ | | **********************************************************************/ | |
| | | | |
| #if !defined(GEOGRAPHICLIB_DMS_HPP) | | #if !defined(GEOGRAPHICLIB_DMS_HPP) | |
|
| #define GEOGRAPHICLIB_DMS_HPP "$Id: DMS.hpp 6954 2011-02-16 14:36:56Z karne
y $" | | #define GEOGRAPHICLIB_DMS_HPP "$Id: 782c1f225888b3c0b4ff19d0b677e549bd36fa1
c $" | |
| | | | |
|
| #include "GeographicLib/Constants.hpp" | | | |
| #include <sstream> | | #include <sstream> | |
| #include <iomanip> | | #include <iomanip> | |
|
| | | #include <GeographicLib/Constants.hpp> | |
| | | | |
| | | #if defined(_MSC_VER) | |
| | | // Squelch warnings about dll vs string | |
| | | #pragma warning (push) | |
| | | #pragma warning (disable: 4251) | |
| | | #endif | |
| | | | |
| namespace GeographicLib { | | namespace GeographicLib { | |
| | | | |
| /** | | /** | |
| * \brief Convert between degrees and %DMS representation. | | * \brief Convert between degrees and %DMS representation. | |
| * | | * | |
| * Parse a string representing degree, minutes, and seconds and return th
e | | * Parse a string representing degree, minutes, and seconds and return th
e | |
| * angle in degrees and format an angle in degrees as degree, minutes, an
d | | * angle in degrees and format an angle in degrees as degree, minutes, an
d | |
|
| * seconds. | | * seconds. In addition, handle NANs and infinities on input and output. | |
| **********************************************************************/ | | **********************************************************************/ | |
|
| class DMS { | | class GEOGRAPHIC_EXPORT DMS { | |
| private: | | private: | |
| typedef Math::real real; | | typedef Math::real real; | |
| static int lookup(const std::string& s, char c) throw() { | | static int lookup(const std::string& s, char c) throw() { | |
| std::string::size_type r = s.find(toupper(c)); | | std::string::size_type r = s.find(toupper(c)); | |
| return r == std::string::npos ? -1 : int(r); | | return r == std::string::npos ? -1 : int(r); | |
| } | | } | |
| template<typename T> static std::string str(T x) { | | template<typename T> static std::string str(T x) { | |
| std::ostringstream s; s << x; return s.str(); | | std::ostringstream s; s << x; return s.str(); | |
| } | | } | |
|
| static const std::string hemispheres; | | static const std::string hemispheres_; | |
| static const std::string signs; | | static const std::string signs_; | |
| static const std::string digits; | | static const std::string digits_; | |
| static const std::string dmsindicators; | | static const std::string dmsindicators_; | |
| static const std::string components[3]; | | static const std::string components_[3]; | |
| static Math::real NumMatch(const std::string& s); | | static Math::real NumMatch(const std::string& s); | |
| DMS(); // Disable constructor | | DMS(); // Disable constructor | |
| | | | |
| public: | | public: | |
| | | | |
| /** | | /** | |
| * Indicator for presence of hemisphere indicator (N/S/E/W) on latitude
s | | * Indicator for presence of hemisphere indicator (N/S/E/W) on latitude
s | |
| * and longitudes. | | * and longitudes. | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| enum flag { | | enum flag { | |
| | | | |
| skipping to change at line 109 | | skipping to change at line 115 | |
| | | | |
| /** | | /** | |
| * Convert a string in DMS to an angle. | | * Convert a string in DMS to an angle. | |
| * | | * | |
| * @param[in] dms string input. | | * @param[in] dms string input. | |
| * @param[out] ind a DMS::flag value indicating the presence of a | | * @param[out] ind a DMS::flag value indicating the presence of a | |
| * hemisphere indicator. | | * hemisphere indicator. | |
| * @return angle (degrees). | | * @return angle (degrees). | |
| * | | * | |
| * Degrees, minutes, and seconds are indicated by the letters d, ', &qu
ot;, | | * Degrees, minutes, and seconds are indicated by the letters d, ', &qu
ot;, | |
|
| * and these components may only be given in this order. Any (but not | | * and these components_ may only be given in this order. Any (but not | |
| all) | | all) | |
| * components may be omitted. The last component indicator may be omit | | * components_ may be omitted. The last component indicator may be omi | |
| ted | | tted | |
| * and is assumed to be tbe next smallest unit (thus 33d10 is interpret
ed | | * and is assumed to be tbe next smallest unit (thus 33d10 is interpret
ed | |
| * as 33d10'). The final component may be a decimal fraction but the | | * as 33d10'). The final component may be a decimal fraction but the | |
|
| * non-final components must be integers. The integer parts of the min | | * non-final components_ must be integers. The integer parts of the mi | |
| utes | | nutes | |
| * and seconds components must be less than 60. A single leading sign | | * and seconds components_ must be less than 60. A single leading sign | |
| is | | is | |
| * permitted. A hemisphere designator (N, E, W, S) may be added to tbe | | * permitted. A hemisphere designator (N, E, W, S) may be added to tbe | |
| * beginning or end of the string. The result is multiplied by the imp
lied | | * beginning or end of the string. The result is multiplied by the imp
lied | |
| * signed of the hemisphere designator (negative for S and W). In addi
tion | | * signed of the hemisphere designator (negative for S and W). In addi
tion | |
| * \e ind is set to DMS::LATITUDE if N or S is present, to DMS::LONGITU
DE | | * \e ind is set to DMS::LATITUDE if N or S is present, to DMS::LONGITU
DE | |
| * if E or W is present, and to DMS::NONE otherwise. Throws an error o
n a | | * if E or W is present, and to DMS::NONE otherwise. Throws an error o
n a | |
| * malformed string. No check is performed on the range of the result. | | * malformed string. No check is performed on the range of the result. | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static Math::real Decode(const std::string& dms, flag& ind); | | static Math::real Decode(const std::string& dms, flag& ind); | |
| | | | |
| /** | | /** | |
| * Convert DMS to an angle. | | * Convert DMS to an angle. | |
| * | | * | |
| * @param[in] d degrees. | | * @param[in] d degrees. | |
| * @param[in] m arc minutes. | | * @param[in] m arc minutes. | |
| * @param[in] s arc seconds. | | * @param[in] s arc seconds. | |
| * @return angle (degrees) | | * @return angle (degrees) | |
| * | | * | |
|
| * This does not propagate the sign on \e d to the other components, so | | * This does not propagate the sign on \e d to the other components_, s
o | |
| * -3d20' would need to be represented as - DMS::Decode(3.0, 20.0) or | | * -3d20' would need to be represented as - DMS::Decode(3.0, 20.0) or | |
| * DMS::Decode(-3.0, -20.0). | | * DMS::Decode(-3.0, -20.0). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static Math::real Decode(real d, real m = 0, real s = 0) throw() | | static Math::real Decode(real d, real m = 0, real s = 0) throw() | |
| { return d + (m + s/real(60))/real(60); } | | { return d + (m + s/real(60))/real(60); } | |
| | | | |
| /** | | /** | |
| * Convert a string to a real number. | | * Convert a string to a real number. | |
| * | | * | |
| * @param[in] str string input. | | * @param[in] str string input. | |
| * @return decoded number. | | * @return decoded number. | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static Math::real Decode(const std::string& str); | | static Math::real Decode(const std::string& str); | |
| | | | |
| /** | | /** | |
|
| * Convert a pait of strings to latitude and longitude. | | * Convert a pair of strings to latitude and longitude. | |
| * | | * | |
| * @param[in] dmsa first string. | | * @param[in] dmsa first string. | |
| * @param[in] dmsb second string. | | * @param[in] dmsb second string. | |
| * @param[out] lat latitude. | | * @param[out] lat latitude. | |
| * @param[out] lon longitude. | | * @param[out] lon longitude. | |
| * | | * | |
| * By default, the \e lat (resp., \e lon) is assigned to the results of | | * By default, the \e lat (resp., \e lon) is assigned to the results of | |
| * decoding \e dmsa (resp., \e dmsb). However this is overridden if ei
ther | | * decoding \e dmsa (resp., \e dmsb). However this is overridden if ei
ther | |
| * \e dmsa or \e dmsb contain a latitude or longitude hemisphere design
ator | | * \e dmsa or \e dmsb contain a latitude or longitude hemisphere design
ator | |
| * (N, S, E, W). Throws an error if the decoded numbers are out of the | | * (N, S, E, W). Throws an error if the decoded numbers are out of the | |
| | | | |
| skipping to change at line 197 | | skipping to change at line 203 | |
| * the range [-180<sup>o</sup>, 180<sup>o</sup>). | | * the range [-180<sup>o</sup>, 180<sup>o</sup>). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static Math::real DecodeAzimuth(const std::string& azistr); | | static Math::real DecodeAzimuth(const std::string& azistr); | |
| | | | |
| /** | | /** | |
| * Convert angle (in degrees) into a DMS string. | | * Convert angle (in degrees) into a DMS string. | |
| * | | * | |
| * @param[in] angle input angle (degrees) | | * @param[in] angle input angle (degrees) | |
| * @param[in] trailing DMS::component value indicating the trailing uni
ts | | * @param[in] trailing DMS::component value indicating the trailing uni
ts | |
| * on the string and this is given as a decimal number if necessary. | | * on the string and this is given as a decimal number if necessary. | |
|
| * @param[in] prec the number of digits after the decimal point for the | | * @param[in] prec the number of digits_ after the decimal point for th
e | |
| * trailing component. | | * trailing component. | |
| * @param[in] ind DMS::flag value indicated additional formatting. | | * @param[in] ind DMS::flag value indicated additional formatting. | |
| * @return formatted string | | * @return formatted string | |
| * | | * | |
| * The interpretation of \e ind is as follows: | | * The interpretation of \e ind is as follows: | |
| * - ind == DMS::NONE, signed result no leading zeros on degrees except
in | | * - ind == DMS::NONE, signed result no leading zeros on degrees except
in | |
| * the units place, e.g., -8d03'. | | * the units place, e.g., -8d03'. | |
| * - ind == DMS::LATITUDE, trailing N or S hemisphere designator, no si
gn, | | * - ind == DMS::LATITUDE, trailing N or S hemisphere designator, no si
gn, | |
|
| * pad degrees to 2 digits, e.g., 08d03'S. | | * pad degrees to 2 digits_, e.g., 08d03'S. | |
| * - ind == DMS::LONGITUDE, trailing E or W hemisphere designator, no | | * - ind == DMS::LONGITUDE, trailing E or W hemisphere designator, no | |
|
| * sign, pad degrees to 3 digits, e.g., 008d03'W. | | * sign, pad degrees to 3 digits_, e.g., 008d03'W. | |
| * - ind == DMS::AZIMUTH, convert to the range [0, 360<sup>o</sup>), no | | * - ind == DMS::AZIMUTH, convert to the range [0, 360<sup>o</sup>), no | |
|
| * sign, pad degrees to 3 digits, , e.g., 351d57'. | | * sign, pad degrees to 3 digits_, , e.g., 351d57'. | |
| * . | | * . | |
|
| * The integer parts of the minutes and seconds components are always g | | * The integer parts of the minutes and seconds components_ are always | |
| iven | | given | |
| * with 2 digits. | | * with 2 digits_. | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static std::string Encode(real angle, component trailing, unsigned prec
, | | static std::string Encode(real angle, component trailing, unsigned prec
, | |
| flag ind = NONE); | | flag ind = NONE); | |
| | | | |
| /** | | /** | |
| * Convert angle into a DMS string selecting the trailing component | | * Convert angle into a DMS string selecting the trailing component | |
| * based on the precision. | | * based on the precision. | |
| * | | * | |
| * @param[in] angle input angle (degrees) | | * @param[in] angle input angle (degrees) | |
| * @param[in] prec the precision relative to 1 degree. | | * @param[in] prec the precision relative to 1 degree. | |
| | | | |
| skipping to change at line 276 | | skipping to change at line 282 | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static void Encode(real ang, real& d, real& m, real& s) throw() { | | static void Encode(real ang, real& d, real& m, real& s) throw() { | |
| d = int(ang); ang = 60 * (ang - d); | | d = int(ang); ang = 60 * (ang - d); | |
| m = int(ang); s = 60 * (ang - m); | | m = int(ang); s = 60 * (ang - m); | |
| } | | } | |
| | | | |
| }; | | }; | |
| | | | |
| } // namespace GeographicLib | | } // namespace GeographicLib | |
| | | | |
|
| | | #if defined(_MSC_VER) | |
| | | #pragma warning (pop) | |
| #endif | | #endif | |
|
| | | | |
| | | #endif // GEOGRAPHICLIB_DMS_HPP | |
| | | | |
End of changes. 17 change blocks. |
| 26 lines changed or deleted | | 34 lines changed or added | |
|
| Geocentric.hpp | | Geocentric.hpp | |
| /** | | /** | |
| * \file Geocentric.hpp | | * \file Geocentric.hpp | |
| * \brief Header for GeographicLib::Geocentric class | | * \brief Header for GeographicLib::Geocentric class | |
| * | | * | |
| * Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co
m> | | * Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co
m> | |
| * and licensed under the LGPL. For more information, see | | * and licensed under the LGPL. For more information, see | |
| * http://geographiclib.sourceforge.net/ | | * http://geographiclib.sourceforge.net/ | |
| **********************************************************************/ | | **********************************************************************/ | |
| | | | |
| #if !defined(GEOGRAPHICLIB_GEOCENTRIC_HPP) | | #if !defined(GEOGRAPHICLIB_GEOCENTRIC_HPP) | |
|
| #define GEOGRAPHICLIB_GEOCENTRIC_HPP "$Id: Geocentric.hpp 6967 2011-02-19 1
5:53:41Z karney $" | | #define GEOGRAPHICLIB_GEOCENTRIC_HPP "$Id: e3512898a482b7eccb92d2ca020e63ae
f5614446 $" | |
| | | | |
|
| #include "GeographicLib/Constants.hpp" | | | |
| #include <vector> | | #include <vector> | |
| #include <algorithm> | | #include <algorithm> | |
|
| | | #include <GeographicLib/Constants.hpp> | |
| | | | |
| namespace GeographicLib { | | namespace GeographicLib { | |
| | | | |
| /** | | /** | |
| * \brief %Geocentric coordinates | | * \brief %Geocentric coordinates | |
| * | | * | |
| * Convert between geodetic coordinates latitude = \e lat, longitude = \e | | * Convert between geodetic coordinates latitude = \e lat, longitude = \e | |
| * lon, height = \e h (measured vertically from the surface of the ellips
oid) | | * lon, height = \e h (measured vertically from the surface of the ellips
oid) | |
| * to geocentric coordinates (\e x, \e y, \e z). The origin of geocentri
c | | * to geocentric coordinates (\e x, \e y, \e z). The origin of geocentri
c | |
| * coordinates is at the center of the earth. The \e z axis goes thru th
e | | * coordinates is at the center of the earth. The \e z axis goes thru th
e | |
| | | | |
| skipping to change at line 41 | | skipping to change at line 41 | |
| * straightforward. For the reverse transformation we use | | * straightforward. For the reverse transformation we use | |
| * - H. Vermeille, | | * - H. Vermeille, | |
| * <a href="http://dx.doi.org/10.1007/s00190-002-0273-6"> Direct | | * <a href="http://dx.doi.org/10.1007/s00190-002-0273-6"> Direct | |
| * transformation from geocentric coordinates to geodetic coordinates</
a>, | | * transformation from geocentric coordinates to geodetic coordinates</
a>, | |
| * J. Geodesy 76, 451–454 (2002). | | * J. Geodesy 76, 451–454 (2002). | |
| * . | | * . | |
| * Several changes have been made to ensure that the method returns accur
ate | | * Several changes have been made to ensure that the method returns accur
ate | |
| * results for all finite inputs (even if \e h is infinite). The changes
are | | * results for all finite inputs (even if \e h is infinite). The changes
are | |
| * described in Appendix B of | | * described in Appendix B of | |
| * - C. F. F. Karney, | | * - C. F. F. Karney, | |
|
| * <a href="http://arxiv.org/abs/1102.1215">Geodesics | | * <a href="http://arxiv.org/abs/1102.1215v1">Geodesics | |
| * on an ellipsoid of revolution</a>, | | * on an ellipsoid of revolution</a>, | |
| * Feb. 2011; | | * Feb. 2011; | |
| * preprint | | * preprint | |
|
| * <a href="http://arxiv.org/abs/1102.1215">arxiv:1102.1215</a>. | | * <a href="http://arxiv.org/abs/1102.1215v1">arxiv:1102.1215v1</a>. | |
| * . | | * . | |
| * See \ref geocentric for more information. | | * See \ref geocentric for more information. | |
| * | | * | |
| * The errors in these routines are close to round-off. Specifically, fo
r | | * The errors in these routines are close to round-off. Specifically, fo
r | |
| * points within 5000 km of the surface of the ellipsoid (either inside o
r | | * points within 5000 km of the surface of the ellipsoid (either inside o
r | |
| * outside the ellipsoid), the error is bounded by 7 nm for the WGS84 | | * outside the ellipsoid), the error is bounded by 7 nm for the WGS84 | |
| * ellipsoid. See \ref geocentric for further information on the errors. | | * ellipsoid. See \ref geocentric for further information on the errors. | |
| **********************************************************************/ | | **********************************************************************/ | |
| | | | |
|
| class Geocentric { | | class GEOGRAPHIC_EXPORT Geocentric { | |
| private: | | private: | |
| typedef Math::real real; | | typedef Math::real real; | |
| friend class LocalCartesian; | | friend class LocalCartesian; | |
|
| static const size_t dim = 3, dim2 = dim * dim; | | static const size_t dim_ = 3; | |
| | | static const size_t dim2_ = dim_ * dim_; | |
| const real _a, _r, _f, _e2, _e2m, _e2a, _e4a, _maxrad; | | const real _a, _r, _f, _e2, _e2m, _e2a, _e4a, _maxrad; | |
|
| static inline real sq(real x) throw() { return x * x; } | | | |
| // Actually this can be static because it doesn't depend on the ellipso
id. | | // Actually this can be static because it doesn't depend on the ellipso
id. | |
| // But let's be more general than that. | | // But let's be more general than that. | |
| void Rotation(real sphi, real cphi, real slam, real clam, | | void Rotation(real sphi, real cphi, real slam, real clam, | |
|
| real M[dim2]) const throw(); | | real M[dim2_]) const throw(); | |
| void IntForward(real lat, real lon, real h, real& x, real& y, real& z, | | void IntForward(real lat, real lon, real h, real& x, real& y, real& z, | |
|
| real M[dim2]) const throw(); | | real M[dim2_]) const throw(); | |
| void IntReverse(real x, real y, real z, real& lat, real& lon, real& h, | | void IntReverse(real x, real y, real z, real& lat, real& lon, real& h, | |
|
| real M[dim2]) const throw(); | | real M[dim2_]) const throw(); | |
| public: | | public: | |
| | | | |
| /** | | /** | |
| * Constructor for a ellipsoid with | | * Constructor for a ellipsoid with | |
| * | | * | |
| * @param[in] a equatorial radius (meters) | | * @param[in] a equatorial radius (meters) | |
| * @param[in] r reciprocal flattening. Setting \e r = 0 implies \e r =
inf | | * @param[in] r reciprocal flattening. Setting \e r = 0 implies \e r =
inf | |
| * or flattening = 0 (i.e., a sphere). Negative \e r indicates a pro
late | | * or flattening = 0 (i.e., a sphere). Negative \e r indicates a pro
late | |
| * ellipsoid. | | * ellipsoid. | |
| * | | * | |
| | | | |
| skipping to change at line 122 | | skipping to change at line 122 | |
| * @param[out] z geocentric coordinate (meters). | | * @param[out] z geocentric coordinate (meters). | |
| * @param[out] M if the length of the vector is 9, fill with the rotati
on | | * @param[out] M if the length of the vector is 9, fill with the rotati
on | |
| * matrix in row-major order. | | * matrix in row-major order. | |
| * | | * | |
| * Pre-multiplying a unit vector in local cartesian coordinates (east, | | * Pre-multiplying a unit vector in local cartesian coordinates (east, | |
| * north, up) by \e M transforms the vector to geocentric coordinates. | | * north, up) by \e M transforms the vector to geocentric coordinates. | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| void Forward(real lat, real lon, real h, real& x, real& y, real& z, | | void Forward(real lat, real lon, real h, real& x, real& y, real& z, | |
| std::vector<real>& M) | | std::vector<real>& M) | |
| const throw() { | | const throw() { | |
|
| real t[dim2]; | | real t[dim2_]; | |
| IntForward(lat, lon, h, x, y, z, t); | | IntForward(lat, lon, h, x, y, z, t); | |
|
| if (M.end() == M.begin() + dim2) | | if (M.end() == M.begin() + dim2_) | |
| copy(t, t + dim2, M.begin()); | | copy(t, t + dim2_, M.begin()); | |
| } | | } | |
| | | | |
| /** | | /** | |
| * Convert from geocentric to geodetic to coordinates. | | * Convert from geocentric to geodetic to coordinates. | |
| * | | * | |
| * @param[in] x geocentric coordinate (meters). | | * @param[in] x geocentric coordinate (meters). | |
| * @param[in] y geocentric coordinate (meters). | | * @param[in] y geocentric coordinate (meters). | |
| * @param[in] z geocentric coordinate (meters). | | * @param[in] z geocentric coordinate (meters). | |
| * @param[out] lat latitude of point (degrees). | | * @param[out] lat latitude of point (degrees). | |
| * @param[out] lon longitude of point (degrees). | | * @param[out] lon longitude of point (degrees). | |
| | | | |
| skipping to change at line 171 | | skipping to change at line 171 | |
| * @param[out] M if the length of the vector is 9, fill with the rotati
on | | * @param[out] M if the length of the vector is 9, fill with the rotati
on | |
| * matrix in row-major order. | | * matrix in row-major order. | |
| * | | * | |
| * Pre-multiplying a unit vector in geocentric coordinates by the trans
pose | | * Pre-multiplying a unit vector in geocentric coordinates by the trans
pose | |
| * of \e M transforms the vector to local cartesian coordinates (east, | | * of \e M transforms the vector to local cartesian coordinates (east, | |
| * north, up). | | * north, up). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| void Reverse(real x, real y, real z, real& lat, real& lon, real& h, | | void Reverse(real x, real y, real z, real& lat, real& lon, real& h, | |
| std::vector<real>& M) | | std::vector<real>& M) | |
| const throw() { | | const throw() { | |
|
| real t[dim2]; | | real t[dim2_]; | |
| IntReverse(x, y, z, lat, lon, h, t); | | IntReverse(x, y, z, lat, lon, h, t); | |
|
| if (M.end() == M.begin() + dim2) | | if (M.end() == M.begin() + dim2_) | |
| copy(t, t + dim2, M.begin()); | | copy(t, t + dim2_, M.begin()); | |
| } | | } | |
| | | | |
| /** \name Inspector functions | | /** \name Inspector functions | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| ///@{ | | ///@{ | |
| /** | | /** | |
| * @return \e a the equatorial radius of the ellipsoid (meters). This
is | | * @return \e a the equatorial radius of the ellipsoid (meters). This
is | |
| * the value used in the constructor. | | * the value used in the constructor. | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| Math::real MajorRadius() const throw() { return _a; } | | Math::real MajorRadius() const throw() { return _a; } | |
| | | | |
| skipping to change at line 202 | | skipping to change at line 202 | |
| ///@} | | ///@} | |
| | | | |
| /** | | /** | |
| * A global instantiation of Geocentric with the parameters for the WGS
84 | | * A global instantiation of Geocentric with the parameters for the WGS
84 | |
| * ellipsoid. | | * ellipsoid. | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static const Geocentric WGS84; | | static const Geocentric WGS84; | |
| }; | | }; | |
| | | | |
| } // namespace GeographicLib | | } // namespace GeographicLib | |
|
| #endif | | | |
| | | #endif // GEOGRAPHICLIB_GEOCENTRIC_HPP | |
| | | | |
End of changes. 16 change blocks. |
| 16 lines changed or deleted | | 16 lines changed or added | |
|
| Geodesic.hpp | | Geodesic.hpp | |
| /** | | /** | |
| * \file Geodesic.hpp | | * \file Geodesic.hpp | |
| * \brief Header for GeographicLib::Geodesic class | | * \brief Header for GeographicLib::Geodesic class | |
| * | | * | |
| * Copyright (c) Charles Karney (2009, 2010, 2011) <charles@karney.com> | | * Copyright (c) Charles Karney (2009, 2010, 2011) <charles@karney.com> | |
| * and licensed under the LGPL. For more information, see | | * and licensed under the LGPL. For more information, see | |
| * http://geographiclib.sourceforge.net/ | | * http://geographiclib.sourceforge.net/ | |
| **********************************************************************/ | | **********************************************************************/ | |
| | | | |
| #if !defined(GEOGRAPHICLIB_GEODESIC_HPP) | | #if !defined(GEOGRAPHICLIB_GEODESIC_HPP) | |
|
| #define GEOGRAPHICLIB_GEODESIC_HPP "$Id: Geodesic.hpp 6950 2011-02-11 04:09
:24Z karney $" | | #define GEOGRAPHICLIB_GEODESIC_HPP "$Id: 052cedecb4ece3fa1d25103ae750a4e05f
ee2126 $" | |
| | | | |
|
| #include "GeographicLib/Constants.hpp" | | #include <GeographicLib/Constants.hpp> | |
| | | | |
| #if !defined(GEOD_ORD) | | #if !defined(GEOD_ORD) | |
| /** | | /** | |
| * The order of the expansions used by Geodesic. | | * The order of the expansions used by Geodesic. | |
| **********************************************************************/ | | **********************************************************************/ | |
| #define GEOD_ORD (GEOGRAPHICLIB_PREC == 1 ? 6 : GEOGRAPHICLIB_PREC == 0 ? 3
: 7) | | #define GEOD_ORD (GEOGRAPHICLIB_PREC == 1 ? 6 : GEOGRAPHICLIB_PREC == 0 ? 3
: 7) | |
| #endif | | #endif | |
| | | | |
| namespace GeographicLib { | | namespace GeographicLib { | |
| | | | |
| | | | |
| skipping to change at line 94 | | skipping to change at line 94 | |
| * | | * | |
| * Overloaded versions of Geodesic::Direct, Geodesic::ArcDirect, and | | * Overloaded versions of Geodesic::Direct, Geodesic::ArcDirect, and | |
| * Geodesic::Inverse allow these quantities to be returned. In addition | | * Geodesic::Inverse allow these quantities to be returned. In addition | |
| * there are general functions Geodesic::GenDirect, and Geodesic::GenInve
rse | | * there are general functions Geodesic::GenDirect, and Geodesic::GenInve
rse | |
| * which allow an arbitrary set of results to be computed. | | * which allow an arbitrary set of results to be computed. | |
| * | | * | |
| * Additional functionality if provided by the GeodesicLine class, which | | * Additional functionality if provided by the GeodesicLine class, which | |
| * allows a sequence of points along a geodesic to be computed. | | * allows a sequence of points along a geodesic to be computed. | |
| * | | * | |
| * The calculations are accurate to better than 15 nm. See Sec. 9 of | | * The calculations are accurate to better than 15 nm. See Sec. 9 of | |
|
| * <a href="http://arxiv.org/abs/1102.1215">arXiv:1102.1215</a> for detai
ls. | | * <a href="http://arxiv.org/abs/1102.1215v1">arXiv:1102.1215v1</a> for d
etails. | |
| * | | * | |
| * The algorithms are described in | | * The algorithms are described in | |
| * - C. F. F. Karney, | | * - C. F. F. Karney, | |
|
| * <a href="http://arxiv.org/abs/1102.1215">Geodesics | | * <a href="http://arxiv.org/abs/1102.1215v1">Geodesics | |
| * on an ellipsoid of revolution</a>, | | * on an ellipsoid of revolution</a>, | |
| * Feb. 2011; | | * Feb. 2011; | |
| * preprint | | * preprint | |
|
| * <a href="http://arxiv.org/abs/1102.1215">arXiv:1102.1215</a>. | | * <a href="http://arxiv.org/abs/1102.1215v1">arXiv:1102.1215v1</a>. | |
| * . | | * . | |
| * For more information on geodesics see \ref geodesic. | | * For more information on geodesics see \ref geodesic. | |
| **********************************************************************/ | | **********************************************************************/ | |
| | | | |
|
| class Geodesic { | | class GEOGRAPHIC_EXPORT Geodesic { | |
| private: | | private: | |
| typedef Math::real real; | | typedef Math::real real; | |
| friend class GeodesicLine; | | friend class GeodesicLine; | |
|
| static const int nA1 = GEOD_ORD, nC1 = GEOD_ORD, nC1p = GEOD_ORD, | | static const int nA1_ = GEOD_ORD; | |
| nA2 = GEOD_ORD, nC2 = GEOD_ORD, | | static const int nC1_ = GEOD_ORD; | |
| nA3 = GEOD_ORD, nA3x = nA3, | | static const int nC1p_ = GEOD_ORD; | |
| nC3 = GEOD_ORD, nC3x = (nC3 * (nC3 - 1)) / 2, | | static const int nA2_ = GEOD_ORD; | |
| nC4 = GEOD_ORD, nC4x = (nC4 * (nC4 + 1)) / 2; | | static const int nC2_ = GEOD_ORD; | |
| static const unsigned maxit = 50; | | static const int nA3_ = GEOD_ORD; | |
| | | static const int nA3x_ = nA3_; | |
| | | static const int nC3_ = GEOD_ORD; | |
| | | static const int nC3x_ = (nC3_ * (nC3_ - 1)) / 2; | |
| | | static const int nC4_ = GEOD_ORD; | |
| | | static const int nC4x_ = (nC4_ * (nC4_ + 1)) / 2; | |
| | | static const unsigned maxit_ = 50; | |
| | | | |
|
| static inline real sq(real x) throw() { return x * x; } | | | |
| void Lengths(real eps, real sig12, | | void Lengths(real eps, real sig12, | |
| real ssig1, real csig1, real ssig2, real csig2, | | real ssig1, real csig1, real ssig2, real csig2, | |
| real cbet1, real cbet2, | | real cbet1, real cbet2, | |
| real& s12s, real& m12a, real& m0, | | real& s12s, real& m12a, real& m0, | |
| bool scalep, real& M12, real& M21, | | bool scalep, real& M12, real& M21, | |
| real tc[], real zc[]) const throw(); | | real tc[], real zc[]) const throw(); | |
| static real Astroid(real R, real z) throw(); | | static real Astroid(real R, real z) throw(); | |
| real InverseStart(real sbet1, real cbet1, real sbet2, real cbet2, | | real InverseStart(real sbet1, real cbet1, real sbet2, real cbet2, | |
| real lam12, | | real lam12, | |
| real& salp1, real& calp1, | | real& salp1, real& calp1, | |
| real& salp2, real& calp2, | | real& salp2, real& calp2, | |
| real C1a[], real C2a[]) const throw(); | | real C1a[], real C2a[]) const throw(); | |
| real Lambda12(real sbet1, real cbet1, real sbet2, real cbet2, | | real Lambda12(real sbet1, real cbet1, real sbet2, real cbet2, | |
| real salp1, real calp1, | | real salp1, real calp1, | |
| real& salp2, real& calp2, real& sig12, | | real& salp2, real& calp2, real& sig12, | |
| real& ssig1, real& csig1, real& ssig2, real& csig2, | | real& ssig1, real& csig1, real& ssig2, real& csig2, | |
| real& eps, real& domg12, bool diffp, real& dlam12, | | real& eps, real& domg12, bool diffp, real& dlam12, | |
| real C1a[], real C2a[], real C3a[]) | | real C1a[], real C2a[], real C3a[]) | |
| const throw(); | | const throw(); | |
| | | | |
|
| static const real eps2, tol0, tol1, tol2, xthresh; | | static const real eps2_; | |
| | | static const real tol0_; | |
| | | static const real tol1_; | |
| | | static const real tol2_; | |
| | | static const real xthresh_; | |
| const real _a, _r, _f, _f1, _e2, _ep2, _n, _b, _c2, _etol2; | | const real _a, _r, _f, _f1, _e2, _ep2, _n, _b, _c2, _etol2; | |
|
| real _A3x[nA3x], _C3x[nC3x], _C4x[nC4x]; | | real _A3x[nA3x_], _C3x[nC3x_], _C4x[nC4x_]; | |
| static real SinCosSeries(bool sinp, | | static real SinCosSeries(bool sinp, | |
| real sinx, real cosx, const real c[], int n) | | real sinx, real cosx, const real c[], int n) | |
| throw(); | | throw(); | |
| static inline real AngNormalize(real x) throw() { | | static inline real AngNormalize(real x) throw() { | |
| // Place angle in [-180, 180). Assumes x is in [-540, 540). | | // Place angle in [-180, 180). Assumes x is in [-540, 540). | |
| return x >= 180 ? x - 360 : x < -180 ? x + 360 : x; | | return x >= 180 ? x - 360 : x < -180 ? x + 360 : x; | |
| } | | } | |
| static inline real AngRound(real x) throw() { | | static inline real AngRound(real x) throw() { | |
| // The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^
57 | | // The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^
57 | |
| // for reals = 0.7 pm on the earth if x is an angle in degrees. (Thi
s | | // for reals = 0.7 pm on the earth if x is an angle in degrees. (Thi
s | |
| | | | |
| skipping to change at line 820 | | skipping to change at line 829 | |
| return s12; | | return s12; | |
| } else { | | } else { | |
| real s12 = s12_a12; | | real s12 = s12_a12; | |
| return Direct(lat1, lon1, azi1, s12, lat2, lon2, azi2, m12); | | return Direct(lat1, lon1, azi1, s12, lat2, lon2, azi2, m12); | |
| } | | } | |
| } | | } | |
| ///@} | | ///@} | |
| }; | | }; | |
| | | | |
| } // namespace GeographicLib | | } // namespace GeographicLib | |
|
| #endif | | | |
| | | #endif // GEOGRAPHICLIB_GEODESIC_HPP | |
| | | | |
End of changes. 11 change blocks. |
| 15 lines changed or deleted | | 24 lines changed or added | |
|
| GeodesicLine.hpp | | GeodesicLine.hpp | |
| /** | | /** | |
| * \file GeodesicLine.hpp | | * \file GeodesicLine.hpp | |
| * \brief Header for GeographicLib::GeodesicLine class | | * \brief Header for GeographicLib::GeodesicLine class | |
| * | | * | |
| * Copyright (c) Charles Karney (2009, 2010, 2011) <charles@karney.com> | | * Copyright (c) Charles Karney (2009, 2010, 2011) <charles@karney.com> | |
| * and licensed under the LGPL. For more information, see | | * and licensed under the LGPL. For more information, see | |
| * http://geographiclib.sourceforge.net/ | | * http://geographiclib.sourceforge.net/ | |
| **********************************************************************/ | | **********************************************************************/ | |
| | | | |
| #if !defined(GEOGRAPHICLIB_GEODESICLINE_HPP) | | #if !defined(GEOGRAPHICLIB_GEODESICLINE_HPP) | |
|
| #define GEOGRAPHICLIB_GEODESICLINE_HPP "$Id: GeodesicLine.hpp 6950 2011-02-
11 04:09:24Z karney $" | | #define GEOGRAPHICLIB_GEODESICLINE_HPP "$Id: 5d443e356d7037207c0298f930948d
0712c2fe3d $" | |
| | | | |
|
| #include "GeographicLib/Constants.hpp" | | #include <GeographicLib/Constants.hpp> | |
| #include "GeographicLib/Geodesic.hpp" | | #include <GeographicLib/Geodesic.hpp> | |
| | | | |
| namespace GeographicLib { | | namespace GeographicLib { | |
| | | | |
| /** | | /** | |
| * \brief A geodesic line. | | * \brief A geodesic line. | |
| * | | * | |
| * GeodesicLine facilitates the determination of a series of points on a | | * GeodesicLine facilitates the determination of a series of points on a | |
| * single geodesic. The starting point (\e lat1, \e lon1) and the azimut
h \e | | * single geodesic. The starting point (\e lat1, \e lon1) and the azimut
h \e | |
| * azi1 are specified in the constructor. GeodesicLine.Position returns
the | | * azi1 are specified in the constructor. GeodesicLine.Position returns
the | |
| * location of point 2 a distance \e s12 along the geodesic. Alternative
ly | | * location of point 2 a distance \e s12 along the geodesic. Alternative
ly | |
| * GeodesicLine.ArcPosition gives the position of point 2 an arc length \
e | | * GeodesicLine.ArcPosition gives the position of point 2 an arc length \
e | |
| * a12 along the geodesic. An example of use of this class is: | | * a12 along the geodesic. An example of use of this class is: | |
| \code | | \code | |
| #include <iostream> | | #include <iostream> | |
| #include <iomanip> | | #include <iomanip> | |
| #include <cmath> | | #include <cmath> | |
|
| #include "GeographicLib/Geodesic.hpp" | | #include <GeographicLib/Geodesic.hpp> | |
| #include "GeographicLib/GeodesicLine.hpp" | | #include <GeographicLib/GeodesicLine.hpp> | |
| | | | |
| int main() { | | int main() { | |
|
| // Print waypoints between JFK and SIN at | | // Print waypoints between JFK and SIN at approximately | |
| // approximately 100km intervals. | | // equal 100km intervals | |
| double | | double | |
| lat1 = 40.640, lon1 = -73.779, // JFK | | lat1 = 40.640, lon1 = -73.779, // JFK | |
|
| lat2 = 1.359, lon2 = 177.486; // SIN | | lat2 = 1.359, lon2 = 177.486; // SIN | |
| const GeographicLib::Geodesic& | | const GeographicLib::Geodesic& | |
| g = GeographicLib::Geodesic::WGS84; | | g = GeographicLib::Geodesic::WGS84; | |
|
| double azi1, azi2, | | double s12, azi1, azi2, | |
| a12 = g.Inverse(lat1, lon1, lat2, lon2, azi1, azi2); | | a12 = g.Inverse(lat1, lon1, lat2, lon2, s12, azi1, azi2); | |
| double ds = 100e3; // Nominal distance between points = 100 km | | double ds = 100e3; // Nominal distance between points = 100 km | |
| int num = std::ceil(a12/ (90 * ds / 10e6)); // 90 deg = 10e6 m | | int num = std::ceil(s12 / ds); // Number of intervals | |
| double da = a12/num; // Arc length between points | | double da = a12/num; // Arc length of each interval | |
| const GeographicLib::GeodesicLine l(g, lat1, lon1, azi1); | | const GeographicLib::GeodesicLine l(g, lat1, lon1, azi1); | |
| std::cout << std::fixed << std::setprecision(3); | | std::cout << std::fixed << std::setprecision(3); | |
| for (int i = 0; i <= num; ++i) { | | for (int i = 0; i <= num; ++i) { | |
| double lat, lon; | | double lat, lon; | |
| l.ArcPosition(i * da, lat, lon); | | l.ArcPosition(i * da, lat, lon); | |
| std::cout << lat << " " << lon << "\n"; | | std::cout << lat << " " << lon << "\n"; | |
| } | | } | |
| } | | } | |
| \endcode | | \endcode | |
| * The default copy constructor and assignment operators work with this | | * The default copy constructor and assignment operators work with this | |
| * class, so that, for example, the previous example could start | | * class, so that, for example, the previous example could start | |
| \code | | \code | |
| GeographicLib::GeodesicLine l; | | GeographicLib::GeodesicLine l; | |
| l = g.Line(lat1, lon1, azi1); | | l = g.Line(lat1, lon1, azi1); | |
| \endcode | | \endcode | |
| * Similarly, a vector can be used to hold GeodesicLine objects. | | * Similarly, a vector can be used to hold GeodesicLine objects. | |
| * | | * | |
| * The calculations are accurate to better than 15 nm. See Sec. 9 of | | * The calculations are accurate to better than 15 nm. See Sec. 9 of | |
|
| * <a href="http://arxiv.org/abs/1102.1215">arXiv:1102.1215</a> for detai
ls. | | * <a href="http://arxiv.org/abs/1102.1215v1">arXiv:1102.1215v1</a> for d
etails. | |
| * | | * | |
| * The algorithms are described in | | * The algorithms are described in | |
| * - C. F. F. Karney, | | * - C. F. F. Karney, | |
|
| * <a href="http://arxiv.org/abs/1102.1215">Geodesics | | * <a href="http://arxiv.org/abs/1102.1215v1">Geodesics | |
| * on an ellipsoid of revolution</a>, | | * on an ellipsoid of revolution</a>, | |
| * Feb. 2011; | | * Feb. 2011; | |
| * preprint | | * preprint | |
|
| * <a href="http://arxiv.org/abs/1102.1215">arXiv:1102.1215</a>. | | * <a href="http://arxiv.org/abs/1102.1215v1">arXiv:1102.1215v1</a>. | |
| * . | | * . | |
| * For more information on geodesics see \ref geodesic. | | * For more information on geodesics see \ref geodesic. | |
| **********************************************************************/ | | **********************************************************************/ | |
| | | | |
|
| class GeodesicLine { | | class GEOGRAPHIC_EXPORT GeodesicLine { | |
| private: | | private: | |
| typedef Math::real real; | | typedef Math::real real; | |
| friend class Geodesic; | | friend class Geodesic; | |
|
| static const int nC1 = Geodesic::nC1, nC1p = Geodesic::nC1p, | | static const int nC1_ = Geodesic::nC1_; | |
| nC2 = Geodesic::nC2, nC3 = Geodesic::nC3, nC4 = Geodesic::nC4; | | static const int nC1p_ = Geodesic::nC1p_; | |
| | | static const int nC2_ = Geodesic::nC2_; | |
| | | static const int nC3_ = Geodesic::nC3_; | |
| | | static const int nC4_ = Geodesic::nC4_; | |
| | | | |
| real _lat1, _lon1, _azi1; | | real _lat1, _lon1, _azi1; | |
| real _a, _r, _b, _c2, _f1, _salp0, _calp0, _k2, | | real _a, _r, _b, _c2, _f1, _salp0, _calp0, _k2, | |
| _salp1, _calp1, _ssig1, _csig1, _stau1, _ctau1, _somg1, _comg1, | | _salp1, _calp1, _ssig1, _csig1, _stau1, _ctau1, _somg1, _comg1, | |
| _A1m1, _A2m1, _A3c, _B11, _B21, _B31, _A4, _B41; | | _A1m1, _A2m1, _A3c, _B11, _B21, _B31, _A4, _B41; | |
| // index zero elements of _C1a, _C1pa, _C2a, _C3a are unused | | // index zero elements of _C1a, _C1pa, _C2a, _C3a are unused | |
|
| real _C1a[nC1 + 1], _C1pa[nC1p + 1], _C2a[nC2 + 1], _C3a[nC3], | | real _C1a[nC1_ + 1], _C1pa[nC1p_ + 1], _C2a[nC2_ + 1], _C3a[nC3_], | |
| _C4a[nC4]; // all the elements of _C4a are used | | _C4a[nC4_]; // all the elements of _C4a are used | |
| bool _areap; | | bool _areap; | |
| unsigned _caps; | | unsigned _caps; | |
| | | | |
|
| static inline real sq(real x) throw() { return x * x; } | | | |
| enum captype { | | enum captype { | |
| CAP_NONE = Geodesic::CAP_NONE, | | CAP_NONE = Geodesic::CAP_NONE, | |
| CAP_C1 = Geodesic::CAP_C1, | | CAP_C1 = Geodesic::CAP_C1, | |
| CAP_C1p = Geodesic::CAP_C1p, | | CAP_C1p = Geodesic::CAP_C1p, | |
| CAP_C2 = Geodesic::CAP_C2, | | CAP_C2 = Geodesic::CAP_C2, | |
| CAP_C3 = Geodesic::CAP_C3, | | CAP_C3 = Geodesic::CAP_C3, | |
| CAP_C4 = Geodesic::CAP_C4, | | CAP_C4 = Geodesic::CAP_C4, | |
| CAP_ALL = Geodesic::CAP_ALL, | | CAP_ALL = Geodesic::CAP_ALL, | |
| OUT_ALL = Geodesic::OUT_ALL, | | OUT_ALL = Geodesic::OUT_ALL, | |
| }; | | }; | |
| | | | |
| skipping to change at line 633 | | skipping to change at line 635 | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| void Scale(real a12, real& M12, real& M21) const throw() { | | void Scale(real a12, real& M12, real& M21) const throw() { | |
| real lat2, lon2, azi2, s12; | | real lat2, lon2, azi2, s12; | |
| ArcPosition(a12, lat2, lon2, azi2, s12, M12, M21); | | ArcPosition(a12, lat2, lon2, azi2, s12, M12, M21); | |
| } | | } | |
| ///@} | | ///@} | |
| | | | |
| }; | | }; | |
| | | | |
| } // namespace GeographicLib | | } // namespace GeographicLib | |
|
| #endif | | | |
| | | #endif // GEOGRAPHICLIB_GEODESICLINE_HPP | |
| | | | |
End of changes. 14 change blocks. |
| 22 lines changed or deleted | | 24 lines changed or added | |
|
| Geoid.hpp | | Geoid.hpp | |
| /** | | /** | |
| * \file Geoid.hpp | | * \file Geoid.hpp | |
| * \brief Header for GeographicLib::Geoid class | | * \brief Header for GeographicLib::Geoid class | |
| * | | * | |
| * Copyright (c) Charles Karney (2009, 2010, 2011) <charles@karney.com> | | * Copyright (c) Charles Karney (2009, 2010, 2011) <charles@karney.com> | |
| * and licensed under the LGPL. For more information, see | | * and licensed under the LGPL. For more information, see | |
| * http://geographiclib.sourceforge.net/ | | * http://geographiclib.sourceforge.net/ | |
| **********************************************************************/ | | **********************************************************************/ | |
| | | | |
| #if !defined(GEOGRAPHICLIB_GEOID_HPP) | | #if !defined(GEOGRAPHICLIB_GEOID_HPP) | |
|
| #define GEOGRAPHICLIB_GEOID_HPP "$Id: Geoid.hpp 6937 2011-02-01 20:17:13Z k
arney $" | | #define GEOGRAPHICLIB_GEOID_HPP "$Id: a7310c17a3b20283af1640db0dc1de40974d9
10c $" | |
| | | | |
|
| #include "GeographicLib/Constants.hpp" | | | |
| #include <string> | | #include <string> | |
| #include <vector> | | #include <vector> | |
| #include <fstream> | | #include <fstream> | |
|
| | | #include <GeographicLib/Constants.hpp> | |
| | | | |
| | | #if defined(_MSC_VER) | |
| | | // Squelch warnings about dll vs vector | |
| | | #pragma warning (push) | |
| | | #pragma warning (disable: 4251) | |
| | | #endif | |
| | | | |
| namespace GeographicLib { | | namespace GeographicLib { | |
| | | | |
| /** | | /** | |
| * \brief Looking up the height of the geoid | | * \brief Looking up the height of the geoid | |
| * | | * | |
| * This class evaluated the height of one of the standard geoids, EGM84, | | * This class evaluated the height of one of the standard geoids, EGM84, | |
| * EGM96, or EGM2008 by bilinear or cubic interpolation into a rectangula
r | | * EGM96, or EGM2008 by bilinear or cubic interpolation into a rectangula
r | |
| * grid of data. These geoid models are documented in | | * grid of data. These geoid models are documented in | |
| * - EGM84: | | * - EGM84: | |
| | | | |
| skipping to change at line 45 | | skipping to change at line 51 | |
| * this class evaluates the height by interpolation inot a grid of | | * this class evaluates the height by interpolation inot a grid of | |
| * precomputed values. | | * precomputed values. | |
| * | | * | |
| * See \ref geoid for details of how to install the data sets, the data | | * See \ref geoid for details of how to install the data sets, the data | |
| * format, estimates of the interpolation errors, and how to use caching. | | * format, estimates of the interpolation errors, and how to use caching. | |
| * | | * | |
| * In addition to returning the geoid height, the gradient of the geoid c
an | | * In addition to returning the geoid height, the gradient of the geoid c
an | |
| * be calculated. The gradient is defined as the rate of change of the g
eoid | | * be calculated. The gradient is defined as the rate of change of the g
eoid | |
| * as a function of position on the ellipsoid. This uses the parameters
for | | * as a function of position on the ellipsoid. This uses the parameters
for | |
| * the WGS84 ellipsoid. The gradient defined in terms of the interpolate
d | | * the WGS84 ellipsoid. The gradient defined in terms of the interpolate
d | |
|
| * heights. | | * heights. As a result of the way that the geoid data is stored, the | |
| | | * calculation of gradients can result in large quantization errors. Thi | |
| | | s is | |
| | | * particularly acute at high latitudes and for the easterly gradient. | |
| * | | * | |
| * This class is typically \e not thread safe in that a single instantiat
ion | | * This class is typically \e not thread safe in that a single instantiat
ion | |
| * cannot be safely used by multiple threads because of the way the objec
t | | * cannot be safely used by multiple threads because of the way the objec
t | |
| * reads the data set and because it maintains a single-cell cache. If | | * reads the data set and because it maintains a single-cell cache. If | |
| * multiple threads need to calculate geoid heights they should all const
ruct | | * multiple threads need to calculate geoid heights they should all const
ruct | |
| * thread-local instantiations. Alternatively, set the optional \e | | * thread-local instantiations. Alternatively, set the optional \e | |
| * threadsafe parameter to true in the constuctor. This causes the | | * threadsafe parameter to true in the constuctor. This causes the | |
| * constructor to read all the data into memory and to turn off the | | * constructor to read all the data into memory and to turn off the | |
| * single-cell caching which results in a Geoid object which \e is thread | | * single-cell caching which results in a Geoid object which \e is thread | |
| * safe. | | * safe. | |
| **********************************************************************/ | | **********************************************************************/ | |
| | | | |
|
| class Geoid { | | class GEOGRAPHIC_EXPORT Geoid { | |
| private: | | private: | |
| typedef Math::real real; | | typedef Math::real real; | |
|
| static const unsigned stencilsize = 12; | | static const unsigned stencilsize_ = 12; | |
| static const unsigned nterms = ((3 + 1) * (3 + 2))/2; // for a cubic fi | | static const unsigned nterms_ = ((3 + 1) * (3 + 2))/2; // for a cubic f | |
| t | | it | |
| static const real c0, c0n, c0s; | | static const real c0_; | |
| static const real c3[stencilsize * nterms]; | | static const real c0n_; | |
| static const real c3n[stencilsize * nterms]; | | static const real c0s_; | |
| static const real c3s[stencilsize * nterms]; | | static const real c3_[stencilsize_ * nterms_]; | |
| | | static const real c3n_[stencilsize_ * nterms_]; | |
| | | static const real c3s_[stencilsize_ * nterms_]; | |
| | | | |
| std::string _name, _dir, _filename; | | std::string _name, _dir, _filename; | |
| const bool _cubic; | | const bool _cubic; | |
| const real _a, _e2, _degree, _eps; | | const real _a, _e2, _degree, _eps; | |
| mutable std::ifstream _file; | | mutable std::ifstream _file; | |
| real _rlonres, _rlatres; | | real _rlonres, _rlatres; | |
| std::string _description, _datetime; | | std::string _description, _datetime; | |
| real _offset, _scale, _maxerror, _rmserror; | | real _offset, _scale, _maxerror, _rmserror; | |
| int _width, _height; | | int _width, _height; | |
| unsigned long long _datastart, _swidth; | | unsigned long long _datastart, _swidth; | |
| bool _threadsafe; | | bool _threadsafe; | |
| // Area cache | | // Area cache | |
| mutable std::vector< std::vector<unsigned short> > _data; | | mutable std::vector< std::vector<unsigned short> > _data; | |
| mutable bool _cache; | | mutable bool _cache; | |
| // NE corner and extent of cache | | // NE corner and extent of cache | |
| mutable int _xoffset, _yoffset, _xsize, _ysize; | | mutable int _xoffset, _yoffset, _xsize, _ysize; | |
| // Cell cache | | // Cell cache | |
| mutable int _ix, _iy; | | mutable int _ix, _iy; | |
| mutable real _v00, _v01, _v10, _v11; | | mutable real _v00, _v01, _v10, _v11; | |
|
| mutable real _t[nterms]; | | mutable real _t[nterms_]; | |
| void filepos(int ix, int iy) const { | | void filepos(int ix, int iy) const { | |
| _file.seekg( | | _file.seekg( | |
| #if !(defined(__GNUC__) && __GNUC__ < 4) | | #if !(defined(__GNUC__) && __GNUC__ < 4) | |
| // g++ 3.x doesn't know about the cast to streamoff. | | // g++ 3.x doesn't know about the cast to streamoff. | |
| std::ios::streamoff | | std::ios::streamoff | |
| #endif | | #endif | |
| (_datastart + 2ULL * (unsigned(iy)*_swidth + unsigned(ix)
))); | | (_datastart + 2ULL * (unsigned(iy)*_swidth + unsigned(ix)
))); | |
| } | | } | |
| real rawval(int ix, int iy) const { | | real rawval(int ix, int iy) const { | |
| if (ix < 0) | | if (ix < 0) | |
| | | | |
| skipping to change at line 450 | | skipping to change at line 460 | |
| static std::string DefaultPath(); | | static std::string DefaultPath(); | |
| | | | |
| /** | | /** | |
| * <b>DEPRECATED</b> Return the value of the environment variable | | * <b>DEPRECATED</b> Return the value of the environment variable | |
| * GEOID_PATH. | | * GEOID_PATH. | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static std::string GeoidPath(); | | static std::string GeoidPath(); | |
| }; | | }; | |
| | | | |
| } // namespace GeographicLib | | } // namespace GeographicLib | |
|
| | | | |
| | | #if defined(_MSC_VER) | |
| | | #pragma warning (pop) | |
| #endif | | #endif | |
|
| | | | |
| | | #endif // GEOGRAPHICLIB_GEOID_HPP | |
| | | | |
End of changes. 9 change blocks. |
| 12 lines changed or deleted | | 26 lines changed or added | |
|
| LambertConformalConic.hpp | | LambertConformalConic.hpp | |
| /** | | /** | |
| * \file LambertConformalConic.hpp | | * \file LambertConformalConic.hpp | |
| * \brief Header for GeographicLib::LambertConformalConic class | | * \brief Header for GeographicLib::LambertConformalConic class | |
| * | | * | |
| * Copyright (c) Charles Karney (2010, 2011) <charles@karney.com> and licen
sed | | * Copyright (c) Charles Karney (2010, 2011) <charles@karney.com> and licen
sed | |
| * under the LGPL. For more information, see | | * under the LGPL. For more information, see | |
| * http://geographiclib.sourceforge.net/ | | * http://geographiclib.sourceforge.net/ | |
| **********************************************************************/ | | **********************************************************************/ | |
| | | | |
| #if !defined(GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_HPP) | | #if !defined(GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_HPP) | |
|
| #define GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_HPP "$Id: LambertConformalConic
.hpp 6937 2011-02-01 20:17:13Z karney $" | | #define GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_HPP "$Id: 2d1efcd0fe69672e8c21b
87133c7c04e8a7b45f6 $" | |
| | | | |
|
| #include "GeographicLib/Constants.hpp" | | | |
| #include <algorithm> | | #include <algorithm> | |
|
| | | #include <GeographicLib/Constants.hpp> | |
| | | | |
| namespace GeographicLib { | | namespace GeographicLib { | |
| | | | |
| /** | | /** | |
| * \brief Lambert Conformal Conic Projection | | * \brief Lambert Conformal Conic Projection | |
| * | | * | |
| * Implementation taken from the report, | | * Implementation taken from the report, | |
| * - J. P. Snyder, | | * - J. P. Snyder, | |
| * <a href="http://pubs.er.usgs.gov/usgspubs/pp/pp1395"> Map Projection
s: A | | * <a href="http://pubs.er.usgs.gov/usgspubs/pp/pp1395"> Map Projection
s: A | |
| * Working Manual</a>, USGS Professional Paper 1395 (1987), | | * Working Manual</a>, USGS Professional Paper 1395 (1987), | |
| * pp. 107–109. | | * pp. 107–109. | |
| * | | * | |
| * This is a implementation of the equations in Snyder except that divide
d | | * This is a implementation of the equations in Snyder except that divide
d | |
| * differences have been used to transform the expressions into ones whic
h | | * differences have been used to transform the expressions into ones whic
h | |
| * may be evaluated accurately and that Newton's method is used to invert
the | | * may be evaluated accurately and that Newton's method is used to invert
the | |
| * projection. In this implementation, the projection correctly becomes
the | | * projection. In this implementation, the projection correctly becomes
the | |
|
| * Mercator projection or the polar sterographic projection when the stan
dard | | * Mercator_ projection or the polar sterographic projection when the sta
ndard | |
| * latitude is the equator or a pole. The accuracy of the projections is | | * latitude is the equator or a pole. The accuracy of the projections is | |
| * about 10 nm. | | * about 10 nm. | |
| * | | * | |
| * The ellipsoid parameters, the standard parallels, and the scale on the | | * The ellipsoid parameters, the standard parallels, and the scale on the | |
| * standard parallels are set in the constructor. Internally, the case w
ith | | * standard parallels are set in the constructor. Internally, the case w
ith | |
| * two standard parallels is converted into a single standard parallel, t
he | | * two standard parallels is converted into a single standard parallel, t
he | |
| * latitude of tangency (also the latitude of minimum scale), with a scal
e | | * latitude of tangency (also the latitude of minimum scale), with a scal
e | |
| * specified on this parallel. This latitude is also used as the latitud
e of | | * specified on this parallel. This latitude is also used as the latitud
e of | |
| * origin which is returned by LambertConformalConic::OriginLatitude. Th
e | | * origin which is returned by LambertConformalConic::OriginLatitude. Th
e | |
| * scale on the latitude of origin is given by | | * scale on the latitude of origin is given by | |
| | | | |
| skipping to change at line 56 | | skipping to change at line 56 | |
| * LambertConformalConic::Reverse functions. There is no provision in th
is | | * LambertConformalConic::Reverse functions. There is no provision in th
is | |
| * class for specifying a false easting or false northing or a different | | * class for specifying a false easting or false northing or a different | |
| * latitude of origin. However these are can be simply included by the | | * latitude of origin. However these are can be simply included by the | |
| * calling function. For example the Pennsylvania South state coordinate | | * calling function. For example the Pennsylvania South state coordinate | |
| * system (<a href="http://www.spatialreference.org/ref/epsg/3364/"> | | * system (<a href="http://www.spatialreference.org/ref/epsg/3364/"> | |
| * EPSG:3364</a>) is obtained by: | | * EPSG:3364</a>) is obtained by: | |
| \code | | \code | |
| const double | | const double | |
| a = GeographicLib::Constants::WGS84_a<double>(), | | a = GeographicLib::Constants::WGS84_a<double>(), | |
| r = 298.257222101, // GRS80 | | r = 298.257222101, // GRS80 | |
|
| lat1 = 39 + 56/60.0, lat1 = 40 + 58/60.0, // standard parallels | | lat1 = 40 + 58/60.0, lat2 = 39 + 56/60.0, // standard parallels | |
| k1 = 1, // scale | | k1 = 1, // scale | |
|
| lat0 = 39 + 20/60.0, lon0 = 75 + 45/60.0, // origin | | lat0 = 39 + 20/60.0, lon0 =-77 - 45/60.0, // origin | |
| fe = 600000, fn = 0; // false easting and northin
g | | fe = 600000, fn = 0; // false easting and northin
g | |
| // Set up basic projection | | // Set up basic projection | |
| const GeographicLib::LambertConformalConic PASouth(a, r, lat1, lat2, k1)
; | | const GeographicLib::LambertConformalConic PASouth(a, r, lat1, lat2, k1)
; | |
| double x0, y0; | | double x0, y0; | |
| { | | { | |
| // Transform origin point | | // Transform origin point | |
| PASouth.Forward(lon0, lat0, lon0, x0, y0); | | PASouth.Forward(lon0, lat0, lon0, x0, y0); | |
| x0 -= fe; y0 -= fn; // Combine result with false origin | | x0 -= fe; y0 -= fn; // Combine result with false origin | |
| } | | } | |
| double lat, lon, x, y; | | double lat, lon, x, y; | |
| | | | |
| skipping to change at line 81 | | skipping to change at line 81 | |
| PASouth.Forward(lon0, lat, lon, x, y); | | PASouth.Forward(lon0, lat, lon, x, y); | |
| x -= x0; y -= y0; | | x -= x0; y -= y0; | |
| std::cout << x << " " << y << "\n"; | | std::cout << x << " " << y << "\n"; | |
| // Sample conversion from PASouth grid to geodetic | | // Sample conversion from PASouth grid to geodetic | |
| std::cin >> x >> y; | | std::cin >> x >> y; | |
| x += x0; y += y0; | | x += x0; y += y0; | |
| PASouth.Reverse(lon0, x, y, lat, lon); | | PASouth.Reverse(lon0, x, y, lat, lon); | |
| std::cout << lat << " " << lon << "\n"; | | std::cout << lat << " " << lon << "\n"; | |
| \endcode | | \endcode | |
| **********************************************************************/ | | **********************************************************************/ | |
|
| class LambertConformalConic { | | class GEOGRAPHIC_EXPORT LambertConformalConic { | |
| private: | | private: | |
| typedef Math::real real; | | typedef Math::real real; | |
| const real _a, _r, _f, _fm, _e2, _e, _e2m; | | const real _a, _r, _f, _fm, _e2, _e, _e2m; | |
| real _sign, _n, _nc, _t0nm1, _scale, _lat0, _k0; | | real _sign, _n, _nc, _t0nm1, _scale, _lat0, _k0; | |
| real _scbet0, _tchi0, _scchi0, _psi0, _nrho0; | | real _scbet0, _tchi0, _scchi0, _psi0, _nrho0; | |
|
| static const real eps, epsx, tol, ahypover; | | static const real eps_; | |
| static const int numit = 5; | | static const real epsx_; | |
| static inline real sq(real x) throw() { return x * x; } | | static const real tol_; | |
| | | static const real ahypover_; | |
| | | static const int numit_ = 5; | |
| static inline real hyp(real x) throw() { return Math::hypot(real(1), x)
; } | | static inline real hyp(real x) throw() { return Math::hypot(real(1), x)
; } | |
| // e * atanh(e * x) = log( ((1 + e*x)/(1 - e*x))^(e/2) ) if f >= 0 | | // e * atanh(e * x) = log( ((1 + e*x)/(1 - e*x))^(e/2) ) if f >= 0 | |
| // - sqrt(-e2) * atan( sqrt(-e2) * x) if f < 0 | | // - sqrt(-e2) * atan( sqrt(-e2) * x) if f < 0 | |
| inline real eatanhe(real x) const throw() { | | inline real eatanhe(real x) const throw() { | |
| return _f >= 0 ? _e * Math::atanh(_e * x) : - _e * std::atan(_e * x); | | return _f >= 0 ? _e * Math::atanh(_e * x) : - _e * std::atan(_e * x); | |
| } | | } | |
| // Divided differences | | // Divided differences | |
| // Definition: Df(x,y) = (f(x)-f(y))/(x-y) | | // Definition: Df(x,y) = (f(x)-f(y))/(x-y) | |
| // See: W. M. Kahan and R. J. Fateman, | | // See: W. M. Kahan and R. J. Fateman, | |
| // Symbolic computation of divided differences, | | // Symbolic computation of divided differences, | |
| // SIGSAM Bull. 33(3), 7-28 (1999) | | // SIGSAM Bull. 33(3), 7-28 (1999) | |
|
| // http://doi.acm.org/10.1145/334714.334716 | | // http://dx.doi.org/10.1145/334714.334716 | |
| // http://www.cs.berkeley.edu/~fateman/papers/divdiff.pdf | | // http://www.cs.berkeley.edu/~fateman/papers/divdiff.pdf | |
| // | | // | |
| // General rules | | // General rules | |
| // h(x) = f(g(x)): Dh(x,y) = Df(g(x),g(y))*Dg(x,y) | | // h(x) = f(g(x)): Dh(x,y) = Df(g(x),g(y))*Dg(x,y) | |
| // h(x) = f(x)*g(x): | | // h(x) = f(x)*g(x): | |
| // Dh(x,y) = Df(x,y)*g(x) + Dg(x,y)*f(y) | | // Dh(x,y) = Df(x,y)*g(x) + Dg(x,y)*f(y) | |
| // = Df(x,y)*g(y) + Dg(x,y)*f(x) | | // = Df(x,y)*g(y) + Dg(x,y)*f(x) | |
| // = Df(x,y)*(g(x)+g(y))/2 + Dg(x,y)*(f(x)+f(y))/2 | | // = Df(x,y)*(g(x)+g(y))/2 + Dg(x,y)*(f(x)+f(y))/2 | |
| // | | // | |
| // hyp(x) = sqrt(1+x^2): Dhyp(x,y) = (x+y)/(hyp(x)+hyp(y)) | | // hyp(x) = sqrt(1+x^2): Dhyp(x,y) = (x+y)/(hyp(x)+hyp(y)) | |
| static inline real Dhyp(real x, real y, real hx, real hy) throw() | | static inline real Dhyp(real x, real y, real hx, real hy) throw() | |
| // hx = hyp(x) | | // hx = hyp(x) | |
| { return (x + y) / (hx + hy); } | | { return (x + y) / (hx + hy); } | |
| // sn(x) = x/sqrt(1+x^2): Dsn(x,y) = (x+y)/((sn(x)+sn(y))*(1+x^2)*(1+y^
2)) | | // sn(x) = x/sqrt(1+x^2): Dsn(x,y) = (x+y)/((sn(x)+sn(y))*(1+x^2)*(1+y^
2)) | |
| static inline real Dsn(real x, real y, real sx, real sy) throw() { | | static inline real Dsn(real x, real y, real sx, real sy) throw() { | |
| // sx = x/hyp(x) | | // sx = x/hyp(x) | |
| real t = x * y; | | real t = x * y; | |
|
| return t > 0 ? (x + y) * sq( (sx * sy)/t ) / (sx + sy) : | | return t > 0 ? (x + y) * Math::sq( (sx * sy)/t ) / (sx + sy) : | |
| (x - y != 0 ? (sx - sy) / (x - y) : 1); | | (x - y != 0 ? (sx - sy) / (x - y) : 1); | |
| } | | } | |
| // Dlog1p(x,y) = log1p((x-y)/(1+y)/(x-y) | | // Dlog1p(x,y) = log1p((x-y)/(1+y)/(x-y) | |
| static inline real Dlog1p(real x, real y) throw() { | | static inline real Dlog1p(real x, real y) throw() { | |
| real t = x - y; if (t < 0) { t = -t; y = x; } | | real t = x - y; if (t < 0) { t = -t; y = x; } | |
| return t != 0 ? Math::log1p(t / (1 + y)) / t : 1 / (1 + x); | | return t != 0 ? Math::log1p(t / (1 + y)) / t : 1 / (1 + x); | |
| } | | } | |
| // Dexp(x,y) = exp((x+y)/2) * 2*sinh((x-y)/2)/(x-y) | | // Dexp(x,y) = exp((x+y)/2) * 2*sinh((x-y)/2)/(x-y) | |
| static inline real Dexp(real x, real y) throw() { | | static inline real Dexp(real x, real y) throw() { | |
| real t = (x - y)/2; | | real t = (x - y)/2; | |
| | | | |
| skipping to change at line 332 | | skipping to change at line 334 | |
| /** | | /** | |
| * @return central scale for the projection. This is the scale on the | | * @return central scale for the projection. This is the scale on the | |
| * latitude of origin. | | * latitude of origin. | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| Math::real CentralScale() const throw() { return _k0; } | | Math::real CentralScale() const throw() { return _k0; } | |
| ///@} | | ///@} | |
| | | | |
| /** | | /** | |
| * A global instantiation of LambertConformalConic with the WGS84 | | * A global instantiation of LambertConformalConic with the WGS84 | |
| * ellipsoid, \e stdlat = 0, and \e k0 = 1. This degenerates to the | | * ellipsoid, \e stdlat = 0, and \e k0 = 1. This degenerates to the | |
|
| * Mercator projection. | | * Mercator_ projection. | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static const LambertConformalConic Mercator; | | static const LambertConformalConic Mercator; | |
| }; | | }; | |
| | | | |
| } // namespace GeographicLib | | } // namespace GeographicLib | |
| | | | |
|
| #endif | | #endif // GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_HPP | |
| | | | |
End of changes. 12 change blocks. |
| 12 lines changed or deleted | | 14 lines changed or added | |
|
| LocalCartesian.hpp | | LocalCartesian.hpp | |
| /** | | /** | |
| * \file LocalCartesian.hpp | | * \file LocalCartesian.hpp | |
| * \brief Header for GeographicLib::LocalCartesian class | | * \brief Header for GeographicLib::LocalCartesian class | |
| * | | * | |
|
| * Copyright (c) Charles Karney (2008, 2009, 2010) <charles@karney.com> | | * Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co
m> | |
| * and licensed under the LGPL. For more information, see | | * and licensed under the LGPL. For more information, see | |
| * http://geographiclib.sourceforge.net/ | | * http://geographiclib.sourceforge.net/ | |
| **********************************************************************/ | | **********************************************************************/ | |
| | | | |
| #if !defined(GEOGRAPHICLIB_LOCALCARTESIAN_HPP) | | #if !defined(GEOGRAPHICLIB_LOCALCARTESIAN_HPP) | |
|
| #define GEOGRAPHICLIB_LOCALCARTESIAN_HPP "$Id: LocalCartesian.hpp 6952 2011
-02-14 20:26:44Z karney $" | | #define GEOGRAPHICLIB_LOCALCARTESIAN_HPP "$Id: f9346d8010f4b5b4e74a65961b47
2c73578dfbfe $" | |
| | | | |
|
| #include "GeographicLib/Geocentric.hpp" | | #include <GeographicLib/Geocentric.hpp> | |
| #include "GeographicLib/Constants.hpp" | | #include <GeographicLib/Constants.hpp> | |
| | | | |
| namespace GeographicLib { | | namespace GeographicLib { | |
| | | | |
| /** | | /** | |
| * \brief Local Cartesian coordinates | | * \brief Local Cartesian coordinates | |
| * | | * | |
| * Convert between geodetic coordinates latitude = \e lat, longitude = \e | | * Convert between geodetic coordinates latitude = \e lat, longitude = \e | |
| * lon, height = \e h (measured vertically from the surface of the ellips
oid) | | * lon, height = \e h (measured vertically from the surface of the ellips
oid) | |
| * to local cartesian coordinates (\e x, \e y, \e z). The origin of loca
l | | * to local cartesian coordinates (\e x, \e y, \e z). The origin of loca
l | |
| * cartesian coordinate system is at \e lat = \e lat0, \e lon = \e lon0,
\e h | | * cartesian coordinate system is at \e lat = \e lat0, \e lon = \e lon0,
\e h | |
| * = \e h0. The \e z axis is normal to the ellipsoid; the \e y axis point
s | | * = \e h0. The \e z axis is normal to the ellipsoid; the \e y axis point
s | |
| * due north. The plane \e z = - \e h0 is tangent to the ellipsoid. | | * due north. The plane \e z = - \e h0 is tangent to the ellipsoid. | |
| * | | * | |
| * The conversions all take place via geocentric coordinates using a | | * The conversions all take place via geocentric coordinates using a | |
| * Geocentric object (by default Geocentric::WGS84). | | * Geocentric object (by default Geocentric::WGS84). | |
| **********************************************************************/ | | **********************************************************************/ | |
| | | | |
|
| class LocalCartesian { | | class GEOGRAPHIC_EXPORT LocalCartesian { | |
| private: | | private: | |
| typedef Math::real real; | | typedef Math::real real; | |
|
| static const size_t dim = 3, dim2 = dim * dim; | | static const size_t dim_ = 3; | |
| | | static const size_t dim2_ = dim_ * dim_; | |
| const Geocentric _earth; | | const Geocentric _earth; | |
| real _lat0, _lon0, _h0; | | real _lat0, _lon0, _h0; | |
|
| real _x0, _y0, _z0, _r[dim2]; | | real _x0, _y0, _z0, _r[dim2_]; | |
| void IntForward(real lat, real lon, real h, real& x, real& y, real& z, | | void IntForward(real lat, real lon, real h, real& x, real& y, real& z, | |
|
| real M[dim2]) const throw(); | | real M[dim2_]) const throw(); | |
| void IntReverse(real x, real y, real z, real& lat, real& lon, real& h, | | void IntReverse(real x, real y, real z, real& lat, real& lon, real& h, | |
|
| real M[dim2]) const throw(); | | real M[dim2_]) const throw(); | |
| void MatrixMultiply(real M[dim2]) const throw(); | | void MatrixMultiply(real M[dim2_]) const throw(); | |
| public: | | public: | |
| | | | |
| /** | | /** | |
| * Constructor setting the origin. | | * Constructor setting the origin. | |
| * | | * | |
| * @param[in] lat0 latitude at origin (degrees). | | * @param[in] lat0 latitude at origin (degrees). | |
| * @param[in] lon0 longitude at origin (degrees). | | * @param[in] lon0 longitude at origin (degrees). | |
| * @param[in] h0 height above ellipsoid at origin (meters); default 0. | | * @param[in] h0 height above ellipsoid at origin (meters); default 0. | |
| * @param[in] earth Geocentric object for the transformation; default | | * @param[in] earth Geocentric object for the transformation; default | |
| * Geocentric::WGS84. | | * Geocentric::WGS84. | |
| | | | |
| skipping to change at line 121 | | skipping to change at line 122 | |
| * @param[out] M if the length of the vector is 9, fill with the rotati
on | | * @param[out] M if the length of the vector is 9, fill with the rotati
on | |
| * matrix in row-major order. | | * matrix in row-major order. | |
| * | | * | |
| * Pre-multiplying a unit vector in local cartesian coordinates at (lat
, | | * Pre-multiplying a unit vector in local cartesian coordinates at (lat
, | |
| * lon, h) by \e M transforms the vector to local cartesian coordinates
at | | * lon, h) by \e M transforms the vector to local cartesian coordinates
at | |
| * (lat0, lon0, h0). | | * (lat0, lon0, h0). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| void Forward(real lat, real lon, real h, real& x, real& y, real& z, | | void Forward(real lat, real lon, real h, real& x, real& y, real& z, | |
| std::vector<real>& M) | | std::vector<real>& M) | |
| const throw() { | | const throw() { | |
|
| real t[dim2]; | | real t[dim2_]; | |
| IntForward(lat, lon, h, x, y, z, t); | | IntForward(lat, lon, h, x, y, z, t); | |
|
| if (M.end() == M.begin() + dim2) | | if (M.end() == M.begin() + dim2_) | |
| copy(t, t + dim2, M.begin()); | | copy(t, t + dim2_, M.begin()); | |
| } | | } | |
| | | | |
| /** | | /** | |
| * Convert from local cartesian to geodetic coordinates. | | * Convert from local cartesian to geodetic coordinates. | |
| * | | * | |
| * @param[in] x local cartesian coordinate (meters). | | * @param[in] x local cartesian coordinate (meters). | |
| * @param[in] y local cartesian coordinate (meters). | | * @param[in] y local cartesian coordinate (meters). | |
| * @param[in] z local cartesian coordinate (meters). | | * @param[in] z local cartesian coordinate (meters). | |
| * @param[out] lat latitude of point (degrees). | | * @param[out] lat latitude of point (degrees). | |
| * @param[out] lon longitude of point (degrees). | | * @param[out] lon longitude of point (degrees). | |
| | | | |
| skipping to change at line 164 | | skipping to change at line 165 | |
| * @param[out] M if the length of the vector is 9, fill with the rotati
on | | * @param[out] M if the length of the vector is 9, fill with the rotati
on | |
| * matrix in row-major order. | | * matrix in row-major order. | |
| * | | * | |
| * Pre-multiplying a unit vector in local cartesian coordinates at (lat
0, | | * Pre-multiplying a unit vector in local cartesian coordinates at (lat
0, | |
| * lon0, h0) by the transpose of \e M transforms the vector to local | | * lon0, h0) by the transpose of \e M transforms the vector to local | |
| * cartesian coordinates at (lat, lon, h). | | * cartesian coordinates at (lat, lon, h). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| void Reverse(real x, real y, real z, real& lat, real& lon, real& h, | | void Reverse(real x, real y, real z, real& lat, real& lon, real& h, | |
| std::vector<real>& M) | | std::vector<real>& M) | |
| const throw() { | | const throw() { | |
|
| real t[dim2]; | | real t[dim2_]; | |
| IntReverse(x, y, z, lat, lon, h, t); | | IntReverse(x, y, z, lat, lon, h, t); | |
|
| if (M.end() == M.begin() + dim2) | | if (M.end() == M.begin() + dim2_) | |
| copy(t, t + dim2, M.begin()); | | copy(t, t + dim2_, M.begin()); | |
| } | | } | |
| | | | |
| /** \name Inspector functions | | /** \name Inspector functions | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| ///@{ | | ///@{ | |
| /** | | /** | |
| * @return latitude of the origin (degrees). | | * @return latitude of the origin (degrees). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| Math::real LatitudeOrigin() const throw() { return _lat0; } | | Math::real LatitudeOrigin() const throw() { return _lat0; } | |
| | | | |
| | | | |
| skipping to change at line 208 | | skipping to change at line 209 | |
| * constructor. A value of 0 is returned for a sphere (infinite inve
rse | | * constructor. A value of 0 is returned for a sphere (infinite inve
rse | |
| * flattening). | | * flattening). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| Math::real InverseFlattening() const throw() | | Math::real InverseFlattening() const throw() | |
| { return _earth.InverseFlattening(); } | | { return _earth.InverseFlattening(); } | |
| ///@} | | ///@} | |
| }; | | }; | |
| | | | |
| } // namespace GeographicLib | | } // namespace GeographicLib | |
| | | | |
|
| #endif | | #endif // GEOGRAPHICLIB_LOCALCARTESIAN_HPP | |
| | | | |
End of changes. 13 change blocks. |
| 16 lines changed or deleted | | 17 lines changed or added | |
|
| MGRS.hpp | | MGRS.hpp | |
| /** | | /** | |
| * \file MGRS.hpp | | * \file MGRS.hpp | |
| * \brief Header for GeographicLib::MGRS class | | * \brief Header for GeographicLib::MGRS class | |
| * | | * | |
|
| * Copyright (c) Charles Karney (2008, 2009, 2010) <charles@karney.com> | | * Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co
m> | |
| * and licensed under the LGPL. For more information, see | | * and licensed under the LGPL. For more information, see | |
| * http://geographiclib.sourceforge.net/ | | * http://geographiclib.sourceforge.net/ | |
| **********************************************************************/ | | **********************************************************************/ | |
| | | | |
| #if !defined(GEOGRAPHICLIB_MGRS_HPP) | | #if !defined(GEOGRAPHICLIB_MGRS_HPP) | |
|
| #define GEOGRAPHICLIB_MGRS_HPP "$Id: MGRS.hpp 6911 2010-12-09 23:13:55Z kar
ney $" | | #define GEOGRAPHICLIB_MGRS_HPP "$Id: c818eb3e77b330b44be1412ea54172a9458a67
0c $" | |
| | | | |
|
| #include "GeographicLib/Constants.hpp" | | | |
| #include "GeographicLib/UTMUPS.hpp" | | | |
| #include <sstream> | | #include <sstream> | |
|
| | | #include <GeographicLib/Constants.hpp> | |
| | | #include <GeographicLib/UTMUPS.hpp> | |
| | | | |
| | | #if defined(_MSC_VER) | |
| | | // Squelch warnings about dll vs string | |
| | | #pragma warning (push) | |
| | | #pragma warning (disable: 4251) | |
| | | #endif | |
| | | | |
| namespace GeographicLib { | | namespace GeographicLib { | |
| | | | |
| /** | | /** | |
| * \brief Convert between UTM/UPS and %MGRS | | * \brief Convert between UTM/UPS and %MGRS | |
| * | | * | |
| * MGRS is defined in Chapter 3 of | | * MGRS is defined in Chapter 3 of | |
| * - J. W. Hager, L. L. Fry, S. S. Jacks, D. R. Hill, | | * - J. W. Hager, L. L. Fry, S. S. Jacks, D. R. Hill, | |
| * <a href="http://earth-info.nga.mil/GandG/publications/tm8358.1/pdf/T
M8358_1.pdf"> | | * <a href="http://earth-info.nga.mil/GandG/publications/tm8358.1/pdf/T
M8358_1.pdf"> | |
| | | | |
| | | | |
| skipping to change at line 54 | | skipping to change at line 60 | |
| * class. | | * class. | |
| * | | * | |
| * The <a href="http://www.nga.mil">NGA</a> software package | | * The <a href="http://www.nga.mil">NGA</a> software package | |
| * <a href="http://earth-info.nga.mil/GandG/geotrans/index.html">geotrans
</a> | | * <a href="http://earth-info.nga.mil/GandG/geotrans/index.html">geotrans
</a> | |
| * also provides conversions to and from MGRS. Version 3.0 (and earlier) | | * also provides conversions to and from MGRS. Version 3.0 (and earlier) | |
| * suffers from some drawbacks: | | * suffers from some drawbacks: | |
| * - Inconsistent rules are used to determine the whether a particular MG
RS | | * - Inconsistent rules are used to determine the whether a particular MG
RS | |
| * coordinate is legal. A more systematic approach is taken here. | | * coordinate is legal. A more systematic approach is taken here. | |
| * - The underlying projections are not very accurately implemented. | | * - The underlying projections are not very accurately implemented. | |
| **********************************************************************/ | | **********************************************************************/ | |
|
| class MGRS { | | class GEOGRAPHIC_EXPORT MGRS { | |
| private: | | private: | |
| typedef Math::real real; | | typedef Math::real real; | |
|
| // The smallest length s.t., 1.0e7 - eps < 1.0e7 (approx 1.9 nm) | | // The smallest length s.t., 1.0e7 - eps_ < 1.0e7 (approx 1.9 nm) | |
| static const real eps; | | static const real eps_; | |
| // The smallest angle s.t., 90 - eps < 90 (approx 50e-12 arcsec) | | // The smallest angle s.t., 90 - eps_ < 90 (approx 50e-12 arcsec) | |
| static const real angeps; | | static const real angeps_; | |
| static const std::string hemispheres; | | static const std::string hemispheres_; | |
| static const std::string utmcols[3]; | | static const std::string utmcols_[3]; | |
| static const std::string utmrow; | | static const std::string utmrow_; | |
| static const std::string upscols[4]; | | static const std::string upscols_[4]; | |
| static const std::string upsrows[2]; | | static const std::string upsrows_[2]; | |
| static const std::string latband; | | static const std::string latband_; | |
| static const std::string upsband; | | static const std::string upsband_; | |
| static const std::string digits; | | static const std::string digits_; | |
| | | | |
| static const int mineasting[4]; | | static const int mineasting_[4]; | |
| static const int maxeasting[4]; | | static const int maxeasting_[4]; | |
| static const int minnorthing[4]; | | static const int minnorthing_[4]; | |
| static const int maxnorthing[4]; | | static const int maxnorthing_[4]; | |
| enum { | | enum { | |
|
| base = 10, | | base_ = 10, | |
| // Top-level tiles are 10^5 m = 100km on a side | | // Top-level tiles are 10^5 m = 100km on a side | |
|
| tilelevel = 5, | | tilelevel_ = 5, | |
| // Period of UTM row letters | | // Period of UTM row letters | |
|
| utmrowperiod = 20, | | utmrowperiod_ = 20, | |
| // Row letters are shifted by 5 for even zones | | // Row letters are shifted by 5 for even zones | |
|
| utmevenrowshift = 5, | | utmevenrowshift_ = 5, | |
| // Maximum precision is um | | // Maximum precision is um | |
|
| maxprec = 5 + 6, | | maxprec_ = 5 + 6, | |
| }; | | }; | |
| static void CheckCoords(bool utmp, bool& northp, real& x, real& y); | | static void CheckCoords(bool utmp, bool& northp, real& x, real& y); | |
| static int lookup(const std::string& s, char c) throw() { | | static int lookup(const std::string& s, char c) throw() { | |
| std::string::size_type r = s.find(toupper(c)); | | std::string::size_type r = s.find(toupper(c)); | |
| return r == std::string::npos ? -1 : int(r); | | return r == std::string::npos ? -1 : int(r); | |
| } | | } | |
| template<typename T> static std::string str(T x) { | | template<typename T> static std::string str(T x) { | |
| std::ostringstream s; s << x; return s.str(); | | std::ostringstream s; s << x; return s.str(); | |
| } | | } | |
| static int UTMRow(int iband, int icol, int irow) throw(); | | static int UTMRow(int iband, int icol, int irow) throw(); | |
| | | | |
| friend class UTMUPS; // UTMUPS::StandardZone calls LatitudeBand | | friend class UTMUPS; // UTMUPS::StandardZone calls LatitudeBand | |
| // Return latitude band number [-10, 10) for the give latitude (degrees
). | | // Return latitude band number [-10, 10) for the give latitude (degrees
). | |
| // The bands are reckoned in include their southern edges. | | // The bands are reckoned in include their southern edges. | |
| static int LatitudeBand(real lat) throw() { | | static int LatitudeBand(real lat) throw() { | |
| int ilat = int(std::floor(lat)); | | int ilat = int(std::floor(lat)); | |
| return (std::max)(-10, (std::min)(9, (ilat + 80)/8 - 10)); | | return (std::max)(-10, (std::min)(9, (ilat + 80)/8 - 10)); | |
| } | | } | |
| // UTMUPS access these enums | | // UTMUPS access these enums | |
| enum { | | enum { | |
|
| tile = 100000, // Size MGRS blocks | | tile_ = 100000, // Size MGRS blocks | |
| minutmcol = 1, | | minutmcol_ = 1, | |
| maxutmcol = 9, | | maxutmcol_ = 9, | |
| minutmSrow = 10, | | minutmSrow_ = 10, | |
| maxutmSrow = 100, // Also used for UTM S false northing | | maxutmSrow_ = 100, // Also used for UTM S false northing | |
| minutmNrow = 0, // Also used for UTM N false northing | | minutmNrow_ = 0, // Also used for UTM N false northing | |
| maxutmNrow = 95, | | maxutmNrow_ = 95, | |
| minupsSind = 8, // These 4 ind's apply to easting and north | | minupsSind_ = 8, // These 4 ind's apply to easting and nort | |
| ing | | hing | |
| maxupsSind = 32, | | maxupsSind_ = 32, | |
| minupsNind = 13, | | minupsNind_ = 13, | |
| maxupsNind = 27, | | maxupsNind_ = 27, | |
| upseasting = 20, // Also used for UPS false northing | | upseasting_ = 20, // Also used for UPS false northing | |
| utmeasting = 5, // UTM false easting | | utmeasting_ = 5, // UTM false easting | |
| // Difference between S hemisphere northing and N hemisphere northing | | // Difference between S hemisphere northing and N hemisphere northing | |
|
| utmNshift = (maxutmSrow - minutmNrow) * tile | | utmNshift_ = (maxutmSrow_ - minutmNrow_) * tile_ | |
| }; | | }; | |
| MGRS(); // Disable constructor | | MGRS(); // Disable constructor | |
| | | | |
| public: | | public: | |
| | | | |
| /** | | /** | |
| * Convert UTM or UPS coordinate to an MGRS coordinate. | | * Convert UTM or UPS coordinate to an MGRS coordinate. | |
| * | | * | |
| * @param[in] zone UTM zone (zero means UPS). | | * @param[in] zone UTM zone (zero means UPS). | |
| * @param[in] northp hemisphere (true means north, false means south). | | * @param[in] northp hemisphere (true means north, false means south). | |
| | | | |
| skipping to change at line 276 | | skipping to change at line 282 | |
| * | | * | |
| * (The WGS84 value is returned because the UTM and UPS projections are | | * (The WGS84 value is returned because the UTM and UPS projections are | |
| * based on this ellipsoid.) | | * based on this ellipsoid.) | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static Math::real InverseFlattening() throw() | | static Math::real InverseFlattening() throw() | |
| { return UTMUPS::InverseFlattening(); } | | { return UTMUPS::InverseFlattening(); } | |
| ///@} | | ///@} | |
| }; | | }; | |
| | | | |
| } // namespace GeographicLib | | } // namespace GeographicLib | |
|
| | | | |
| | | #if defined(_MSC_VER) | |
| | | #pragma warning (pop) | |
| #endif | | #endif | |
|
| | | | |
| | | #endif // GEOGRAPHICLIB_MGRS_HPP | |
| | | | |
End of changes. 15 change blocks. |
| 42 lines changed or deleted | | 51 lines changed or added | |
|
| OSGB.hpp | | OSGB.hpp | |
| /** | | /** | |
| * \file OSGB.hpp | | * \file OSGB.hpp | |
| * \brief Header for GeographicLib::OSGB class | | * \brief Header for GeographicLib::OSGB class | |
| * | | * | |
|
| * Copyright (c) Charles Karney (2010) <charles@karney.com> and licensed un | | * Copyright (c) Charles Karney (2010, 2011) <charles@karney.com> and licen | |
| der | | sed | |
| * the LGPL. For more information, see http://geographiclib.sourceforge.ne | | * under the LGPL. For more information, see | |
| t/ | | * http://geographiclib.sourceforge.net/ | |
| **********************************************************************/ | | **********************************************************************/ | |
| | | | |
| #if !defined(GEOGRAPHICLIB_OSGB_HPP) | | #if !defined(GEOGRAPHICLIB_OSGB_HPP) | |
|
| #define GEOGRAPHICLIB_OSGB_HPP "$Id: OSGB.hpp 6911 2010-12-09 23:13:55Z kar
ney $" | | #define GEOGRAPHICLIB_OSGB_HPP "$Id: 35c0f9c3cf61ce5fae409218463cce0c5558a7
e0 $" | |
| | | | |
|
| #include "GeographicLib/Constants.hpp" | | | |
| #include "GeographicLib/TransverseMercator.hpp" | | | |
| #include <string> | | #include <string> | |
| #include <sstream> | | #include <sstream> | |
|
| | | #include <GeographicLib/Constants.hpp> | |
| | | #include <GeographicLib/TransverseMercator.hpp> | |
| | | | |
| | | #if defined(_MSC_VER) | |
| | | // Squelch warnings about dll vs string | |
| | | #pragma warning (push) | |
| | | #pragma warning (disable: 4251) | |
| | | #endif | |
| | | | |
| namespace GeographicLib { | | namespace GeographicLib { | |
| | | | |
| /** | | /** | |
| * \brief Ordnance Survey grid system for Great Britain | | * \brief Ordnance Survey grid system for Great Britain | |
| * | | * | |
| * The class implements the coordinate system used by the Ordnance Survey
for | | * The class implements the coordinate system used by the Ordnance Survey
for | |
| * maps of Great Britain and conversions to the grid reference system. | | * maps of Great Britain and conversions to the grid reference system. | |
| * | | * | |
| * See | | * See | |
| * - <a href="http://www.ordnancesurvey.co.uk/oswebsite/gps/docs/A_Guide_
to_Coordinate_Systems_in_Great_Britain.pdf"> | | * - <a href="http://www.ordnancesurvey.co.uk/oswebsite/gps/docs/A_Guide_
to_Coordinate_Systems_in_Great_Britain.pdf"> | |
| * A guide to coordinate systems in Great Britain</a> | | * A guide to coordinate systems in Great Britain</a> | |
| * - <a href="http://www.ordnancesurvey.co.uk/oswebsite/gps/information/c
oordinatesystemsinfo/guidetonationalgrid/page1.html"> | | * - <a href="http://www.ordnancesurvey.co.uk/oswebsite/gps/information/c
oordinatesystemsinfo/guidetonationalgrid/page1.html"> | |
| * Guide to National Grid</a> | | * Guide to National Grid</a> | |
| * | | * | |
| * \b WARNING: the latitudes and longitudes for the Ordnance Survey grid | | * \b WARNING: the latitudes and longitudes for the Ordnance Survey grid | |
| * system do not use the WGS84 datum. Do not use the values returned by
this | | * system do not use the WGS84 datum. Do not use the values returned by
this | |
| * class in the UTMUPS, MGRS, or Geoid classes without first converting t
he | | * class in the UTMUPS, MGRS, or Geoid classes without first converting t
he | |
| * datum (and vice versa). | | * datum (and vice versa). | |
| **********************************************************************/ | | **********************************************************************/ | |
|
| class OSGB { | | class GEOGRAPHIC_EXPORT OSGB { | |
| private: | | private: | |
| typedef Math::real real; | | typedef Math::real real; | |
|
| static const std::string letters; | | static const std::string letters_; | |
| static const std::string digits; | | static const std::string digits_; | |
| static const TransverseMercator OSGBTM; | | static const TransverseMercator OSGBTM_; | |
| static const real northoffset; | | static const real northoffset_; | |
| enum { | | enum { | |
|
| base = 10, | | base_ = 10, | |
| tile = 100000, | | tile_ = 100000, | |
| tilelevel = 5, | | tilelevel_ = 5, | |
| tilegrid = 5, | | tilegrid_ = 5, | |
| tileoffx = 2 * tilegrid, | | tileoffx_ = 2 * tilegrid_, | |
| tileoffy = 1 * tilegrid, | | tileoffy_ = 1 * tilegrid_, | |
| minx = - tileoffx * tile, | | minx_ = - tileoffx_ * tile_, | |
| miny = - tileoffy * tile, | | miny_ = - tileoffy_ * tile_, | |
| maxx = (tilegrid*tilegrid - tileoffx) * tile, | | maxx_ = (tilegrid_*tilegrid_ - tileoffx_) * tile_, | |
| maxy = (tilegrid*tilegrid - tileoffy) * tile, | | maxy_ = (tilegrid_*tilegrid_ - tileoffy_) * tile_, | |
| // Maximum precision is um | | // Maximum precision is um | |
|
| maxprec = 5 + 6, | | maxprec_ = 5 + 6, | |
| }; | | }; | |
| static real computenorthoffset() throw(); | | static real computenorthoffset() throw(); | |
| static void CheckCoords(real x, real y); | | static void CheckCoords(real x, real y); | |
| static int lookup(const std::string& s, char c) throw() { | | static int lookup(const std::string& s, char c) throw() { | |
| std::string::size_type r = s.find(toupper(c)); | | std::string::size_type r = s.find(toupper(c)); | |
| return r == std::string::npos ? -1 : int(r); | | return r == std::string::npos ? -1 : int(r); | |
| } | | } | |
| template<typename T> static std::string str(T x) { | | template<typename T> static std::string str(T x) { | |
| std::ostringstream s; s << x; return s.str(); | | std::ostringstream s; s << x; return s.str(); | |
| } | | } | |
| | | | |
| skipping to change at line 85 | | skipping to change at line 92 | |
| * @param[out] x easting of point (meters). | | * @param[out] x easting of point (meters). | |
| * @param[out] y northing of point (meters). | | * @param[out] y northing of point (meters). | |
| * @param[out] gamma meridian convergence at point (degrees). | | * @param[out] gamma meridian convergence at point (degrees). | |
| * @param[out] k scale of projection at point. | | * @param[out] k scale of projection at point. | |
| * | | * | |
| * \e lat should be in the range [-90, 90]; \e lon and \e lon0 should b
e in | | * \e lat should be in the range [-90, 90]; \e lon and \e lon0 should b
e in | |
| * the range [-180, 360]. | | * the range [-180, 360]. | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static void Forward(real lat, real lon, | | static void Forward(real lat, real lon, | |
| real& x, real& y, real& gamma, real& k) throw() { | | real& x, real& y, real& gamma, real& k) throw() { | |
|
| OSGBTM.Forward(OriginLongitude(), lat, lon, x, y, gamma, k); | | OSGBTM_.Forward(OriginLongitude(), lat, lon, x, y, gamma, k); | |
| x += FalseEasting(); | | x += FalseEasting(); | |
|
| y += northoffset; | | y += northoffset_; | |
| } | | } | |
| | | | |
| /** | | /** | |
| * Reverse projection, from OSGB coordinates to geographic. | | * Reverse projection, from OSGB coordinates to geographic. | |
| * | | * | |
| * @param[in] x easting of point (meters). | | * @param[in] x easting of point (meters). | |
| * @param[in] y northing of point (meters). | | * @param[in] y northing of point (meters). | |
| * @param[out] lat latitude of point (degrees). | | * @param[out] lat latitude of point (degrees). | |
| * @param[out] lon longitude of point (degrees). | | * @param[out] lon longitude of point (degrees). | |
| * @param[out] gamma meridian convergence at point (degrees). | | * @param[out] gamma meridian convergence at point (degrees). | |
| * @param[out] k scale of projection at point. | | * @param[out] k scale of projection at point. | |
| * | | * | |
| * The value of \e lon returned is in the range [-180, 180). | | * The value of \e lon returned is in the range [-180, 180). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| | | | |
| static void Reverse(real x, real y, | | static void Reverse(real x, real y, | |
| real& lat, real& lon, real& gamma, real& k) throw()
{ | | real& lat, real& lon, real& gamma, real& k) throw()
{ | |
| x -= FalseEasting(); | | x -= FalseEasting(); | |
|
| y -= northoffset; | | y -= northoffset_; | |
| OSGBTM.Reverse(OriginLongitude(), x, y, lat, lon, gamma, k); | | OSGBTM_.Reverse(OriginLongitude(), x, y, lat, lon, gamma, k); | |
| } | | } | |
| | | | |
| /** | | /** | |
| * OSGB::Forward without returning the convergence and scale. | | * OSGB::Forward without returning the convergence and scale. | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static void Forward(real lat, real lon, real& x, real& y) throw() { | | static void Forward(real lat, real lon, real& x, real& y) throw() { | |
| real gamma, k; | | real gamma, k; | |
| Forward(lat, lon, x, y, gamma, k); | | Forward(lat, lon, x, y, gamma, k); | |
| } | | } | |
| | | | |
| | | | |
| skipping to change at line 162 | | skipping to change at line 169 | |
| /** | | /** | |
| * Convert OSGB coordinates to a grid reference. | | * Convert OSGB coordinates to a grid reference. | |
| * | | * | |
| * @param[in] gridref National Grid reference. | | * @param[in] gridref National Grid reference. | |
| * @param[out] x easting of point (meters). | | * @param[out] x easting of point (meters). | |
| * @param[out] y northing of point (meters). | | * @param[out] y northing of point (meters). | |
| * @param[out] prec precision relative to 100 km. | | * @param[out] prec precision relative to 100 km. | |
| * @param[in] centerp if true (default), return center of the grid squa
re, | | * @param[in] centerp if true (default), return center of the grid squa
re, | |
| * else return SW (lower left) corner. | | * else return SW (lower left) corner. | |
| * | | * | |
|
| * The grid reference must be of the form: two letters (not including I | | * The grid reference must be of the form: two letters_ (not including | |
| ) | | I) | |
| * followed by an even number of digits (up to 22). | | * followed by an even number of digits_ (up to 22). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static void GridReference(const std::string& gridref, | | static void GridReference(const std::string& gridref, | |
| real& x, real& y, int& prec, | | real& x, real& y, int& prec, | |
| bool centerp = true); | | bool centerp = true); | |
| | | | |
| /** \name Inspector functions | | /** \name Inspector functions | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| ///@{ | | ///@{ | |
| /** | | /** | |
| * @return \e a the equatorial radius of the Airy 1830 ellipsoid (meter
s). | | * @return \e a the equatorial radius of the Airy 1830 ellipsoid (meter
s). | |
| * | | * | |
| * This is 20923713 ft converted to meters using the rule 1 ft = | | * This is 20923713 ft converted to meters using the rule 1 ft = | |
| * 10^(9.48401603-10) m. (The Airy 1830 value is returned because the
OSGB | | * 10^(9.48401603-10) m. (The Airy 1830 value is returned because the
OSGB | |
| * projection is based on this ellipsoid.) | | * projection is based on this ellipsoid.) | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static Math::real MajorRadius() throw() | | static Math::real MajorRadius() throw() | |
|
| { return real(6377563.396032066440602120008397385037L); } | | // result is about 6377563.3960320664406 m | |
| | | { return real(20923713) * std::pow(real(10), real(0.48401603L) - 1); } | |
| | | | |
| /** | | /** | |
| * @return \e r the inverse flattening of the Airy 1830 ellipsoid. | | * @return \e r the inverse flattening of the Airy 1830 ellipsoid. | |
| * | | * | |
| * For the Airy 1830 ellipsoid, \e a = 20923713 ft and \e b = 20853810
ft; | | * For the Airy 1830 ellipsoid, \e a = 20923713 ft and \e b = 20853810
ft; | |
| * thus the inverse flattening = 20923713/(20923713 - 20853810) = | | * thus the inverse flattening = 20923713/(20923713 - 20853810) = | |
|
| * 2324857/7767 = 299.32496459... (The Airy 1830 value is returned bec | | * 299.32496459... (The Airy 1830 value is returned because the OSGB | |
| ause | | * projection is based on this ellipsoid.) | |
| * the OSGB projection is based on this ellipsoid.) | | | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static Math::real InverseFlattening() throw() | | static Math::real InverseFlattening() throw() | |
|
| { return real(2324857)/real(7767); } | | { return real(20923713) / real(20923713 - 20853810); } | |
| | | | |
| /** | | /** | |
| * @return \e k0 central scale for the OSGB projection (0.9996012717). | | * @return \e k0 central scale for the OSGB projection (0.9996012717). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static Math::real CentralScale() throw() | | static Math::real CentralScale() throw() | |
| { return real(0.9996012717L); } | | { return real(0.9996012717L); } | |
| | | | |
| /** | | /** | |
| * @return latitude of the origin for the OSGB projection (49 degrees). | | * @return latitude of the origin for the OSGB projection (49 degrees). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| | | | |
| skipping to change at line 223 | | skipping to change at line 231 | |
| | | | |
| /** | | /** | |
| * @return false easting the OSGB projection (400000 meters). | | * @return false easting the OSGB projection (400000 meters). | |
| **********************************************************************
/ | | **********************************************************************
/ | |
| static Math::real FalseEasting() throw() { return real(400000); } | | static Math::real FalseEasting() throw() { return real(400000); } | |
| ///@} | | ///@} | |
| | | | |
| }; | | }; | |
| | | | |
| } // namespace GeographicLib | | } // namespace GeographicLib | |
|
| | | | |
| | | #if defined(_MSC_VER) | |
| | | #pragma warning (pop) | |
| | | #endif | |
| | | | |
| #endif | | #endif | |
| | | | |
End of changes. 16 change blocks. |
| 35 lines changed or deleted | | 46 lines changed or added | |
|