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


 AzimuthalEquidistant.hpp   AzimuthalEquidistant.hpp 
/** /**
* \file AzimuthalEquidistant.hpp * \file AzimuthalEquidistant.hpp
* \brief Header for GeographicLib::AzimuthalEquidistant class * \brief Header for GeographicLib::AzimuthalEquidistant class
* *
* Copyright (c) Charles Karney (2009, 2010) <charles@karney.com> * Copyright (c) Charles Karney (2009, 2010, 2011) <charles@karney.com> and
* and licensed under the LGPL. For more information, see * licensed under the LGPL. For more information, see
* http://geographiclib.sourceforge.net/ * http://geographiclib.sourceforge.net/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_AZIMUTHALEQUIDISTANT_HPP) #if !defined(GEOGRAPHICLIB_AZIMUTHALEQUIDISTANT_HPP)
#define GEOGRAPHICLIB_AZIMUTHALEQUIDISTANT_HPP "$Id: AzimuthalEquidistant.h pp 6867 2010-09-11 13:04:26Z karney $" #define GEOGRAPHICLIB_AZIMUTHALEQUIDISTANT_HPP "$Id: 526db2f041dd6d71d05acd 60c7429ad2835ab79e $"
#include "GeographicLib/Geodesic.hpp" #include <GeographicLib/Geodesic.hpp>
#include "GeographicLib/Constants.hpp" #include <GeographicLib/Constants.hpp>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief Azimuthal Equidistant Projection. * \brief Azimuthal Equidistant Projection.
* *
* Azimuthal equidistant projection centered at an arbitrary position on the * Azimuthal equidistant projection centered at an arbitrary position on the
* ellipsoid. For a point in projected space (\e x, \e y), the geodesic * ellipsoid. For a point in projected space (\e x, \e y), the geodesic
* distance from the center position is hypot(\e x, \e y) and the azimuth of * distance from the center position is hypot(\e x, \e y) and the azimuth of
* the geodesic from the center point is atan2(\e x, \e y). The Forward and * the geodesic from the center point is atan2(\e x, \e y). The Forward and
* Reverse methods also return the azimuth \e azi of the geodesic at (\e x, * Reverse methods also return the azimuth \e azi of the geodesic at (\e x,
* \e y) and reciprocal scale \e rk in the azimuthal direction which, * \e y) and reciprocal scale \e rk in the azimuthal direction which,
* together with the basic properties of the projection, serve to specify * together with the basic properties of the projection, serve to specify
* completely the local affine transformation between geographic and * completely the local affine transformation between geographic and
* projected coordinates. * projected coordinates.
* *
* The conversions all take place using a Geodesic object (by default * The conversions all take place using a Geodesic object (by default
* Geodesic::WGS84). For more information on geodesics see \ref geodesic . * Geodesic::WGS84). For more information on geodesics see \ref geodesic .
**********************************************************************/ **********************************************************************/
class AzimuthalEquidistant { class GEOGRAPHIC_EXPORT AzimuthalEquidistant {
private: private:
typedef Math::real real; typedef Math::real real;
const Geodesic _earth; const Geodesic _earth;
static const real eps; static const real eps_;
public: public:
/** /**
* Constructor for AzimuthalEquidistant. * Constructor for AzimuthalEquidistant.
* *
* @param[in] earth the Geodesic object to use for geodesic calculation s. * @param[in] earth the Geodesic object to use for geodesic calculation s.
* By default this uses the WGS84 ellipsoid. * By default this uses the WGS84 ellipsoid.
********************************************************************** / ********************************************************************** /
explicit AzimuthalEquidistant(const Geodesic& earth = Geodesic::WGS84) explicit AzimuthalEquidistant(const Geodesic& earth = Geodesic::WGS84)
throw() : _earth(earth) {} throw() : _earth(earth) {}
skipping to change at line 135 skipping to change at line 135
* value inherited from the Geodesic object used in the constructor. A * value inherited from the Geodesic object used in the constructor. A
* value of 0 is returned for a sphere (infinite inverse flattening). * value of 0 is returned for a sphere (infinite inverse flattening).
********************************************************************** / ********************************************************************** /
Math::real InverseFlattening() const throw() Math::real InverseFlattening() const throw()
{ return _earth.InverseFlattening(); } { return _earth.InverseFlattening(); }
///@} ///@}
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif #endif // GEOGRAPHICLIB_AZIMUTHALEQUIDISTANT_HPP
 End of changes. 6 change blocks. 
7 lines changed or deleted 7 lines changed or added


 CassiniSoldner.hpp   CassiniSoldner.hpp 
/** /**
* \file CassiniSoldner.hpp * \file CassiniSoldner.hpp
* \brief Header for GeographicLib::CassiniSoldner class * \brief Header for GeographicLib::CassiniSoldner class
* *
* Copyright (c) Charles Karney (2009, 2010) <charles@karney.com> * Copyright (c) Charles Karney (2009, 2010, 2011) <charles@karney.com> and
* and licensed under the LGPL. For more information, see * licensed under the LGPL. For more information, see
* http://geographiclib.sourceforge.net/ * http://geographiclib.sourceforge.net/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_CASSINISOLDNER_HPP) #if !defined(GEOGRAPHICLIB_CASSINISOLDNER_HPP)
#define GEOGRAPHICLIB_CASSINISOLDNER_HPP "$Id: CassiniSoldner.hpp 6867 2010 -09-11 13:04:26Z karney $" #define GEOGRAPHICLIB_CASSINISOLDNER_HPP "$Id: 7d36c41fb37aafe5f0f467755684 85be6f2d6502 $"
#include "GeographicLib/Geodesic.hpp" #include <GeographicLib/Geodesic.hpp>
#include "GeographicLib/GeodesicLine.hpp" #include <GeographicLib/GeodesicLine.hpp>
#include "GeographicLib/Constants.hpp" #include <GeographicLib/Constants.hpp>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief Cassini-Soldner Projection. * \brief Cassini-Soldner Projection.
* *
* Cassini-Soldner projection centered at an arbitrary position, \e lat0, \e * Cassini-Soldner projection centered at an arbitrary position, \e lat0, \e
* lon0, on the ellipsoid. This projection is a transverse cylindrical * lon0, on the ellipsoid. This projection is a transverse cylindrical
* equidistant projection. The projection from (\e lat, \e lon) to easti ng * equidistant projection. The projection from (\e lat, \e lon) to easti ng
* and northing (\e x, \e y) is defined by geodesics as follows. Go nort h * and northing (\e x, \e y) is defined by geodesics as follows. Go nort h
skipping to change at line 62 skipping to change at line 62
* Geodesic::WGS84). For more information on geodesics see \ref geodesic . * Geodesic::WGS84). For more information on geodesics see \ref geodesic .
* The determination of (\e lat1, \e lon1) in the forward projection is b y * The determination of (\e lat1, \e lon1) in the forward projection is b y
* solving the inverse geodesic problem for (\e lat, \e lon) and its twin * solving the inverse geodesic problem for (\e lat, \e lon) and its twin
* obtained by reflection in the meridional plane. The scale is found by * obtained by reflection in the meridional plane. The scale is found by
* determining where two neighboring geodesics intersecting the central * determining where two neighboring geodesics intersecting the central
* meridan at \e lat1 and \e lat1 + \e dlat1 intersect and taking the rat io * meridan at \e lat1 and \e lat1 + \e dlat1 intersect and taking the rat io
* of the reduced lengths for the two geodesics between that point and, * of the reduced lengths for the two geodesics between that point and,
* respectively, (\e lat1, \e lon1) and (\e lat, \e lon). * respectively, (\e lat1, \e lon1) and (\e lat, \e lon).
**********************************************************************/ **********************************************************************/
class CassiniSoldner { class GEOGRAPHIC_EXPORT CassiniSoldner {
private: private:
typedef Math::real real; typedef Math::real real;
const Geodesic _earth; const Geodesic _earth;
GeodesicLine _meridian; GeodesicLine _meridian;
real _sbet0, _cbet0; real _sbet0, _cbet0;
static const real eps1, eps2; static const real eps1_;
static const unsigned maxit = 10; static const real eps2_;
static const unsigned maxit_ = 10;
static inline real sq(real x) throw() { return x * x; }
// The following private helper functions are copied from Geodesic. // The following private helper functions are copied from Geodesic.
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
// is about 1000 times more resolution than we get with angles around 90 // is about 1000 times more resolution than we get with angles around 90
// degrees.) We use this to avoid having to deal with near singular // degrees.) We use this to avoid having to deal with near singular
skipping to change at line 228 skipping to change at line 228
* value inherited from the Geodesic object used in the constructor. A * value inherited from the Geodesic object used in the constructor. A
* value of 0 is returned for a sphere (infinite inverse flattening). * value of 0 is returned for a sphere (infinite inverse flattening).
********************************************************************** / ********************************************************************** /
Math::real InverseFlattening() const throw() Math::real InverseFlattening() const throw()
{ return _earth.InverseFlattening(); } { return _earth.InverseFlattening(); }
///@} ///@}
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif #endif // GEOGRAPHICLIB_CASSINISOLDNER_HPP
 End of changes. 7 change blocks. 
10 lines changed or deleted 10 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


 EllipticFunction.hpp   EllipticFunction.hpp 
/** /**
* \file EllipticFunction.hpp * \file EllipticFunction.hpp
* \brief Header for GeographicLib::EllipticFunction class * \brief Header for GeographicLib::EllipticFunction class
* *
* Copyright (c) Charles Karney (2008, 2009) <charles@karney.com> * Copyright (c) Charles Karney (2008, 2009, 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_ELLIPTICFUNCTION_HPP) #if !defined(GEOGRAPHICLIB_ELLIPTICFUNCTION_HPP)
#define GEOGRAPHICLIB_ELLIPTICFUNCTION_HPP "$Id: EllipticFunction.hpp 6866 2010-09-11 02:15:29Z karney $" #define GEOGRAPHICLIB_ELLIPTICFUNCTION_HPP "$Id: 8e5a2b0b4722672f75bb33acb8 233a78f00d8839 $"
#include "GeographicLib/Constants.hpp" #include <GeographicLib/Constants.hpp>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief Elliptic functions needed for TransverseMercatorExact * \brief Elliptic functions needed for TransverseMercatorExact
* *
* This provides the subset of elliptic functions needed for * This provides the subset of elliptic functions needed for
* TransverseMercatorExact. For a given ellipsoid, only parameters \e * TransverseMercatorExact. For a given ellipsoid, only parameters \e
* e<sup>2</sup> and 1 - \e e<sup>2</sup> are needed. This class taken t he * e<sup>2</sup> and 1 - \e e<sup>2</sup> are needed. This class taken t he
* parameter as a constructor parameters and caches the values of the * parameter as a constructor parameters and caches the values of the
skipping to change at line 42 skipping to change at line 42
* . * .
* The computation of the Jacobi elliptic functions uses the algorithm gi ven * The computation of the Jacobi elliptic functions uses the algorithm gi ven
* in * in
* - R. Bulirsch, * - R. Bulirsch,
* <a href="http://dx.doi.org/10.1007/BF01397975"> Numerical Calculatio n of * <a href="http://dx.doi.org/10.1007/BF01397975"> Numerical Calculatio n of
* Elliptic Integrals and Elliptic Functions</a>, Numericshe Mathematik 7, * Elliptic Integrals and Elliptic Functions</a>, Numericshe Mathematik 7,
* 78&ndash;90 (1965). * 78&ndash;90 (1965).
* . * .
* The notation follows Abramowitz and Stegun, Chapters 16 and 17. * The notation follows Abramowitz and Stegun, Chapters 16 and 17.
**********************************************************************/ **********************************************************************/
class EllipticFunction { class GEOGRAPHIC_EXPORT EllipticFunction {
private: private:
typedef Math::real real; typedef Math::real real;
static const real tol, tolRF, tolRD, tolRG0, tolJAC, tolJAC1; static const real tol_;
enum { num = 10 }; // Max depth required for sncndn. Probably 5 is eno static const real tolRF_;
ugh. static const real tolRD_;
static const real tolRG0_;
static const real tolJAC_;
static const real tolJAC1_;
enum { num_ = 10 }; // Max depth required for sncndn. Probably 5 is en
ough.
static real RF(real x, real y, real z) throw(); static real RF(real x, real y, real z) throw();
static real RD(real x, real y, real z) throw(); static real RD(real x, real y, real z) throw();
static real RG0(real x, real y) throw(); static real RG0(real x, real y) throw();
const real _m, _m1; const real _m, _m1;
mutable bool _init; mutable bool _init;
mutable real _kc, _ec, _kec; mutable real _kc, _ec, _kec;
bool Init() const throw(); bool Init() const throw();
public: public:
/** /**
skipping to change at line 126 skipping to change at line 131
* *
* Instead of specifying the ampltiude \e phi, we provide \e sn = sin(\ e * Instead of specifying the ampltiude \e phi, we provide \e sn = sin(\ e
* phi), \e cn = cos(\e phi), \e dn = sqrt(1 - \e m sin<sup>2</sup>(\e * phi), \e cn = cos(\e phi), \e dn = sqrt(1 - \e m sin<sup>2</sup>(\e
* phi)). * phi)).
********************************************************************** / ********************************************************************** /
Math::real E(real sn, real cn, real dn) const throw(); Math::real E(real sn, real cn, real dn) const throw();
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif #endif // GEOGRAPHICLIB_ELLIPTICFUNCTION_HPP
 End of changes. 6 change blocks. 
7 lines changed or deleted 12 lines changed or added


 GeoCoords.hpp   GeoCoords.hpp 
/** /**
* \file GeoCoords.hpp * \file GeoCoords.hpp
* \brief Header for GeographicLib::GeoCoords class * \brief Header for GeographicLib::GeoCoords class
* *
* Copyright (c) Charles Karney (2008, 2009, 2010) <charles@karney.com> * Copyright (c) Charles Karney (2008, 2009, 2010) <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/
**********************************************************************/ **********************************************************************/
#ifndef GEOGRAPHICLIB_GEOCOORDS_HPP #ifndef GEOGRAPHICLIB_GEOCOORDS_HPP
#define GEOGRAPHICLIB_GEOCOORDS_HPP "$Id: GeoCoords.hpp 6876 2010-10-18 13: 47:55Z karney $" #define GEOGRAPHICLIB_GEOCOORDS_HPP "$Id: 43489a5a6c17b341fd293b326d0340e71 9df01e9 $"
#include "GeographicLib/UTMUPS.hpp" #include <GeographicLib/UTMUPS.hpp>
#include "GeographicLib/Constants.hpp" #include <GeographicLib/Constants.hpp>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief Conversion between geographic coordinates * \brief Conversion between geographic coordinates
* *
* This class stores a geographic position which may be set via the * This class stores a geographic position which may be set via the
* constructors or Reset via * constructors or Reset via
* - latitude and longitude * - latitude and longitude
* - UTM or UPS coordinates * - UTM or UPS coordinates
skipping to change at line 43 skipping to change at line 43
* zone. A method SetAltZone is provided to set the alternate UPS/UTM zo ne. * zone. A method SetAltZone is provided to set the alternate UPS/UTM zo ne.
* *
* Methods are provided to return the geographic coordinates, the input U TM * Methods are provided to return the geographic coordinates, the input U TM
* or UPS coordinates (and associated meridian convergence and scale), or * or UPS coordinates (and associated meridian convergence and scale), or
* alternate UTM or UPS coordinates (and their associated meridian * alternate UTM or UPS coordinates (and their associated meridian
* convergence and scale). * convergence and scale).
* *
* Once the input string has been parsed, you can print the result out in any * Once the input string has been parsed, you can print the result out in any
* of the formats, decimal degrees, degrees minutes seconds, MGRS, UTM/UP S. * of the formats, decimal degrees, degrees minutes seconds, MGRS, UTM/UP S.
**********************************************************************/ **********************************************************************/
class GeoCoords { class GEOGRAPHIC_EXPORT GeoCoords {
private: private:
typedef Math::real real; typedef Math::real real;
real _lat, _long, _easting, _northing, _gamma, _k; real _lat, _long, _easting, _northing, _gamma, _k;
bool _northp; bool _northp;
int _zone; // See UTMUPS::zonespec int _zone; // See UTMUPS::zonespec
mutable real _alt_easting, _alt_northing, _alt_gamma, _alt_k; mutable real _alt_easting, _alt_northing, _alt_gamma, _alt_k;
mutable int _alt_zone; mutable int _alt_zone;
void CopyToAlt() const throw() { void CopyToAlt() const throw() {
_alt_easting = _easting; _alt_easting = _easting;
skipping to change at line 432 skipping to change at line 432
* *
* (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.)
********************************************************************** / ********************************************************************** /
Math::real InverseFlattening() const throw() Math::real InverseFlattening() const throw()
{ return UTMUPS::InverseFlattening(); } { return UTMUPS::InverseFlattening(); }
///@} ///@}
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif // GEOGRAPHICLIB_GEOCOORDS_HPP #endif // GEOGRAPHICLIB_GEOCOORDS_HPP
 End of changes. 4 change blocks. 
4 lines changed or deleted 5 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&ndash;454 (2002). * J. Geodesy 76, 451&ndash;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


 Gnomonic.hpp   Gnomonic.hpp 
/** /**
* \file Gnomonic.hpp * \file Gnomonic.hpp
* \brief Header for GeographicLib::Gnomonic class * \brief Header for GeographicLib::Gnomonic 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_GNOMONIC_HPP) #if !defined(GEOGRAPHICLIB_GNOMONIC_HPP)
#define GEOGRAPHICLIB_GNOMONIC_HPP "$Id: Gnomonic.hpp 6968 2011-02-19 15:58 :56Z karney $" #define GEOGRAPHICLIB_GNOMONIC_HPP "$Id: bf3b7f9f87a997fd135064f9b3233ef267 eb8394 $"
#include "GeographicLib/Geodesic.hpp" #include <GeographicLib/Geodesic.hpp>
#include "GeographicLib/GeodesicLine.hpp" #include <GeographicLib/GeodesicLine.hpp>
#include "GeographicLib/Constants.hpp" #include <GeographicLib/Constants.hpp>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief %Gnomonic Projection. * \brief %Gnomonic Projection.
* *
* %Gnomonic projection centered at an arbitrary position \e C on the * %Gnomonic projection centered at an arbitrary position \e C on the
* ellipsoid. This projection is derived in Section 13 of * ellipsoid. This projection is derived in Section 13 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>.
* . * .
* The projection of \e P is defined as follows: compute the * The projection of \e P is defined as follows: compute the
* geodesic line from \e C to \e P; compute the reduced length \e m12, * geodesic line from \e C to \e P; compute the reduced length \e m12,
* geodesic scale \e M12, and \e rho = \e m12/\e M12; finally \e x = \e r ho * geodesic scale \e M12, and \e rho = \e m12/\e M12; finally \e x = \e r ho
* sin \e azi1; \e y = \e rho cos \e azi1, where \e azi1 is the azimuth o f * sin \e azi1; \e y = \e rho cos \e azi1, where \e azi1 is the azimuth o f
* the geodesic at \e C. The Gnomonic::Forward and Gnomonic::Reverse met hods * the geodesic at \e C. The Gnomonic::Forward and Gnomonic::Reverse met hods
* also return the azimuth \e azi of the geodesic at \e P and reciprocal * also return the azimuth \e azi of the geodesic at \e P and reciprocal
* scale \e rk in the azimuthal direction. The scale in the radial direc tion * scale \e rk in the azimuthal direction. The scale in the radial direc tion
* if 1/\e rk<sup>2</sup>. * if 1/\e rk<sup>2</sup>.
* *
skipping to change at line 90 skipping to change at line 90
* generalizes the spherical great circle to a normal section. This wa s * generalizes the spherical great circle to a normal section. This wa s
* proposed by I. G. Letoval'tsev, Generalization of the %Gnomonic * proposed by I. G. Letoval'tsev, Generalization of the %Gnomonic
* Projection for a Spheroid and the Principal Geodetic Problems Involv ed * Projection for a Spheroid and the Principal Geodetic Problems Involv ed
* in the Alignment of Surface Routes, Geodesy and Aerophotography (5), * in the Alignment of Surface Routes, Geodesy and Aerophotography (5),
* 271-274 (1963). * 271-274 (1963).
* - The projection given here. This causes geodesics close to the cente r * - The projection given here. This causes geodesics close to the cente r
* point to appear as straight lines in the projection; i.e., it * point to appear as straight lines in the projection; i.e., it
* generalizes the spherical great circle to a geodesic. * generalizes the spherical great circle to a geodesic.
**********************************************************************/ **********************************************************************/
class Gnomonic { class GEOGRAPHIC_EXPORT Gnomonic {
private: private:
typedef Math::real real; typedef Math::real real;
const Geodesic _earth; const Geodesic _earth;
real _a, _f; real _a, _f;
static const real eps0, eps; static const real eps0_;
static const int numit = 5; static const real eps_;
static const int numit_ = 5;
public: public:
/** /**
* Constructor for Gnomonic. * Constructor for Gnomonic.
* *
* @param[in] earth the Geodesic object to use for geodesic calculation s. * @param[in] earth the Geodesic object to use for geodesic calculation s.
* By default this uses the WGS84 ellipsoid. * By default this uses the WGS84 ellipsoid.
********************************************************************** / ********************************************************************** /
explicit Gnomonic(const Geodesic& earth = Geodesic::WGS84) explicit Gnomonic(const Geodesic& earth = Geodesic::WGS84)
throw() throw()
skipping to change at line 203 skipping to change at line 204
* value inherited from the Geodesic object used in the constructor. A * value inherited from the Geodesic object used in the constructor. A
* value of 0 is returned for a sphere (infinite inverse flattening). * value of 0 is returned for a sphere (infinite inverse flattening).
********************************************************************** / ********************************************************************** /
Math::real InverseFlattening() const throw() Math::real InverseFlattening() const throw()
{ return _earth.InverseFlattening(); } { return _earth.InverseFlattening(); }
///@} ///@}
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif #endif // GEOGRAPHICLIB_GNOMONIC_HPP
 End of changes. 7 change blocks. 
9 lines changed or deleted 10 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&ndash;109. * pp. 107&ndash;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


 PolarStereographic.hpp   PolarStereographic.hpp 
/** /**
* \file PolarStereographic.hpp * \file PolarStereographic.hpp
* \brief Header for GeographicLib::PolarStereographic class * \brief Header for GeographicLib::PolarStereographic 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_POLARSTEREOGRAPHIC_HPP) #if !defined(GEOGRAPHICLIB_POLARSTEREOGRAPHIC_HPP)
#define GEOGRAPHICLIB_POLARSTEREOGRAPHIC_HPP "$Id: PolarStereographic.hpp 6 906 2010-12-02 22:10:56Z karney $" #define GEOGRAPHICLIB_POLARSTEREOGRAPHIC_HPP "$Id: 75b815b23a06168b724d4bc6 779bad53029ede58 $"
#include "GeographicLib/Constants.hpp" #include <GeographicLib/Constants.hpp>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief Polar Stereographic Projection * \brief Polar Stereographic 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. 160&ndash;163. * pp. 160&ndash;163.
* *
* This is a straightforward implementation of the equations in Snyder ex cept * This is a straightforward implementation of the equations in Snyder ex cept
* that Newton's method is used to invert the projection. * that Newton's method is used to invert the projection.
**********************************************************************/ **********************************************************************/
class PolarStereographic { class GEOGRAPHIC_EXPORT PolarStereographic {
private: private:
typedef Math::real real; typedef Math::real real;
// _Cx used to be _C but g++ 3.4 has a macro of that name // _Cx used to be _C but g++ 3.4 has a macro of that name
const real _a, _r, _f, _e2, _e, _e2m, _Cx, _c; const real _a, _r, _f, _e2, _e, _e2m, _Cx, _c;
real _k0; real _k0;
static const real tol, overflow; static const real tol_;
static const int numit = 5; static const real overflow_;
static inline real sq(real x) throw() { return x * x; } static const int numit_ = 5;
// Return e * atanh(e * x) for f >= 0, else return // Return e * atanh(e * x) for f >= 0, else return
// - sqrt(-e2) * atan( sqrt(-e2) * x) for f < 0 // - sqrt(-e2) * atan( sqrt(-e2) * x) for 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);
} }
public: public:
/** /**
* Constructor for a ellipsoid with * Constructor for a ellipsoid with
* *
skipping to change at line 160 skipping to change at line 160
/** /**
* A global instantiation of PolarStereographic with the WGS84 ellipsoi d * A global instantiation of PolarStereographic with the WGS84 ellipsoi d
* and the UPS scale factor. However, unlike UPS, no false easting or * and the UPS scale factor. However, unlike UPS, no false easting or
* northing is added. * northing is added.
********************************************************************** / ********************************************************************** /
static const PolarStereographic UPS; static const PolarStereographic UPS;
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif #endif // GEOGRAPHICLIB_POLARSTEREOGRAPHIC_HPP
 End of changes. 6 change blocks. 
7 lines changed or deleted 7 lines changed or added


 TransverseMercator.hpp   TransverseMercator.hpp 
/** /**
* \file TransverseMercator.hpp * \file TransverseMercator.hpp
* \brief Header for GeographicLib::TransverseMercator class * \brief Header for GeographicLib::TransverseMercator 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_TRANSVERSEMERCATOR_HPP) #if !defined(GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP)
#define GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP "$Id: TransverseMercator.hpp 6 950 2011-02-11 04:09:24Z karney $" #define GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP "$Id: efc1fd9413c24c01583c80fd 99435b84f0be1597 $"
#include "GeographicLib/Constants.hpp" #include <GeographicLib/Constants.hpp>
#if !defined(TM_TX_MAXPOW) #if !defined(TM_TX_MAXPOW)
/** /**
* The order of the series approximation used in TransverseMercator. * The order of the series approximation used in TransverseMercator.
* TM_TX_MAXPOW can be set to any integer in [4, 8]. * TM_TX_MAXPOW can be set to any integer in [4, 8].
**********************************************************************/ **********************************************************************/
#define TM_TX_MAXPOW \ #define TM_TX_MAXPOW \
(GEOGRAPHICLIB_PREC == 1 ? 6 : GEOGRAPHICLIB_PREC == 0 ? 4 : 8) (GEOGRAPHICLIB_PREC == 1 ? 6 : GEOGRAPHICLIB_PREC == 0 ? 4 : 8)
#endif #endif
skipping to change at line 71 skipping to change at line 71
* the UTM projections. A more complicated example is the British Nation al * the UTM projections. A more complicated example is the British Nation al
* Grid (<a href="http://www.spatialreference.org/ref/epsg/7405/"> * Grid (<a href="http://www.spatialreference.org/ref/epsg/7405/">
* EPSG:7405</a>) which requires the use of a latitude of origin. This i s * EPSG:7405</a>) which requires the use of a latitude of origin. This i s
* implemented by the GeographicLib::OSGB class. * implemented by the GeographicLib::OSGB class.
* *
* See TransverseMercator.cpp for more information on the implementation. * See TransverseMercator.cpp for more information on the implementation.
* *
* See \ref transversemercator for a discussion of this projection. * See \ref transversemercator for a discussion of this projection.
**********************************************************************/ **********************************************************************/
class TransverseMercator { class GEOGRAPHIC_EXPORT TransverseMercator {
private: private:
typedef Math::real real; typedef Math::real real;
static const int maxpow = TM_TX_MAXPOW; static const int maxpow_ = TM_TX_MAXPOW;
static const real tol, overflow; static const real tol_;
static const int numit = 5; static const real overflow_;
static const int numit_ = 5;
const real _a, _r, _f, _k0, _e2, _e, _e2m, _c, _n; const real _a, _r, _f, _k0, _e2, _e, _e2m, _c, _n;
// _alp[0] and _bet[0] unused // _alp[0] and _bet[0] unused
real _a1, _b1, _alp[maxpow + 1], _bet[maxpow + 1]; real _a1, _b1, _alp[maxpow_ + 1], _bet[maxpow_ + 1];
static inline real sq(real x) throw() { return x * x; }
// tan(x) for x in [-pi/2, pi/2] ensuring that the sign is right // tan(x) for x in [-pi/2, pi/2] ensuring that the sign is right
static inline real tanx(real x) throw() { static inline real tanx(real x) throw() {
real t = std::tan(x); real t = std::tan(x);
// Write the tests this way to ensure that tanx(NaN()) is NaN() // Write the tests this way to ensure that tanx(NaN()) is NaN()
return x >= 0 ? (!(t < 0) ? t : overflow) : (!(t >= 0) ? t : -overflo w); return x >= 0 ? (!(t < 0) ? t : overflow_) : (!(t >= 0) ? t : -overfl ow_);
} }
// Return e * atanh(e * x) for f >= 0, else return // Return e * atanh(e * x) for f >= 0, else return
// - sqrt(-e2) * atan( sqrt(-e2) * x) for f < 0 // - sqrt(-e2) * atan( sqrt(-e2) * x) for 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);
} }
public: public:
/** /**
* Constructor for a ellipsoid with * Constructor for a ellipsoid with
skipping to change at line 193 skipping to change at line 193
/** /**
* A global instantiation of TransverseMercator with the WGS84 ellipsoi d * A global instantiation of TransverseMercator with the WGS84 ellipsoi d
* and the UTM scale factor. However, unlike UTM, no false easting or * and the UTM scale factor. However, unlike UTM, no false easting or
* northing is added. * northing is added.
********************************************************************** / ********************************************************************** /
static const TransverseMercator UTM; static const TransverseMercator UTM;
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif #endif // GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP
 End of changes. 7 change blocks. 
9 lines changed or deleted 9 lines changed or added


 TransverseMercatorExact.hpp   TransverseMercatorExact.hpp 
/** /**
* \file TransverseMercatorExact.hpp * \file TransverseMercatorExact.hpp
* \brief Header for GeographicLib::TransverseMercatorExact class * \brief Header for GeographicLib::TransverseMercatorExact 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_TRANSVERSEMERCATOREXACT_HPP) #if !defined(GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_HPP)
#define GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_HPP "$Id: TransverseMercatorE xact.hpp 6950 2011-02-11 04:09:24Z karney $" #define GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_HPP "$Id: 06681bb71d3ef73c0f4 6224eed80ca77fbeb6fc8 $"
#include "GeographicLib/Constants.hpp" #include <GeographicLib/Constants.hpp>
#include "GeographicLib/EllipticFunction.hpp" #include <GeographicLib/EllipticFunction.hpp>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief An exact implementation of the Transverse Mercator Projection * \brief An exact implementation of the Transverse Mercator Projection
* *
* Implementation of the Transverse Mercator Projection given in * Implementation of the Transverse Mercator Projection given in
* - L. P. Lee, * - L. P. Lee,
* <a href="http://dx.doi.org/10.3138/X687-1574-4325-WM62"> Conformal * <a href="http://dx.doi.org/10.3138/X687-1574-4325-WM62"> Conformal
* Projections Based On Jacobian Elliptic Functions</a>, Part V of * Projections Based On Jacobian Elliptic Functions</a>, Part V of
skipping to change at line 67 skipping to change at line 67
* and TransverseMercatorExact::Reverse functions. The latitude of origi n is * and TransverseMercatorExact::Reverse functions. The latitude of origi n is
* taken to be the equator. See the documentation on TransverseMercator for * taken to be the equator. See the documentation on TransverseMercator for
* how to include a false easting, false northing, or a latitude of origi n. * how to include a false easting, false northing, or a latitude of origi n.
* *
* See TransverseMercatorExact.cpp for more information on the * See TransverseMercatorExact.cpp for more information on the
* implementation. * implementation.
* *
* See \ref transversemercator for a discussion of this projection. * See \ref transversemercator for a discussion of this projection.
**********************************************************************/ **********************************************************************/
class TransverseMercatorExact { class GEOGRAPHIC_EXPORT TransverseMercatorExact {
private: private:
typedef Math::real real; typedef Math::real real;
static const real tol, tol1, tol2, taytol, overflow; static const real tol_;
static const int numit = 10; static const real tol1_;
static const real tol2_;
static const real taytol_;
static const real overflow_;
static const int numit_ = 10;
const real _a, _r, _f, _k0, _mu, _mv, _e, _ep2; const real _a, _r, _f, _k0, _mu, _mv, _e, _ep2;
const bool _extendp; const bool _extendp;
const EllipticFunction _Eu, _Ev; const EllipticFunction _Eu, _Ev;
static inline real sq(real x) throw() { return x * x; }
// tan(x) for x in [-pi/2, pi/2] ensuring that the sign is right // tan(x) for x in [-pi/2, pi/2] ensuring that the sign is right
static inline real tanx(real x) throw() { static inline real tanx(real x) throw() {
real t = std::tan(x); real t = std::tan(x);
// Write the tests this way to ensure that tanx(NaN()) is NaN() // Write the tests this way to ensure that tanx(NaN()) is NaN()
return x >= 0 ? (!(t < 0) ? t : overflow) : (!(t >= 0) ? t : -overflo w); return x >= 0 ? (!(t < 0) ? t : overflow_) : (!(t >= 0) ? t : -overfl ow_);
} }
real taup(real tau) const throw(); real taup(real tau) const throw();
real taupinv(real taup) const throw(); real taupinv(real taup) const throw();
void zeta(real u, real snu, real cnu, real dnu, void zeta(real u, real snu, real cnu, real dnu,
real v, real snv, real cnv, real dnv, real v, real snv, real cnv, real dnv,
real& taup, real& lam) const throw(); real& taup, real& lam) const throw();
void dwdzeta(real u, real snu, real cnu, real dnu, void dwdzeta(real u, real snu, real cnu, real dnu,
skipping to change at line 250 skipping to change at line 253
/** /**
* A global instantiation of TransverseMercatorExact with the WGS84 * A global instantiation of TransverseMercatorExact with the WGS84
* ellipsoid and the UTM scale factor. However, unlike UTM, no false * ellipsoid and the UTM scale factor. However, unlike UTM, no false
* easting or northing is added. * easting or northing is added.
********************************************************************** / ********************************************************************** /
static const TransverseMercatorExact UTM; static const TransverseMercatorExact UTM;
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif #endif // GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_HPP
 End of changes. 7 change blocks. 
8 lines changed or deleted 11 lines changed or added


 UTMUPS.hpp   UTMUPS.hpp 
/** /**
* \file UTMUPS.hpp * \file UTMUPS.hpp
* \brief Header for GeographicLib::UTMUPS class * \brief Header for GeographicLib::UTMUPS 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_UTMUPS_HPP) #if !defined(GEOGRAPHICLIB_UTMUPS_HPP)
#define GEOGRAPHICLIB_UTMUPS_HPP "$Id: UTMUPS.hpp 6937 2011-02-01 20:17:13Z karney $" #define GEOGRAPHICLIB_UTMUPS_HPP "$Id: 280b3ae1d12d4ed1f8a429abf99be36a0546 c309 $"
#include "GeographicLib/Constants.hpp"
#include <sstream> #include <sstream>
#include <GeographicLib/Constants.hpp>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief Convert between Geographic coordinates and UTM/UPS * \brief Convert between Geographic coordinates and UTM/UPS
* *
* UTM and UPS are defined * UTM and UPS are defined
* - J. W. Hager, J. F. Behensky, and B. W. Drew, * - J. W. Hager, J. F. Behensky, and B. W. Drew,
* <a href="http://earth-info.nga.mil/GandG/publications/tm8358.2/TM835 8_2.pdf"> * <a href="http://earth-info.nga.mil/GandG/publications/tm8358.2/TM835 8_2.pdf">
* The Universal Grids: Universal Transverse Mercator (UTM) and Univers al * The Universal Grids: Universal Transverse Mercator (UTM) and Univers al
skipping to change at line 59 skipping to change at line 59
* generous overlap between UTM and UPS and between UTM zones. * generous overlap between UTM and UPS and between UTM zones.
* *
* 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 UTM and UPS. Version 2.4.2 (and * also provides conversions to and from UTM and UPS. Version 2.4.2 (and
* earlier) suffers from some drawbacks: * earlier) suffers from some drawbacks:
* - Inconsistent rules are used to determine the whether a particular UT M or * - Inconsistent rules are used to determine the whether a particular UT M or
* UPS coordinate is legal. A more systematic approach is taken here. * UPS 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 UTMUPS { class GEOGRAPHIC_EXPORT UTMUPS {
private: private:
typedef Math::real real; typedef Math::real real;
static const real falseeasting[4]; static const real falseeasting_[4];
static const real falsenorthing[4]; static const real falsenorthing_[4];
static const real mineasting[4]; static const real mineasting_[4];
static const real maxeasting[4]; static const real maxeasting_[4];
static const real minnorthing[4]; static const real minnorthing_[4];
static const real maxnorthing[4]; static const real maxnorthing_[4];
static real CentralMeridian(int zone) throw() static real CentralMeridian(int zone) throw()
{ return real(6 * zone - 183); } { return real(6 * zone - 183); }
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 void CheckLatLon(real lat, real lon); static void CheckLatLon(real lat, real lon);
// Throw an error if easting or northing are outside standard ranges. If // Throw an error if easting or northing are outside standard ranges. If
// throwp = false, return bool instead. // throwp = false, return bool instead.
static bool CheckCoords(bool utmp, bool northp, real x, real y, static bool CheckCoords(bool utmp, bool northp, real x, real y,
bool msgrlimits = false, bool throwp = true); bool msgrlimits = false, bool throwp = true);
skipping to change at line 313 skipping to change at line 313
* *
* (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 Constants::WGS84_r<real>(); } { return Constants::WGS84_r<real>(); }
///@} ///@}
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif
#endif // GEOGRAPHICLIB_UTMUPS_HPP
 End of changes. 6 change blocks. 
9 lines changed or deleted 9 lines changed or added

This html diff was produced by rfcdiff 1.41. The latest version is available from http://tools.ietf.org/tools/rfcdiff/