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) <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/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_CONSTANTS_HPP) #if !defined(GEOGRAPHICLIB_CONSTANTS_HPP)
#define GEOGRAPHICLIB_CONSTANTS_HPP "$Id: Constants.hpp 6867 2010-09-11 13: 04:26Z karney $" #define GEOGRAPHICLIB_CONSTANTS_HPP "$Id: Constants.hpp 6891 2010-11-16 16: 10:53Z karney $"
/** /**
* A simple compile-time assert. This is designed to be compatible with th e * A simple compile-time assert. This is designed to be compatible with th e
* C++0X static_assert. * C++0X static_assert.
**********************************************************************/ **********************************************************************/
#if !defined(STATIC_ASSERT) #if !defined(STATIC_ASSERT)
#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
#if defined(__GNUC__) #if defined(__GNUC__)
skipping to change at line 45 skipping to change at line 45
/** /**
* 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 75 skipping to change at line 80
**********************************************************************/ **********************************************************************/
class Math { class 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)
/**
* 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.
**********************************************************************
/
typedef long double extended;
#else
typedef double extended;
#endif
#if GEOGRAPHICLIB_PREC == 1 #if GEOGRAPHICLIB_PREC == 1
/** /**
* The real type for %GeographicLib. Nearly all the testing has been do ne * The real type for %GeographicLib. Nearly all the testing has been do ne
* with \e real = double. However, the algorithms should also work wit h * with \e real = double. However, the algorithms should also work wit h
* float and long double (where available). * float and long double (where available).
********************************************************************** / ********************************************************************** /
typedef double real; typedef double real;
#elif GEOGRAPHICLIB_PREC == 0 #elif GEOGRAPHICLIB_PREC == 0
typedef float real; typedef float real;
#elif GEOGRAPHICLIB_PREC == 2 #elif GEOGRAPHICLIB_PREC == 2
typedef long double real; typedef extended real;
#else #else
typedef double real; typedef double real;
#endif #endif
/**
* @return \e pi
**********************************************************************
/
static inline real pi() throw()
// good for about 168-bit accuracy
{ return real(3.1415926535897932384626433832795028841971693993751L); }
/**
* @return the number of radians in a degree.
**********************************************************************
/
static inline real degree() throw() { return pi() / 180; }
#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>).
********************************************************************** / ********************************************************************** /
static inline real hypot(real x, real y) throw() { static inline real hypot(real x, real y) throw() {
x = std::abs(x); x = std::abs(x);
skipping to change at line 277 skipping to change at line 304
/** /**
* 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.
********************************************************************** / ********************************************************************** /
static inline real NaN() throw() { static inline real NaN() throw() {
return std::numeric_limits<real>::has_quiet_NaN ? return std::numeric_limits<real>::has_quiet_NaN ?
std::numeric_limits<real>::quiet_NaN() : std::numeric_limits<real>::quiet_NaN() :
std::numeric_limits<real>::max(); std::numeric_limits<real>::max();
} }
/**
* @return \e pi in extended precision
**********************************************************************
/
static inline extended epi() throw()
// good for about 168-bit accuracy
{ return extended(3.1415926535897932384626433832795028841971693993751L)
; }
/**
* @return the number of radians in a degree in extended precision.
**********************************************************************
/
static inline extended edegree() throw() { return epi() / 180; }
}; };
/** /**
* \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 Constants {
private: private:
typedef Math::real real; typedef Math::real real;
Constants(); // Disable constructor Constants(); // Disable constructor
public: public:
/** /**
* @return \e pi * @return \e pi. This duplicates Math::pi().
********************************************************************** / ********************************************************************** /
static inline Math::real pi() throw() static inline Math::real pi() throw() { return Math::pi(); }
// good for about 123-bit accuracy
{ return real(3.141592653589793238462643383279502884L); }
/** /**
* @return the number of radians in a degree. * @return the number of radians in a degree. This duplicates
* Math::degree().
********************************************************************** / ********************************************************************** /
static inline Math::real degree() throw() { return pi() / 180; } static inline Math::real degree() throw() { return Math::degree(); }
/** /**
* @return the number of radians in an arcminute. * @return the number of radians in an arcminute.
********************************************************************** / ********************************************************************** /
static inline Math::real arcminute() throw() { return degree() / 60; } static inline Math::real arcminute() throw() { return Math::degree() / 60; }
/** /**
* @return the number of radians in an arcsecond. * @return the number of radians in an arcsecond.
********************************************************************** / ********************************************************************** /
static inline Math::real arcsecond() throw() { return arcminute() / 60; static inline Math::real arcsecond() throw()
} { return Math::degree() / 3600; }
/** \name Ellipsoid parameters /** \name Ellipsoid parameters
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* @return equatorial radius of WGS84 ellipsoid * @return the equatorial radius of WGS84 ellipsoid
********************************************************************** / ********************************************************************** /
static inline Math::real WGS84_a() throw() { return 6378137 * meter(); } static inline Math::real WGS84_a() throw() { return 6378137 * meter(); }
/** /**
* @return reciprocal flattening of WGS84 ellipsoid * @return the reciprocal flattening of WGS84 ellipsoid
********************************************************************** / ********************************************************************** /
static inline Math::real WGS84_r() throw() { return real(298.257223563L ); } static inline Math::real WGS84_r() throw() { return real(298.257223563L ); }
/** /**
* @return central scale factor for UTM * @return the central scale factor for UTM
********************************************************************** / ********************************************************************** /
static inline Math::real UTM_k0() throw() {return real(0.9996L); } static inline Math::real UTM_k0() throw() {return real(0.9996L); }
/** /**
* @return central scale factor for UPS * @return the central scale factor for UPS
********************************************************************** / ********************************************************************** /
static inline Math::real UPS_k0() throw() { return real(0.994L); } static inline Math::real UPS_k0() throw() { return real(0.994L); }
///@} ///@}
/** \name SI units /** \name SI units
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* @return the number of meters in a meter. * @return the number of meters in a meter.
* *
* This is unity, but this lets the internal system of units be changed if * This is unity, but this lets the internal system of units be changed if
* necessary. * necessary.
********************************************************************** / ********************************************************************** /
static inline Math::real meter() throw() { return real(1); } static inline Math::real meter() throw() { return real(1); }
/** /**
* @return the number of meters in a kilometer. * @return the number of meters in a kilometer.
********************************************************************** / ********************************************************************** /
static inline Math::real kilometer() throw() { return 1000 * meter(); } static inline Math::real kilometer() throw() { return 1000 * meter(); }
///@}
/** /**
* @return the number of meters in a nautical mile (approximately 1 arc * @return the number of meters in a nautical mile (approximately 1 arc
* minute) * minute)
********************************************************************** / ********************************************************************** /
static inline Math::real nauticalmile() throw() { return 1852 * meter() ; } static inline Math::real nauticalmile() throw() { return 1852 * meter() ; }
///@}
/** \name Anachronistic British units /** \name Anachronistic British units
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* @return the number of meters in an international foot. * @return the number of meters in an international foot.
********************************************************************** / ********************************************************************** /
static inline Math::real foot() throw() static inline Math::real foot() throw()
{ return real(0.0254L) * 12 * meter(); } { return real(0.0254L) * 12 * meter(); }
/** /**
 End of changes. 18 change blocks. 
17 lines changed or deleted 61 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 6867 2010-09-11 13: 04:26Z karney $" #define GEOGRAPHICLIB_GEOCOORDS_HPP "$Id: GeoCoords.hpp 6876 2010-10-18 13: 47:55Z karney $"
#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
skipping to change at line 415 skipping to change at line 415
********************************************************************** / ********************************************************************** /
std::string AltUTMUPSRepresentation(int prec = 0) const; std::string AltUTMUPSRepresentation(int prec = 0) const;
///@} ///@}
/** \name Inspector functions /** \name Inspector functions
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* @return \e a the equatorial radius of the WGS84 ellipsoid (meters). * @return \e a the equatorial radius of the WGS84 ellipsoid (meters).
* *
* (The WGS84 values is returned because the UTM and UPS projections ar e * (The WGS84 value is returned because the UTM and UPS projections are
* based on this ellipsoid.) * based on this ellipsoid.)
********************************************************************** / ********************************************************************** /
Math::real MajorRadius() const throw() { return UTMUPS::MajorRadius(); } Math::real MajorRadius() const throw() { return UTMUPS::MajorRadius(); }
/** /**
* @return \e r the inverse flattening of the WGS84 ellipsoid. * @return \e r the inverse flattening of the WGS84 ellipsoid.
* *
* (The WGS84 values is returned because the UTM and UPS projections ar e * (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. 3 change blocks. 
3 lines changed or deleted 3 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) <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/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_GEOCENTRIC_HPP) #if !defined(GEOGRAPHICLIB_GEOCENTRIC_HPP)
#define GEOGRAPHICLIB_GEOCENTRIC_HPP "$Id: Geocentric.hpp 6867 2010-09-11 1 3:04:26Z karney $" #define GEOGRAPHICLIB_GEOCENTRIC_HPP "$Id: Geocentric.hpp 6876 2010-10-18 1 3:47:55Z karney $"
#include "GeographicLib/Constants.hpp" #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)
skipping to change at line 124 skipping to change at line 124
* value used in the constructor. A value of 0 is returned for a sph ere * value used in the constructor. A value of 0 is returned for a sph ere
* (infinite inverse flattening). * (infinite inverse flattening).
********************************************************************** / ********************************************************************** /
Math::real InverseFlattening() const throw() { return _r; } Math::real InverseFlattening() const throw() { return _r; }
///@} ///@}
/** /**
* 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.
********************************************************************** / ********************************************************************** /
const static Geocentric WGS84; static const Geocentric WGS84;
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif #endif
 End of changes. 2 change blocks. 
2 lines changed or deleted 2 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) <charles@karney.com> * Copyright (c) Charles Karney (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/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_GEODESIC_HPP) #if !defined(GEOGRAPHICLIB_GEODESIC_HPP)
#define GEOGRAPHICLIB_GEODESIC_HPP "$Id: Geodesic.hpp 6867 2010-09-11 13:04 :26Z karney $" #define GEOGRAPHICLIB_GEODESIC_HPP "$Id: Geodesic.hpp 6876 2010-10-18 13:47 :55Z karney $"
#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
skipping to change at line 126 skipping to change at line 126
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, 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, tol0, tol1, tol2, 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() {
skipping to change at line 765 skipping to change at line 765
* (infinite inverse flattening). * (infinite inverse flattening).
********************************************************************** / ********************************************************************** /
Math::real InverseFlattening() const throw() { return _r; } Math::real InverseFlattening() const throw() { return _r; }
/** /**
* @return total area of ellipsoid in meters<sup>2</sup>. The area of a * @return total area of ellipsoid in meters<sup>2</sup>. The area of a
* polygon encircling a pole can be found by adding * polygon encircling a pole can be found by adding
* Geodesic::EllipsoidArea()/2 to the sum of \e S12 for each side of the * Geodesic::EllipsoidArea()/2 to the sum of \e S12 for each side of the
* polygon. * polygon.
********************************************************************** / ********************************************************************** /
Math::real EllipsoidArea() const throw() { Math::real EllipsoidArea() const throw() { return 4 * Math::pi() * _c2;
return 4 * Constants::pi() * _c2; }
}
///@} ///@}
/** /**
* A global instantiation of Geodesic with the parameters for the WGS84 * A global instantiation of Geodesic with the parameters for the WGS84
* ellipsoid. * ellipsoid.
********************************************************************** / ********************************************************************** /
const static Geodesic WGS84; static const Geodesic WGS84;
/** \name Deprecated function. /** \name Deprecated function.
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* <b>DEPRECATED</b> Perform the direct geodesic calculation. Given a * <b>DEPRECATED</b> Perform the direct geodesic calculation. Given a
* latitude, \e lat1, longitude, \e lon1, and azimuth \e azi1 (degrees) for * latitude, \e lat1, longitude, \e lon1, and azimuth \e azi1 (degrees) for
* point 1 and a range, \e s12 (meters) from point 1 to point 2, return the * point 1 and a range, \e s12 (meters) from point 1 to point 2, return the
* latitude, \e lat2, longitude, \e lon2, and forward azimuth, \e azi2 * latitude, \e lat2, longitude, \e lon2, and forward azimuth, \e azi2
* (degrees) for point 2 and the reduced length \e m12 (meters). If ei ther * (degrees) for point 2 and the reduced length \e m12 (meters). If ei ther
 End of changes. 4 change blocks. 
6 lines changed or deleted 5 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) <charles@karney.com> * Copyright (c) Charles Karney (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/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_GEODESICLINE_HPP) #if !defined(GEOGRAPHICLIB_GEODESICLINE_HPP)
#define GEOGRAPHICLIB_GEODESICLINE_HPP "$Id: GeodesicLine.hpp 6867 2010-09- 11 13:04:26Z karney $" #define GEOGRAPHICLIB_GEODESICLINE_HPP "$Id: GeodesicLine.hpp 6875 2010-10- 02 19:31:54Z karney $"
#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
skipping to change at line 540 skipping to change at line 540
/** /**
* @return \e azi1 the azimuth (degrees) of the geodesic line at point 1. * @return \e azi1 the azimuth (degrees) of the geodesic line at point 1.
********************************************************************** / ********************************************************************** /
Math::real Azimuth() const throw() { return Init() ? _azi1 : 0; } Math::real Azimuth() const throw() { return Init() ? _azi1 : 0; }
/** /**
* @return \e azi0 the azimuth (degrees) of the geodesic line as it cro sses * @return \e azi0 the azimuth (degrees) of the geodesic line as it cro sses
* the equator in a northward direction. * the equator in a northward direction.
********************************************************************** / ********************************************************************** /
Math::real EquatorialAzimuth() const throw() { Math::real EquatorialAzimuth() const throw() {
return Init() ? atan2(_salp0, _calp0) / Constants::degree() : 0; return Init() ? atan2(_salp0, _calp0) / Math::degree() : 0;
} }
/** /**
* @return \e a1 the arc length (degrees) between the northward equator ial * @return \e a1 the arc length (degrees) between the northward equator ial
* crossing and point 1. * crossing and point 1.
********************************************************************** / ********************************************************************** /
Math::real EquatorialArc() const throw() { Math::real EquatorialArc() const throw() {
return Init() ? atan2(_ssig1, _csig1) / Constants::degree() : 0; return Init() ? atan2(_ssig1, _csig1) / Math::degree() : 0;
} }
/** /**
* @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 inherited from the Geodesic object used in the construct or. * the value inherited from the Geodesic object used in the construct or.
********************************************************************** / ********************************************************************** /
Math::real MajorRadius() const throw() { return Init() ? _a : 0; } Math::real MajorRadius() const throw() { return Init() ? _a : 0; }
/** /**
* @return \e r the inverse flattening of the ellipsoid. This is the * @return \e r the inverse flattening of the ellipsoid. This is the
 End of changes. 3 change blocks. 
3 lines changed or deleted 3 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) <charles@karney.com> * Copyright (c) Charles Karney (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/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_GEOID_HPP) #if !defined(GEOGRAPHICLIB_GEOID_HPP)
#define GEOGRAPHICLIB_GEOID_HPP "$Id: Geoid.hpp 6868 2010-09-12 10:27:11Z k arney $" #define GEOGRAPHICLIB_GEOID_HPP "$Id: Geoid.hpp 6888 2010-11-13 15:19:42Z k arney $"
#include "GeographicLib/Constants.hpp" #include "GeographicLib/Constants.hpp"
#include <string>
#include <vector> #include <vector>
#include <fstream> #include <fstream>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief Computing the height of the geoid * \brief Computing 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
skipping to change at line 46 skipping to change at line 47
* *
* 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.
* *
* This class is \e not thread safe in that a single instantiation cannot * This class is typically \e not thread safe in that a single instantiat
be ion
* safely used by multiple threads. If multiple threads need to calculat * cannot be safely used by multiple threads because of the way the objec
e t
* geoid heights they should all construct thread-local instantiations. * 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
* thread-local instantiations. Alternatively, set the optional \e
* threadsafe parameter to true in the constuctor. This causes 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
* safe.
**********************************************************************/ **********************************************************************/
class Geoid { class 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 t static const unsigned nterms = ((3 + 1) * (3 + 2))/2; // for a cubic fi t
static const real c0, c0n, c0s; static const real c0, c0n, c0s;
static const real c3[stencilsize * nterms]; static const real c3[stencilsize * nterms];
static const real c3n[stencilsize * nterms]; static const real c3n[stencilsize * nterms];
skipping to change at line 70 skipping to change at line 77
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;
// 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 {
// g++ 3.x doesn't know about the cast to std::ios::streamoff. Just _file.seekg(
// remove the cast in this case. #if !(defined(__GNUC__) && __GNUC__ < 4)
_file.seekg(std::ios::streamoff(_datastart + // g++ 3.x doesn't know about the cast to streamoff.
2ULL * (unsigned(iy) * _swidth + std::ios::streamoff
unsigned(ix)))); #endif
(_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)
ix += _width; ix += _width;
else if (ix >= _width) else if (ix >= _width)
ix -= _width; ix -= _width;
if (_cache && iy >= _yoffset && iy < _yoffset + _ysize && if (_cache && iy >= _yoffset && iy < _yoffset + _ysize &&
((ix >= _xoffset && ix < _xoffset + _xsize) || ((ix >= _xoffset && ix < _xoffset + _xsize) ||
(ix + _width >= _xoffset && ix + _width < _xoffset + _xsize))) { (ix + _width >= _xoffset && ix + _width < _xoffset + _xsize))) {
return real(_data[iy - _yoffset] return real(_data[iy - _yoffset]
skipping to change at line 118 skipping to change at line 127
catch (const std::exception& e) { catch (const std::exception& e) {
// throw GeographicErr("Error reading " + _filename + ": " // throw GeographicErr("Error reading " + _filename + ": "
// + e.what()); // + e.what());
// triggers complaints about the "binary '+'" under Visual Studio . // triggers complaints about the "binary '+'" under Visual Studio .
// So use '+=' instead. // So use '+=' instead.
std::string err("Error reading "); std::string err("Error reading ");
err += _filename; err += _filename;
err += ": "; err += ": ";
err += e.what(); err += e.what();
throw GeographicErr(err); throw GeographicErr(err);
} }
} }
} }
real height(real lat, real lon, bool gradp, real height(real lat, real lon, bool gradp,
real& grade, real& gradn) const; real& grade, real& gradn) const;
Geoid(const Geoid&); // copy constructor not allowed Geoid(const Geoid&); // copy constructor not allowed
Geoid& operator=(const Geoid&); // copy assignment not allowed Geoid& operator=(const Geoid&); // copy assignment not allowed
public: public:
/** /**
skipping to change at line 157 skipping to change at line 165
}; };
/** \name Setting up the geoid /** \name Setting up the geoid
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* Construct a Geoid. * Construct a Geoid.
* *
* @param[in] name the name of the geoid. * @param[in] name the name of the geoid.
* @param[in] path (optional) directory for data file. * @param[in] path (optional) directory for data file.
* @param[in] cubic interpolation method; false means bilinear, true (t * @param[in] cubic (optional) interpolation method; false means biline
he ar,
* default) means cubic. * true (the default) means cubic.
* @param[in] threadsafe (optional), if true, construct a thread safe
* object. The default is false
* *
* The data file is formed by appending ".pgm" to the name. If \e path is * The data file is formed by appending ".pgm" to the name. If \e path is
* specified (and is non-empty), then the file is loaded from directory , \e * specified (and is non-empty), then the file is loaded from directory , \e
* path. Otherwise the path is given by the GEOID_PATH environment * path. Otherwise the path is given by the GEOID_PATH environment
* variable. If that is undefined, a compile-time default path is used * variable. If that is undefined, a compile-time default path is used
* (/usr/local/share/GeographicLib/geoids on non-Windows systems and * (/usr/local/share/GeographicLib/geoids on non-Windows systems and
* C:/cygwin/usr/local/share/GeographicLib/geoids on Windows systems). * C:/cygwin/usr/local/share/GeographicLib/geoids on Windows systems).
* This may throw an error because the file does not exist, is unreadab * This may throw an exception because the file does not exist, is
le, * unreadable, or is corrupt. If the \e threadsafe parameter is true,
* or is corrupt. the
* data set is read into memory (which this may also cause an exception
to
* be thrown), the data file is closed, and single-cell caching is turn
ed
* off; this results in a Geoid object which \e is thread safe.
********************************************************************** / ********************************************************************** /
explicit Geoid(const std::string& name, const std::string& path = "", explicit Geoid(const std::string& name, const std::string& path = "",
bool cubic = true); bool cubic = true, bool threadsafe = false);
/** /**
* Set up a cache. * Set up a cache.
* *
* @param[in] south latitude (degrees) of the south edge of the cached area. * @param[in] south latitude (degrees) of the south edge of the cached area.
* @param[in] west longitude (degrees) of the west edge of the cached a rea. * @param[in] west longitude (degrees) of the west edge of the cached a rea.
* @param[in] north latitude (degrees) of the north edge of the cached area. * @param[in] north latitude (degrees) of the north edge of the cached area.
* @param[in] east longitude (degrees) of the east edge of the cached a rea. * @param[in] east longitude (degrees) of the east edge of the cached a rea.
* *
* Cache the data for the specified rectangular area . \e east is alwa ys * Cache the data for the specified rectangular area . \e east is alwa ys
* interpreted as being east of \e west, if necessary by adding * interpreted as being east of \e west, if necessary by adding
* 360<sup>o</sup> to its value. This may throw an error because of * 360<sup>o</sup> to its value. This may throw an error because of
* insufficent memory or because of an error reading the data from the * insufficent memory or because of an error reading the data from the
* file. In this case, you can catch the error and either do nothing ( you * file. In this case, you can catch the error and either do nothing ( you
* will have no cache in this case) or try again with a smaller area. \e * will have no cache in this case) or try again with a smaller area. \e
* south and \e north should be in the range [-90, 90]; \e west and \e east * south and \e north should be in the range [-90, 90]; \e west and \e east
* should be in the range [-180, 360]. * should be in the range [-180, 360]. An exception is thrown if this
* routine is called on a thread safe Geoid.
********************************************************************** / ********************************************************************** /
void CacheArea(real south, real west, real north, real east) const; void CacheArea(real south, real west, real north, real east) const;
/** /**
* Cache all the data. On most computers, this is fast for data sets w ith * Cache all the data. On most computers, this is fast for data sets w ith
* grid resolution of 5' or coarser. For a 1' grid, the required RAM i s * grid resolution of 5' or coarser. For a 1' grid, the required RAM i s
* 450MB; a 2.5' grid needs 72MB; and a 5' grid needs 18MB. This may t hrow * 450MB; a 2.5' grid needs 72MB; and a 5' grid needs 18MB. This may t hrow
* an error because of insufficent memory or because of an error readin g * an error because of insufficent memory or because of an error readin g
* the data from the file. In this case, you can catch the error and * the data from the file. In this case, you can catch the error and
* either do nothing (you will have no cache in this case) or try using * either do nothing (you will have no cache in this case) or try using
* Geoid::CacheArea on a specific area. * Geoid::CacheArea on a specific area. An exception is thrown if this
* routine is called on a thread safe Geoid.
********************************************************************** / ********************************************************************** /
void CacheAll() const { CacheArea(real(-90), real(0), void CacheAll() const { CacheArea(real(-90), real(0),
real(90), real(360)); } real(90), real(360)); }
/** /**
* Clear the cache. This never throws an error. * Clear the cache. This never throws an error. (This does nothing wi
th a
* thread safe Geoid.)
********************************************************************** / ********************************************************************** /
void CacheClear() const throw(); void CacheClear() const throw();
///@} ///@}
/** \name Compute geoid heights /** \name Compute geoid heights
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* Compute the geoid height at a point * Compute the geoid height at a point
skipping to change at line 340 skipping to change at line 356
/** /**
* @return scale (meters). * @return scale (meters).
* *
* This in used in converting from the pixel values in the data file to * This in used in converting from the pixel values in the data file to
* geoid heights. * geoid heights.
********************************************************************** / ********************************************************************** /
Math::real Scale() const throw() { return _scale; } Math::real Scale() const throw() { return _scale; }
/** /**
* @return true if the object is constructed to be thread safe.
**********************************************************************
/
bool ThreadSafe() const throw() { return _threadsafe; }
/**
* @return true if a data cache is active. * @return true if a data cache is active.
********************************************************************** / ********************************************************************** /
bool Cache() const throw() { return _cache; } bool Cache() const throw() { return _cache; }
/** /**
* @return west edge of the cached area; the cache includes this edge. * @return west edge of the cached area; the cache includes this edge.
********************************************************************** / ********************************************************************** /
Math::real CacheWest() const throw() { Math::real CacheWest() const throw() {
return _cache ? ((_xoffset + (_xsize == _width ? 0 : _cubic) return _cache ? ((_xoffset + (_xsize == _width ? 0 : _cubic)
+ _width/2) % _width - _width/2) / _rlonres : + _width/2) % _width - _width/2) / _rlonres :
skipping to change at line 381 skipping to change at line 402
* @return south edge of the cached area; the cache excludes this edge * @return south edge of the cached area; the cache excludes this edge
* unless it's the south pole. * unless it's the south pole.
********************************************************************** / ********************************************************************** /
Math::real CacheSouth() const throw() { Math::real CacheSouth() const throw() {
return _cache ? 90 - ( _yoffset + _ysize - 1 - _cubic) / _rlatres : 0 ; return _cache ? 90 - ( _yoffset + _ysize - 1 - _cubic) / _rlatres : 0 ;
} }
/** /**
* @return \e a the equatorial radius of the WGS84 ellipsoid (meters). * @return \e a the equatorial radius of the WGS84 ellipsoid (meters).
* *
* (The WGS84 values is returned because the supported geoid models are all * (The WGS84 value is returned because the supported geoid models are all
* based on this ellipsoid.) * based on this ellipsoid.)
********************************************************************** / ********************************************************************** /
Math::real MajorRadius() const throw() { return Constants::WGS84_a(); } Math::real MajorRadius() const throw() { return Constants::WGS84_a(); }
/** /**
* @return \e r the inverse flattening of the WGS84 ellipsoid. * @return \e r the inverse flattening of the WGS84 ellipsoid.
* *
* (The WGS84 values is returned because the supported geoid models are all * (The WGS84 value is returned because the supported geoid models are all
* based on this ellipsoid.) * based on this ellipsoid.)
********************************************************************** / ********************************************************************** /
Math::real InverseFlattening() const throw() Math::real InverseFlattening() const throw()
{ return Constants::WGS84_r(); } { return Constants::WGS84_r(); }
///@} ///@}
/** /**
* Return the compile-time default path for the geoid data files. * Return the compile-time default path for the geoid data files.
********************************************************************** / ********************************************************************** /
static std::string DefaultPath(); static std::string DefaultPath();
 End of changes. 15 change blocks. 
24 lines changed or deleted 51 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 (2009, 2010) <charles@karney.com> * Copyright (c) Charles Karney (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/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_HPP) #if !defined(GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_HPP)
#define GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_HPP "$Id: LambertConformalConic .hpp 6867 2010-09-11 13:04:26Z karney $" #define GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_HPP "$Id: LambertConformalConic .hpp 6898 2010-11-19 19:21:49Z karney $"
#include "GeographicLib/Constants.hpp" #include "GeographicLib/Constants.hpp"
#include <algorithm> #include <algorithm>
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 straightforward implementation of the equations in Snyder ex * This is a implementation of the equations in Snyder except that divide
cept d
* that Newton's method is used to invert the projection and that several * differences have been used to transform the expressions into ones whic
of h
* the formulas are modified so that the projection correctly degenerates * may be evaluated accurately and that Newton's method is used to invert
to the
* the Mercator projection or the polar sterographic projection. * projection. In this implementation, the projection correctly becomes
the
* Mercator projection or the polar sterographic projection when the stan
dard
* latitude is the equator or a pole. The accuracy of the projections is
* 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. If the standard parall * standard parallels are set in the constructor. Internally, the case w
els ith
* are both at a single pole the projection becomes the polar stereograph * two standard parallels is converted into a single standard parallel, t
ic he
* projection (compare with the PolarStereographic class). If the standa * latitude of tangency (also the latitude of minimum scale). This latit
rd ude
* parallels are symmetric around equator, the projection becomes the * is also used as the latitude of origin which is returned by
* Mercator projection. The central meridian (which is a trivial shift o * LambertConformalConic::OriginLatitude. The scale on the latitude of
f * origin is given by LambertConformalConic::CentralScale. The case with
* the longitude) is specified as the \e lon0 argument of the two
* distinct standard parallels where one is a pole is singular and is
* disallowed. The central meridian (which is a trivial shift
* of the longitude) is specified as the \e lon0 argument of the
* LambertConformalConic::Forward and LambertConformalConic::Reverse * LambertConformalConic::Forward and LambertConformalConic::Reverse
* functions. The latitude of origin is taken to be the latitude of tang * functions. There is no provision in this class for specifying a false
ency * easting or false northing or a different latitude of origin. However
* which lies between the standard parallels which is given by * these are can be simply included by the calling function. For example
* LambertConformalConic::OriginLatitude. There is no provision in this the
* class for specifying a false easting or false northing or a different * Pennsylvania South state coordinate system
* latitude of origin. However these are can be simply included by the * (<a href="http://www.spatialreference.org/ref/epsg/3364/"> EPSG:3364</
* calling function. For example the Pennsylvania South state coordinate a>)
* system (<a href="http://www.spatialreference.org/ref/epsg/3364/"> * is obtained by:
* EPSG:3364</a>) is obtained by:
\code \code
const double const double
a = GeographicLib::Constants::WGS84_a(), r = 298.257222101, // GRS80 a = GeographicLib::Constants::WGS84_a(), r = 298.257222101, // GRS80
lat1 = 39 + 56/60.0, lat1 = 40 + 58/60.0, // standard parallels lat1 = 39 + 56/60.0, lat1 = 40 + 58/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 = 75 + 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;
skipping to change at line 79 skipping to change at line 83
// 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 LambertConformalConic {
private: private:
typedef Math::real real; typedef Math::real real;
const real _a, _r, _f, _e2, _e, _e2m; const real _a, _r, _f, _fm, _e2, _e, _e2m;
real _sign, _n, _nc, _lt0, _t0n, _t0nm1, _scale, _lat0, _k0; real _sign, _n, _nc, _t0nm1, _scale, _lat0, _k0;
static const real eps, eps2, epsx, tol, ahypover; real _scbet0, _tchi0, _scchi0, _psi0, _nrho0;
static const real eps, epsx, tol, ahypover;
static const int numit = 5; static const int numit = 5;
static inline real sq(real x) throw() { return x * x; } static inline real sq(real x) throw() { return x * 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);
} }
inline real mf(real sphi, real cphi) const throw() { // Divided differences
return cphi/std::sqrt(1 - _e2 * sq(sphi)); // Snyder's m, p 108, eq 1 // Definition: Df(x,y) = (f(x)-f(y))/(x-y)
4-15 // See: W. M. Kahan and R. J. Fateman,
// Symbolic computation of divided differences,
// SIGSAM Bull. 33(3), 7-28 (1999)
// http://doi.acm.org/10.1145/334714.334716
// http://www.cs.berkeley.edu/~fateman/papers/divdiff.pdf
//
// General rules
// h(x) = f(g(x)): Dh(x,y) = Df(g(x),g(y))*Dg(x,y)
// 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
//
// 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()
// hx = hyp(x)
{ 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))
static inline real Dsn(real x, real y, real sx, real sy) throw() {
// sx = x/hyp(x)
real t = x * y;
return t > 0 ? (x+y) * sq( (sx*sy)/t ) / (sx+sy) :
(x-y != 0 ? (sx-sy) / (x-y) : 1);
} }
inline real tf(real sphi, real cphi) const throw() { // Dlog1p(x,y) = log1p((x-y)/(1+y)/(x-y)
// Snyder's t, p 108, eq 15-9a static inline real Dlog1p(real x, real y) throw() {
// First factor is sqrt((1 - sphi) / (1 + sphi)) real t = x - y; if (t < 0) { t = -t; y = x; }
// Second factor is ((1 + e * sphi)/(1 - e * sphi)) ^ (e/2) return t != 0 ? Math::log1p(t/(1+y)) / t : 1/(1+x);
return (sphi >= 0 ? cphi / (1 + sphi) : (1 - sphi) / cphi) *
std::exp(eatanhe(sphi));
} }
inline real logtf(real sphi, real cphi) const throw() { // Dexp(x,y) = exp((x+y)/2) * 2*sinh((x-y)/2)/(x-y)
// Return log(t). This retains relative accuracy near the equator wh static inline real Dexp(real x, real y) throw() {
ere real t = (x - y)/2;
// t -> 1. This is the negative of the standard Mercator northing. return (t != 0 ? sinh(t)/t : real(1)) * exp((x + y)/2);
Note
// that log( sqrt((1 - sin(phi)) / (1 + sin(phi))) ) = -asinh(tan(phi
))
return -Math::asinh(sphi/std::max(epsx, cphi)) + eatanhe(sphi);
} }
inline real logmtf(real sphi) const throw() { // Dsinh(x,y) = 2*sinh((x-y)/2)/(x-y) * cosh((x+y)/2)
// Return log(m/t). This allows the cancellation of cphi in m and t. // cosh((x+y)/2) = (c+sinh(x)*sinh(y)/c)/2
return std::log((1 + sphi)/std::sqrt(1 - _e2 * sq(sphi))) - // c=sqrt((1+cosh(x))*(1+cosh(y)))
eatanhe(sphi); // cosh((x+y)/2) = sqrt( (sinh(x)*sinh(y) + cosh(x)*cosh(y) + 1)/2 )
static inline real Dsinh(real x, real y, real sx, real sy, real cx, rea
l cy)
// sx = sinh(x), cx = cosh(x)
throw() {
// real t = (x - y)/2, c = sqrt((1 + cx) * (1 + cy));
// return (t != 0 ? sinh(t)/t : real(1)) * (c + sx * sy / c) /2;
real t = (x - y)/2;
return (t != 0 ? sinh(t)/t : real(1)) * sqrt((sx * sy + cx * cy + 1)
/2);
}
// Dasinh(x,y) = asinh((x-y)*(x+y)/(x*sqrt(1+y^2)+y*sqrt(1+x^2)))/(x-y)
// = asinh((x*sqrt(1+y^2)-y*sqrt(1+x^2)))/(x-y)
static inline real Dasinh(real x, real y, real hx, real hy) throw() {
// hx = hyp(x)
real t = x - y;
return t != 0 ?
Math::asinh(x*y > 0 ? t * (x+y) / (x*hy + y*hx) : x*hy - y*hx) / t
:
1/hx;
}
// Deatanhe(x,y) = eatanhe((x-y)/(1-e^2*x*y))/(x-y)
inline real Deatanhe(real x, real y) const throw() {
real t = x - y, d = (1 - _e2 * x * y);
return t != 0 ? eatanhe(t / d) / t : _e2 / d;
} }
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();
public: public:
/** /**
* Constructor with a single standard parallel. * Constructor with a single standard parallel.
* *
* @param[in] a equatorial radius of ellipsoid (meters) * @param[in] a equatorial radius of ellipsoid (meters)
* @param[in] r reciprocal flattening of ellipsoid. Setting \e r = 0 * @param[in] r reciprocal flattening of ellipsoid. Setting \e r = 0
* implies \e r = inf or flattening = 0 (i.e., a sphere). Negative \ e r * implies \e r = inf or flattening = 0 (i.e., a sphere). Negative \ e r
skipping to change at line 160 skipping to change at line 205
* @param[in] r reciprocal flattening of ellipsoid. Setting \e r = 0 * @param[in] r reciprocal flattening of ellipsoid. Setting \e r = 0
* implies \e r = inf or flattening = 0 (i.e., a sphere). Negative \ e r * implies \e r = inf or flattening = 0 (i.e., a sphere). Negative \ e r
* indicates a prolate ellipsoid. * indicates a prolate ellipsoid.
* @param[in] sinlat1 sine of first standard parallel. * @param[in] sinlat1 sine of first standard parallel.
* @param[in] coslat1 cosine of first standard parallel. * @param[in] coslat1 cosine of first standard parallel.
* @param[in] sinlat2 sine of second standard parallel. * @param[in] sinlat2 sine of second standard parallel.
* @param[in] coslat2 cosine of second standard parallel. * @param[in] coslat2 cosine of second standard parallel.
* @param[in] k1 scale on the standard parallels. * @param[in] k1 scale on the standard parallels.
* *
* This allows parallels close to the poles to be specified accurately. * This allows parallels close to the poles to be specified accurately.
* This routine computes the latitude of origin and the scale at this
* latitude. In the case where \e lat1 and \e lat2 are different, the
* errors in this routines are as follows: if abs(\e lat2 - \e lat1) <=
* 160<sup>o</sup> and max(abs(\e lat1), abs(\e lat2)) <=
* 89.9999<sup>o</sup>, then the error in the latitude of origin is les
s
* than 4e-14d and the relative error in the scale is less than 7e-15.
********************************************************************** / ********************************************************************** /
LambertConformalConic(real a, real r, LambertConformalConic(real a, real r,
real sinlat1, real coslat1, real sinlat1, real coslat1,
real sinlat2, real coslat2, real sinlat2, real coslat2,
real k1); real k1);
/** /**
* Set the scale for the projection. * Set the scale for the projection.
* *
* @param[in] lat (degrees). * @param[in] lat (degrees).
skipping to change at line 191 skipping to change at line 242
* @param[in] lat latitude of point (degrees). * @param[in] lat latitude of point (degrees).
* @param[in] lon longitude of point (degrees). * @param[in] lon longitude of point (degrees).
* @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.
* *
* The latitude origin is given by LambertConformalConic::LatitudeOrigi n(). * The latitude origin is given by LambertConformalConic::LatitudeOrigi n().
* No false easting or northing is added and \e lat should be in the ra nge * No false easting or northing is added and \e lat should be in the ra nge
* [-90, 90]; \e lon and \e lon0 should be in the range [-180, 360]. T he * [-90, 90]; \e lon and \e lon0 should be in the range [-180, 360]. T he
* values of \e x and \e y returned for points which project to infinit * error in the projection is less than about 10 nm (true distance) and
y the
* (i.e., one or both of the poles) will be large but finite. The valu * errors in the meridian convergence and scale are consistent with thi
e of s.
* \e k returned for the poles in infinite (unless \e lat equals the * The values of \e x and \e y returned for points which project to
* latitude origin). * infinity (i.e., one or both of the poles) will be large but finite.
********************************************************************** / ********************************************************************** /
void Forward(real lon0, real lat, real lon, void Forward(real lon0, real lat, real lon,
real& x, real& y, real& gamma, real& k) const throw(); real& x, real& y, real& gamma, real& k) const throw();
/** /**
* Reverse projection, from Lambert conformal conic to geographic. * Reverse projection, from Lambert conformal conic to geographic.
* *
* @param[in] lon0 central meridian longitude (degrees). * @param[in] lon0 central meridian longitude (degrees).
* @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 latitude origin is given by LambertConformalConic::LatitudeOrigi n(). * The latitude origin is given by LambertConformalConic::LatitudeOrigi n().
* No false easting or northing is added. \e lon0 should be in the ran ge * No false easting or northing is added. \e lon0 should be in the ran ge
* [-180, 360]. The value of \e lon returned is in the range [-180, 18 0). * [-180, 360]. The value of \e lon returned is in the range [-180, 18 0).
* The error in the projection is less than about 10 nm (true distance)
and
* the errors in the meridian convergence and scale are consistent with
* this.
********************************************************************** / ********************************************************************** /
void Reverse(real lon0, real x, real y, void Reverse(real lon0, real x, real y,
real& lat, real& lon, real& gamma, real& k) const throw(); real& lat, real& lon, real& gamma, real& k) const throw();
/** /**
* LambertConformalConic::Forward without returning the convergence and * LambertConformalConic::Forward without returning the convergence and
* scale. * scale.
********************************************************************** / ********************************************************************** /
void Forward(real lon0, real lat, real lon, void Forward(real lon0, real lat, real lon,
real& x, real& y) const throw() { real& x, real& y) const throw() {
skipping to change at line 274 skipping to change at line 328
* 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 ellip soid * A global instantiation of LambertConformalConic with the WGS84 ellip soid
* and the UPS scale factor and \e stdlat = 0. This degenerates to the * and the UPS scale factor and \e stdlat = 0. This degenerates to the
* Mercator projection. * Mercator projection.
********************************************************************** / ********************************************************************** /
const static LambertConformalConic Mercator; static const LambertConformalConic Mercator;
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif #endif
 End of changes. 14 change blocks. 
58 lines changed or deleted 117 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) <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_MGRS_HPP) #if !defined(GEOGRAPHICLIB_MGRS_HPP)
#define GEOGRAPHICLIB_MGRS_HPP "$Id: MGRS.hpp 6867 2010-09-11 13:04:26Z kar ney $" #define GEOGRAPHICLIB_MGRS_HPP "$Id: MGRS.hpp 6876 2010-10-18 13:47:55Z kar ney $"
#include "GeographicLib/Constants.hpp" #include "GeographicLib/Constants.hpp"
#include "GeographicLib/UTMUPS.hpp" #include "GeographicLib/UTMUPS.hpp"
#include <sstream> #include <sstream>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief Convert between UTM/UPS and %MGRS * \brief Convert between UTM/UPS and %MGRS
* *
skipping to change at line 43 skipping to change at line 43
* - Forward followed by Reverse and vice versa is approximately the * - Forward followed by Reverse and vice versa is approximately the
* identity. (This is affected in predictable ways by errors in * identity. (This is affected in predictable ways by errors in
* determining the latitude band and by loss of precision in the MGRS * determining the latitude band and by loss of precision in the MGRS
* coordinates.) * coordinates.)
* - All MGRS coordinates truncate to legal 100km blocks. All MGRS * - All MGRS coordinates truncate to legal 100km blocks. All MGRS
* coordinates with a legal 100km block prefix are legal (even though t he * coordinates with a legal 100km block prefix are legal (even though t he
* latitude band letter may now belong to a neighboring band). * latitude band letter may now belong to a neighboring band).
* - The range of UTM/UPS coordinates allowed for conversion to MGRS * - The range of UTM/UPS coordinates allowed for conversion to MGRS
* coordinates is the maximum consistent with staying within the letter * coordinates is the maximum consistent with staying within the letter
* ranges of the MGRS scheme. * ranges of the MGRS scheme.
* - All the transformations are implemented as static methods in the MGR
S
* 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 MGRS {
skipping to change at line 81 skipping to change at line 83
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();
skipping to change at line 125 skipping to change at line 127
// 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).
* @param[in] x (meters) * @param[in] x (meters).
* @param[in] y (meters) * @param[in] y (meters).
* @param[in] prec precision relative to 100 km * @param[in] prec precision relative to 100 km.
* @param[out] mgrs MGRS string. * @param[out] mgrs MGRS string.
* *
* \e prec specifies the precision of the MSGRS string as follows: * \e prec specifies the precision of the MSGRS string as follows:
* - prec = 0 (min), 100km * - prec = 0 (min), 100km
* - prec = 1, 10km * - prec = 1, 10km
* - prec = 2, 1km * - prec = 2, 1km
* - prec = 3, 100m * - prec = 3, 100m
* - prec = 4, 10m * - prec = 4, 10m
* - prec = 5, 1m * - prec = 5, 1m
* - prec = 6, 0.1m * - prec = 6, 0.1m
skipping to change at line 219 skipping to change at line 221
/** /**
* Convert a MGRS coordinate to UTM or UPS coordinates. * Convert a MGRS coordinate to UTM or UPS coordinates.
* *
* @param[in] mgrs MGRS string. * @param[in] mgrs MGRS string.
* @param[out] zone UTM zone (zero means UPS). * @param[out] zone UTM zone (zero means UPS).
* @param[out] northp hemisphere (true means north, false means south). * @param[out] northp hemisphere (true means north, false means south).
* @param[out] x (meters). * @param[out] x (meters).
* @param[out] y (meters). * @param[out] y (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 MGRS squar e, * @param[in] centerp if true (default), return center of the MGRS squa re,
* else return SW (lower left) corner. * else return SW (lower left) corner.
* *
* All conversions from MGRS to UTM/UPS are permitted provided the MGRS * All conversions from MGRS to UTM/UPS are permitted provided the MGRS
* coordinate is a possible result of a conversion in the other directi on. * coordinate is a possible result of a conversion in the other directi on.
* (The leading 0 may be dropped from an input MGRS coordinate for UTM * (The leading 0 may be dropped from an input MGRS coordinate for UTM
* zones 1&ndash;9.) In addition, MGRS coordinates with a neighboring * zones 1&ndash;9.) In addition, MGRS coordinates with a neighboring
* latitude band letter are permitted provided that some portion of the * latitude band letter are permitted provided that some portion of the
* 100km block is within the given latitude band. Thus * 100km block is within the given latitude band. Thus
* - 38VLS and 38WLS are allowed (latitude 64N intersects the square * - 38VLS and 38WLS are allowed (latitude 64N intersects the square
* 38[VW]LS); but 38VMS is not permitted (all of 38VMS is north of 64N) * 38[VW]LS); but 38VMS is not permitted (all of 38VMS is north of 64N)
skipping to change at line 254 skipping to change at line 256
static void Reverse(const std::string& mgrs, static void Reverse(const std::string& mgrs,
int& zone, bool& northp, real& x, real& y, int& zone, bool& northp, real& x, real& y,
int& prec, bool centerp = true); int& prec, bool centerp = true);
/** \name Inspector functions /** \name Inspector functions
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* @return \e a the equatorial radius of the WGS84 ellipsoid (meters). * @return \e a the equatorial radius of the WGS84 ellipsoid (meters).
* *
* (The WGS84 values is returned because the UTM and UPS projections ar e * (The WGS84 value is returned because the UTM and UPS projections are
* based on this ellipsoid.) * based on this ellipsoid.)
********************************************************************** / ********************************************************************** /
static Math::real MajorRadius() throw() { return UTMUPS::MajorRadius() ; } static Math::real MajorRadius() throw() { return UTMUPS::MajorRadius(); }
/** /**
* @return \e r the inverse flattening of the WGS84 ellipsoid. * @return \e r the inverse flattening of the WGS84 ellipsoid.
* *
* (The WGS84 values is returned because the UTM and UPS projections ar e * (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
#endif #endif
 End of changes. 8 change blocks. 
11 lines changed or deleted 14 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) <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_POLARSTEREOGRAPHIC_HPP) #if !defined(GEOGRAPHICLIB_POLARSTEREOGRAPHIC_HPP)
#define GEOGRAPHICLIB_POLARSTEREOGRAPHIC_HPP "$Id: PolarStereographic.hpp 6 867 2010-09-11 13:04:26Z karney $" #define GEOGRAPHICLIB_POLARSTEREOGRAPHIC_HPP "$Id: PolarStereographic.hpp 6 876 2010-10-18 13:47:55Z karney $"
#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,
skipping to change at line 154 skipping to change at line 154
* PolarStereographic::SetScale. * PolarStereographic::SetScale.
********************************************************************** / ********************************************************************** /
Math::real CentralScale() const throw() { return _k0; } Math::real CentralScale() const throw() { return _k0; }
///@} ///@}
/** /**
* 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.
********************************************************************** / ********************************************************************** /
const static PolarStereographic UPS; static const PolarStereographic UPS;
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif #endif
 End of changes. 2 change blocks. 
2 lines changed or deleted 2 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) <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/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP) #if !defined(GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP)
#define GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP "$Id: TransverseMercator.hpp 6 867 2010-09-11 13:04:26Z karney $" #define GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP "$Id: TransverseMercator.hpp 6 876 2010-10-18 13:47:55Z karney $"
#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)
skipping to change at line 57 skipping to change at line 57
* The central meridian (which is a trivial shift of the longitude) is * The central meridian (which is a trivial shift of the longitude) is
* specified as the \e lon0 argument of the TransverseMercator::Forward a nd * specified as the \e lon0 argument of the TransverseMercator::Forward a nd
* TransverseMercator::Reverse functions. The latitude of origin is take n to * TransverseMercator::Reverse functions. The latitude of origin is take n to
* be the equator. There is no provision in this class for specifying a * be the equator. There is no provision in this class for specifying a
* false easting or false northing or a different latitude of origin. * false easting or false northing or a different latitude of origin.
* However these are can be simply included by the calling funtcion. For * However these are can be simply included by the calling funtcion. For
* example, the UTMUPS class applies the false easting and false northing for * example, the UTMUPS class applies the false easting and false northing for
* 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
* accommodated by (constants from * implemented by the GeographicLib::OSGB class.
* <a href="http://www.ordnancesurvey.co.uk/oswebsite/gps/information/coo
rdinatesystemsinfo/guidecontents/guidea.html">
* A guide to coordinate systems in Great Britain</a>):
\code
const double
a = 6377563.396, b = 6356256.910, r = a/(a - b), // Airy 1830 ellipsoi
d
k0 = 0.9996012717, lat0 = 49, lon0 = -2, // central scale and origin
fe = 400000, fn = -100000; // false easting and northing
// Set up basic projection
const GeographicLib::TransverseMercator OSGB(a, r, k0);
double x0, y0;
{
// Transform origin point
OSGB.Forward(lon0, lat0, lon0, x0, y0);
x0 -= fe; y0 -= fn; // Combine result with false origin
}
double lat, lon, x, y;
// Sample conversion from geodetic to OSGB grid
std::cin >> lat >> lon;
OSGB.Forward(lon0, lat, lon, x, y);
x -= x0; y -= y0;
std::cout << x << " " << y << "\n";
// Sample conversion from OSGB grid to geodetic
std::cin >> x >> y;
x += x0; y += y0;
OSGB.Reverse(lon0, x, y, lat, lon);
std::cout << lat << " " << lon << "\n";
\endcode
* *
* 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 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;
skipping to change at line 207 skipping to change at line 180
* k0 used in the constructor and is the scale on the central meridia n. * k0 used in the constructor and is the scale on the central meridia n.
********************************************************************** / ********************************************************************** /
Math::real CentralScale() const throw() { return _k0; } Math::real CentralScale() const throw() { return _k0; }
///@} ///@}
/** /**
* 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.
********************************************************************** / ********************************************************************** /
const static TransverseMercator UTM; static const TransverseMercator UTM;
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif #endif
 End of changes. 3 change blocks. 
32 lines changed or deleted 3 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) <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/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_HPP) #if !defined(GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_HPP)
#define GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_HPP "$Id: TransverseMercatorE xact.hpp 6867 2010-09-11 13:04:26Z karney $" #define GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_HPP "$Id: TransverseMercatorE xact.hpp 6876 2010-10-18 13:47:55Z karney $"
#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
skipping to change at line 236 skipping to change at line 236
* k0 used in the constructor and is the scale on the central meridia n. * k0 used in the constructor and is the scale on the central meridia n.
********************************************************************** / ********************************************************************** /
Math::real CentralScale() const throw() { return _k0; } Math::real CentralScale() const throw() { return _k0; }
///@} ///@}
/** /**
* 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.
********************************************************************** / ********************************************************************** /
const static TransverseMercatorExact UTM; static const TransverseMercatorExact UTM;
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif #endif
 End of changes. 2 change blocks. 
2 lines changed or deleted 2 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) <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/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_UTMUPS_HPP) #if !defined(GEOGRAPHICLIB_UTMUPS_HPP)
#define GEOGRAPHICLIB_UTMUPS_HPP "$Id: UTMUPS.hpp 6867 2010-09-11 13:04:26Z karney $" #define GEOGRAPHICLIB_UTMUPS_HPP "$Id: UTMUPS.hpp 6888 2010-11-13 15:19:42Z karney $"
#include "GeographicLib/Constants.hpp" #include "GeographicLib/Constants.hpp"
#include <sstream> #include <sstream>
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
skipping to change at line 288 skipping to change at line 288
* (10<sup>7</sup>). * (10<sup>7</sup>).
********************************************************************** / ********************************************************************** /
static Math::real UTMShift() throw(); static Math::real UTMShift() throw();
/** \name Inspector functions /** \name Inspector functions
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* @return \e a the equatorial radius of the WGS84 ellipsoid (meters). * @return \e a the equatorial radius of the WGS84 ellipsoid (meters).
* *
* (The WGS84 values is returned because the UTM and UPS projections ar e * (The WGS84 value is returned because the UTM and UPS projections are
* based on this ellipsoid.) * based on this ellipsoid.)
********************************************************************** / ********************************************************************** /
static Math::real MajorRadius() throw() { return Constants::WGS84_a(); } static Math::real MajorRadius() throw() { return Constants::WGS84_a(); }
/** /**
* @return \e r the inverse flattening of the WGS84 ellipsoid. * @return \e r the inverse flattening of the WGS84 ellipsoid.
* *
* (The WGS84 values is returned because the UTM and UPS projections ar e * (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(); } { return Constants::WGS84_r(); }
///@} ///@}
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif #endif
 End of changes. 3 change blocks. 
3 lines changed or deleted 3 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/