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) <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_AZIMUTHALEQUIDISTANT_HPP) #if !defined(GEOGRAPHICLIB_AZIMUTHALEQUIDISTANT_HPP)
#define GEOGRAPHICLIB_AZIMUTHALEQUIDISTANT_HPP "$Id: AzimuthalEquidistant.h pp 6778 2010-01-02 21:29:34Z karney $" #define GEOGRAPHICLIB_AZIMUTHALEQUIDISTANT_HPP "$Id: AzimuthalEquidistant.h pp 6867 2010-09-11 13:04:26Z karney $"
#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 GeographicLib::Geodesic object * The conversions all take place using a Geodesic object (by default
(by * Geodesic::WGS84). For more information on geodesics see \ref geodesic
* default GeographicLib::Geodesic::WGS84). For more information on .
* geodesics see \ref geodesic.
**********************************************************************/ **********************************************************************/
class AzimuthalEquidistant { class 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 setting the Geodesic object to * Constructor for AzimuthalEquidistant.
use *
* for geodesic calculations. By default this uses the WGS84 ellipsoid * @param[in] earth the Geodesic object to use for geodesic calculation
. s.
* 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) {}
/** /**
* Convert from latitude \e lat (degrees) and longitude \e lon (degrees * Forward projection, from geographic to azimuthal equidistant.
) to *
* azimuthal equidistant easting \e x (meters) and northing \e y (meter * @param[in] lat0 latitude of center point of projection (degrees).
s). * @param[in] lon0 longitude of center point of projection (degrees).
* The center of the projection is at latitude \e lat0 (degrees) and * @param[in] lat latitude of point (degrees).
* longitude \e lon0 (degrees). Also return the azimuth \e azi (degree * @param[in] lon longitude of point (degrees).
s) * @param[out] x easting of point (meters).
* and the reciprocal of the azimuthal scale \e rk. \e lat0 and \e lat * @param[out] y northing of point (meters).
* should be in the range [-90, 90] and \e lon0 and \e lon should be in * @param[out] azi azimuth of geodesic at point (degrees).
the * @param[out] rk reciprocal of azimuthal scale at point.
* range [-180, 360]. The scale of the projection is 1 in the "radial" *
* direction, \e azi clockwise from true north, and is 1/\e rk in the * \e lat0 and \e lat should be in the range [-90, 90] and \e lon0 and
* direction perpendicular to this. A call to Forward followed by a ca \e
ll * lon should be in the range [-180, 360]. The scale of the projection
* to Reverse will return the original (\e lat, \e lon) (to within is
* roundoff). * 1 in the "radial" direction, \e azi clockwise from true north, and i
s
* 1/\e rk in the direction perpendicular to this. A call to Forward
* followed by a call to Reverse will return the original (\e lat, \e l
on)
* (to within roundoff).
********************************************************************** / ********************************************************************** /
void Forward(real lat0, real lon0, real lat, real lon, void Forward(real lat0, real lon0, real lat, real lon,
real& x, real& y, real& azi, real& rk) const throw(); real& x, real& y, real& azi, real& rk) const throw();
/** /**
* Convert from azimuthal equidistant easting \e x (meters) and northin * Reverse projection, from azimuthal equidistant to geographic.
g \e *
* y (meters) to latitude \e lat (degrees) and longitude \e lon (degree * @param[in] lat0 latitude of center point of projection (degrees).
s). * @param[in] lon0 longitude of center point of projection (degrees).
* The center of the projection is at latitude \e lat0 (degrees) and * @param[in] x easting of point (meters).
* longitude \e lon0 (degrees). Also return the azimuth \e azi (degree * @param[in] y northing of point (meters).
s) * @param[out] lat latitude of point (degrees).
* and the reciprocal of the azimuthal scale \e rk. \e lat0 should be * @param[out] lon longitude of point (degrees).
in * @param[out] azi azimuth of geodesic at point (degrees).
* the range [-90, 90] and \e lon0 should be in the range [-180, 360]. * @param[out] rk reciprocal of azimuthal scale at point.
\e *
* lat will be in the range [-90, 90] and \e lon will be in the range * \e lat0 should be in the range [-90, 90] and \e lon0 should be in th
* [-180, 180). The scale of the projection is 1 in the "radial" e
* direction, \e azi clockwise from true north, and is 1/\e rk in the * range [-180, 360]. \e lat will be in the range [-90, 90] and \e lon
* direction perpendicular to this. A call to Reverse followed by a ca * will be in the range [-180, 180). The scale of the projection is 1
ll in
* to Forward will return the original (\e x, \e y) (to roundoff) only * the "radial" direction, \e azi clockwise from true north, and is 1/\
if e rk
* the geodesic to (\e x, \e y) is a shortest path. * in the direction perpendicular to this. A call to Reverse followed
by a
* call to Forward will return the original (\e x, \e y) (to roundoff)
only
* if the geodesic to (\e x, \e y) is a shortest path.
********************************************************************** / ********************************************************************** /
void Reverse(real lat0, real lon0, real x, real y, void Reverse(real lat0, real lon0, real x, real y,
real& lat, real& lon, real& azi, real& rk) const throw(); real& lat, real& lon, real& azi, real& rk) const throw();
/** /**
* The major radius of the ellipsoid (meters). This is that value of \ * AzimuthalEquidistant::Forward without returning the azimuth and scal
e a e.
* inherited from the Geodesic object used in the constructor. **********************************************************************
/
void Forward(real lat0, real lon0, real lat, real lon,
real& x, real& y) const throw() {
real azi, rk;
Forward(lat0, lon0, lat, lon, x, y, azi, rk);
}
/**
* AzimuthalEquidistant::Reverse without returning the azimuth and scal
e.
**********************************************************************
/
void Reverse(real lat0, real lon0, real x, real y,
real& lat, real& lon) const throw() {
real azi, rk;
Reverse(lat0, lon0, x, y, lat, lon, azi, rk);
}
/** \name Inspector functions
**********************************************************************
/
///@{
/**
* @return \e a the equatorial radius of the ellipsoid (meters). This
is
* the value inherited from the Geodesic object used in the construct
or.
********************************************************************** / ********************************************************************** /
Math::real MajorRadius() const throw() { return _earth.MajorRadius(); } Math::real MajorRadius() const throw() { return _earth.MajorRadius(); }
/** /**
* The inverse flattening of the ellipsoid. This is that value of \e r * @return \e r the inverse flattening of the ellipsoid. This is the
* inherited from the Geodesic object used in the constructor. A value * value inherited from the Geodesic object used in the constructor.
of A
* 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
 End of changes. 8 change blocks. 
51 lines changed or deleted 88 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) <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_CASSINISOLDNER_HPP) #if !defined(GEOGRAPHICLIB_CASSINISOLDNER_HPP)
#define GEOGRAPHICLIB_CASSINISOLDNER_HPP "$Id: CassiniSoldner.hpp 6827 2010 -05-20 19:56:18Z karney $" #define GEOGRAPHICLIB_CASSINISOLDNER_HPP "$Id: CassiniSoldner.hpp 6867 2010 -09-11 13:04:26Z karney $"
#include "GeographicLib/Geodesic.hpp" #include "GeographicLib/Geodesic.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
skipping to change at line 50 skipping to change at line 51
* for which there may be two equally short routes for either leg of the * for which there may be two equally short routes for either leg of the
* path. * path.
* *
* Because of the properties of geodesics, the (\e x, \e y) grid is * Because of the properties of geodesics, the (\e x, \e y) grid is
* orthogonal. The scale in the easting direction is unity. The scale, \e * orthogonal. The scale in the easting direction is unity. The scale, \e
* k, in the northing direction is unity on the central meridian and * k, in the northing direction is unity on the central meridian and
* increases away from the central meridian. The projection routines ret urn * increases away from the central meridian. The projection routines ret urn
* \e azi, the true bearing of the easting direction, and \e rk = 1/\e k, the * \e azi, the true bearing of the easting direction, and \e rk = 1/\e k, the
* reciprocal of the scale in the northing direction. * reciprocal of the scale in the northing direction.
* *
* The conversions all take place using a GeographicLib::Geodesic object * The conversions all take place using a Geodesic object (by default
(by * Geodesic::WGS84). For more information on geodesics see \ref geodesic
* default GeographicLib::Geodesic::WGS84). For more information on .
* geodesics see \ref geodesic. The determination of (\e lat1, \e lon1) * The determination of (\e lat1, \e lon1) in the forward projection is b
in y
* the forward projection is by solving the inverse geodesic problem for * solving the inverse geodesic problem for (\e lat, \e lon) and its twin
(\e * obtained by reflection in the meridional plane. The scale is found by
* lat, \e lon) and its twin obtained by reflection in the meridional pla * determining where two neighboring geodesics intersecting the central
ne. * meridan at \e lat1 and \e lat1 + \e dlat1 intersect and taking the rat
* The scale is found by determining where two neighboring geodesics io
* intersecting the central meridan at \e lat1 and \e lat1 + \e dlat1 * of the reduced lengths for the two geodesics between that point and,
* intersect and taking the ratio of the reduced lengths for the two * respectively, (\e lat1, \e lon1) and (\e lat, \e lon).
* geodesics between that point and, respectively, (\e lat1, \e lon1) and
(\e
* lat, \e lon).
**********************************************************************/ **********************************************************************/
class CassiniSoldner { class 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, eps2;
static const unsigned maxit = 10; static const unsigned maxit = 10;
skipping to change at line 97 skipping to change at line 97
return x < 0 ? -y : y; return x < 0 ? -y : y;
} }
static inline void SinCosNorm(real& sinx, real& cosx) throw() { static inline void SinCosNorm(real& sinx, real& cosx) throw() {
real r = Math::hypot(sinx, cosx); real r = Math::hypot(sinx, cosx);
sinx /= r; sinx /= r;
cosx /= r; cosx /= r;
} }
public: public:
/** /**
* Constructor for CassiniSoldner setting the Geodesic object, \e earth * Constructor for CassiniSoldner.
, to *
* use for geodesic calculations. By default this uses the WGS84 * @param[in] earth the Geodesic object to use for geodesic calculation
* ellipsoid. This constructor makes an "uninitialized" object. Call s.
* Reset to set the central latitude and longuitude, prior to calling * By default this uses the WGS84 ellipsoid.
* Forward and Reverse. *
* This constructor makes an "uninitialized" object. Call Reset to set
the
* central latitude and longuitude, prior to calling Forward and Revers
e.
********************************************************************** / ********************************************************************** /
explicit CassiniSoldner(const Geodesic& earth = Geodesic::WGS84) throw( ) explicit CassiniSoldner(const Geodesic& earth = Geodesic::WGS84) throw( )
: _earth(earth) {} : _earth(earth) {}
/** /**
* Constructor for CassiniSoldner setting the center point, \e lat0, \e * Constructor for CassiniSoldner specifying a center point.
* lon0 (degrees) of the projection and the Geodesic object, \e earth, *
to * @param[in] lat0 latitude of center point of projection (degrees).
* use for geodesic calculations. By default this uses the WGS84 * @param[in] lon0 longitude of center point of projection (degrees).
* ellipsoid. * @param[in] earth the Geodesic object to use for geodesic calculation
s.
* By default this uses the WGS84 ellipsoid.
*
* \e lat0 should be in the range [-90, 90] and \e lon0 should be in th
e
* range [-180, 360].
********************************************************************** / ********************************************************************** /
CassiniSoldner(real lat0, real lon0, CassiniSoldner(real lat0, real lon0,
const Geodesic& earth = Geodesic::WGS84) throw() const Geodesic& earth = Geodesic::WGS84) throw()
: _earth(earth) { : _earth(earth) {
Reset(lat0, lon0); Reset(lat0, lon0);
} }
/** /**
* Set the central latititude to \e lat0 and central longitude to \e lo * Set the central point of the projection
n0 *
* (degrees). \e lat0 should be in the range [-90, 90] and \e lon0 sho * @param[in] lat0 latitude of center point of projection (degrees).
uld * @param[in] lon0 longitude of center point of projection (degrees).
* be in the range [-180, 360]. *
* \e lat0 should be in the range [-90, 90] and \e lon0 should be in th
e
* range [-180, 360].
********************************************************************** / ********************************************************************** /
void Reset(real lat0, real lon0) throw(); void Reset(real lat0, real lon0) throw();
/** /**
* Convert from latitude \e lat (degrees) and longitude \e lon (degrees * Forward projection, from geographic to Cassini-Soldner.
) to *
* Cassini-Soldner easting \e x (meters) and northing \e y (meters). A * @param[in] lat latitude of point (degrees).
lso * @param[in] lon longitude of point (degrees).
* return the azimuth of the easting direction \e azi (degrees) and the * @param[out] x easting of point (meters).
* reciprocal of the northing scale \e rk. \e lat should be in the ran * @param[out] y northing of point (meters).
ge * @param[out] azi azimuth of easting direction at point (degrees).
* [-90, 90] and \e lon should be in the range [-180, 360]. A call to * @param[out] rk reciprocal of azimuthal northing scale at point.
* Forward followed by a call to Reverse will return the original (\e l *
at, * \e lat should be in the range [-90, 90] and \e lon should be in the
* \e lon) (to within roundoff). The routine does nothing if the origi * range [-180, 360]. A call to Forward followed by a call to Reverse
n will
* has not been set. * return the original (\e lat, \e lon) (to within roundoff). The rout
ine
* does nothing if the origin has not been set.
********************************************************************** / ********************************************************************** /
void Forward(real lat, real lon, void Forward(real lat, real lon,
real& x, real& y, real& azi, real& rk) const throw(); real& x, real& y, real& azi, real& rk) const throw();
/** /**
* Convert from Cassini-Soldner easting \e x (meters) and northing \e y * Reverse projection, from Cassini-Soldner to geographic.
* (meters) to latitude \e lat (degrees) and longitude \e lon (degrees) *
. * @param[in] x easting of point (meters).
* Also return the azimuth of the easting direction \e azi (degrees) an * @param[in] y northing of point (meters).
d * @param[out] lat latitude of point (degrees).
* the reciprocal of the northing scale \e rk. A call to Reverse follo * @param[out] lon longitude of point (degrees).
wed * @param[out] azi azimuth of easting direction at point (degrees).
* by a call to Forward will return the original (\e x, \e y) (to withi * @param[out] rk reciprocal of azimuthal northing scale at point.
n *
* roundoff), provided that \e x and \e y are sufficiently small not to * A call to Reverse followed by a call to Forward will return the orig
* "wrap around" the earth. The routine does nothing if the origin has inal
not * (\e x, \e y) (to within roundoff), provided that \e x and \e y are
* been set. * sufficiently small not to "wrap around" the earth. The routine does
* nothing if the origin has not been set.
********************************************************************** / ********************************************************************** /
void Reverse(real x, real y, void Reverse(real x, real y,
real& lat, real& lon, real& azi, real& rk) const throw(); real& lat, real& lon, real& azi, real& rk) const throw();
/** /**
* Has this object been initialized with an origin? * CassiniSoldner::Forward without returning the azimuth and scale.
**********************************************************************
/
void Forward(real lat, real lon,
real& x, real& y) const throw() {
real azi, rk;
Forward(lat, lon, x, y, azi, rk);
}
/**
* CassiniSoldner::Reverse without returning the azimuth and scale.
**********************************************************************
/
void Reverse(real x, real y,
real& lat, real& lon) const throw() {
real azi, rk;
Reverse(x, y, lat, lon, azi, rk);
}
/** \name Inspector functions
**********************************************************************
/
///@{
/**
* @return true if the object has been initialized.
********************************************************************** / ********************************************************************** /
bool Init() const throw() { return _meridian.Init(); } bool Init() const throw() { return _meridian.Init(); }
/** /**
* Return the latitude of the origin (degrees). * @return \e lat0 the latitude of origin (degrees).
********************************************************************** / ********************************************************************** /
Math::real LatitudeOrigin() const throw() Math::real LatitudeOrigin() const throw()
{ return _meridian.Latitude(); } { return _meridian.Latitude(); }
/** /**
* Return the longitude of the origin (degrees). * @return \e lon0 the longitude of origin (degrees).
********************************************************************** / ********************************************************************** /
Math::real LongitudeOrigin() const throw() Math::real LongitudeOrigin() const throw()
{ return _meridian.Longitude(); } { return _meridian.Longitude(); }
/** /**
* The major radius of the ellipsoid (meters). This is that value of \ * @return \e a the equatorial radius of the ellipsoid (meters). This
e a is
* inherited from the Geodesic object used in the constructor. * the value inherited from the Geodesic object used in the construct
or.
********************************************************************** / ********************************************************************** /
Math::real MajorRadius() const throw() { return _earth.MajorRadius(); } Math::real MajorRadius() const throw() { return _earth.MajorRadius(); }
/** /**
* The inverse flattening of the ellipsoid. This is that value of \e r * @return \e r the inverse flattening of the ellipsoid. This is the
* inherited from the Geodesic object used in the constructor. A value * value inherited from the Geodesic object used in the constructor.
of A
* 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
 End of changes. 14 change blocks. 
68 lines changed or deleted 108 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) <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 6842 2010-07-16 20: 51:06Z karney $" #define GEOGRAPHICLIB_CONSTANTS_HPP "$Id: Constants.hpp 6867 2010-09-11 13: 04:26Z 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 92 skipping to change at line 92
#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 long double real;
#else #else
typedef double real; typedef double real;
#endif #endif
#if defined(DOXYGEN) #if defined(DOXYGEN)
/** /**
* Return sqrt(\e x<sup>2</sup> + \e y<sup>2</sup>) avoiding underflow * The hypotenuse function avoiding underflow and overflow.
and *
* overflow. * @param[in] x
* @param[in] y
* @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);
y = std::abs(y); y = std::abs(y);
real real
a = std::max(x, y), a = std::max(x, y),
b = std::min(x, y) / a; b = std::min(x, y) / a;
return a * std::sqrt(1 + b * b); return a * std::sqrt(1 + b * b);
} }
#elif defined(_MSC_VER) #elif defined(_MSC_VER)
skipping to change at line 126 skipping to change at line 129
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(__NO_LONG_DOUBLE_MATH)
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) #if defined(DOXYGEN) || defined(_MSC_VER)
/** /**
* Return 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.
*
* @param[in] x
* @return exp(\e x) - 1.
********************************************************************** / ********************************************************************** /
static inline real expm1(real x) throw() { static inline real expm1(real x) throw() {
volatile real volatile real
y = std::exp(x), y = std::exp(x),
z = y - 1; z = y - 1;
// The reasoning here is similar to that for log1p. The expression // The reasoning here is similar to that for log1p. The expression
// mathematically reduces to exp(x) - 1, and the factor z/log(y) = (y - // mathematically reduces to exp(x) - 1, and the factor z/log(y) = (y -
// 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);
skipping to change at line 151 skipping to change at line 157
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(__NO_LONG_DOUBLE_MATH)
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) #if defined(DOXYGEN) || defined(_MSC_VER)
/** /**
* Return log(\e x + 1) accurate near \e x = 0. This is taken See * log(\e x + 1) accurate near \e x = 0. This is taken See
* D. Goldberg, * D. Goldberg,
* <a href="http://docs.sun.com/source/806-3568/ncg_goldberg.html"> Wha t * <a href="http://docs.sun.com/source/806-3568/ncg_goldberg.html"> Wha t
* every computer scientist should know about floating-point arithmetic </a> * every computer scientist should know about floating-point arithmetic </a>
* (1991), Theorem 4. See also, Higham (op. cit.), Answer to Problem 1 .5, * (1991), Theorem 4. See also, Higham (op. cit.), Answer to Problem 1 .5,
* p 528. * p 528.
*
* @param[in] x
* @return log(\e x + 1).
********************************************************************** / ********************************************************************** /
static inline real log1p(real x) throw() { static inline real log1p(real x) throw() {
volatile real volatile real
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
// approx x, thus log(y)/z (which is nearly constant near z = 0) retu rns // approx x, thus log(y)/z (which is nearly constant near z = 0) retu rns
// 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;
skipping to change at line 179 skipping to change at line 188
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(__NO_LONG_DOUBLE_MATH)
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) #if defined(DOXYGEN) || defined(_MSC_VER)
/** /**
* Return asinh(\e x). This is defined in terms of Math::log1p(\e x) i * The inverse hyperbolic sine function. This is defined in terms of
n * Math::log1p(\e x) in order to maintain accuracy near \e x = 0. In
* order to maintain accuracy near \e x = 0. In addition, the odd pari * addition, the odd parity of the function is enforced.
ty *
* of the function is enforced. * @param[in] x
* @return asinh(\e x).
********************************************************************** / ********************************************************************** /
static inline real asinh(real x) throw() { static inline real asinh(real x) throw() {
real y = std::abs(x); // Enforce odd parity real y = std::abs(x); // Enforce odd parity
y = log1p(y * (1 + y/(hypot(real(1), y) + 1))); y = log1p(y * (1 + y/(hypot(real(1), y) + 1)));
return x < 0 ? -y : y; return x < 0 ? -y : y;
} }
#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(__NO_LONG_DOUBLE_MATH)
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) #if defined(DOXYGEN) || defined(_MSC_VER)
/** /**
* Return atanh(\e x). This is defined in terms of Math::log1p(\e x) i * The inverse hyperbolic tangent function. This is defined in terms o
n f
* order to maintain accuracy near \e x = 0. In addition, the odd pari * Math::log1p(\e x) in order to maintain accuracy near \e x = 0. In
ty * addition, the odd parity of the function is enforced.
* of the function is enforced. *
* @param[in] x
* @return atanh(\e x).
********************************************************************** / ********************************************************************** /
static inline real atanh(real x) throw() { static inline real atanh(real x) throw() {
real y = std::abs(x); // Enforce odd parity real 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;
} }
#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(__NO_LONG_DOUBLE_MATH)
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) #if defined(DOXYGEN) || defined(_MSC_VER)
/** /**
* Return the real cube root of \e x. * The cube root function.
*
* @param[in] x
* @return the real cube root of \e x.
********************************************************************** / ********************************************************************** /
static inline real cbrt(real x) throw() { static inline real cbrt(real x) throw() {
real y = std::pow(std::abs(x), 1/real(3)); // Return the real cube ro ot real y = std::pow(std::abs(x), 1/real(3)); // Return the real cube ro ot
return x < 0 ? -y : y; return x < 0 ? -y : y;
} }
#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(__NO_LONG_DOUBLE_MATH)
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
/** /**
* Return true if number is finite, false if NaN or infinite. * Test for finiteness.
*
* @param[in] x
* @return true if number is finite, false if NaN or infinite.
********************************************************************** / ********************************************************************** /
static inline bool isfinite(real x) throw() { static inline bool isfinite(real x) throw() {
#if defined(DOXYGEN) #if defined(DOXYGEN)
return std::abs(x) <= std::numeric_limits<real>::max(); return std::abs(x) <= std::numeric_limits<real>::max();
#elif defined(_MSC_VER) #elif defined(_MSC_VER)
return _finite(x) != 0; return _finite(x) != 0;
#else #else
return std::isfinite(x) != 0; return std::isfinite(x) != 0;
#endif #endif
} }
/** /**
* Return a NaN if available, otherwise return the max real. * The NaN (not a number)
*
* @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();
} }
}; };
/** /**
* \brief %Constants needed by %GeographicLib * \brief %Constants needed by %GeographicLib
skipping to change at line 269 skipping to change at line 292
* 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:
/** /**
* pi * @return \e pi
********************************************************************** / ********************************************************************** /
static inline Math::real pi() throw() static inline Math::real pi() throw()
// good for about 123-bit accuracy // good for about 123-bit accuracy
{ return real(3.141592653589793238462643383279502884L); } { return real(3.141592653589793238462643383279502884L); }
/** /**
* Factor to convert from degrees to radians * @return the number of radians in a degree.
********************************************************************** / ********************************************************************** /
static inline Math::real degree() throw() { return pi() / 180; } static inline Math::real degree() throw() { return pi() / 180; }
/** /**
* Factor to convert from minutes to radians * @return the number of radians in an arcminute.
********************************************************************** / ********************************************************************** /
static inline Math::real arcminute() throw() { return degree() / 60; } static inline Math::real arcminute() throw() { return degree() / 60; }
/** /**
* Factor to convert from seconds to radians * @return the number of radians in an arcsecond.
********************************************************************** / ********************************************************************** /
static inline Math::real arcsecond() throw() { return arcminute() / 60; } static inline Math::real arcsecond() throw() { return arcminute() / 60; }
/** \name Ellipsoid parameters /** \name Ellipsoid parameters
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* Major radius of WGS84 ellipsoid * @return 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(); }
/** /**
* Reciprocal flattening of WGS84 ellipsoid * @return 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 ); }
/** /**
* Central scale factor for UTM * @return 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); }
/** /**
* Central scale factor for UPS * @return 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
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* Factor to convert from meters to meters (i.e., 1, but this lets the * @return the number of meters in a meter.
* internal system of units be changed if necessary). *
* This is unity, but this lets the internal system of units be changed
if
* necessary.
********************************************************************** / ********************************************************************** /
static inline Math::real meter() throw() { return real(1); } static inline Math::real meter() throw() { return real(1); }
/** /**
* Factor to convert from kilometers to meters. * @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(); }
///@} ///@}
/** /**
* Factor to convert from nautical miles (approximately 1 arc minute) t * @return the number of meters in a nautical mile (approximately 1 arc
o * minute)
* meters.
********************************************************************** / ********************************************************************** /
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
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* Factor to convert from international feet to meters. * @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(); }
/** /**
* Factor to convert from yards to meters. * @return the number of meters in a yard.
********************************************************************** / ********************************************************************** /
static inline Math::real yard() throw() { return 3 * foot(); } static inline Math::real yard() throw() { return 3 * foot(); }
/** /**
* Factor to convert from fathoms to meters. * @return the number of meters in a fathom.
********************************************************************** / ********************************************************************** /
static inline Math::real fathom() throw() { return 2 * yard(); } static inline Math::real fathom() throw() { return 2 * yard(); }
/** /**
* Factor to convert from chains to meters. * @return the number of meters in a chain.
********************************************************************** / ********************************************************************** /
static inline Math::real chain() throw() { return 22 * yard(); } static inline Math::real chain() throw() { return 22 * yard(); }
/** /**
* Factor to convert from furlongs to meters. * @return the number of meters in a furlong.
********************************************************************** / ********************************************************************** /
static inline Math::real furlong() throw() { return 10 * chain(); } static inline Math::real furlong() throw() { return 10 * chain(); }
/** /**
* Factor to convert from statute miles to meters. * @return the number of meters in a statute mile.
********************************************************************** / ********************************************************************** /
static inline Math::real mile() throw() { return 8 * furlong(); } static inline Math::real mile() throw() { return 8 * furlong(); }
///@} ///@}
/** \name Anachronistic US units /** \name Anachronistic US units
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* Factor to convert from US survery feet to meters. * @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(); } { return real(1200) / real(3937) * meter(); }
///@} ///@}
}; };
/** /**
* \brief %Exception handling for %GeographicLib * \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 off 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 takes a string message, \e msg, which is accessible in t * Constructor
he *
* catch clause, via what(). * @param[in] msg a string message, which is accessible in the catch
* 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
 End of changes. 30 change blocks. 
43 lines changed or deleted 65 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) <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_DMS_HPP) #if !defined(GEOGRAPHICLIB_DMS_HPP)
#define GEOGRAPHICLIB_DMS_HPP "$Id: DMS.hpp 6827 2010-05-20 19:56:18Z karne y $" #define GEOGRAPHICLIB_DMS_HPP "$Id: DMS.hpp 6866 2010-09-11 02:15:29Z karne y $"
#include "GeographicLib/Constants.hpp" #include "GeographicLib/Constants.hpp"
#include <sstream> #include <sstream>
#include <iomanip> #include <iomanip>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief Convert between degrees and %DMS representation. * \brief Convert between degrees and %DMS representation.
* *
skipping to change at line 47 skipping to change at line 47
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];
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. DMS::AZIMUTH is used in Encode to indicate output i * and longitudes.
n
* [000, 360) with no letter indicator. DMS::NUMBER is used in Encode
to
* indicate the output of a plain number.
********************************************************************** / ********************************************************************** /
enum flag { NONE = 0, LATITUDE = 1, LONGITUDE = 2, AZIMUTH = 3, enum flag {
NUMBER = 4, }; /**
* No indicator present.
* @hideinitializer
********************************************************************
**/
NONE = 0,
/**
* Latitude indicator (N/S) present.
* @hideinitializer
********************************************************************
**/
LATITUDE = 1,
/**
* Longitude indicator (E/W) present.
* @hideinitializer
********************************************************************
**/
LONGITUDE = 2,
/**
* Used in Encode to indicate output of an azimuth in [000, 360) with
no
* letter indicator.
* @hideinitializer
********************************************************************
**/
AZIMUTH = 3,
/**
* Used in Encode to indicate output of a plain number.
* @hideinitializer
********************************************************************
**/
NUMBER = 4,
};
/** /**
* Indicator for trailing units on an angle. * Indicator for trailing units on an angle.
********************************************************************** / ********************************************************************** /
enum component { DEGREE = 0, MINUTE = 1, SECOND = 2 }; enum component {
/**
/** * Trailing unit is degrees.
* Read a string \e dms in DMS format and return the resulting angle in * @hideinitializer
* degrees. Degrees, minutes, and seconds are indicated by the letters ********************************************************************
d, **/
* ', &quot;, and these components may only be given in this order. An DEGREE = 0,
y /**
* (but not all) components may be omitted. The last component indicat * Trailing unit is arc minutes.
or * @hideinitializer
* may be omitted and is assumed to be tbe next smallest unit (thus 33d ********************************************************************
10 **/
* is interpreted as 33d10'). The final component may be a decimal MINUTE = 1,
* fraction but the non-final components must be integers. The integer /**
* parts of the minutes and seconds components must be less than 60. A * Trailing unit is arc seconds.
* single leading sign is permitted. A hemisphere designator (N, E, W, * @hideinitializer
S) ********************************************************************
* may be added to tbe beginning or end of the string. The result is **/
* multiplied by the implied signed of the hemisphere designator (negat SECOND = 2,
ive };
* for S and W). In addition \e flag is used to indicate whether such
a /**
* designator was found and whether it implies that the angle is a lati * Convert a string in DMS to an angle.
tude *
* (N or S) or longitude (E or W). Throws an error on a malformed stri * @param[in] dms string input.
ng. * @param[out] ind a DMS::flag value indicating the presence of a
* No check is performed on the range of the result. * hemisphere indicator.
* @return angle (degrees).
*
* 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
all)
* components may be omitted. The last component indicator may be omit
ted
* 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
* non-final components must be integers. The integer parts of the min
utes
* and seconds components must be less than 60. A single leading sign
is
* 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
* 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
* 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.
********************************************************************** / ********************************************************************** /
static Math::real Decode(const std::string& dms, flag& ind); static Math::real Decode(const std::string& dms, flag& ind);
/** /**
* Convert real degrees \e d, minutes \e m, and seconds \e s, to degree * Convert DMS to an angle.
s. *
* @param[in] d degrees.
* @param[in] m arc minutes.
* @param[in] s arc seconds.
* @return angle (degress)
*
* 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, so
* -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 \e str to a real number. * Convert a string to a real number.
*
* @param[in] str string input.
* @return decoded number.
********************************************************************** / ********************************************************************** /
static Math::real Decode(const std::string& str); static Math::real Decode(const std::string& str);
/** /**
* Convert two strings \e dmsa and \e dmsb to a latitude, \e lat, and * Convert a pait of strings to latitude and longitude.
* longitude, \e lon. By default, the \e lat (resp., \e lon) is assign *
ed * @param[in] dmsa first string.
* to the results of decoding \e dmsa (resp., \e dmsb). However this i * @param[in] dmsb second string.
s * @param[out] lat latitude.
* overridden if either \e dmsa or \e dmsb contain a latitude or longit * @param[out] lon longitude.
ude *
* hemisphere designator (N, S, E, W). Throws an error if the decoded * By default, the \e lat (resp., \e lon) is assigned to the results of
* numbers are out of the ranges [-90<sup>o</sup>, 90<sup>o</sup>] for * decoding \e dmsa (resp., \e dmsb). However this is overridden if ei
* latitude and [-180<sup>o</sup>, 360<sup>o</sup>] for longitude and, ther
in * \e dmsa or \e dmsb contain a latitude or longitude hemisphere design
* which case \e lat and \e lon are unchanged. Finally the longitude i ator
s * (N, S, E, W). Throws an error if the decoded numbers are out of the
* reduced to the range [-180<sup>o</sup>, 180<sup>o</sup>). * ranges [-90<sup>o</sup>, 90<sup>o</sup>] for latitude and
* [-180<sup>o</sup>, 360<sup>o</sup>] for longitude and, in which case
\e
* lat and \e lon are unchanged. Finally the longitude is reduced to t
he
* range [-180<sup>o</sup>, 180<sup>o</sup>).
********************************************************************** / ********************************************************************** /
static void DecodeLatLon(const std::string& dmsa, const std::string& dm sb, static void DecodeLatLon(const std::string& dmsa, const std::string& dm sb,
real& lat, real& lon); real& lat, real& lon);
/** /**
* Convert a string \e angstr to an angle in degrees. No hemisphere * Convert a string to an angle in degrees.
* designator is allowed and no check is done on the range of the resul *
t. * @param[in] angstr input string.
* @return angle (degrees)
*
* No hemisphere designator is allowed and no check is done on the rang
e of
* the result.
********************************************************************** / ********************************************************************** /
static Math::real DecodeAngle(const std::string& angstr); static Math::real DecodeAngle(const std::string& angstr);
/** /**
* Convert a string \e azistr to an azimuth in degrees. A hemisphere * Convert a string to an azimuth in degrees.
* designator E/W can be used; the result is multiplied by -1 if W is *
* present. Throws an error if the result is out of the range * @param[in] azistr input string.
* @return azimuth (degrees)
*
* A hemisphere designator E/W can be used; the result is multiplied by
-1
* if W is present. Throws an error if the result is out of the range
* [-180<sup>o</sup>, 360<sup>o</sup>]. Finally the azimuth is reduced to * [-180<sup>o</sup>, 360<sup>o</sup>]. Finally the azimuth is reduced to
* 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 \e angle (in degrees) into a DMS string. \e trailing indica * Convert angle (in degrees) into a DMS string.
tes *
* the least significant component of the string (and this component is * @param[in] angle input angle (degrees)
* given as a decimal number if necessary). \e prec indicates the numb * @param[in] trailing DMS::component value indicating the trailing uni
er ts
* of digits after the decimal point for the trailing component. \e fl * on the string and this is given as a decimal number if necessary.
ag * @param[in] prec the number of digits after the decimal point for the
* indicates additional formating as follows * trailing component.
* - flag == DMS::NONE, signed result no leading zeros on degrees excep * @param[in] ind DMS::flag value indicated additional formatting.
t in * @return formatting string
*
* The interpretation of \e ind is as follows:
* - ind == DMS::NONE, signed result no leading zeros on degrees except
in
* the units place, e.g., -8d03'. * the units place, e.g., -8d03'.
* - flag == DMS::LATITUDE, trailing N or S hemisphere designator, no s ign, * - ind == DMS::LATITUDE, trailing N or S hemisphere designator, no si gn,
* pad degress to 2 digits, e.g., 08d03'S. * pad degress to 2 digits, e.g., 08d03'S.
* - flag == DMS::LONGITUDE, trailing E or W hemisphere designator, no * - ind == DMS::LONGITUDE, trailing E or W hemisphere designator, no
* sign, pad degress to 3 digits, e.g., 008d03'W. * sign, pad degress to 3 digits, e.g., 008d03'W.
* - flag == DMS::AZIMUTH, convert to the range [0, 360<sup>o</sup>), n o * - 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 iven * The integer parts of the minutes and seconds components are always g iven
* 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 \e angle into a DMS string selecting the trailing component * Convert angle into a DMS string selecting the trailing component
* based on \e prec. \e prec indicates the precision relative to 1 deg * based on the precision.
ree, *
* e.g., \e prec = 3 gives a result accurate to 0.1' and \e prec = 4 gi * @param[in] angle input angle (degrees)
ves * @param[in] prec the precision relative to 1 degree.
* a result accurate to 1&quot;. If \e ind is DMS::NUMBER, then merely * @param[in] ind DMS::flag value indicated additional formatting.
* format \e angle as a number in fixed format with precision \e prec. * @return formatting string
*
* \e prec indicates the precision relative to 1 degree, e.g., \e prec
= 3
* gives a result accurate to 0.1' and \e prec = 4 gives a result accur
ate
* to 1&quot;. \e ind is interpreted as in DMS::Encode with the additi
onal
* facility at DMS::NUMBER treats \e angle a number in fixed format wit
h
* precision \e prec.
********************************************************************** / ********************************************************************** /
static std::string Encode(real angle, unsigned prec, flag ind = NONE) { static std::string Encode(real angle, unsigned prec, flag ind = NONE) {
if (ind == NUMBER) { if (ind == NUMBER) {
std::ostringstream s; std::ostringstream s;
s << std::fixed << std::setprecision(prec) << angle; s << std::fixed << std::setprecision(prec) << angle;
return s.str(); return s.str();
} else } else
return Encode(angle, return Encode(angle,
prec < 2 ? DEGREE : (prec < 4 ? MINUTE : SECOND), prec < 2 ? DEGREE : (prec < 4 ? MINUTE : SECOND),
prec < 2 ? prec : (prec < 4 ? prec - 2 : prec - 4), prec < 2 ? prec : (prec < 4 ? prec - 2 : prec - 4),
ind); ind);
} }
/** /**
* Split angle, \e ang, into degrees, \e d, and minutes \e m. * Split angle into degrees and mintues
*
* @param[in] ang angle (degrees)
* @param[out] d degrees (an integer returned as a real)
* @param[out] m arc minutes.
********************************************************************** / ********************************************************************** /
static void Encode(real ang, real& d, real& m) throw() { static void Encode(real ang, real& d, real& m) throw() {
d = int(ang); m = 60 * (ang - d); d = int(ang); m = 60 * (ang - d);
} }
/** /**
* Split angle, \e ang, into degrees, \e d, minutes, \e m, and seconds * Split angle into degrees and mintues and seconds.
\e *
* s. * @param[in] ang angle (degrees)
* @param[out] d degrees (an integer returned as a real)
* @param[out] m arc minutes (an integer returned as a real)
* @param[out] s arc seconds.
********************************************************************** / ********************************************************************** /
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
 End of changes. 16 change blocks. 
82 lines changed or deleted 179 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) <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 6838 2010-06-22 21:26:37Z karney $" #define GEOGRAPHICLIB_ELLIPTICFUNCTION_HPP "$Id: EllipticFunction.hpp 6866 2010-09-11 02:15:29Z karney $"
#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
skipping to change at line 57 skipping to change at line 57
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:
/** /**
* Constructor with parameter \e m which must lie in [0, 1]. (No check * Constructor.
ing *
* is done.) * @param[in] m the parameter which must lie in [0, 1]. (No checking
* is done.)
********************************************************************** / ********************************************************************** /
explicit EllipticFunction(real m) throw(); explicit EllipticFunction(real m) throw();
/** /**
* The parameter \e m. * @return the parameter \e m.
********************************************************************** / ********************************************************************** /
Math::real m() const throw() { return _m; } Math::real m() const throw() { return _m; }
/** /**
* The complementary parameter \e m' = (1 - \e m). * @return the complementary parameter \e m' = (1 - \e m).
********************************************************************** / ********************************************************************** /
Math::real m1() const throw() { return _m1; } Math::real m1() const throw() { return _m1; }
/** /**
* The complete integral of first kind, \e K(\e m). * @return the complete integral of first kind, \e K(\e m).
********************************************************************** / ********************************************************************** /
Math::real K() const throw() { _init || Init(); return _kc; } Math::real K() const throw() { _init || Init(); return _kc; }
/** /**
* The complete integral of second kind, \e E(\e m). * @return the complete integral of second kind, \e E(\e m).
********************************************************************** / ********************************************************************** /
Math::real E() const throw() { _init || Init(); return _ec; } Math::real E() const throw() { _init || Init(); return _ec; }
/** /**
* The difference \e K(\e m) - \e E(\e m) (which can be computed direct * @return the difference \e K(\e m) - \e E(\e m) (which can be compute
ly). d
* directly).
********************************************************************** / ********************************************************************** /
Math::real KE() const throw() { _init || Init(); return _kec; } Math::real KE() const throw() { _init || Init(); return _kec; }
/** /**
* The Jacobi elliptic functions sn(<i>x</i>|<i>m</i>), * The Jacobi elliptic functions.
* cn(<i>x</i>|<i>m</i>), and dn(<i>x</i>|<i>m</i>) with argument \e x. *
* The results are returned in \e sn, \e cn, and \e dn. * @param[in] x the argument.
* @param[out] sn sn(<i>x</i>|<i>m</i>).
* @param[out] cn cn(<i>x</i>|<i>m</i>).
* @param[out] dn dn(<i>x</i>|<i>m</i>).
********************************************************************** / ********************************************************************** /
void sncndn(real x, real& sn, real& cn, real& dn) const throw(); void sncndn(real x, real& sn, real& cn, real& dn) const throw();
/** /**
* The incomplete integral of the second kind = int sqrt(1 - \e m * The incomplete integral of the second kind.
* sin<sup>2</sup>(\e phi)) \e dphi. *
* @param[in] phi
* @return int sqrt(1 - \e m sin<sup>2</sup>(\e phi)) \e dphi.
********************************************************************** / ********************************************************************** /
Math::real E(real phi) const throw(); Math::real E(real phi) const throw();
/** /**
* The incomplete integral of the second kind = int dn(\e w)<sup>2</sup * The incomplete integral of the second kind in terms of Jacobi ellipt
> \e ic
* dw (A+S 17.2.10). Instead of specifying the ampltiude \e phi, we * functions
* provide \e sn = sin(\e phi), \e cn = cos(\e phi), \e dn = sqrt(1 - \ *
e m * @param[in] sn
* sin<sup>2</sup>(\e phi)). * @param[in] cn
* @param[in] dn
* @return int dn(\e w)<sup>2</sup> \e dw (A+S 17.2.10).
*
* 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)).
********************************************************************** / ********************************************************************** /
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
 End of changes. 10 change blocks. 
21 lines changed or deleted 35 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 6785 2010-01-05 22: 15:42Z karney $" #define GEOGRAPHICLIB_GEOCOORDS_HPP "$Id: GeoCoords.hpp 6867 2010-09-11 13: 04:26Z 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 64 skipping to change at line 64
_alt_northing = _northing; _alt_northing = _northing;
_alt_gamma = _gamma; _alt_gamma = _gamma;
_alt_k = _k; _alt_k = _k;
_alt_zone = _zone; _alt_zone = _zone;
} }
void UTMUPSString(int zone, real easting, real northing, void UTMUPSString(int zone, real easting, real northing,
int prec, std::string& utm) const; int prec, std::string& utm) const;
void FixHemisphere(); void FixHemisphere();
public: public:
/** \name Initializing the GeoCoords object
**********************************************************************
/
///@{
/** /**
* The default contructor is equivalent to \e latitude = 90<sup>o</sup> , \e * The default contructor is equivalent to \e latitude = 90<sup>o</sup> , \e
* longitude = 0<sup>o</sup>. * longitude = 0<sup>o</sup>.
********************************************************************** / ********************************************************************** /
GeoCoords() throw() GeoCoords() throw()
// This is the N pole // This is the N pole
: _lat(90) : _lat(90)
, _long(0) , _long(0)
, _easting(2000000) , _easting(2000000)
, _northing(2000000) , _northing(2000000)
, _northp(true) , _northp(true)
, _zone(0) , _zone(0)
{ CopyToAlt(); } { CopyToAlt(); }
/** /**
* Construct from a string.
*
* @param[in] s 1-element, 2-element, or 3-element string representatio
n of
* the position.
* @param[in] centerp governs the interpretation of MGRS coordinates (s
ee
* below).
*
* Parse as a string and interpret it as a geographic position. The in put * Parse as a string and interpret it as a geographic position. The in put
* string is broken into space (or comma) separated pieces and Basic * string is broken into space (or comma) separated pieces and Basic
* decision on which format is based on number of components * decision on which format is based on number of components
* -# MGRS * -# MGRS
* -# "Lat Long" or "Long Lat" * -# "Lat Long" or "Long Lat"
* -# "Zone Easting Northing" or "Easting Northing Zone" * -# "Zone Easting Northing" or "Easting Northing Zone"
* *
* The following inputs are approximately the same (Ar Ramadi Bridge, I raq) * The following inputs are approximately the same (Ar Ramadi Bridge, I raq)
* - Latitude and Longitude * - Latitude and Longitude
* - 33.44 43.27 * - 33.44 43.27
skipping to change at line 147 skipping to change at line 157
* . * .
* Otherwise, the "south-west" corner of the square is used, i.e., * Otherwise, the "south-west" corner of the square is used, i.e.,
* - 38SMB = 38N 400000 3600000 * - 38SMB = 38N 400000 3600000
* - 38SMB4484 = 38N 444000 3684000 * - 38SMB4484 = 38N 444000 3684000
* - 38SMB44148470 = 38N 444140 3684700 * - 38SMB44148470 = 38N 444140 3684700
********************************************************************** / ********************************************************************** /
explicit GeoCoords(const std::string& s, bool centerp = true) explicit GeoCoords(const std::string& s, bool centerp = true)
{ Reset(s, centerp); } { Reset(s, centerp); }
/** /**
* Specify the location in terms of \e latitude (degrees) and \e longit * Construct from geographic coordinates.
ude *
* (degrees). Use \e zone to force the UTM/UPS representation to use a * @param[in] latitude (degrees).
* specified zone using the rules given in UTMUPS::zonespec. * @param[in] longitude (degrees).
* @param[in] zone if specified, force the UTM/UPS representation to us
e a
* specified zone using the rules given in UTMUPS::zonespec.
********************************************************************** / ********************************************************************** /
GeoCoords(real latitude, real longitude, int zone = UTMUPS::STANDARD) { GeoCoords(real latitude, real longitude, int zone = UTMUPS::STANDARD) {
Reset(latitude, longitude, zone); Reset(latitude, longitude, zone);
} }
/** /**
* Specify the location in terms of UPS/UPS \e zone (zero means UPS), * Construct from UTM/UPS coordinates.
* hemisphere \e northp (false means south, true means north), \e easti *
ng * @param[in] zone UTM zone (zero means UPS).
* (meters) and \e northing (meters). * @param[in] northp hemisphere (true means north, false means south).
* @param[in] easting (meters).
* @param[in] northing (meters).
********************************************************************** / ********************************************************************** /
GeoCoords(int zone, bool northp, real easting, real northing) { GeoCoords(int zone, bool northp, real easting, real northing) {
Reset(zone, northp, easting, northing); Reset(zone, northp, easting, northing);
} }
/** /**
* Reset the location as a 1-element, 2-element, or 3-element string. * Reset the location from a string. See
See * GeoCoords(const std::string& s, bool centerp).
* GeoCoords(const string& s, bool centerp).
********************************************************************** / ********************************************************************** /
void Reset(const std::string& s, bool centerp = true); void Reset(const std::string& s, bool centerp = true);
/** /**
* Reset the location in terms of \e latitude and \e longitude. See * Reset the location in terms of geographic coordinates. See
* GeoCoords(real latitude, real longitude, int zone). * GeoCoords(real latitude, real longitude, int zone).
********************************************************************** / ********************************************************************** /
void Reset(real latitude, real longitude, int zone = UTMUPS::STANDARD) { void Reset(real latitude, real longitude, int zone = UTMUPS::STANDARD) {
UTMUPS::Forward(latitude, longitude, UTMUPS::Forward(latitude, longitude,
_zone, _northp, _easting, _northing, _gamma, _k, _zone, _northp, _easting, _northing, _gamma, _k,
zone); zone);
_lat = latitude; _lat = latitude;
_long = longitude; _long = longitude;
if (_long >= 180) if (_long >= 180)
_long -= 360; _long -= 360;
CopyToAlt(); CopyToAlt();
} }
/** /**
* Reset the location in terms of UPS/UPS \e zone, hemisphere \e northp * Reset the location in terms of UPS/UPS coordinates. See
, \e * GeoCoords(int zone, bool northp, real easting, real northing).
* easting, and \e northing. See GeoCoords(int zone, bool northp,
* real easting, real northing).
********************************************************************** / ********************************************************************** /
void Reset(int zone, bool northp, real easting, real northing) { void Reset(int zone, bool northp, real easting, real northing) {
UTMUPS::Reverse(zone, northp, easting, northing, UTMUPS::Reverse(zone, northp, easting, northing,
_lat, _long, _gamma, _k); _lat, _long, _gamma, _k);
_zone = zone; _zone = zone;
_northp = northp; _northp = northp;
_easting = easting; _easting = easting;
_northing = northing; _northing = northing;
FixHemisphere(); FixHemisphere();
CopyToAlt(); CopyToAlt();
} }
///@}
/** \name Querying the GeoCoords object
**********************************************************************
/
///@{
/** /**
* Return latitude (degrees) * @return latitude (degrees)
********************************************************************** / ********************************************************************** /
Math::real Latitude() const throw() { return _lat; } Math::real Latitude() const throw() { return _lat; }
/** /**
* Return longitude (degrees) * @return longitude (degrees)
********************************************************************** / ********************************************************************** /
Math::real Longitude() const throw() { return _long; } Math::real Longitude() const throw() { return _long; }
/** /**
* Return easting (meters) * @return easting (meters)
********************************************************************** / ********************************************************************** /
Math::real Easting() const throw() { return _easting; } Math::real Easting() const throw() { return _easting; }
/** /**
* Return northing (meters) * @return northing (meters)
********************************************************************** / ********************************************************************** /
Math::real Northing() const throw() { return _northing; } Math::real Northing() const throw() { return _northing; }
/** /**
* Return meridian convergence (degrees) for the UTM/UPS projection. * @return meridian convergence (degrees) for the UTM/UPS projection.
********************************************************************** / ********************************************************************** /
Math::real Convergence() const throw() { return _gamma; } Math::real Convergence() const throw() { return _gamma; }
/** /**
* Return scale for the UTM/UPS projection. * @return scale for the UTM/UPS projection.
********************************************************************** / ********************************************************************** /
Math::real Scale() const throw() { return _k; } Math::real Scale() const throw() { return _k; }
/** /**
* Return hemisphere (false means south, true means north). * @return hemisphere (false means south, true means north).
********************************************************************** / ********************************************************************** /
bool Northp() const throw() { return _northp; } bool Northp() const throw() { return _northp; }
/** /**
* Return hemisphere letter N or S. * @return hemisphere letter N or S.
********************************************************************** / ********************************************************************** /
char Hemisphere() const throw() { return _northp ? 'N' : 'S'; } char Hemisphere() const throw() { return _northp ? 'N' : 'S'; }
/** /**
* Return the zone corresponding to the input (return 0 for UPS). * @return the zone corresponding to the input (return 0 for UPS).
********************************************************************** / ********************************************************************** /
int Zone() const throw() { return _zone; } int Zone() const throw() { return _zone; }
///@}
/** \name Setting and querying the alternate zone
**********************************************************************
/
///@{
/** /**
* Use zone number, \e zone, for the alternate representation. See * Specify alternate zone number.
* UTMUPS::zonespec for more information on the interpretation of \e zo *
ne. * @param[in] zone zone number for the alternate representation.
* Note that \e zone == UTMUPS::STANDARD (the default) use the standard *
UPS * See UTMUPS::zonespec for more information on the interpretation of \
* or UTM zone, UTMUPS::MATCH does nothing retaining the existing alter e
nate * zone. Note that \e zone == UTMUPS::STANDARD (the default) use the
* representation. Before this is called the alternate zone is the inp * standard UPS or UTM zone, UTMUPS::MATCH does nothing retaining the
ut * existing alternate representation. Before this is called the altern
* zone. ate
* zone is the input zone.
********************************************************************** / ********************************************************************** /
void SetAltZone(int zone = UTMUPS::STANDARD) const { void SetAltZone(int zone = UTMUPS::STANDARD) const {
if (zone == UTMUPS::MATCH) if (zone == UTMUPS::MATCH)
return; return;
zone = UTMUPS::StandardZone(_lat, _long, zone); zone = UTMUPS::StandardZone(_lat, _long, zone);
if (zone == _zone) if (zone == _zone)
CopyToAlt(); CopyToAlt();
else { else {
bool northp; bool northp;
UTMUPS::Forward(_lat, _long, UTMUPS::Forward(_lat, _long,
_alt_zone, northp, _alt_zone, northp,
_alt_easting, _alt_northing, _alt_gamma, _alt_k, _alt_easting, _alt_northing, _alt_gamma, _alt_k,
zone); zone);
} }
} }
/** /**
* Returns the current alternate zone (return 0 for UPS). * @return current alternate zone (return 0 for UPS).
********************************************************************** / ********************************************************************** /
int AltZone() const throw() { return _alt_zone; } int AltZone() const throw() { return _alt_zone; }
/** /**
* Return easting (meters) for alternate zone. * @return easting (meters) for alternate zone.
********************************************************************** / ********************************************************************** /
Math::real AltEasting() const throw() { return _alt_easting; } Math::real AltEasting() const throw() { return _alt_easting; }
/** /**
* Return northing (meters) for alternate zone. * @return northing (meters) for alternate zone.
********************************************************************** / ********************************************************************** /
Math::real AltNorthing() const throw() { return _alt_northing; } Math::real AltNorthing() const throw() { return _alt_northing; }
/** /**
* Return meridian convergence (degrees) for altermate zone. * @return meridian convergence (degrees) for altermate zone.
********************************************************************** / ********************************************************************** /
Math::real AltConvergence() const throw() { return _alt_gamma; } Math::real AltConvergence() const throw() { return _alt_gamma; }
/** /**
* Return scale for altermate zone. * @return scale for altermate zone.
********************************************************************** / ********************************************************************** /
Math::real AltScale() const throw() { return _alt_k; } Math::real AltScale() const throw() { return _alt_k; }
///@}
/** \name String representations of the GeoCoords object
**********************************************************************
/
///@{
/** /**
* Return string with latitude and longitude as signed decimal degrees. * String representation with latitude and longitude as signed decimal
* Precision \e prec specifies accuracy of representation as follows: * degrees.
*
* @param[in] prec precision (relative to about 1m).
* @return decimal latitude/longitude string representation.
*
* Precision specifies accuracy of representation as follows:
* - prec = -5 (min), 1d * - prec = -5 (min), 1d
* - prec = 0, 10<sup>-5</sup>d (about 1m) * - prec = 0, 10<sup>-5</sup>d (about 1m)
* - prec = 3, 10<sup>-8</sup>d * - prec = 3, 10<sup>-8</sup>d
* - prec = 9 (max), 10<sup>-14</sup>d * - prec = 9 (max), 10<sup>-14</sup>d
********************************************************************** / ********************************************************************** /
std::string GeoRepresentation(int prec = 0) const; std::string GeoRepresentation(int prec = 0) const;
/** /**
* Return string with latitude and longitude as degrees, minutes, secon * String representation with latitude and longitude as degrees, minute
ds, s,
* and hemisphere. Precision \e prec specifies accuracy of representat * seconds, and hemisphere.
ion *
* as follows: * @param[in] prec precision (relative to about 1m)
* @return DMS latitude/longitude string representation.
*
* Precision specifies accuracy of representation as follows:
* - prec = -5 (min), 1d * - prec = -5 (min), 1d
* - prec = -4, 0.1d * - prec = -4, 0.1d
* - prec = -3, 1' * - prec = -3, 1'
* - prec = -2, 0.1' * - prec = -2, 0.1'
* - prec = -1, 1" * - prec = -1, 1"
* - prec = 0, 0.1" (about 3m) * - prec = 0, 0.1" (about 3m)
* - prec = 1, 0.01" * - prec = 1, 0.01"
* - prec = 10 (max), 10<sup>-11</sup>" * - prec = 10 (max), 10<sup>-11</sup>"
********************************************************************** / ********************************************************************** /
std::string DMSRepresentation(int prec = 0) const; std::string DMSRepresentation(int prec = 0) const;
/** /**
* Return MGRS string. This gives the coordinates of the enclosing gri * MGRS string.
d *
* square with size given by the precision \e prec. Thus 38N 444180 * @param[in] prec precision (relative to about 1m).
* 3684790 converted to a MGRS coordinate at precision -2 (100m) is * @return MGRS string.
* 38SMB441847 and not 38SMB442848. Precision \e prec specifies the *
* precision of the MSGRS string as follows: * This gives the coordinates of the enclosing grid square with size gi
ven
* by the precision. Thus 38N 444180 3684790 converted to a MGRS
* coordinate at precision -2 (100m) is 38SMB441847 and not 38SMB442848
.
* \e prec specifies the precision of the MSGRS string as follows:
* - prec = -5 (min), 100km * - prec = -5 (min), 100km
* - prec = -4, 10km * - prec = -4, 10km
* - prec = -3, 1km * - prec = -3, 1km
* - prec = -2, 100m * - prec = -2, 100m
* - prec = -1, 10m * - prec = -1, 10m
* - prec = 0, 1m * - prec = 0, 1m
* - prec = 1, 0.1m * - prec = 1, 0.1m
* - prec = 6 (max), 1um * - prec = 6 (max), 1um
********************************************************************** / ********************************************************************** /
std::string MGRSRepresentation(int prec = 0) const; std::string MGRSRepresentation(int prec = 0) const;
/** /**
* Return string consisting of UTM/UPS zone designator, easting, and * UTM/UPS string.
* northing, Precision \e prec specifies accuracy of representation *
* as follows: * @param[in] prec precision (relative to about 1m)
* @return UTM/UPS string representation: zone designator, easting, and
* northing.
*
* Precision specifies accuracy of representation as follows:
* - prec = -5 (min), 100km * - prec = -5 (min), 100km
* - prec = -3, 1km * - prec = -3, 1km
* - prec = 0, 1m * - prec = 0, 1m
* - prec = 3, 1mm * - prec = 3, 1mm
* - prec = 6, 1um * - prec = 6, 1um
* - prec = 9 (max), 1nm * - prec = 9 (max), 1nm
********************************************************************** / ********************************************************************** /
std::string UTMUPSRepresentation(int prec = 0) const; std::string UTMUPSRepresentation(int prec = 0) const;
/** /**
* Return MGRS string using alternative zone. See MGRSRepresentation f * MGRS string for the alternate zone. See GeoCoords::MGRSRepresentati
or on.
* the interpretation of \e prec.
********************************************************************** / ********************************************************************** /
std::string AltMGRSRepresentation(int prec = 0) const; std::string AltMGRSRepresentation(int prec = 0) const;
/** /**
* Return string consisting of alternate UTM/UPS zone designator, easti * UTM/UPS string for the alternate zone. See
ng, * GeoCoords::UTMUPSRepresentation.
* and northing. See UTMUPSRepresentation for the interpretation of \e
* prec.
********************************************************************** / ********************************************************************** /
std::string AltUTMUPSRepresentation(int prec = 0) const; std::string AltUTMUPSRepresentation(int prec = 0) const;
///@}
/** \name Inspector functions
**********************************************************************
/
///@{
/** /**
* The major radius of the ellipsoid (meters). This is the value for t * @return \e a the equatorial radius of the WGS84 ellipsoid (meters).
he *
* WGS84 ellipsoid because the UTM and UPS projections are based on thi * (The WGS84 values is returned because the UTM and UPS projections ar
s e
* ellipsoid. * based on this ellipsoid.)
********************************************************************** / ********************************************************************** /
Math::real MajorRadius() const throw() { return UTMUPS::MajorRadius(); } Math::real MajorRadius() const throw() { return UTMUPS::MajorRadius(); }
/** /**
* The inverse flattening of the ellipsoid. This is the value for the * @return \e r the inverse flattening of the WGS84 ellipsoid.
* WGS84 ellipsoid because the UTM and UPS projections are based on thi *
s * (The WGS84 values is returned because the UTM and UPS projections ar
* ellipsoid. e
* 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. 39 change blocks. 
73 lines changed or deleted 126 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 6827 2010-05-20 1 9:56:18Z karney $" #define GEOGRAPHICLIB_GEOCENTRIC_HPP "$Id: Geocentric.hpp 6867 2010-09-11 1 3:04:26Z 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 53 skipping to change at line 53
**********************************************************************/ **********************************************************************/
class Geocentric { class Geocentric {
private: private:
typedef Math::real real; typedef Math::real real;
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; } static inline real sq(real x) throw() { return x * x; }
public: public:
/** /**
* Constructor for a ellipsoid radius \e a (meters) and reciprocal * Constructor for a ellipsoid with
* flattening \e r. Setting \e r = 0 implies \e r = inf or flattening *
= 0 * @param[in] a equatorial radius (meters)
* (i.e., a sphere). Negative \e r indicates a prolate spheroid. An * @param[in] r reciprocal flattening. Setting \e r = 0 implies \e r =
* exception is thrown if \e a is not positive. inf
* or flattening = 0 (i.e., a sphere). Negative \e r indicates a pro
late
* ellipsoid.
*
* An exception is thrown if either of the axes of the ellipsoid is
* non-positive.
********************************************************************** / ********************************************************************** /
Geocentric(real a, real r); Geocentric(real a, real r);
/** /**
* Convert from geodetic coordinates \e lat, \e lon (degrees), \e h * Convert from geodetic to geocentric coordinates.
* (meters) to geocentric coordinates \e x, \e y, \e z (meters). \e la *
t * @param[in] lat latitude of point (degrees).
* should be in the range [-90, 90]; \e lon and \e lon0 should be in th * @param[in] lon longitude of point (degrees).
e * @param[in] h height of point above the ellipsoid (meters).
* range [-180, 360]. * @param[out] x geocentric coordinate (meters).
* @param[out] y geocentric coordinate (meters).
* @param[out] z geocentric coordinate (meters).
*
* \e lat should be in the range [-90, 90]; \e lon and \e lon0 should b
e in
* the range [-180, 360].
********************************************************************** / ********************************************************************** /
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)
const throw(); const throw();
/** /**
* Convert from geocentric coordinates \e x, \e y, \e z (meters) to * Convert from geocentric to geodetic to coordinates.
* geodetic \e lat, \e lon (degrees), \e h (meters). In general there *
are * @param[in] x geocentric coordinate (meters).
* multiple solutions and the result which minimizes the absolute value * @param[in] y geocentric coordinate (meters).
of * @param[in] z geocentric coordinate (meters).
* @param[out] lat latitude of point (degrees).
* @param[out] lon longitude of point (degrees).
* @param[out] h height of point above the ellipsoid (meters).
*
* In general there are multiple solutions and the result which maximiz
es
* \e h is returned. If there are still multiple solutions with differ ent * \e h is returned. If there are still multiple solutions with differ ent
* latitutes (applies only if \e z = 0), then the solution with \e lat > 0 * latitutes (applies only if \e z = 0), then the solution with \e lat > 0
* is returned. If there are still multiple solutions with different * is returned. If there are still multiple solutions with different
* longitudes (applies only if \e x = \e y = 0) then \e lon = 0 is * longitudes (applies only if \e x = \e y = 0) then \e lon = 0 is
* returned. The value of \e h returned satisfies \e h >= - \e a (1 - \e * returned. The value of \e h returned satisfies \e h >= - \e a (1 - \e
* e<sup>2</sup>) / sqrt(1 - \e e<sup>2</sup> sin<sup>2</sup>\e lat). The * e<sup>2</sup>) / sqrt(1 - \e e<sup>2</sup> sin<sup>2</sup>\e lat). The
* value of \e lon returned is in the range [-180, 180). * value of \e lon returned is in the range [-180, 180).
********************************************************************** / ********************************************************************** /
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)
const throw(); const throw();
/** \name Inspector functions
**********************************************************************
/
///@{
/** /**
* The major radius of the ellipsoid (meters). This is that value of \ * @return \e a the equatorial radius of the ellipsoid (meters). This
e a is
* used in the constructor. * the value used in the constructor.
********************************************************************** / ********************************************************************** /
Math::real MajorRadius() const throw() { return _a; } Math::real MajorRadius() const throw() { return _a; }
/** /**
* The inverse flattening of the ellipsoid. This is that value of \e r * @return \e r the inverse flattening of the ellipsoid. This is the
* used in the constructor. A value of 0 is returned for a sphere * value used in the constructor. A value of 0 is returned for a sph
* (infinite inverse flattening). ere
* (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; const static Geocentric WGS84;
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif #endif
 End of changes. 8 change blocks. 
23 lines changed or deleted 47 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 6827 2010-05-20 19:56 :18Z karney $" #define GEOGRAPHICLIB_GEODESIC_HPP "$Id: Geodesic.hpp 6867 2010-09-11 13:04 :26Z 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 36 skipping to change at line 36
/** /**
* \brief %Geodesic calculations * \brief %Geodesic calculations
* *
* The shortest path between two points on a ellipsoid at (\e lat1, \e lo n1) * The shortest path between two points on a ellipsoid at (\e lat1, \e lo n1)
* and (\e lat2, \e lon2) is called the geodesic. Its length is \e s12 a nd * and (\e lat2, \e lon2) is called the geodesic. Its length is \e s12 a nd
* the geodesic from point 1 to point 2 has azimuths \e azi1 and \e azi2 at * the geodesic from point 1 to point 2 has azimuths \e azi1 and \e azi2 at
* the two end points. (The azimuth is the heading measured clockwise fr om * the two end points. (The azimuth is the heading measured clockwise fr om
* north. \e azi2 is the "forward" azimuth, i.e., the heading that takes you * north. \e azi2 is the "forward" azimuth, i.e., the heading that takes you
* beyond point 2 not back to point 1.) * beyond point 2 not back to point 1.)
* *
* If we fix the first point and increase \e s12 by \e ds12, then the sec
ond
* point is displaced \e ds12 in the direction \e azi2. Similarly we
* increase \e azi1 by \e dazi1 (radians), the the second point is displa
ced
* \e m12 \e dazi1 in the direction \e azi2 + 90<sup>o</sup>. The quanti
ty
* \e m12 is called the "reduced length" and is symmetric under interchan
ge
* of the two points. On a flat surface, he have \e m12 = \e s12. The r
atio
* \e s12/\e m12 gives the azimuthal scale for an azimuthal equidistant
* projection.
*
* Given \e lat1, \e lon1, \e azi1, and \e s12, we can determine \e lat2, \e * Given \e lat1, \e lon1, \e azi1, and \e s12, we can determine \e lat2, \e
* lon2, \e azi2, \e m12. This is the \e direct geodesic problem. (If \ * lon2, and \e azi2. This is the \e direct geodesic problem and its
e * solution is given by the function Geodesic::Direct. (If \e s12 is
* s12 is sufficiently large that the geodesic wraps more than halfway ar * sufficiently large that the geodesic wraps more than halfway around th
ound e
* the earth, there will be another geodesic between the points with a * earth, there will be another geodesic between the points with a smalle
* smaller \e s12.) r \e
* s12.)
* *
* Given \e lat1, \e lon1, \e lat2, and \e lon2, we can determine \e azi1 , \e * Given \e lat1, \e lon1, \e lat2, and \e lon2, we can determine \e azi1 , \e
* azi2, \e s12, \e m12. This is the \e inverse geodesic problem. Usual * azi2, and \e s12. This is the \e inverse geodesic problem, whose solu
ly, tion
* the solution to the inverse problem is unique. In cases where there a * is given by Geodesic::Inverse. Usually, the solution to the inverse
re * problem is unique. In cases where there are muliple solutions (all wi
* muliple solutions (all with the same \e s12, of course), all the solut th
ions * the same \e s12, of course), all the solutions can be easily generated
* can be easily generated once a particular solution is provided. * once a particular solution is provided.
*
* The standard way of specifying the direct problem is the specify the
* distance \e s12 to the second point. However it is sometimes useful
* instead to specify the the arc length \e a12 (in degrees) on the auxil
iary
* sphere. This is a mathematical construct used in solving the geodesic
* problems. The solution of the direct problem in this form is provide
by
* Geodesic::ArcDirect. An arc length in excess of 180<sup>o</sup> indic
ates
* that the geodesic is not a shortest path. In addition, the arc length
* between an equatorial crossing and the next extremum of latitude for a
* geodesic is 90<sup>o</sup>.
*
* This class can also calculate several other quantities related to
* geodesics. These are:
* - <i>reduced length</i>. If we fix the first point and increase \e az
i1
* by \e dazi1 (radians), the the second point is displaced \e m12 \e d
azi1
* in the direction \e azi2 + 90<sup>o</sup>. The quantity \e m12 is
* called the "reduced length" and is symmetric under interchange of th
e
* two points. On a flat surface, we have \e m12 = \e s12. The ratio
\e
* s12/\e m12 gives the azimuthal scale for an azimuthal equidistant
* projection.
* - <i>geodesic scale</i>. Consider a reference geodesic and a second
* geodesic parallel to this one at point 1 and separated by a small
* distance \e dt. The separation of the two geodesics at point 2 is \
e
* M12 \e dt where \e M12 is called the "geodesic scale". \e M21 is
* defined similarly (with the geodesics being parallel at point 2). O
n a
* flat surface, we have \e M12 = \e M21 = 1. The quantity 1/\e M12 gi
ves
* the scale of the Cassini-Soldner projection.
* - <i>area</i>. Consider the quadrilateral bounded by the following li
nes:
* the geodesic from point 1 to point 2, the meridian from point 2 to t
he
* equator, the equator from \e lon2 to \e lon1, the meridian from the
* equator to point 1. The area of this quadrilateral is represented b
y \e
* S12 with a clockwise traversal of the perimeter counting as a positi
ve
* area and it can be used to compute the area of any simple geodesic
* polygon.
*
* Overloaded versions of Geodesic::Direct, Geodesic::ArcDirect, and
* Geodesic::Inverse allow these quantities to be returned. In addition
* there are general functions Geodesic::GenDirect, and Geodesic::GenInve
rse
* which allow an arbitrary set of results to be computed.
* *
* As an alternative to using distance to measure \e s12, the class can a * Additional functionality if provided by the GeodesicLine class, which
lso * allows a sequence of points along a geodesic to be computed.
* use the arc length \e a12 (in degrees) on the auxiliary sphere. This
is a
* mathematical construct used in solving the geodesic problems. However
, an
* arc length in excess of 180<sup>o</sup> indicates that the geodesic is
not
* a shortest path. In addition, the arc length between an equatorial
* crossing and the next extremum of latitude for a geodesic is
* 90<sup>o</sup>.
* *
* The calculations are accurate to better than 15 nm. (See \ref geoderr ors * The calculations are accurate to better than 15 nm. (See \ref geoderr ors
* for details.) * for details.)
* *
* For more information on geodesics see \ref geodesic. * For more information on geodesics see \ref geodesic.
**********************************************************************/ **********************************************************************/
class Geodesic { class 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, nC1 = GEOD_ORD, nC1p = GEOD_ORD,
nA2 = GEOD_ORD, nC2 = GEOD_ORD, nA3 = GEOD_ORD, nC3 = GEOD_ORD; nA2 = GEOD_ORD, nC2 = GEOD_ORD,
nA3 = GEOD_ORD, nA3x = nA3,
nC3 = GEOD_ORD, nC3x = (nC3 * (nC3 - 1)) / 2,
nC4 = GEOD_ORD, nC4x = (nC4 * (nC4 + 1)) / 2;
static const unsigned maxit = 50; static const unsigned maxit = 50;
static inline real sq(real x) throw() { return x * x; } 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,
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, bool diffp, real& dlam12, real& eps, 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, _etol2; const real _a, _r, _f, _f1, _e2, _ep2, _n, _b, _c2, _etol2;
static real SinSeries(real sinx, real cosx, const real c[], int n) real _A3x[nA3x], _C3x[nC3x], _C4x[nC4x];
static real SinCosSeries(bool sinp,
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
// 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 131 skipping to change at line 164
sinx /= r; sinx /= r;
cosx /= r; cosx /= r;
} }
// These are Maxima generated functions to provide series approximation s to // These are Maxima generated functions to provide series approximation s to
// the integrals for the ellipsoidal geodesic. // the integrals for the ellipsoidal geodesic.
static real A1m1f(real eps) throw(); static real A1m1f(real eps) throw();
static void C1f(real eps, real c[]) throw(); static void C1f(real eps, real c[]) throw();
static void C1pf(real eps, real c[]) throw(); static void C1pf(real eps, real c[]) throw();
static real A2m1f(real eps) throw(); static real A2m1f(real eps) throw();
static void C2f(real eps, real c[]) throw(); static void C2f(real eps, real c[]) throw();
static real A3f(real f, real eps) throw(); void A3coeff() throw();
static void C3f(real f, real eps, real c[]) throw(); real A3f(real eps) const throw();
void C3coeff() throw();
void C3f(real eps, real c[]) const throw();
void C4coeff() throw();
void C4f(real k2, real c[]) const throw();
enum captype {
CAP_NONE = 0U,
CAP_C1 = 1U<<0,
CAP_C1p = 1U<<1,
CAP_C2 = 1U<<2,
CAP_C3 = 1U<<3,
CAP_C4 = 1U<<4,
CAP_ALL = 0x1FU,
OUT_ALL = 0x7F80U,
};
public: public:
/** /**
* Constructor for a ellipsoid with equatorial radius \e a (meters) and * Bit masks for what calculations to do. These masks do double duty.
* reciprocal flattening \e r. Setting \e r = 0 implies \e r = inf or * They signify to the GeodesicLine::GeodesicLine constructor and to
* flattening = 0 (i.e., a sphere). Negative \e r indicates a prolate * Geodesic::Line what capabilities should be included in the GeodesicL
* ellipsoid. An exception is thrown if \e a is not positive. ine
* object. They also specify which results to return in the general
* routines Geodesic::GenDirect and Geodesic::GenInverse routines.
* GeodesicLine::mask is a duplication of this enum.
**********************************************************************
/
enum mask {
/**
* No capabilities, no output.
* @hideinitializer
********************************************************************
**/
NONE = 0U,
/**
* Calculate latitude \e lat2. (It's not necessary to include this a
s a
* capability to GeodesicLine because this is included by default.)
* @hideinitializer
********************************************************************
**/
LATITUDE = 1U<<7 | CAP_NONE,
/**
* Calculate longitude \e lon2.
* @hideinitializer
********************************************************************
**/
LONGITUDE = 1U<<8 | CAP_C3,
/**
* Calculate azimuths \e azi1 and \e azi2. (It's not necessary to
* include this as a capability to GeodesicLine because this is inclu
ded
* by default.)
* @hideinitializer
********************************************************************
**/
AZIMUTH = 1U<<9 | CAP_NONE,
/**
* Calculate distance \e s12.
* @hideinitializer
********************************************************************
**/
DISTANCE = 1U<<10 | CAP_C1,
/**
* Allow distance \e s12 to be used as input in the direct geodesic
* problem.
* @hideinitializer
********************************************************************
**/
DISTANCE_IN = 1U<<11 | CAP_C1 | CAP_C1p,
/**
* Calculate reduced length \e m12.
* @hideinitializer
********************************************************************
**/
REDUCEDLENGTH = 1U<<12 | CAP_C1 | CAP_C2,
/**
* Calculate geodesic scales \e M12 and \e M21.
* @hideinitializer
********************************************************************
**/
GEODESICSCALE = 1U<<13 | CAP_C1 | CAP_C2,
/**
* Calculate area \e S12.
* @hideinitializer
********************************************************************
**/
AREA = 1U<<14 | CAP_C4,
/**
* All capabilities. Calculate everything.
* @hideinitializer
********************************************************************
**/
ALL = OUT_ALL| CAP_ALL,
};
/** \name Constructor
**********************************************************************
/
///@{
/**
* Constructor for a ellipsoid with
*
* @param[in] a equatorial radius (meters)
* @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
* ellipsoid.
*
* An exception is thrown if either of the axes of the ellipsoid is
* non-positive.
********************************************************************** / ********************************************************************** /
Geodesic(real a, real r); Geodesic(real a, real r);
///@}
/** \name Direct geodesic problem specified in terms of distance.
**********************************************************************
/
///@{
/** /**
* Perform the direct geodesic calculation. Given a latitude, \e lat1, * Perform the direct geodesic calculation where the length of the geod
* longitude, \e lon1, and azimuth \e azi1 (degrees) for point 1 and a esic
* range, \e s12 (meters) from point 1 to point 2, return the latitude, * is specify in terms of distance.
\e *
* lat2, longitude, \e lon2, and forward azimuth, \e azi2 (degrees) for * @param[in] lat1 latitude of point 1 (degrees).
* point 2 and the reduced length \e m12 (meters). If \e arcmode (defa * @param[in] lon1 longitude of point 1 (degrees).
ult * @param[in] azi1 azimuth at point 1 (degrees).
* false) is set to true, \e s12 is interpreted as the arc length \e a1 * @param[in] s12 distance between point 1 and point 2 (meters); it can
2 be
* (degrees) on the auxiliary sphere. An arc length greater that 180 * signed.
* degrees results in a geodesic which is not a shortest path. For a * @param[out] lat2 latitude of point 2 (degrees).
* @param[out] lon2 longitude of point 2 (degrees).
* @param[out] azi2 (forward) azimuth at point 2 (degrees).
* @param[out] m12 reduced length of geodesic (meters).
* @param[out] M12 geodesic scale of point 2 relative to point 1
* (dimensionless).
* @param[out] M21 geodesic scale of point 1 relative to point 2
* (dimensionless).
* @param[out] S12 area under the geodesic (meters<sup>2</sup>).
* @return \e a12 arc length of between point 1 and point 2 (degrees).
*
* If either point is at a pole, the azimuth is defined by keeping the
* longitude fixed and writing \e lat = 90 - \e eps or -90 + \e eps and
* taking the limit \e eps -> 0 from above. An arc length greater that
180
* degrees signifies a geodesic which is not a shortest path. (For a
* prolate ellipsoid, an additional condition is necessary for a shorte st * prolate ellipsoid, an additional condition is necessary for a shorte st
* path: the longitudinal extent must not exceed of 180 degrees. Retur * path: the longitudinal extent must not exceed of 180 degrees.)
ned *
* value is the arc length \e a12 (degrees) if \e arcmode is false, * The following functions are overloaded versions of Geodesic::Direct
* otherwise it is the distance \e s12 (meters) * which omit some of the output parameters. Note, however, that the a
rc
* length is always computed and returned as the function value.
********************************************************************** / ********************************************************************** /
Math::real Direct(real lat1, real lon1, real azi1, real s12, Math::real Direct(real lat1, real lon1, real azi1, real s12,
real& lat2, real& lon2, real& azi2, real& m12, real& lat2, real& lon2, real& azi2,
bool arcmode = false) const throw(); real& m12, real& M12, real& M21, real& S12)
const throw() {
real t;
return GenDirect(lat1, lon1, azi1, false, s12,
LATITUDE | LONGITUDE | AZIMUTH |
REDUCEDLENGTH | GEODESICSCALE | AREA,
lat2, lon2, azi2, t, m12, M12, M21, S12);
}
/** /**
* Perform the inverse geodesic calculation. Given a latitude, \e lat1 * See the documentation for Geodesic::Direct.
,
* longitude, \e lon1, for point 1 and a latitude, \e lat2, longitude,
\e
* lon2, for point 2 (all in degrees), return the geodesic distance, \e
s12
* (meters), the forward azimuths, \e azi1 and \e azi2 (degrees) at poi
nts
* 1 and 2, and the reduced length \e m12 (meters). Returned value is
the
* arc length \e a12 (degrees) on the auxiliary sphere. The routine us
es
* an iterative method. If the method fails to converge, the negative
of
* the distances (\e s12, \e m12, and \e a12) and reverse of the azimut
hs
* are returned. This is not expected to happen with ellipsoidal model
s of
* the earth. Please report all cases where this occurs.
********************************************************************** / ********************************************************************** /
Math::real Inverse(real lat1, real lon1, real lat2, real lon2, Math::real Direct(real lat1, real lon1, real azi1, real s12,
real& s12, real& azi1, real& azi2, real& m12) real& lat2, real& lon2)
const throw(); const throw() {
real t;
return GenDirect(lat1, lon1, azi1, false, s12,
LATITUDE | LONGITUDE,
lat2, lon2, t, t, t, t, t, t);
}
/** /**
* Set up to do a series of ranges. This returns a GeodesicLine object * See the documentation for Geodesic::Direct.
* with point 1 given by latitude, \e lat1, longitude, \e lon1, and azi
muth
* \e azi1 (degrees). Calls to GeodesicLine::Position return the
* position and azimuth for point 2 a specified distance away. Using
* GeodesicLine::Position is approximately 2.1 faster than calling
* Geodesic::Direct.
********************************************************************** / ********************************************************************** /
GeodesicLine Line(real lat1, real lon1, real azi1) Math::real Direct(real lat1, real lon1, real azi1, real s12,
const throw(); real& lat2, real& lon2, real& azi2)
const throw() {
real t;
return GenDirect(lat1, lon1, azi1, false, s12,
LATITUDE | LONGITUDE | AZIMUTH,
lat2, lon2, azi2, t, t, t, t, t);
}
/** /**
* The major radius of the ellipsoid (meters). This is that value of \ * See the documentation for Geodesic::Direct.
e a
* used in the constructor.
********************************************************************** / ********************************************************************** /
Math::real MajorRadius() const throw() { return _a; } Math::real Direct(real lat1, real lon1, real azi1, real s12,
real& lat2, real& lon2, real& azi2, real& m12)
const throw() {
real t;
return GenDirect(lat1, lon1, azi1, false, s12,
LATITUDE | LONGITUDE | AZIMUTH | REDUCEDLENGTH,
lat2, lon2, azi2, t, m12, t, t, t);
}
/** /**
* The inverse flattening of the ellipsoid. This is that value of \e r * See the documentation for Geodesic::Direct.
* used in the constructor. A value of 0 is returned for a sphere
* (infinite inverse flattening).
********************************************************************** / ********************************************************************** /
Math::real InverseFlattening() const throw() { return _r; } Math::real Direct(real lat1, real lon1, real azi1, real s12,
real& lat2, real& lon2, real& azi2,
real& M12, real& M21)
const throw() {
real t;
return GenDirect(lat1, lon1, azi1, false, s12,
LATITUDE | LONGITUDE | AZIMUTH | GEODESICSCALE,
lat2, lon2, azi2, t, t, M12, M21, t);
}
/** /**
* A global instantiation of Geodesic with the parameters for the WGS84 * See the documentation for Geodesic::Direct.
* ellipsoid.
********************************************************************** / ********************************************************************** /
const static Geodesic WGS84; Math::real Direct(real lat1, real lon1, real azi1, real s12,
real& lat2, real& lon2, real& azi2,
real& m12, real& M12, real& M21)
const throw() {
real t;
return GenDirect(lat1, lon1, azi1, false, s12,
LATITUDE | LONGITUDE | AZIMUTH |
REDUCEDLENGTH | GEODESICSCALE,
lat2, lon2, azi2, t, m12, M12, M21, t);
}
///@}
}; /** \name Direct geodesic problem specified in terms of arc length.
**********************************************************************
/
///@{
/**
* Perform the direct geodesic calculation where the length of the geod
esic
* is specify in terms of arc length.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] azi1 azimuth at point 1 (degrees).
* @param[in] a12 arc length between point 1 and point 2 (degrees); it
can
* be signed.
* @param[out] lat2 latitude of point 2 (degrees).
* @param[out] lon2 longitude of point 2 (degrees).
* @param[out] azi2 (forward) azimuth at point 2 (degrees).
* @param[out] s12 distance between point 1 and point 2 (meters).
* @param[out] m12 reduced length of geodesic (meters).
* @param[out] M12 geodesic scale of point 2 relative to point 1
* (dimensionless).
* @param[out] M21 geodesic scale of point 1 relative to point 2
* (dimensionless).
* @param[out] S12 area under the geodesic (meters<sup>2</sup>).
*
* If either point is at a pole, the azimuth is defined by keeping the
* longitude fixed and writing \e lat = 90 - \e eps or -90 + \e eps and
* taking the limit \e eps -> 0 from above. An arc length greater that
180
* degrees signifies a geodesic which is not a shortest path. (For a
* prolate ellipsoid, an additional condition is necessary for a shorte
st
* path: the longitudinal extent must not exceed of 180 degrees.)
*
* The following functions are overloaded versions of Geodesic::Direct
* which omit some of the output parameters.
**********************************************************************
/
void ArcDirect(real lat1, real lon1, real azi1, real a12,
real& lat2, real& lon2, real& azi2, real& s12,
real& m12, real& M12, real& M21, real& S12)
const throw() {
GenDirect(lat1, lon1, azi1, true, a12,
LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
REDUCEDLENGTH | GEODESICSCALE | AREA,
lat2, lon2, azi2, s12, m12, M12, M21, S12);
}
/** /**
* \brief A geodesic line. * See the documentation for Geodesic::ArcDirect.
* **********************************************************************
* GeodesicLine facilitates the determination of a series of points on a /
* single geodesic. Geodesic.Line returns a GeodesicLine object with the void ArcDirect(real lat1, real lon1, real azi1, real a12,
* geodesic defined by by \e lat1, \e lon1, and \e azi1. real& lat2, real& lon2) const throw() {
* GeodesicLine.Position returns the \e lat2, \e lon2, \e azi2, and \e m1 real t;
2 GenDirect(lat1, lon1, azi1, true, a12,
* given \e s12. An example of use of this class is: LATITUDE | LONGITUDE,
\verbatim lat2, lon2, t, t, t, t, t, t);
// Print positions on a geodesic going through latitude 30, }
// longitude 10 at azimuth 80. Points at intervals of 10km
// in the range [-1000km, 1000km] are given.
GeographicLib::GeodesicLine
line(GeographicLib::Geodesic::WGS84.Line(30.0, 10.0, 80.0));
double step = 10e3;
for (int s = -100; s <= 100; ++s) {
double lat2, lon2, azi2, m12;
double s12 = s * step;
line.Position(s12, lat2, lon2, azi2, m12);
cout << s12 << " " << lat2 << " " << lon2 << " "
<< azi2 << " " << m12 << "\n";
}
\endverbatim
* The default copy constructor and assignment operators work with this
* class, so that, for example, the previous example could start
\verbatim
GeographicLib::GeodesicLine line;
line = GeographicLib::Geodesic::WGS84.Line(30.0, 10.0, 80.0);
...
\endverbatim
* Similarly, a vector can be used to hold GeodesicLine objects.
*
* The calculations are accurate to better than 12 nm. (See \ref geoderr
ors
* for details.)
**********************************************************************/
class GeodesicLine { /**
private: * See the documentation for Geodesic::ArcDirect.
typedef Math::real real; **********************************************************************
friend class Geodesic; /
static const int nC1 = Geodesic::nC1, nC1p = Geodesic::nC1p, void ArcDirect(real lat1, real lon1, real azi1, real a12,
nC2 = Geodesic::nC2, nC3 = Geodesic::nC3; real& lat2, real& lon2, real& azi2) const throw() {
real t;
real _lat1, _lon1, _azi1; GenDirect(lat1, lon1, azi1, true, a12,
real _a, _r, _b, _f1, _salp0, _calp0, _k2, LATITUDE | LONGITUDE | AZIMUTH,
_ssig1, _csig1, _stau1, _ctau1, _somg1, _comg1, lat2, lon2, azi2, t, t, t, t, t);
_A1m1, _A2m1, _A3c, _B11, _B21, _B31; }
// index zero elements of these arrays are unused
real _C1a[nC1 + 1], _C1pa[nC1p + 1], _C2a[nC2 + 1], _C3a[nC3];
static inline real sq(real x) throw() { return x * x; } /**
GeodesicLine(const Geodesic& g, real lat1, real lon1, real azi1) * See the documentation for Geodesic::ArcDirect.
throw(); **********************************************************************
public: /
void ArcDirect(real lat1, real lon1, real azi1, real a12,
real& lat2, real& lon2, real& azi2, real& s12)
const throw() {
real t;
GenDirect(lat1, lon1, azi1, true, a12,
LATITUDE | LONGITUDE | AZIMUTH | DISTANCE,
lat2, lon2, azi2, s12, t, t, t, t);
}
/** /**
* A default constructor. If GeodesicLine::Position is called on the * See the documentation for Geodesic::ArcDirect.
* resulting object, it returns immediately (without doing any
* calculations). The object should be set with a call to Geodesic::Li
ne.
* Use Init() to test whether object is still in this uninitialized sta
te.
**********************************************************************
/
GeodesicLine() throw() : _b(0) {};
/**
* Return the latitude, \e lat2, longitude, \e lon2, and forward azimut
h,
* \e azi2 (degrees) of the point 2 which is a distance, \e s12 (in
* meters), from point 1. Also return the reduced length \e m12 (meter
s).
* \e s12 can be signed. If \e arcmode (default false) is set to true,
\e
* s12 is interpreted as the arc length \e a12 (in degrees) on the
* auxiliary sphere. Returned value is the arc length \e a12 (degrees)
if
* \e arcmode is false, otherwise it is the distance \e s12 (meters).
********************************************************************** / ********************************************************************** /
Math::real Position(real s12, real& lat2, real& lon2, void ArcDirect(real lat1, real lon1, real azi1, real a12,
real& azi2, real &m12, bool arcmode = false) real& lat2, real& lon2, real& azi2,
const throw(); real& s12, real& m12) const throw() {
real t;
GenDirect(lat1, lon1, azi1, true, a12,
LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
REDUCEDLENGTH,
lat2, lon2, azi2, s12, m12, t, t, t);
}
/** /**
* Return the scale of the geodesic line extending an arc length \e a12 * See the documentation for Geodesic::ArcDirect.
* (degrees) from point 1 to point 2. \e M12 (a number) measures the
* convergence of initially parallel geodesics. It is defined by the
* following construction: starting at point 1 proceed at azimuth \e az
i1 +
* 90<sup>o</sup> a small distance \e dt; turn -90<sup>o</sup> and proc
eed
* a distance \e s12 (\e not the arc length \e a12); the distance to po
int
* 2 is given by \e M12 \e dt. \e M21 is defined analogously.
********************************************************************** / ********************************************************************** /
void Scale(real a12, real& M12, real& M21) const throw(); void ArcDirect(real lat1, real lon1, real azi1, real a12,
real& lat2, real& lon2, real& azi2, real& s12,
real& M12, real& M21) const throw() {
real t;
GenDirect(lat1, lon1, azi1, true, a12,
LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
GEODESICSCALE,
lat2, lon2, azi2, s12, t, M12, M21, t);
}
/** /**
* Has this object been initialized so that Position can be called? * See the documentation for Geodesic::ArcDirect.
********************************************************************** / ********************************************************************** /
bool Init() const throw() { return _b > 0; } void ArcDirect(real lat1, real lon1, real azi1, real a12,
real& lat2, real& lon2, real& azi2, real& s12,
real& m12, real& M12, real& M21) const throw() {
real t;
GenDirect(lat1, lon1, azi1, true, a12,
LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
REDUCEDLENGTH | GEODESICSCALE,
lat2, lon2, azi2, s12, m12, M12, M21, t);
}
///@}
/** \name General version of the direct geodesic solution.
**********************************************************************
/
///@{
/** /**
* Return the latitude of point 1 (degrees). * The general direct geodesic calculation. Geodesic::Direct and
* Geodesic::ArcDirect are defined in terms of this function.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] azi1 azimuth at point 1 (degrees).
* @param[in] arcmode boolean flag determining the meaning of the secon
d
* parameter.
* @param[in] s12_a12 if \e arcmode is false, this is the distance betw
een
* point 1 and point 2 (meters); otherwise it is the arc length betwe
en
* point 1 and point 2 (degrees); it can be signed.
* @param[in] outmask a bitor'ed combination of Geodesic::mask values
* specifying which of the following parameters should be set.
* @param[out] lat2 latitude of point 2 (degrees).
* @param[out] lon2 longitude of point 2 (degrees).
* @param[out] azi2 (forward) azimuth at point 2 (degrees).
* @param[out] s12 distance between point 1 and point 2 (meters).
* @param[out] m12 reduced length of geodesic (meters).
* @param[out] M12 geodesic scale of point 2 relative to point 1
* (dimensionless).
* @param[out] M21 geodesic scale of point 1 relative to point 2
* (dimensionless).
* @param[out] S12 area under the geodesic (meters<sup>2</sup>).
* @return \e a12 arc length of between point 1 and point 2 (degrees).
*
* The Geodesic::mask values possible for \e outmask are
* - \e outmask |= Geodesic::LATITUDE for the latitude \e lat2.
* - \e outmask |= Geodesic::LONGITUDE for the latitude \e lon2.
* - \e outmask |= Geodesic::AZIMUTH for the latitude \e azi2.
* - \e outmask |= Geodesic::DISTANCE for the distance \e s12.
* - \e outmask |= Geodesic::REDUCEDLENGTH for the reduced length \e
* m12.
* - \e outmask |= Geodesic::GEODESICSCALE for the geodesic scales \e
* M12 and \e M21.
* - \e outmask |= Geodesic::AREA for the area \e S12.
* .
* The function value \e a12 is always computed and returned and this
* equals \e s12_a12 is \e arcmode is true. If \e outmask includes
* Geodesic::DISTANCE and \e arcmode is false, then \e s12 = \e s12_a12
.
* It is not necessary to include Geodesic::DISTANCE_IN in \e outmask;
this
* is automatically included is \e arcmode is false.
**********************************************************************
/
Math::real GenDirect(real lat1, real lon1, real azi1,
bool arcmode, real s12_a12, unsigned outmask,
real& lat2, real& lon2, real& azi2,
real& s12, real& m12, real& M12, real& M21,
real& S12) const throw();
///@}
/** \name Inverse geodesic problem.
**********************************************************************
/
///@{
/**
* Perform the inverse geodesic calculation.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] lat2 latitude of point 2 (degrees).
* @param[in] lon2 longitude of point 2 (degrees).
* @param[out] s12 distance between point 1 and point 2 (meters).
* @param[out] azi1 azimuth at point 1 (degrees).
* @param[out] azi2 (forward) azimuth at point 1 (degrees).
* @param[out] m12 reduced length of geodesic (meters).
* @param[out] M12 geodesic scale of point 2 relative to point 1
* (dimensionless).
* @param[out] M21 geodesic scale of point 1 relative to point 2
* (dimensionless).
* @param[out] S12 area under the geodesic (meters<sup>2</sup>).
* @return \e a12 arc length of between point 1 and point 2 (degrees).
*
* If either point is at a pole, the azimuth is defined by keeping the
* longitude fixed and writing \e lat = 90 - \e eps or -90 + \e eps and
* taking the limit \e eps -> 0 from above. If the routine fails to
* converge, then all the requested outputs are set to Math::NaN(). Th
is
* is not expected to happen with ellipsoidal models of the earth; plea
se
* report all cases where this occurs.
*
* The following functions are overloaded versions of Geodesic::Inverse
* which omit some of the output parameters. Note, however, that the a
rc
* length is always computed and returned as the function value.
********************************************************************** / ********************************************************************** /
Math::real Latitude() const throw() { return Init() ? _lat1 : 0; } Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
real& s12, real& azi1, real& azi2, real& m12,
real& M12, real& M21, real& S12) const throw() {
return GenInverse(lat1, lon1, lat2, lon2,
DISTANCE | AZIMUTH |
REDUCEDLENGTH | GEODESICSCALE | AREA,
s12, azi1, azi2, m12, M12, M21, S12);
}
/** /**
* Return the longitude of point 1 (degrees). * See the documentation for Geodesic::Inverse.
********************************************************************** / ********************************************************************** /
Math::real Longitude() const throw() { return Init() ? _lon1 : 0; } Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
real& s12) const throw() {
real t;
return GenInverse(lat1, lon1, lat2, lon2,
DISTANCE,
s12, t, t, t, t, t, t);
}
/** /**
* Return the azimuth (degrees) of the geodesic line as it passes throu * See the documentation for Geodesic::Inverse.
gh
* point 1.
********************************************************************** / ********************************************************************** /
Math::real Azimuth() const throw() { return Init() ? _azi1 : 0; } Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
real& azi1, real& azi2) const throw() {
real t;
return GenInverse(lat1, lon1, lat2, lon2,
AZIMUTH,
t, azi1, azi2, t, t, t, t);
}
/** /**
* Return the azimuth (degrees) of the geodesic line as it crosses the * See the documentation for Geodesic::Inverse.
* equator in a northward direction.
********************************************************************** / ********************************************************************** /
Math::real EquatorialAzimuth() const throw() { Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
return Init() ? atan2(_salp0, _calp0) / Constants::degree() : 0; real& s12, real& azi1, real& azi2)
const throw() {
real t;
return GenInverse(lat1, lon1, lat2, lon2,
DISTANCE | AZIMUTH,
s12, azi1, azi2, t, t, t, t);
} }
/** /**
* Return the arc length (degrees) between the northward equatorial * See the documentation for Geodesic::Inverse.
* crossing and point 1.
********************************************************************** / ********************************************************************** /
Math::real EquatorialArc() const throw() { Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
return Init() ? atan2(_ssig1, _csig1) / Constants::degree() : 0; real& s12, real& azi1, real& azi2, real& m12)
const throw() {
real t;
return GenInverse(lat1, lon1, lat2, lon2,
DISTANCE | AZIMUTH | REDUCEDLENGTH,
s12, azi1, azi2, m12, t, t, t);
} }
/** /**
* The major radius of the ellipsoid (meters). This is that value of \ * See the documentation for Geodesic::Inverse.
e a
* inherited from the Geodesic object used in the constructor.
********************************************************************** / ********************************************************************** /
Math::real MajorRadius() const throw() { return Init() ? _a : 0; } Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
real& s12, real& azi1, real& azi2,
real& M12, real& M21) const throw() {
real t;
return GenInverse(lat1, lon1, lat2, lon2,
DISTANCE | AZIMUTH | GEODESICSCALE,
s12, azi1, azi2, t, M12, M21, t);
}
/** /**
* The inverse flattening of the ellipsoid. This is that value of \e r * See the documentation for Geodesic::Inverse.
* inherited from the Geodesic object used in the constructor. A value
of
* 0 is returned for a sphere (infinite inverse flattening).
********************************************************************** / ********************************************************************** /
Math::real InverseFlattening() const throw() { return Init() ? _r : 0; Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
} real& s12, real& azi1, real& azi2, real& m12,
real& M12, real& M21) const throw() {
real t;
return GenInverse(lat1, lon1, lat2, lon2,
DISTANCE | AZIMUTH |
REDUCEDLENGTH | GEODESICSCALE,
s12, azi1, azi2, m12, M12, M21, t);
}
///@}
/** \name General version of inverse geodesic solution.
**********************************************************************
/
///@{
/**
* The general inverse geodesic calculation. Geodesic::Inverse is defi
ned
* in terms of this function.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] lat2 latitude of point 2 (degrees).
* @param[in] lon2 longitude of point 2 (degrees).
* @param[in] outmask a bitor'ed combination of Geodesic::mask values
* specifying which of the following parameters should be set.
* @param[out] s12 distance between point 1 and point 2 (meters).
* @param[out] azi1 azimuth at point 1 (degrees).
* @param[out] azi2 (forward) azimuth at point 1 (degrees).
* @param[out] m12 reduced length of geodesic (meters).
* @param[out] M12 geodesic scale of point 2 relative to point 1
* (dimensionless).
* @param[out] M21 geodesic scale of point 1 relative to point 2
* (dimensionless).
* @param[out] S12 area under the geodesic (meters<sup>2</sup>).
* @return \e a12 arc length of between point 1 and point 2 (degrees).
*
* The Geodesic::mask values possible for \e outmask are
* - \e outmask |= Geodesic::DISTANCE for the distance \e s12.
* - \e outmask |= Geodesic::AZIMUTH for the latitude \e azi2.
* - \e outmask |= Geodesic::REDUCEDLENGTH for the reduced length \e
* m12.
* - \e outmask |= Geodesic::GEODESICSCALE for the geodesic scales \e
* M12 and \e M21.
* - \e outmask |= Geodesic::AREA for the area \e S12.
* .
* The arc length is always computed and returned as the function value
.
**********************************************************************
/
Math::real GenInverse(real lat1, real lon1, real lat2, real lon2,
unsigned outmask,
real& s12, real& azi1, real& azi2,
real& m12, real& M12, real& M21, real& S12)
const throw();
///@}
/** \name Interface to GeodesicLine.
**********************************************************************
/
///@{
/**
* Set up to do a series of ranges.
*
* @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees).
* @param[in] azi1 azimuth at point 1 (degrees).
* @param[in] caps bitor'ed combination of Geodesic::mask values
* specifying the capabilities the GeodesicLine object should possess
,
* i.e., which quantities can be returned in calls to
* GeodesicLib::Position.
*
* The Geodesic::mask values are
* - \e caps |= Geodesic::LATITUDE for the latitude \e lat2; this is
* added automatically
* - \e caps |= Geodesic::LONGITUDE for the latitude \e lon2
* - \e caps |= Geodesic::AZIMUTH for the latitude \e azi2; this is
* added automatically
* - \e caps |= Geodesic::DISTANCE for the distance \e s12
* - \e caps |= Geodesic::REDUCEDLENGTH for the reduced length \e m12
* - \e caps |= Geodesic::GEODESICSCALE for the geodesic scales \e M12
* and \e M21
* - \e caps |= Geodesic::AREA for the area \e S12
* - \e caps |= Geodesic::DISTANCE_IN permits the length of the
* geodesic to be given in terms of \e s12; without this capability t
he
* length can only be specified in terms of arc length.
* .
* The default value of \e caps is Geodesic::ALL which turns on all the
* capabilities.
*
* If the point is at a pole, the azimuth is defined by keeping the \e
lon1
* fixed and writing \e lat1 = 90 - \e eps or -90 + \e eps and taking t
he
* limit \e eps -> 0 from above.
**********************************************************************
/
GeodesicLine Line(real lat1, real lon1, real azi1, unsigned caps = ALL)
const throw();
///@}
/** \name Inspector functions.
**********************************************************************
/
///@{
/**
* @return \e a the equatorial radius of the ellipsoid (meters). This
is
* the value used in the constructor.
**********************************************************************
/
Math::real MajorRadius() const throw() { return _a; }
/**
* @return \e r the inverse flattening of the ellipsoid. This is the
* value used in the constructor. A value of 0 is returned for a sph
ere
* (infinite inverse flattening).
**********************************************************************
/
Math::real InverseFlattening() const throw() { return _r; }
/**
* @return total area of ellipsoid in meters<sup>2</sup>. The area of
a
* polygon encircling a pole can be found by adding
* Geodesic::EllipsoidArea()/2 to the sum of \e S12 for each side of
the
* polygon.
**********************************************************************
/
Math::real EllipsoidArea() const throw() {
return 4 * Constants::pi() * _c2;
}
///@}
/**
* A global instantiation of Geodesic with the parameters for the WGS84
* ellipsoid.
**********************************************************************
/
const static Geodesic WGS84;
/** \name Deprecated function.
**********************************************************************
/
///@{
/**
* <b>DEPRECATED</b> Perform the direct geodesic calculation. Given a
* 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
* 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
* point is at a pole, the azimuth is defined by keeping the longitude
* fixed and writing \e lat = 90 - \e eps or -90 + \e eps and taking th
e
* limit \e eps -> 0 from above. If \e arcmode (default false) is set
to
* true, \e s12 is interpreted as the arc length \e a12 (degrees) on th
e
* auxiliary sphere. An arc length greater that 180 degrees results in
a
* geodesic which is not a shortest path. For a prolate ellipsoid, an
* additional condition is necessary for a shortest path: the longitudi
nal
* extent must not exceed of 180 degrees. Returned value is the arc le
ngth
* \e a12 (degrees) if \e arcmode is false, otherwise it is the distanc
e \e
* s12 (meters).
**********************************************************************
/
Math::real Direct(real lat1, real lon1, real azi1, real s12_a12,
real& lat2, real& lon2, real& azi2, real& m12,
bool arcmode) const throw() {
if (arcmode) {
real a12 = s12_a12, s12;
ArcDirect(lat1, lon1, azi1, a12, lat2, lon2, azi2, s12, m12);
return s12;
} else {
real s12 = s12_a12;
return Direct(lat1, lon1, azi1, s12, lat2, lon2, azi2, m12);
}
}
///@}
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif #endif
 End of changes. 49 change blocks. 
227 lines changed or deleted 737 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 6785 2010-01-05 22:15:42Z k arney $" #define GEOGRAPHICLIB_GEOID_HPP "$Id: Geoid.hpp 6868 2010-09-12 10:27:11Z k arney $"
#include "GeographicLib/Constants.hpp" #include "GeographicLib/Constants.hpp"
#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
* *
skipping to change at line 80 skipping to change at line 80
// 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
// remove the cast in this case.
_file.seekg(std::ios::streamoff(_datastart + _file.seekg(std::ios::streamoff(_datastart +
2ULL * (unsigned(iy) * _swidth + 2ULL * (unsigned(iy) * _swidth +
unsigned(ix)))); 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 &&
skipping to change at line 127 skipping to change at line 129
} }
} }
} }
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:
/** /**
* Create a Geoid loading the data for geoid \e name. The data file is * Flags indicating conversions between heights above the geoid and hei
* formed by appending ".pgm" to the name. If \e path is specified (an ghts
d is * above the ellipsoid.
* non-empty), then the file is loaded from directory, \e path. Otherw **********************************************************************
ise /
* the path is given by the GEOID_PATH environment variable. If that i enum convertflag {
s /**
* undefined, a compile-time default path is used * The multiplier for converting from heights above the geoid to heig
hts
* above the ellipsoid.
********************************************************************
**/
ELLIPSOIDTOGEOID = -1,
/**
* No conversion.
********************************************************************
**/
NONE = 0,
/**
* The multiplier for converting from heights above the ellipsoid to
* heights above the geoid.
********************************************************************
**/
GEOIDTOELLIPSOID = 1,
};
/** \name Setting up the geoid
**********************************************************************
/
///@{
/**
* Construct a Geoid.
*
* @param[in] name the name of the geoid.
* @param[in] path (optional) directory for data file.
* @param[in] cubic interpolation method; false means bilinear, true (t
he
* default) means cubic.
*
* 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
* path. Otherwise the path is given by the GEOID_PATH environment
* 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).
The * This may throw an error because the file does not exist, is unreadab
* final \e cubic argument specifies whether to use bilinear (\e cubic le,
= * or is corrupt.
* false) or cubic (\e cubic = true, the default) interpolation. This
may
* throw an error because the file does not exist, is unreadable, or is
* corrupt.
********************************************************************** / ********************************************************************** /
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);
/** /**
* Cache the data for the rectangular area defined by the four argument * Set up a cache.
s \e *
* south, \e west, \e north, \e east (all in degrees). \e east is alwa * @param[in] south latitude (degrees) of the south edge of the cached
ys area.
* @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] 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
* 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].
********************************************************************** / ********************************************************************** /
void CacheArea(real south, real west, real north, real east) const; void CacheArea(real south, real west, real north, real east) const;
skipping to change at line 172 skipping to change at line 208
* Geoid::CacheArea on a specific area. * Geoid::CacheArea on a specific area.
********************************************************************** / ********************************************************************** /
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.
********************************************************************** / ********************************************************************** /
void CacheClear() const throw(); void CacheClear() const throw();
///@}
/** \name Compute geoid heights
**********************************************************************
/
///@{
/** /**
* Return the geoid height in meters for latitude \e lat (in [-90, 90]) * Compute the geoid height at a point
and *
* longitude \e lon (in [-180,360]), both in degrees. This may throw a * @param[in] lat latitude of the point (degrees).
n * @param[in] lon longitude of the point (degrees).
* error because of an error reading data from disk. However, it will * @return geoid height (meters).
not *
* throw if (\e lat, \e lon) is within a successfully cached area. * The latitude should be in [-90, 90] and longitude should bein
* [-180,360]. This may throw an error because of an error reading dat
a
* from disk. However, it will not throw if (\e lat, \e lon) is within
a
* successfully cached area.
********************************************************************** / ********************************************************************** /
Math::real operator()(real lat, real lon) const { Math::real operator()(real lat, real lon) const {
real gradn, grade; real gradn, grade;
return height(lat, lon, false, gradn, grade); return height(lat, lon, false, gradn, grade);
} }
/** /**
* Return the geoid height in meters for latitude \e lat (in [-90, 90]) * Compute the geoid height and gradient at a point
and *
* longitude \e lon (in [-180,360]), both in degrees. In addition comp * @param[in] lat latitude of the point (degrees).
ute * @param[in] lon longitude of the point (degrees).
* the gradient of the geoid height in the northerly \e gradn and easte * @param[out] gradn northerly gradient (dimensionless).
rly * @param[out] grade easterly gradient (dimensionless).
* \e grade directions. This may throw an error because of an error * @return geoid height (meters).
* reading data from disk. However, it will not throw if (\e lat, \e l *
on) * The latitude should be in [-90, 90] and longitude should be in [-180
* is within a successfully cached area. ,
* 360]. This may throw an error because of an error reading data from
* disk. However, it will not throw if (\e lat, \e lon) is within a
* successfully cached area.
********************************************************************** / ********************************************************************** /
Math::real operator()(real lat, real lon, real& gradn, real& grade) con st { Math::real operator()(real lat, real lon, real& gradn, real& grade) con st {
return height(lat, lon, true, gradn, grade); return height(lat, lon, true, gradn, grade);
} }
/** /**
* Return the geoid description if available in the data file. If abse * Convert a height above the geoid to a height above the ellipsoid and
nt, * vice versa.
* return "NONE". *
* @param[in] lat latitude of the point (degrees).
* @param[in] lon longitude of the point (degrees).
* @param[in] h height of the point (degrees).
* @param[in] d a Geoid::convertflag specifying the direction of the
* conversion; Geoid::GEOIDTOELLIPSOID means convert a height above t
he
* geoid to a height above the ellipsoid; Geoid::ELLIPSOIDTOGEOID mea
ns
* convert a height above the ellipsoid to a height above the geoid.
* @return converted height (meters).
**********************************************************************
/
Math::real ConvertHeight(real lat, real lon, real h,
convertflag d) const {
real gradn, grade;
return h + real(d) * height(lat, lon, true, gradn, grade);
}
///@}
/** \name Inspector functions
**********************************************************************
/
///@{
/**
* @return geoid description, if available, in the data file; if
* absent, return "NONE".
********************************************************************** / ********************************************************************** /
const std::string& Description() const throw() { return _description; } const std::string& Description() const throw() { return _description; }
/** /**
* Return the date of the data file. If absent, return "UNKNOWN". * @return date of the data file; if absent, return "UNKNOWN".
********************************************************************** / ********************************************************************** /
const std::string& DateTime() const throw() { return _datetime; } const std::string& DateTime() const throw() { return _datetime; }
/** /**
* Return the full file name used to load the geoid data. * @return full file name used to load the geoid data.
********************************************************************** / ********************************************************************** /
const std::string& GeoidFile() const throw() { return _filename; } const std::string& GeoidFile() const throw() { return _filename; }
/** /**
* Return the "name" used to load the geoid data (from the first argume * @return "name" used to load the geoid data (from the first argument
nt of
* of the constructor). * the constructor).
********************************************************************** / ********************************************************************** /
const std::string& GeoidName() const throw() { return _name; } const std::string& GeoidName() const throw() { return _name; }
/** /**
* Return the directory used to load the geoid data. * @return directory used to load the geoid data.
********************************************************************** / ********************************************************************** /
const std::string& GeoidDirectory() const throw() { return _dir; } const std::string& GeoidDirectory() const throw() { return _dir; }
/** /**
* Return the interpolation method (cubic or bilinear). * @return interpolation method ("cubic" or "bilinear").
********************************************************************** / ********************************************************************** /
const std::string Interpolation() const const std::string Interpolation() const
{ return std::string(_cubic ? "cubic" : "bilinear"); } { return std::string(_cubic ? "cubic" : "bilinear"); }
/** /**
* Return a estimate of the maximum interpolation and quantization erro * @return estimate of the maximum interpolation and quantization error
r * (meters).
* (meters). This relies on the value being stored in the data file. *
If * This relies on the value being stored in the data file. If the valu
* the value is absent, return -1. e is
* absent, return -1.
********************************************************************** / ********************************************************************** /
Math::real MaxError() const throw() { return _maxerror; } Math::real MaxError() const throw() { return _maxerror; }
/** /**
* Return a estimate of the RMS interpolation and quantization error * @return estimate of the RMS interpolation and quantization error
* (meters). This relies on the value being stored in the data file. * (meters).
If *
* the value is absent, return -1. * This relies on the value being stored in the data file. If the valu
e is
* absent, return -1.
********************************************************************** / ********************************************************************** /
Math::real RMSError() const throw() { return _rmserror; } Math::real RMSError() const throw() { return _rmserror; }
/** /**
* Return offset (meters) for converting pixel values to geoid heights. * @return offset (meters).
*
* This in used in converting from the pixel values in the data file to
* geoid heights.
********************************************************************** / ********************************************************************** /
Math::real Offset() const throw() { return _offset; } Math::real Offset() const throw() { return _offset; }
/** /**
* Return scale (meters) for converting pixel values to geoid * @return scale (meters).
* heights. *
* This in used in converting from the pixel values in the data file to
* geoid heights.
********************************************************************** / ********************************************************************** /
Math::real Scale() const throw() { return _scale; } Math::real Scale() const throw() { return _scale; }
/** /**
* Is a data cache active? * @return true if a data cache is active.
********************************************************************** / ********************************************************************** /
bool Cache() const throw() { return _cache; } bool Cache() const throw() { return _cache; }
/** /**
* Return the west edge of the cached area. The cache includes this ed ge. * @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 :
0; 0;
} }
/** /**
* Return the east edge of the cached area. The cache excludes this ed ge. * @return east edge of the cached area; the cache excludes this edge.
********************************************************************** / ********************************************************************** /
Math::real CacheEast() const throw() { Math::real CacheEast() const throw() {
return _cache ? return _cache ?
CacheWest() + CacheWest() +
(_xsize - (_xsize == _width ? 0 : 1 + 2 * _cubic)) / _rlonres : (_xsize - (_xsize == _width ? 0 : 1 + 2 * _cubic)) / _rlonres :
0; 0;
} }
/** /**
* Return the north edge of the cached area. The cache includes this e dge. * @return north edge of the cached area; the cache includes this edge.
********************************************************************** / ********************************************************************** /
Math::real CacheNorth() const throw() { Math::real CacheNorth() const throw() {
return _cache ? 90 - (_yoffset + _cubic) / _rlatres : 0; return _cache ? 90 - (_yoffset + _cubic) / _rlatres : 0;
} }
/** /**
* Return the south edge of the cached area. The cache excludes this e * @return south edge of the cached area; the cache excludes this edge
dge * 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 ;
} }
/** /**
* The major radius of the ellipsoid (meters). This is the value for t * @return \e a the equatorial radius of the WGS84 ellipsoid (meters).
he *
* WGS84 ellipsoid because the supported geoid models are all based on * (The WGS84 values is returned because the supported geoid models are
this all
* ellipsoid. * based on this ellipsoid.)
********************************************************************** / ********************************************************************** /
Math::real MajorRadius() const throw() { return Constants::WGS84_a(); } Math::real MajorRadius() const throw() { return Constants::WGS84_a(); }
/** /**
* The inverse flattening of the ellipsoid. This is the value for the * @return \e r the inverse flattening of the WGS84 ellipsoid.
* WGS84 ellipsoid because the supported geoid models are all based on *
this * (The WGS84 values is returned because the supported geoid models are
* ellipsoid. all
* 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();
/** /**
* Return the value of the environment variable GEOID_PATH. * Return the value of the environment variable GEOID_PATH.
********************************************************************** / ********************************************************************** /
static std::string GeoidPath(); static std::string GeoidPath();
 End of changes. 27 change blocks. 
76 lines changed or deleted 171 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 (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_GNOMONIC_HPP) #if !defined(GEOGRAPHICLIB_GNOMONIC_HPP)
#define GEOGRAPHICLIB_GNOMONIC_HPP "$Id: Gnomonic.hpp 6845 2010-07-21 10:23 :46Z karney $" #define GEOGRAPHICLIB_GNOMONIC_HPP "$Id: Gnomonic.hpp 6867 2010-09-11 13:04 :26Z karney $"
#include "GeographicLib/Geodesic.hpp" #include "GeographicLib/Geodesic.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. The projection of \e P is defined as follows: compute the * ellipsoid. 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,
skipping to change at line 38 skipping to change at line 39
* azimuthal direction. The scale in the radial direction if 1/\e * azimuthal direction. The scale in the radial direction if 1/\e
* rk<sup>2</sup>. * rk<sup>2</sup>.
* *
* For a sphere, \e rho is reduces to \e a tan(\e s12/\e a), where \e s12 is * For a sphere, \e rho is reduces to \e a tan(\e s12/\e a), where \e s12 is
* the length of the geodesic from \e C to \e P, and the gnomonic project ion * the length of the geodesic from \e C to \e P, and the gnomonic project ion
* has the property that all geodesics appear as straight lines. For an * has the property that all geodesics appear as straight lines. For an
* ellipsoid, this property holds only for geodesics interesting the cent ers. * ellipsoid, this property holds only for geodesics interesting the cent ers.
* However geodesic segments close to the center are approximately straig ht. * However geodesic segments close to the center are approximately straig ht.
* *
* Consider a geodesic segment of length \e l. Let \e T be the point on the * Consider a geodesic segment of length \e l. Let \e T be the point on the
* geodesic closest to \e C the center of the projection and \e t be the * geodesic (extended if necessary) closest to \e C the center of the
* distance \e CT. To lowest order, the maximum deviation (as a true * projection and \e t be the distance \e CT. To lowest order, the maxim
* distance) of the corresponding gnomonic line segment (i.e., with the s um
ame * deviation (as a true distance) of the corresponding gnomonic line segm
* end points) from the geodesic is<br> ent
* (i.e., with the same end points) from the geodesic is<br>
* <br> * <br>
* (\e K(T) - \e K(C)) \e l<sup>2</sup> \e t / 32.<br> * (\e K(T) - \e K(C)) \e l<sup>2</sup> \e t / 32.<br>
* <br> * <br>
* where \e K is the Gaussian curvature. * where \e K is the Gaussian curvature.
* *
* This result applies for any surface. For an allipsoid of revoluation, * This result applies for any surface. For an allipsoid of revolution,
* consider all geodesics whose end points are within a distance \e r of \e * consider all geodesics whose end points are within a distance \e r of \e
* C. For a given \e r, the deviation is maximum when the latitude of \e C * C. For a given \e r, the deviation is maximum when the latitude of \e C
* is 45<sup>o</sup>, when endpoints are a distance \e r away, and when t heir * is 45<sup>o</sup>, when endpoints are a distance \e r away, and when t heir
* azimuths from the center are +/- 45<sup>o</sup> or +/- 135<sup>o</sup> . * azimuths from the center are +/- 45<sup>o</sup> or +/- 135<sup>o</sup> .
* To lowest order in \e r and the flattening \e f, the deviation is \e f * To lowest order in \e r and the flattening \e f, the deviation is \e f
* (\e r/2\e a)<sup>3</sup> \e r. * (\e r/2\e a)<sup>3</sup> \e r.
* *
* The conversions all take place using a GeographicLib::Geodesic object * The conversions all take place using a Geodesic object (by default
(by * Geodesic::WGS84). For more information on geodesics see \ref geodesic
* default GeographicLib::Geodesic::WGS84). For more information on .
* geodesics see \ref geodesic.
* *
* <b>CAUTION:</b> The definition of this projection for a sphere is * <b>CAUTION:</b> The definition of this projection for a sphere is
* standard. However, there is no standard for how it should be extended to * standard. However, there is no standard for how it should be extended to
* an ellipsoid. The choices are: * an ellipsoid. The choices are:
* - Declare that the projection is undefined for an ellipsoid. * - Declare that the projection is undefined for an ellipsoid.
* - Project to a tangent plane from the center of the ellipsoid. This * - Project to a tangent plane from the center of the ellipsoid. This
* causes great ellipses to appear as straight lines in the projection; * causes great ellipses to appear as straight lines in the projection;
* i.e., it generalizes the spherical great circle to a great ellipse. * i.e., it generalizes the spherical great circle to a great ellipse.
* This was proposed by Roy Williams, Geometry of Navigation (Horwood, * This was proposed by Roy Williams, Geometry of Navigation (Horwood,
* Chichester, 1998). * Chichester, 1998).
skipping to change at line 89 skipping to change at line 89
* 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 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, eps;
static const int numit = 10; static const int numit = 5;
public: public:
/** /**
* Constructor for Gnomonic setting the Geodesic object to use * Constructor for Gnomonic.
* for geodesic calculations. By default this uses the WGS84 ellipsoid *
. * @param[in] earth the Geodesic object to use for geodesic calculation
s.
* By default this uses the WGS84 ellipsoid.
********************************************************************** / ********************************************************************** /
explicit Gnomonic(const Geodesic& earth = Geodesic::WGS84) explicit Gnomonic(const Geodesic& earth = Geodesic::WGS84)
throw() throw()
: _earth(earth) : _earth(earth)
, _a(_earth.MajorRadius()) , _a(_earth.MajorRadius())
, _f(_earth.InverseFlattening() ? , _f(_earth.InverseFlattening() ?
1/std::abs(_earth.InverseFlattening()) : 0) 1/std::abs(_earth.InverseFlattening()) : 0)
{} {}
/** /**
* Convert from latitude \e lat (degrees) and longitude \e lon (degrees * Forward projection, from geographic to gnomonic.
) to *
* gnomonic easting \e x (meters) and northing \e y (meters). The cent * @param[in] lat0 latitude of center point of projection (degrees).
er * @param[in] lon0 longitude of center point of projection (degrees).
* of the projection is at latitude \e lat0 (degrees) and longitude \e * @param[in] lat latitude of point (degrees).
lon0 * @param[in] lon longitude of point (degrees).
* (degrees). Also return the azimuth \e azi (degrees) and the recipro * @param[out] x easting of point (meters).
cal * @param[out] y northing of point (meters).
* of the azimuthal scale \e rk. \e lat0 and \e lat should be in the r * @param[out] azi azimuth of geodesic at point (degrees).
ange * @param[out] rk reciprocal of azimuthal scale at point.
* [-90, 90] and \e lon0 and \e lon should be in the range [-180, 360]. *
* The scale of the projection is 1/\e rk<sup>2</sup> in the "radial" * \e lat0 and \e lat should be in the range [-90, 90] and \e lon0 and
* direction, \e azi clockwise from true north, and is 1/\e rk in the \e
* direction perpendicular to this. If the point lies "over the horizo * lon should be in the range [-180, 360]. The scale of the projection
n", is
* i.e., if \e rk <= 0, then NaNs are returned for \e x and \e y (the * 1/\e rk<sup>2</sup> in the "radial" direction, \e azi clockwise from
* correct values are returned for \e azi and \e rk). A call to Forwar * true north, and is 1/\e rk in the direction perpendicular to this.
d If
* followed by a call to Reverse will return the original (\e lat, \e l * the point lies "over the horizon", i.e., if \e rk <= 0, then NaNs ar
on) e
* (to within roundoff) provided the point in not over the horizon. * returned for \e x and \e y (the correct values are returned for \e a
zi
* and \e rk). A call to Forward followed by a call to Reverse will re
turn
* the original (\e lat, \e lon) (to within roundoff) provided the poin
t in
* not over the horizon.
********************************************************************** / ********************************************************************** /
void Forward(real lat0, real lon0, real lat, real lon, void Forward(real lat0, real lon0, real lat, real lon,
real& x, real& y, real& azi, real& rk) const throw(); real& x, real& y, real& azi, real& rk) const throw();
/** /**
* Convert from gnomonic easting \e x (meters) and northing \e y (meter * Reverse projection, from gnomonic to geographic.
s) *
* to latitude \e lat (degrees) and longitude \e lon (degrees). The ce * @param[in] lat0 latitude of center point of projection (degrees).
nter * @param[in] lon0 longitude of center point of projection (degrees).
* of the projection is at latitude \e lat0 (degrees) and longitude \e * @param[in] x easting of point (meters).
lon0 * @param[in] y northing of point (meters).
* (degrees). Also return the azimuth \e azi (degrees) and the recipro * @param[out] lat latitude of point (degrees).
cal * @param[out] lon longitude of point (degrees).
* of the azimuthal scale \e rk. \e lat0 should be in the range [-90, * @param[out] azi azimuth of geodesic at point (degrees).
90] * @param[out] rk reciprocal of azimuthal scale at point.
* and \e lon0 should be in the range [-180, 360]. \e lat will be in t *
he * \e lat0 should be in the range [-90, 90] and \e lon0 should be in th
* range [-90, 90] and \e lon will be in the range [-180, 180). The sc e
ale * range [-180, 360]. \e lat will be in the range [-90, 90] and \e lon
* of the projection is 1/\e rk<sup>2</sup> in the "radial" direction, * will be in the range [-180, 180). The scale of the projection is 1/
\e \e
* azi clockwise from true north, and is 1/\e rk in the direction * rk<sup>2</sup> in the "radial" direction, \e azi clockwise from true
* perpendicular to this. Even though all inputs should return a valid * north, and is 1/\e rk in the direction perpendicular to this. Even
\e * though all inputs should return a valid \e lat and \e lon, it's poss
* lat and \e lon, it's possible that the procedure fails to converge f ible
or * that the procedure fails to converge for very large \e x or \e y; in
* very large \e x or \e y; in this case NaNs are returned for all the * this case NaNs are returned for all the output arguments. A call to
* output arguments. A call to Reverse followed by a call to Forward w * Reverse followed by a call to Forward will return the original (\e x
ill , \e
* return the original (\e x, \e y) (to roundoff). * y) (to roundoff).
********************************************************************** / ********************************************************************** /
void Reverse(real lat0, real lon0, real x, real y, void Reverse(real lat0, real lon0, real x, real y,
real& lat, real& lon, real& azi, real& rk) const throw(); real& lat, real& lon, real& azi, real& rk) const throw();
/** /**
* The major radius of the ellipsoid (meters). This is that value of \ * Gnomonic::Forward without returning the azimuth and scale.
e a **********************************************************************
* inherited from the Geodesic object used in the constructor. /
void Forward(real lat0, real lon0, real lat, real lon,
real& x, real& y) const throw() {
real azi, rk;
Forward(lat0, lon0, lat, lon, x, y, azi, rk);
}
/**
* Gnomonic::Reverse without returning the azimuth and scale.
**********************************************************************
/
void Reverse(real lat0, real lon0, real x, real y,
real& lat, real& lon) const throw() {
real azi, rk;
Reverse(lat0, lon0, x, y, lat, lon, azi, rk);
}
/** \name Inspector functions
**********************************************************************
/
///@{
/**
* @return \e a the equatorial radius of the ellipsoid (meters). This
is
* the value inherited from the Geodesic object used in the construct
or.
********************************************************************** / ********************************************************************** /
Math::real MajorRadius() const throw() { return _earth.MajorRadius(); } Math::real MajorRadius() const throw() { return _earth.MajorRadius(); }
/** /**
* The inverse flattening of the ellipsoid. This is that value of \e r * @return \e r the inverse flattening of the ellipsoid. This is the
* inherited from the Geodesic object used in the constructor. A value * value inherited from the Geodesic object used in the constructor.
of A
* 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
 End of changes. 12 change blocks. 
68 lines changed or deleted 103 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 6816 2010-02-05 21:03:10Z karney $" #define GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_HPP "$Id: LambertConformalConic .hpp 6867 2010-09-11 13:04:26Z 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,
skipping to change at line 35 skipping to change at line 35
* pp. 107&ndash;109. * pp. 107&ndash;109.
* *
* 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 and that several of * that Newton's method is used to invert the projection and that several of
* the formulas are modified so that the projection correctly degenerates to * the formulas are modified so that the projection correctly degenerates to
* the Mercator projection or the polar sterographic projection. * the Mercator projection or the polar sterographic projection.
* *
* 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 els * standard parallels are set in the constructor. If the standard parall els
* are both at a single pole the projection becomes the polar stereograph ic * are both at a single pole the projection becomes the polar stereograph ic
* projection. If the standard parallels are symmetric around equator, t * projection (compare with the PolarStereographic class). If the standa
he rd
* projection becomes the Mercator projection. The central meridian (whi * parallels are symmetric around equator, the projection becomes the
ch * Mercator projection. The central meridian (which is a trivial shift o
* is a trivial shift of the longitude) is specified as the \e lon0 argum f
ent * the longitude) is specified as the \e lon0 argument of the
* of the LambertConformalConic::Forward and LambertConformalConic::Rever * LambertConformalConic::Forward and LambertConformalConic::Reverse
se
* functions. The latitude of origin is taken to be the latitude of tang ency * functions. The latitude of origin is taken to be the latitude of tang ency
* which lies between the standard parallels which is given by * which lies between the standard parallels which is given by
* LambertConformalConic::OriginLatitude. There is no provision in this * LambertConformalConic::OriginLatitude. There is no provision in this
* 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:
\verbatim \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;
{ {
double gamma, k;
// Transform origin point // Transform origin point
PASouth.Forward(lon0, lat0, lon0, x0, y0, gamma, k); 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, gamma, k; double lat, lon, x, y;
// Sample conversion from geodetic to PASouth grid // Sample conversion from geodetic to PASouth grid
std::cin >> lat >> lon; std::cin >> lat >> lon;
PASouth.Forward(lon0, lat, lon, x, y, gamma, k); 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, gamma, k); PASouth.Reverse(lon0, x, y, lat, lon);
std::cout << lat << " " << lon << "\n"; std::cout << lat << " " << lon << "\n";
\endverbatim \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, _e2, _e, _e2m;
real _sign, _n, _nc, _lt0, _t0n, _t0nm1, _scale, _lat0, _k0; real _sign, _n, _nc, _lt0, _t0n, _t0nm1, _scale, _lat0, _k0;
static const real eps, eps2, epsx, tol, ahypover; static const real eps, eps2, 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; }
// 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
skipping to change at line 114 skipping to change at line 114
} }
inline real logmtf(real sphi) const throw() { inline real logmtf(real sphi) const throw() {
// Return log(m/t). This allows the cancellation of cphi in m and t. // Return log(m/t). This allows the cancellation of cphi in m and t.
return std::log((1 + sphi)/std::sqrt(1 - _e2 * sq(sphi))) - return std::log((1 + sphi)/std::sqrt(1 - _e2 * sq(sphi))) -
eatanhe(sphi); eatanhe(sphi);
} }
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 for a ellipsoid radius \e a (meters), reciprocal flatten * Constructor with a single standard parallel.
ing *
* \e r, standard parallel (the circle of tangency) \e stdlat (degrees) * @param[in] a equatorial radius of ellipsoid (meters)
, * @param[in] r reciprocal flattening of ellipsoid. Setting \e r = 0
* and scale on the standard parallel \e k0. Setting \e r = 0 implies * implies \e r = inf or flattening = 0 (i.e., a sphere). Negative \
\e r e r
* = inf or flattening = 0 (i.e., a sphere). An exception is thrown if * indicates a prolate ellipsoid.
\e a * @param[in] stdlat standard parallel (degrees), the circle of tangenc
* or \e k0 is not positive or if \e stdlat is not in the range [-90, 9 y.
0]. * @param[in] k0 scale on standard parallel.
*
* An exception is thrown if \e a or \e k0 is not positive or if \e std
lat
* is not in the range [-90, 90].
********************************************************************** / ********************************************************************** /
LambertConformalConic(real a, real r, real stdlat, real k0); LambertConformalConic(real a, real r, real stdlat, real k0);
/** /**
* Constructor for a ellipsoid radius \e a (meters), reciprocal flatten * Constructor with two standard parallels.
ing *
* \e r, standard parallels \e stdlat1 (degrees) and \e stdlat2 (degree * @param[in] a equatorial radius of ellipsoid (meters)
s), * @param[in] r reciprocal flattening of ellipsoid. Setting \e r = 0
* and the scale on the standard parallels \e k1. Setting \e r = 0 imp * implies \e r = inf or flattening = 0 (i.e., a sphere). Negative \
lies e r
* \e r = inf or flattening = 0 (i.e., a sphere). An exception is throw * indicates a prolate ellipsoid.
n if * @param[in] stdlat1 first standard parallel (degrees).
* \e a or \e k0 is not positive or if \e stdlat1 or \e stdlat2 is not * @param[in] stdlat2 second standard parallel (degrees).
in * @param[in] k1 scale on the standard parallels.
* the range [-90, 90]. In addition, if either \e stdlat1 or \e stdlat *
2 is * An exception is thrown if \e a or \e k0 is not positive or if \e std
* a pole, then an exception is thrown if \e stdlat1 is not equal \e lat1
* stdlat2 * or \e stdlat2 is not in the range [-90, 90]. In addition, if either
\e
* stdlat1 or \e stdlat2 is a pole, then an exception is thrown if \e
* stdlat1 is not equal \e stdlat2
********************************************************************** / ********************************************************************** /
LambertConformalConic(real a, real r, real stdlat1, real stdlat2, real k1); LambertConformalConic(real a, real r, real stdlat1, real stdlat2, real k1);
/** /**
* An alternative constructor for 2 standard parallels where the parall * Constructor with two standard parallels specified by sines and cosin
els es.
* are given by their sines and cosines. This allows parallels close t *
o * @param[in] a equatorial radius of ellipsoid (meters)
* the poles to be specified accurately. * @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
* indicates a prolate ellipsoid.
* @param[in] sinlat1 sine of first standard parallel.
* @param[in] coslat1 cosine of first standard parallel.
* @param[in] sinlat2 sine of second standard parallel.
* @param[in] coslat2 cosine of second standard parallel.
* @param[in] k1 scale on the standard parallels.
*
* This allows parallels close to the poles to be specified accurately.
********************************************************************** / ********************************************************************** /
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);
/** /**
* Alter the scale for the projection so that on latitude \e lat, the s * Set the scale for the projection.
cale *
* is \e k (default 1). The allows a "latitude of true scale" to be * @param[in] lat (degrees).
* specified. An exception is thrown if \e k is not positive. * @param[in] k scale at latitude \e lat (default 1).
*
* This allows a "latitude of true scale" to be specified. An exceptio
n is
* thrown if \e k is not positive.
********************************************************************** / ********************************************************************** /
void SetScale(real lat, real k = real(1)); void SetScale(real lat, real k = real(1));
/** /**
* Convert from latitude \e lat (degrees) and longitude \e lon (degrees * Forward projection, from geographic to Lambert conformal conic.
) to *
* Lambert conformal conic easting \e x (meters) and northing \e y * @param[in] lon0 central meridian longitude (degrees).
* (meters). The central meridian is given by \e lon0 (degrees) and th * @param[in] lat latitude of point (degrees).
e * @param[in] lon longitude of point (degrees).
* latitude origin is given by LambertConformalConic::LatitudeOrigin(). * @param[out] x easting of point (meters).
* Also return the meridian convergence \e gamma (degrees) and the scal * @param[out] y northing of point (meters).
e \e * @param[out] gamma meridian convergence at point (degrees).
* k. No false easting or northing is added and \e lat should be in th * @param[out] k scale of projection at point.
e *
* range [-90, 90]; \e lon and \e lon0 should be in the range [-180, 36 * The latitude origin is given by LambertConformalConic::LatitudeOrigi
0]. n().
* The values of \e x and \e y returned for points which project to * No false easting or northing is added and \e lat should be in the ra
* infinity (i.e., one or both of the poles) will be large but finite. nge
The * [-90, 90]; \e lon and \e lon0 should be in the range [-180, 360]. T
* value of \e k returned for the poles in infinite (unless \e lat equa he
ls * values of \e x and \e y returned for points which project to infinit
* the latitude origin). y
* (i.e., one or both of the poles) will be large but finite. The valu
e of
* \e k returned for the poles in infinite (unless \e lat equals the
* latitude origin).
********************************************************************** / ********************************************************************** /
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();
/** /**
* Convert from Lambert conformal conic easting \e x (meters) and north * Reverse projection, from Lambert conformal conic to geographic.
ing *
* \e y (meters) to latitude \e lat (degrees) and longitude \e lon * @param[in] lon0 central meridian longitude (degrees).
* (degrees) . The central meridian is given by \e lon0 (degrees) and * @param[in] x easting of point (meters).
the * @param[in] y northing of point (meters).
* latitude origin is given by LambertConformalConic::LatitudeOrigin(). * @param[out] lat latitude of point (degrees).
* Also return the meridian convergence \e gamma (degrees) and the scal * @param[out] lon longitude of point (degrees).
e \e * @param[out] gamma meridian convergence at point (degrees).
* k. No false easting or northing is added. \e lon0 should be in the * @param[out] k scale of projection at point.
* range [-180, 360]. The value of \e lon returned is in the range [-1 *
80, * The latitude origin is given by LambertConformalConic::LatitudeOrigi
* 180). n().
* 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).
********************************************************************** / ********************************************************************** /
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();
/** /**
* The major radius of the ellipsoid (meters). This is that value of \ * LambertConformalConic::Forward without returning the convergence and
e a * scale.
* used in the constructor. **********************************************************************
/
void Forward(real lon0, real lat, real lon,
real& x, real& y) const throw() {
real gamma, k;
Forward(lon0, lat, lon, x, y, gamma, k);
}
/**
* LambertConformalConic::Reverse without returning the convergence and
* scale.
**********************************************************************
/
void Reverse(real lon0, real x, real y,
real& lat, real& lon) const throw() {
real gamma, k;
Reverse(lon0, x, y, lat, lon, gamma, k);
}
/** \name Inspector functions
**********************************************************************
/
///@{
/**
* @return \e a the equatorial radius of the ellipsoid (meters). This
is
* the value used in the constructor.
********************************************************************** / ********************************************************************** /
Math::real MajorRadius() const throw() { return _a; } Math::real MajorRadius() const throw() { return _a; }
/** /**
* The inverse flattening of the ellipsoid. This is that value of \e r * @return \e r the inverse flattening of the ellipsoid. This is the
* used in the constructor. A value of 0 is returned for a sphere * value used in the constructor. A value of 0 is returned for a sph
* (infinite inverse flattening). ere
* (infinite inverse flattening).
********************************************************************** / ********************************************************************** /
Math::real InverseFlattening() const throw() { return _r; } Math::real InverseFlattening() const throw() { return _r; }
/** /**
* The latitude of the origin for the projection (degrees). This is th * @return latitude of the origin for the projection (degrees).
at *
* latitude of minimum scale and equals the \e stdlat in the 1-parallel * This is the latitude of minimum scale and equals the \e stdlat in th
* constructor and lies between \e stdlat1 and \e stdlat2 in the 2-para e
llel * 1-parallel constructor and lies between \e stdlat1 and \e stdlat2 in
* constructors. the
* 2-parallel constructors.
********************************************************************** / ********************************************************************** /
Math::real OriginLatitude() const throw() { return _lat0; } Math::real OriginLatitude() const throw() { return _lat0; }
/** /**
* The central scale for the projection. This is that 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 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; const static LambertConformalConic Mercator;
}; };
} // namespace GeographicLib } // namespace GeographicLib
 End of changes. 20 change blocks. 
93 lines changed or deleted 149 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) <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_LOCALCARTESIAN_HPP) #if !defined(GEOGRAPHICLIB_LOCALCARTESIAN_HPP)
#define GEOGRAPHICLIB_LOCALCARTESIAN_HPP "$Id: LocalCartesian.hpp 6785 2010 -01-05 22:15:42Z karney $" #define GEOGRAPHICLIB_LOCALCARTESIAN_HPP "$Id: LocalCartesian.hpp 6867 2010 -09-11 13:04:26Z karney $"
#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
* GeographicLib::Geocentric object (by default * Geocentric object (by default Geocentric::WGS84).
* GeographicLib::Geocentric::WGS84).
**********************************************************************/ **********************************************************************/
class LocalCartesian { class LocalCartesian {
private: private:
typedef Math::real real; typedef Math::real real;
const Geocentric _earth; const Geocentric _earth;
real _lat0, _lon0, _h0; real _lat0, _lon0, _h0;
real _x0, _y0, _z0, real _x0, _y0, _z0,
_rxx, _rxy, _rxz, _rxx, _rxy, _rxz,
_ryx, _ryy, _ryz, _ryx, _ryy, _ryz,
_rzx, _rzy, _rzz; _rzx, _rzy, _rzz;
public: public:
/** /**
* Constructor setting the origin to latitude = \e lat0, longitude = \e * Constructor setting the origin.
* lon0 (degrees), height = \e h0 (meters). The optional \e earth argu *
ment * @param[in] lat0 latitude at origin (degrees).
* (default Geocentric::WGS84) specifies the Geocentric object to use f * @param[in] lon0 longitude at origin (degrees).
or * @param[in] h0 height above ellipsoid at origin (meters); default 0.
* the transformation. * @param[in] earth Geocentric object for the transformation; default
* Geocentric::WGS84.
********************************************************************** / ********************************************************************** /
LocalCartesian(real lat0, real lon0, real h0 = 0, LocalCartesian(real lat0, real lon0, real h0 = 0,
const Geocentric& earth = Geocentric::WGS84) throw() const Geocentric& earth = Geocentric::WGS84) throw()
: _earth(earth) : _earth(earth)
{ Reset(lat0, lon0, h0); } { Reset(lat0, lon0, h0); }
/** /**
* Default constructor sets the origin to \e lat0 = 0, \e lon0 = 0, \e * Default constructor.
h0 = *
* 0. The optional \e earth argument (default Geocentric::WGS84) speci * @param[in] earth Geocentric object for the transformation; default
fies * Geocentric::WGS84.
* the Geocentric object to use for the transformation. *
* Sets \e lat0 = 0, \e lon0 = 0, \e h0 = 0.
********************************************************************** / ********************************************************************** /
explicit LocalCartesian(const Geocentric& earth = Geocentric::WGS84) explicit LocalCartesian(const Geocentric& earth = Geocentric::WGS84)
throw() throw()
: _earth(earth) : _earth(earth)
{ Reset(real(0), real(0), real(0)); } { Reset(real(0), real(0), real(0)); }
/** /**
* Change the origin to latitude = \e lat0, longitude = \e lon0 (degree * Reset the origin.
s), *
* height = \e h0 (meters). * @param[in] lat0 latitude at origin (degrees).
* @param[in] lon0 longitude at origin (degrees).
* @param[in] h0 height above ellipsoid at origin (meters); default 0.
********************************************************************** / ********************************************************************** /
void Reset(real lat0, real lon0, real h0 = 0) void Reset(real lat0, real lon0, real h0 = 0)
throw(); throw();
/** /**
* Convert from geodetic coordinates \e lat, \e lon (degrees), \e h * Convert from geodetic to local cartesian coordinates.
* (meters) to local cartesian coordinates \e x, \e y, \e z (meters). *
\e * @param[in] lat latitude of point (degrees).
* lat should be in the range [-90, 90]; \e lon and \e lon0 should be i * @param[in] lon longitude of point (degrees).
n * @param[in] h height of point above the ellipsoid (meters).
* @param[out] x local cartesian coordinate (meters).
* @param[out] y local cartesian coordinate (meters).
* @param[out] z local cartesian coordinate (meters).
*
* \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].
********************************************************************** / ********************************************************************** /
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)
const throw(); const throw();
/** /**
* Convert from local cartesian \e x, \e y, \e z (meters) to geodetic * Convert from local cartesian to geodetic to coordinates.
* coordinates \e lat, \e lon (degrees), \e h (meters). The value of \ *
e * @param[in] x local cartesian coordinate (meters).
* lon returned is in the range [-180, 180). * @param[in] y local cartesian coordinate (meters).
* @param[in] z local cartesian coordinate (meters).
* @param[out] lat latitude of point (degrees).
* @param[out] lon longitude of point (degrees).
* @param[out] h height of point above the ellipsoid (meters).
*
* The value of \e lon returned is in the range [-180, 180).
********************************************************************** / ********************************************************************** /
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)
const throw(); const throw();
/** \name Inspector functions
**********************************************************************
/
///@{
/** /**
* Return the 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; }
/** /**
* Return the longitude of the origin (degrees). * @return longitude of the origin (degrees).
********************************************************************** / ********************************************************************** /
Math::real LongitudeOrigin() const throw() { return _lon0; } Math::real LongitudeOrigin() const throw() { return _lon0; }
/** /**
* Return the height of the origin (meters). * @return height of the origin (meters).
********************************************************************** / ********************************************************************** /
Math::real HeightOrigin() const throw() { return _h0; } Math::real HeightOrigin() const throw() { return _h0; }
/** /**
* The major radius of the ellipsoid (meters). This is that value of \ * @return \e a the equatorial radius of the ellipsoid (meters). This
e a is
* inherited from the Geocentric object used in the constructor. * the value of \e a inherited from the Geocentric object used in the
* constructor.
********************************************************************** / ********************************************************************** /
Math::real MajorRadius() const throw() { return _earth.MajorRadius(); } Math::real MajorRadius() const throw() { return _earth.MajorRadius(); }
/** /**
* The inverse flattening of the ellipsoid. This is that value of \e r * @return \e r the inverse flattening of the ellipsoid. This is the
* inherited from the Geocentric object used in the constructor. A val * value of \e r inherited from the Geocentric object used in the
ue * constructor. A value of 0 is returned for a sphere (infinite inve
* of 0 is returned for a sphere (infinite inverse flattening). rse
* flattening).
********************************************************************** / ********************************************************************** /
Math::real InverseFlattening() const throw() Math::real InverseFlattening() const throw()
{ return _earth.InverseFlattening(); } { return _earth.InverseFlattening(); }
///@}
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif #endif
 End of changes. 14 change blocks. 
36 lines changed or deleted 58 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 6843 2010-07-18 21:18:17Z kar ney $" #define GEOGRAPHICLIB_MGRS_HPP "$Id: MGRS.hpp 6867 2010-09-11 13:04:26Z 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 123 skipping to change at line 123
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. \e zone and \e * Convert UTM or UPS coordinate to an MGRS coordinate.
* northp give input zone (with \e zone = 0 indicating UPS) and hemisph *
ere, * @param[in] zone UTM zone (zero means UPS)
* \e x and \e y are the easting and northing (meters). \e prec indica * @param[in] northp hemisphere (true means north, false means south)
tes * @param[in] x (meters)
* the desired precision with \e prec = 0 (the minimum) meaning 100 km, * @param[in] y (meters)
\e * @param[in] prec precision relative to 100 km
* prec = 5 meaning 1 m, and \e prec == 11 (the maximum) meaning 1 um. * @param[out] mgrs MGRS string.
*
* \e prec specifies the precision of the MSGRS string as follows:
* - prec = 0 (min), 100km
* - prec = 1, 10km
* - prec = 2, 1km
* - prec = 3, 100m
* - prec = 4, 10m
* - prec = 5, 1m
* - prec = 6, 0.1m
* - prec = 11 (max), 1um
* *
* UTM eastings are allowed to be in the range [100 km, 900 km], northi ngs * UTM eastings are allowed to be in the range [100 km, 900 km], northi ngs
* are allowed to be in in [0 km, 9500 km] for the northern hemisphere and * are allowed to be in in [0 km, 9500 km] for the northern hemisphere and
* in [1000 km, 10000 km] for the southern hemisphere. (However UTM * in [1000 km, 10000 km] for the southern hemisphere. (However UTM
* northings can be continued across the equator. So the actual limits on * northings can be continued across the equator. So the actual limits on
* the northings are [-9000 km, 9500 km] for the "northern" hemisphere and * the northings are [-9000 km, 9500 km] for the "northern" hemisphere and
* [1000 km, 19500 km] for the "southern" hemisphere.) * [1000 km, 19500 km] for the "southern" hemisphere.)
* *
* UPS eastings/northings are allowed to be in the range [1300 km, 2700 km] * UPS eastings/northings are allowed to be in the range [1300 km, 2700 km]
* in the northern hemisphere and in [800 km, 3200 km] in the southern * in the northern hemisphere and in [800 km, 3200 km] in the southern
skipping to change at line 179 skipping to change at line 192
* roundoff. * roundoff.
* *
* Return the result via a reference argument to avoid the overhead of * Return the result via a reference argument to avoid the overhead of
* allocating a potentially large number of small strings. If an error is * allocating a potentially large number of small strings. If an error is
* thrown, then \e mgrs is unchanged. * thrown, then \e mgrs is unchanged.
********************************************************************** / ********************************************************************** /
static void Forward(int zone, bool northp, real x, real y, static void Forward(int zone, bool northp, real x, real y,
int prec, std::string& mgrs); int prec, std::string& mgrs);
/** /**
* Convert UTM or UPS coordinates to an MGRS coordinate in case that * Convert UTM or UPS coordinate to an MGRS coordinate when the latitud
* latitude is already known. The latitude is ignored for \e zone = 0 e is
* (UPS); otherwise the latitude is used to determine the latitude band * known.
and *
* this is checked for consistency using the same tests as Reverse. * @param[in] zone UTM zone (zero means UPS).
* @param[in] northp hemisphere (true means north, false means south).
* @param[in] x (meters).
* @param[in] y (meters).
* @param[in] lat latitude (degrees).
* @param[in] prec precision relative to 100 km.
* @param[out] mgrs MGRS string.
*
* The latitude is ignored for \e zone = 0 (UPS); otherwise the latitud
e is
* used to determine the latitude band and this is checked for consiste
ncy
* using the same tests as Reverse.
********************************************************************** / ********************************************************************** /
static void Forward(int zone, bool northp, real x, real y, real lat, static void Forward(int zone, bool northp, real x, real y, real lat,
int prec, std::string& mgrs); int prec, std::string& mgrs);
/** /**
* Convert a MGRS coordinate to UTM or UPS coordinates returning zone \ * Convert a MGRS coordinate to UTM or UPS coordinates.
e *
* zone, hemisphere \e northp, easting \e x (meters), northing \e y * @param[in] mgrs MGRS string.
* (meters) . Also return the precision of the MGRS string (see Forwar * @param[out] zone UTM zone (zero means UPS).
d). * @param[out] northp hemisphere (true means north, false means south).
* If \e centerp = true (default), return center of the MGRS square, el * @param[out] x (meters).
se * @param[out] y (meters).
* return SW (lower left) corner. * @param[out] prec precision relative to 100 km.
* @param[in] centerp if true(default), return center of the MGRS squar
e,
* 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)
* - 38MPE and 38NPF are permitted (they straddle the equator); but 3 8NPE * - 38MPE and 38NPF are permitted (they straddle the equator); but 3 8NPE
skipping to change at line 220 skipping to change at line 248
* from MGRS coordinate. The conversion is exact for prec in [0, 5]. With * from MGRS coordinate. The conversion is exact for prec in [0, 5]. With
* centerp = true the conversion from MGRS to geographic and back is * centerp = true the conversion from MGRS to geographic and back is
* stable. This is not assured if \e centerp = false. * stable. This is not assured if \e centerp = false.
* *
* If an error is thrown, then the arguments are unchanged. * If an error is thrown, then the arguments are unchanged.
********************************************************************** / ********************************************************************** /
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
**********************************************************************
/
///@{
/** /**
* The major radius of the ellipsoid (meters). This is the value for t * @return \e a the equatorial radius of the WGS84 ellipsoid (meters).
he *
* WGS84 ellipsoid because MGRS coordinates are based on this ellipsoid * (The WGS84 values is returned because the UTM and UPS projections ar
. e
* based on this ellipsoid.)
********************************************************************** / ********************************************************************** /
static Math::real MajorRadius() throw() { return UTMUPS::MajorRadius() ; } static Math::real MajorRadius() throw() { return UTMUPS::MajorRadius() ; }
/** /**
* The inverse flattening of the ellipsoid. This is the value for the * @return \e r the inverse flattening of the WGS84 ellipsoid.
* WGS84 ellipsoid because MGRS coordinates are based on this ellipsoid *
. * (The WGS84 values is returned because the UTM and UPS projections ar
e
* 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. 
29 lines changed or deleted 62 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 841 2010-07-11 20:46:32Z karney $" #define GEOGRAPHICLIB_POLARSTEREOGRAPHIC_HPP "$Id: PolarStereographic.hpp 6 867 2010-09-11 13:04:26Z 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,
* <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 PolarStereographic {
private: private:
typedef Math::real real; typedef Math::real real;
const real _a, _r, _f, _e2, _e, _e2m, _C, _c; // _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;
real _k0; real _k0;
static const real tol, overflow; static const real tol, overflow;
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; }
// 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 radius \e a (meters), reciprocal flatten * Constructor for a ellipsoid with
ing *
* \e r, and central scale factor \e k0. Setting \e r <= 0 implies \e * @param[in] a equatorial radius (meters)
r = * @param[in] r reciprocal flattening. Setting \e r = 0 implies \e r =
* inf or flattening = 0 (i.e., a sphere). An exception is thrown if \ inf
e a * or flattening = 0 (i.e., a sphere). Negative \e r indicates a pro
* or \e k0 is not positive. late
* ellipsoid.
* @param[in] k0 central scale factor.
*
* An exception is thrown if either of the axes of the ellipsoid is
* not positive \e a or if \e k0 is not positive.
********************************************************************** / ********************************************************************** /
PolarStereographic(real a, real r, real k0); PolarStereographic(real a, real r, real k0);
/** /**
* Alter the scale for the projection so that on latitude \e lat, the s * Set the scale for the projection.
cale *
* is \e k (default 1). \e lat is interpreted in the context a project * @param[in] lat (degrees).
ion * @param[in] k scale at latitude \e lat (default 1).
* centered at the north pole. The allows a "latitude of true scale" t *
o be * This allows a "latitude of true scale" to be specified. An exceptio
* specified. An exception is thrown if \e k is not positive. n is
* thrown if \e k is not positive.
********************************************************************** / ********************************************************************** /
void SetScale(real lat, real k = real(1)); void SetScale(real lat, real k = real(1));
/** /**
* Convert from latitude \e lat (degrees) and longitude \e lon (degrees * Forward projection, from geographic to polar stereographic.
) to *
* polar stereographic easting \e x (meters) and northing \e y (meters) * @param[in] northp the pole which is the center of projection (true m
. eans
* The projection is about the pole given by \e northp (false means sou * north, false means south).
th, * @param[in] lat latitude of point (degrees).
* true means north). Also return the meridian convergence \e gamma * @param[in] lon longitude of point (degrees).
* (degrees) and the scale \e k. No false easting or northing is added * @param[out] x easting of point (meters).
. * @param[out] y northing of point (meters).
* \e lat should be in the range (-90, 90] for \e northp = true and in * @param[out] gamma meridian convergence at point (degrees).
the * @param[out] k scale of projection at point.
* range [-90, 90) for \e northp = false; \e lon should be in the range *
* [-180, 360]. * No false easting or northing is added. \e lat should be in the rang
e
* (-90, 90] for \e northp = true and in the range [-90, 90) for \e nor
thp
* = false; \e lon should be in the range [-180, 360].
********************************************************************** / ********************************************************************** /
void Forward(bool northp, real lat, real lon, void Forward(bool northp, real lat, real lon,
real& x, real& y, real& gamma, real& k) const throw(); real& x, real& y, real& gamma, real& k) const throw();
/** /**
* Convert from polar stereogrphic easting \e x (meters) and northing \ * Reverse projection, from polar stereographic to geographic.
e y *
* (meters) to latitude \e lat (degrees) and longitude \e lon (degrees) * @param[in] northp the pole which is the center of projection (true m
. eans
* The hemisphere is given by \e northp (false means south, true means * north, false means south).
* north). Also return the meridian convergence \e gamma (degrees) and * @param[in] x easting of point (meters).
the * @param[in] y northing of point (meters).
* scale \e k. No false easting or northing is added. The value of \e * @param[out] lat latitude of point (degrees).
lon * @param[out] lon longitude of point (degrees).
* returned is in the range [-180, 180). * @param[out] gamma meridian convergence at point (degrees).
* @param[out] k scale of projection at point.
*
* No false easting or northing is added. The value of \e lon returned
is
* in the range [-180, 180).
********************************************************************** / ********************************************************************** /
void Reverse(bool northp, real x, real y, void Reverse(bool northp, real x, real y,
real& lat, real& lon, real& gamma, real& k) const throw(); real& lat, real& lon, real& gamma, real& k) const throw();
/** /**
* The major radius of the ellipsoid (meters). This is that value of \ * PolarStereographic::Forward without returning the convergence and sc
e a ale.
* used in the constructor. **********************************************************************
/
void Forward(bool northp, real lat, real lon,
real& x, real& y) const throw() {
real gamma, k;
Forward(northp, lat, lon, x, y, gamma, k);
}
/**
* PolarStereographic::Reverse without returning the convergence and sc
ale.
**********************************************************************
/
void Reverse(bool northp, real x, real y,
real& lat, real& lon) const throw() {
real gamma, k;
Reverse(northp, x, y, lat, lon, gamma, k);
}
/** \name Inspector functions
**********************************************************************
/
///@{
/**
* @return \e a the equatorial radius of the ellipsoid (meters). This
is
* the value used in the constructor.
********************************************************************** / ********************************************************************** /
Math::real MajorRadius() const throw() { return _a; } Math::real MajorRadius() const throw() { return _a; }
/** /**
* The inverse flattening of the ellipsoid. This is that value of \e r * @return \e r the inverse flattening of the ellipsoid. This is the
* used in the constructor. A value of 0 is returned for a sphere * value used in the constructor. A value of 0 is returned for a sph
* (infinite inverse flattening). ere
* (infinite inverse flattening).
********************************************************************** / ********************************************************************** /
Math::real InverseFlattening() const throw() { return _r; } Math::real InverseFlattening() const throw() { return _r; }
/** /**
* The central scale for the projection. This is that value of \e k0 u * The central scale for the projection. This is the value of \e k0 us
sed ed
* in the constructor and is the scale at the pole. * in the constructor and is the scale at the pole unless overridden by
* 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; const static PolarStereographic UPS;
}; };
} // namespace GeographicLib } // namespace GeographicLib
 End of changes. 10 change blocks. 
48 lines changed or deleted 93 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 843 2010-07-18 21:18:17Z karney $" #define GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP "$Id: TransverseMercator.hpp 6 867 2010-09-11 13:04:26Z 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 * The order of the series approximation used in TransverseMercator.
* GeographicLib::TransverseMercator. TM_TX_MAXPOW can be set to any integ * TM_TX_MAXPOW can be set to any integer in [4, 8].
er
* 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
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief Transverse Mercator Projection * \brief Transverse Mercator Projection
* *
skipping to change at line 43 skipping to change at line 42
* <a href="http://dx.doi.org/10.2312/GFZ.b103-krueger28"> Konforme * <a href="http://dx.doi.org/10.2312/GFZ.b103-krueger28"> Konforme
* Abbildung des Erdellipsoids in der Ebene</a> (Conformal mapping of the * Abbildung des Erdellipsoids in der Ebene</a> (Conformal mapping of the
* ellipsoidal earth to the plane), Royal Prussian Geodetic Institute, New * ellipsoidal earth to the plane), Royal Prussian Geodetic Institute, New
* Series 52, 172 pp. (1912). * Series 52, 172 pp. (1912).
* *
* Kr&uuml;ger's method has been extended from 4th to 6th order. The max imum * Kr&uuml;ger's method has been extended from 4th to 6th order. The max imum
* errors is 5 nm (ground distance) for all positions within 35 degrees o f * errors is 5 nm (ground distance) for all positions within 35 degrees o f
* the central meridian. The error in the convergence is 2e-15&quot; and the * the central meridian. The error in the convergence is 2e-15&quot; and the
* relative error in the scale is 6e-12%%. (See \ref tmerrors for the we asel * relative error in the scale is 6e-12%%. (See \ref tmerrors for the we asel
* words.) The speed penalty in going to 6th order is only about 1%. * words.) The speed penalty in going to 6th order is only about 1%.
* GeographicLib::TransverseMercatorExact is an alternative implementatio * TransverseMercatorExact is an alternative implementation of the projec
n of tion
* the projection using exact formulas which yield accurate (to 8 nm) res * using exact formulas which yield accurate (to 8 nm) results over the
ults * entire ellipsoid.
* over the entire ellipsoid.
* *
* The ellipsoid parameters and the central scale are set in the construc tor. * The ellipsoid parameters and the central scale are set in the construc tor.
* 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 * EPSG:7405</a>) which requires the use of a latitude of origin. This i
s accommodated s
* by (constants from * accommodated by (constants from
* <a href="http://www.ordnancesurvey.co.uk/oswebsite/gps/information/coo rdinatesystemsinfo/guidecontents/guidea.html"> * <a href="http://www.ordnancesurvey.co.uk/oswebsite/gps/information/coo rdinatesystemsinfo/guidecontents/guidea.html">
* A guide to coordinate systems in Great Britain</a>): * A guide to coordinate systems in Great Britain</a>):
\verbatim \code
const double const double
a = 6377563.396, b = 6356256.910, r = a/(a - b), // Airy 1830 ellipsoi d 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 k0 = 0.9996012717, lat0 = 49, lon0 = -2, // central scale and origin
fe = 400000, fn = -100000; // false easting and northing fe = 400000, fn = -100000; // false easting and northing
// Set up basic projection // Set up basic projection
const GeographicLib::TransverseMercator OSGB(a, r, k0); const GeographicLib::TransverseMercator OSGB(a, r, k0);
double x0, y0; double x0, y0;
{ {
double gamma, k;
// Transform origin point // Transform origin point
OSGB.Forward(lon0, lat0, lon0, x0, y0, gamma, k); OSGB.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, gamma, k; double lat, lon, x, y;
// Sample conversion from geodetic to OSGB grid // Sample conversion from geodetic to OSGB grid
std::cin >> lat >> lon; std::cin >> lat >> lon;
OSGB.Forward(lon0, lat, lon, x, y, gamma, k); OSGB.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 OSGB grid to geodetic // Sample conversion from OSGB grid to geodetic
std::cin >> x >> y; std::cin >> x >> y;
x += x0; y += y0; x += x0; y += y0;
OSGB.Reverse(lon0, x, y, lat, lon, gamma, k); OSGB.Reverse(lon0, x, y, lat, lon);
std::cout << lat << " " << lon << "\n"; std::cout << lat << " " << lon << "\n";
\endverbatim \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 116 skipping to change at line 114
return x >= 0 ? (t >= 0 ? t : overflow) : (t < 0 ? t : -overflow); return x >= 0 ? (t >= 0 ? t : overflow) : (t < 0 ? t : -overflow);
} }
// 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 radius \e a (meters), reciprocal flatten * Constructor for a ellipsoid with
ing *
* \e r, and central scale factor \e k0. Setting \e r = 0 implies \e r * @param[in] a equatorial radius (meters)
= * @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 inf
* prolate spheroid. An exception is thrown if \e a or \e k0 is * or flattening = 0 (i.e., a sphere). Negative \e r indicates a pro
* not positive. late
* ellipsoid.
* @param[in] k0 central scale factor.
*
* An exception is thrown if either of the axes of the ellipsoid or \e
k0
* is not positive.
********************************************************************** / ********************************************************************** /
TransverseMercator(real a, real r, real k0); TransverseMercator(real a, real r, real k0);
/** /**
* Convert from latitude \e lat (degrees) and longitude \e lon (degrees * Forward projection, from geographic to transverse Mercator.
) to *
* transverse Mercator easting \e x (meters) and northing \e y (meters) * @param[in] lon0 central meridian of the projection (degrees).
. * @param[in] lat latitude of point (degrees).
* The central meridian of the transformation is \e lon0 (degrees). Al * @param[in] lon longitude of point (degrees).
so * @param[out] x easting of point (meters).
* return the meridian convergence \e gamma (degrees) and the scale \e * @param[out] y northing of point (meters).
k. * @param[out] gamma meridian convergence at point (degrees).
* @param[out] k scale of projection at point.
*
* No false easting or northing is added. \e lat should be in the range * No false easting or northing is added. \e lat should be in the range
* [-90, 90]; \e lon and \e lon0 should be in the range [-180, 360]. * [-90, 90]; \e lon and \e lon0 should be in the range [-180, 360].
********************************************************************** / ********************************************************************** /
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();
/** /**
* Convert from transverse Mercator easting \e x (meters) and northing * Reverse projection, from transverse Mercator to geographic.
\e y *
* (meters) to latitude \e lat (degrees) and longitude \e lon (degrees) * @param[in] lon0 central meridian of the projection (degrees).
. * @param[in] x easting of point (meters).
* The central meridian of the transformation is \e lon0 (degrees). Al * @param[in] y northing of point (meters).
so * @param[out] lat latitude of point (degrees).
* return the meridian convergence \e gamma (degrees) and the scale \e * @param[out] lon longitude of point (degrees).
k. * @param[out] gamma meridian convergence at point (degrees).
* @param[out] k scale of projection at point.
*
* 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).
********************************************************************** / ********************************************************************** /
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();
/** /**
* The major radius of the ellipsoid (meters). This is that value of \ * TransverseMercator::Forward without returning the convergence and sc
e a ale.
* used in the constructor. **********************************************************************
/
void Forward(real lon0, real lat, real lon,
real& x, real& y) const throw() {
real gamma, k;
Forward(lon0, lat, lon, x, y, gamma, k);
}
/**
* TransverseMercator::Reverse without returning the convergence and sc
ale.
**********************************************************************
/
void Reverse(real lon0, real x, real y,
real& lat, real& lon) const throw() {
real gamma, k;
Reverse(lon0, x, y, lat, lon, gamma, k);
}
/** \name Inspector functions
**********************************************************************
/
///@{
/**
* @return \e a the equatorial radius of the ellipsoid (meters). This
is
* the value used in the constructor.
********************************************************************** / ********************************************************************** /
Math::real MajorRadius() const throw() { return _a; } Math::real MajorRadius() const throw() { return _a; }
/** /**
* The inverse flattening of the ellipsoid. This is that value of \e r * @return \e r the inverse flattening of the ellipsoid. This is the
* used in the constructor. A value of 0 is returned for a sphere * value used in the constructor. A value of 0 is returned for a sph
* (infinite inverse flattening). ere
* (infinite inverse flattening).
********************************************************************** / ********************************************************************** /
Math::real InverseFlattening() const throw() { return _r; } Math::real InverseFlattening() const throw() { return _r; }
/** /**
* The central scale for the projection. This is that value of \e k0 u * @return \e k0 central scale for the projection. This is the value o
sed f \e
* in the constructor and is the scale on the central meridian. * 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; const static TransverseMercator UTM;
}; };
} // namespace GeographicLib } // namespace GeographicLib
 End of changes. 18 change blocks. 
52 lines changed or deleted 87 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 6835 2010-06-15 21:37:16Z karney $" #define GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_HPP "$Id: TransverseMercatorE xact.hpp 6867 2010-09-11 13:04:26Z 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 44 skipping to change at line 44
* nm (ground distance) for the forward and reverse transformations. The * nm (ground distance) for the forward and reverse transformations. The
* error in the convergence is 2e-15&quot;, the relative error in the sca le * error in the convergence is 2e-15&quot;, the relative error in the sca le
* is 7e-12%%. (See \ref tmerrors for the weasel words.) The method is * is 7e-12%%. (See \ref tmerrors for the weasel words.) The method is
* "exact" in the sense that the errors are close to the round-off limit and * "exact" in the sense that the errors are close to the round-off limit and
* that no changes are needed in the algorithms for them to be used with * that no changes are needed in the algorithms for them to be used with
* reals of a higher precision. Thus the errors using long double (with a * reals of a higher precision. Thus the errors using long double (with a
* 64-bit fraction) are about 2000 times smaller than using double (with a * 64-bit fraction) are about 2000 times smaller than using double (with a
* 53-bit fraction). * 53-bit fraction).
* *
* This algorithm is about 4.5 times slower than the 6th-order Kr&uuml;ge r * This algorithm is about 4.5 times slower than the 6th-order Kr&uuml;ge r
* method, GeographicLib::TransverseMercator, taking about 11 us for a * method, TransverseMercator, taking about 11 us for a combined forward
* combined forward and reverse projection on a 2.66 GHz Intel machine (g and
++, * reverse projection on a 2.66 GHz Intel machine (g++, version 4.3.0, -O
* version 4.3.0, -O3). 3).
* *
* The ellipsoid parameters and the central scale are set in the construc tor. * The ellipsoid parameters and the central scale are set in the construc tor.
* 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 TransverseMercatorExact::Forw ard * specified as the \e lon0 argument of the TransverseMercatorExact::Forw ard
* 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 * taken to be the equator. See the documentation on TransverseMercator
* GeographicLib::TransverseMercator for how to include a false easting, for
* false northing, or a latitude of origin. * 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 TransverseMercatorExact {
private: private:
typedef Math::real real; typedef Math::real real;
skipping to change at line 110 skipping to change at line 108
void sigmainv(real xi, real eta, real& u, real& v) const throw(); void sigmainv(real xi, real eta, real& u, real& v) const throw();
void Scale(real tau, real lam, void Scale(real tau, real lam,
real snu, real cnu, real dnu, real snu, real cnu, real dnu,
real snv, real cnv, real dnv, real snv, real cnv, real dnv,
real& gamma, real& k) const throw(); real& gamma, real& k) const throw();
public: public:
/** /**
* Constructor for a ellipsoid radius \e a (meters), reciprocal flatten * Constructor for a ellipsoid with
ing *
* \e r, and central scale factor \e k0. The transverse Mercator * @param[in] a equatorial radius (meters)
* projection has a branch point singularity at \e lat = 0 and \e lon - * @param[in] r reciprocal flattening.
\e * @param[in] k0 central scale factor.
* lon0 = 90 (1 - \e e) or (for TransverseMercatorExact::UTM) x = 18381 * @param[in] extendp use extended domain.
km, *
* y = 0m. The \e extendp argument governs where the branch cut is pla * The transverse Mercator projection has a branch point singularity at
ced. \e
* With \e extendp = false, the "standard" convention is followed, name * lat = 0 and \e lon - \e lon0 = 90 (1 - \e e) or (for
ly * TransverseMercatorExact::UTM) x = 18381 km, y = 0m. The \e extendp
* the cut is placed along x > 18381 km, y = 0m. Forward can be called * argument governs where the branch cut is placed. With \e extendp =
* with any \e lat and \e lon then produces the transformation shown in * false, the "standard" convention is followed, namely the cut is plac
* Lee, Fig 46. Reverse analytically continues this in the +/- \e x ed
* direction. As a consequence, Reverse may map multiple points to the * along x > 18381 km, y = 0m. Forward can be called with any \e lat a
* same geographic location; for example, for TransverseMercatorExact:: nd
UTM, * \e lon then produces the transformation shown in Lee, Fig 46. Rever
* \e x = 22051449.037349 m, \e y = -7131237.022729 m and \e x = se
* 29735142.378357 m, \e y = 4235043.607933 m both map to \e lat = -2 d * analytically continues this in the +/- \e x direction. As a
eg, * consequence, Reverse may map multiple points to the same geographic
* \e lon = 88 deg. * location; for example, for TransverseMercatorExact::UTM, \e x =
* 22051449.037349 m, \e y = -7131237.022729 m and \e x = 29735142.3783
57
* m, \e y = 4235043.607933 m both map to \e lat = -2 deg, \e lon = 88
deg.
* *
* With \e extendp = true, the branch cut is moved to the lower left * With \e extendp = true, the branch cut is moved to the lower left
* quadrant. The various symmetries of the transverse Mercator project ion * quadrant. The various symmetries of the transverse Mercator project ion
* can be used to explore the projection on any sheet. In this mode th e * can be used to explore the projection on any sheet. In this mode th e
* domains of \e lat, \e lon, \e x, and \e y are restricted to * domains of \e lat, \e lon, \e x, and \e y are restricted to
* - the union of * - the union of
* - \e lat in [0, 90] and \e lon - \e lon0 in [0, 90] * - \e lat in [0, 90] and \e lon - \e lon0 in [0, 90]
* - \e lat in (-90, 0] and \e lon - \e lon0 in [90 (1 - \e e), 90] * - \e lat in (-90, 0] and \e lon - \e lon0 in [90 (1 - \e e), 90]
* - the union of * - the union of
* - \e x/(\e k0 \e a) in [0, inf) and * - \e x/(\e k0 \e a) in [0, inf) and
skipping to change at line 145 skipping to change at line 148
* - \e x/(\e k0 \e a) in [K(1 - \e e^2) - E(1 - \e e^2), inf) and * - \e x/(\e k0 \e a) in [K(1 - \e e^2) - E(1 - \e e^2), inf) and
* \e y/(\e k0 \e a) in (-inf, 0] * \e y/(\e k0 \e a) in (-inf, 0]
* . * .
* See \ref extend for a full discussion of the treatment of the branch * See \ref extend for a full discussion of the treatment of the branch
* cut. * cut.
* *
* The method will work for all ellipsoids used in terrestial geodesy. The * The method will work for all ellipsoids used in terrestial geodesy. The
* method cannot be applied directly to the case of a sphere (\e r = in f) * method cannot be applied directly to the case of a sphere (\e r = in f)
* because some the constants characterizing this method diverge in tha t * because some the constants characterizing this method diverge in tha t
* limit, and in practise, \e r should be smaller than about * limit, and in practise, \e r should be smaller than about
* 1/std::numeric_limits<real>::epsilon(). However, * 1/numeric_limits< real >::%epsilon(). However, TransverseMercator
* GeographicLib::TransverseMercator treats the sphere exactly. An * treats the sphere exactly. An exception is thrown if either axis of
* exception is thrown if \e a, \e r, or \e k0 is not positive. the
* ellipsoid or \e k0 is not positive or if \e r < 1.
********************************************************************** / ********************************************************************** /
TransverseMercatorExact(real a, real r, real k0, bool extendp = false); TransverseMercatorExact(real a, real r, real k0, bool extendp = false);
/** /**
* Convert from latitude \e lat (degrees) and longitude \e lon (degrees * Forward projection, from geographic to transverse Mercator.
) to *
* transverse Mercator easting \e x (meters) and northing \e y (meters) * @param[in] lon0 central meridian of the projection (degrees).
. * @param[in] lat latitude of point (degrees).
* The central meridian of the transformation is \e lon0 (degrees). Al * @param[in] lon longitude of point (degrees).
so * @param[out] x easting of point (meters).
* return the meridian convergence \e gamma (degrees) and the scale \e * @param[out] y northing of point (meters).
k. * @param[out] gamma meridian convergence at point (degrees).
* No false easting or northing is added. \e lat should be in the rang * @param[out] k scale of projection at point.
e *
* No false easting or northing is added. \e lat should be in the range
* [-90, 90]; \e lon and \e lon0 should be in the range [-180, 360]. * [-90, 90]; \e lon and \e lon0 should be in the range [-180, 360].
********************************************************************** / ********************************************************************** /
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();
/** /**
* Convert from transverse Mercator easting \e x (meters) and northing * Reverse projection, from transverse Mercator to geographic.
\e y *
* (meters) to latitude \e lat (degrees) and longitude \e lon (degrees) * @param[in] lon0 central meridian of the projection (degrees).
. * @param[in] x easting of point (meters).
* The central meridian of the transformation is \e lon0 (degrees). Al * @param[in] y northing of point (meters).
so * @param[out] lat latitude of point (degrees).
* return the meridian convergence \e gamma (degrees) and the scale \e * @param[out] lon longitude of point (degrees).
k. * @param[out] gamma meridian convergence at point (degrees).
* @param[out] k scale of projection at point.
*
* 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).
********************************************************************** / ********************************************************************** /
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();
/** /**
* The major radius of the ellipsoid (meters). This is that value of \ * TransverseMercatorExact::Forward without returning the convergence a
e a nd
* used in the constructor. * scale.
**********************************************************************
/
void Forward(real lon0, real lat, real lon,
real& x, real& y) const throw() {
real gamma, k;
Forward(lon0, lat, lon, x, y, gamma, k);
}
/**
* TransverseMercatorExact::Reverse without returning the convergence a
nd
* scale.
**********************************************************************
/
void Reverse(real lon0, real x, real y,
real& lat, real& lon) const throw() {
real gamma, k;
Reverse(lon0, x, y, lat, lon, gamma, k);
}
/** \name Inspector functions
**********************************************************************
/
///@{
/**
* @return \e a the equatorial radius of the ellipsoid (meters). This
is
* the value used in the constructor.
********************************************************************** / ********************************************************************** /
Math::real MajorRadius() const throw() { return _a; } Math::real MajorRadius() const throw() { return _a; }
/** /**
* The inverse flattening of the ellipsoid. This is that value of \e r * @return \e r the inverse flattening of the ellipsoid. This is the
* used in the constructor. * value used in the constructor. A value of 0 is returned for a sph
ere
* (infinite inverse flattening).
********************************************************************** / ********************************************************************** /
Math::real InverseFlattening() const throw() { return _r; } Math::real InverseFlattening() const throw() { return _r; }
/** /**
* The central scale for the projection. This is that value of \e k0 u * @return \e k0 central scale for the projection. This is the value o
sed f \e
* in the constructor and is the scale on the central meridian. * 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; const static TransverseMercatorExact UTM;
}; };
} // namespace GeographicLib } // namespace GeographicLib
 End of changes. 11 change blocks. 
58 lines changed or deleted 99 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 6805 2010-01-28 21:18:14Z karney $" #define GEOGRAPHICLIB_UTMUPS_HPP "$Id: UTMUPS.hpp 6867 2010-09-11 13:04:26Z 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 139 skipping to change at line 139
* The smallest UTM zone number. * The smallest UTM zone number.
******************************************************************** **/ ******************************************************************** **/
MINUTMZONE = 1, MINUTMZONE = 1,
/** /**
* The largest UTM zone number. * The largest UTM zone number.
******************************************************************** **/ ******************************************************************** **/
MAXUTMZONE = 60, MAXUTMZONE = 60,
/** /**
* The largest physical zone number. * The largest physical zone number.
******************************************************************** **/ ******************************************************************** **/
MAXZONE = 60 }; MAXZONE = 60,
};
/** /**
* Return the standard zone for latitude \e lat (degrees) and longitude * The standard zone.
\e *
* lon (degrees); see UTMUPS::STANDARD. This is exact. If the optiona * @param[in] lat latitude (degrees).
l * @param[in] lon longitude (degrees).
* argument \e setzone is given then use that zone if it is non-negativ * @param[in] setzone zone override (optional)
e, *
* otherwise apply the rules given in UTMUPS::zonespec. Throws an erro * This is exact. If the optional argument \e setzone is given then us
r if e
* \e setzone is outsize the range [UTMUPS::MINPSEUDOZONE, UTMUPS::MAXZ * that zone if it is non-negative, otherwise apply the rules given in
ONE] * UTMUPS::zonespec. Throws an error if \e setzone is outsize the rang
* = [-3, 60]. e
* [UTMUPS::MINPSEUDOZONE, UTMUPS::MAXZONE] = [-3, 60].
********************************************************************** / ********************************************************************** /
static int StandardZone(real lat, real lon, int setzone = STANDARD); static int StandardZone(real lat, real lon, int setzone = STANDARD);
/** /**
* Convert geographic coordinates to UTM or UPS coordinate. Given lati * Forward projection, from geographic to UTM/UPS.
tude *
* \e lat (degrees), and longitude \e lon (degrees), return \e zone (ze * @param[in] lat latitude of point (degrees).
ro * @param[in] lon longitude of point (degrees).
* indicates UPS), hemisphere \e northp (false means south, true means * @param[out] zone the UTM zone (zero means UPS).
* north), easting \e x (meters), and northing \e y (meters). The pref * @param[out] northp hemisphere of location (true means northern, fals
ered e
* zone for the result can be specified with \e setzone, see * means southern).
* @param[out] x easting of point (meters).
* @param[out] y northing of point (meters).
* @param[out] gamma meridian convergence at point (degrees).
* @param[out] k scale of projection at point.
* @param[in] setzone zone override.
* @param[in] mgrslimits if true enforce the stricted MGRS limits on th
e
* coordinates (default = false).
*
* The prefered zone for the result can be specified with \e setzone, s
ee
* UTMUPS::StandardZone. Throw error if the resulting easting or north ing * UTMUPS::StandardZone. Throw error if the resulting easting or north ing
* is outside the allowed range (see Reverse), in which case the argume nts * is outside the allowed range (see Reverse), in which case the argume nts
* are unchanged. If \e mgrslimits == true, then use the stricter MGRS * are unchanged. This also returns meridian convergence \e gamma
* limits (see Reverse). This also returns meridian convergence \e gam
ma
* (degrees) and scale \e k. The accuracy of the conversion is about 5 nm. * (degrees) and scale \e k. The accuracy of the conversion is about 5 nm.
********************************************************************** / ********************************************************************** /
static void Forward(real lat, real lon, static void Forward(real lat, real lon,
int& zone, bool& northp, real& x, real& y, int& zone, bool& northp, real& x, real& y,
real& gamma, real& k, real& gamma, real& k,
int setzone = STANDARD, bool mgrslimits = false); int setzone = STANDARD, bool mgrslimits = false);
/** /**
* Convert UTM or UPS coordinate to geographic coordinates . Given zon * Reverse projection, from UTM/UPS to geographic.
e \e *
* zone (\e zone == UTMUPS::UPS, 0, indicates UPS), hemisphere \e * @param[in] zone the UTM zone (zero means UPS).
* northp (false means south, true means north), easting \e x (meters), * @param[in] northp hemisphere of location (true means northern, false
and * means southern).
* northing \e y (meters), return latitude \e lat (degrees) and longitu * @param[in] x easting of point (meters).
de * @param[in] y northing of point (meters).
* \e lon (degrees). Throw error if easting or northing is outside the * @param[out] lat latitude of point (degrees).
* allowed range (see below), in which case the arguments are unchanged * @param[out] lon longitude of point (degrees).
. * @param[out] gamma meridian convergence at point (degrees).
* This also returns meridian convergence \e gamma (degrees) and scale * @param[out] k scale of projection at point.
\e * @param[in] mgrslimits if true enforce the stricted MGRS limits on th
* k. The accuracy of the conversion is about 5nm. e
* coordinates (default = false).
*
* Throw error if easting or northing is outside the allowed range (see
* below), in which case the arguments are unchanged. The accuracy of
the
* conversion is about 5nm.
* *
* UTM eastings are allowed to be in the range [0km, 1000km], northings are * UTM eastings are allowed to be in the range [0km, 1000km], northings are
* allowed to be in in [0km, 9600km] for the northern hemisphere and in * allowed to be in in [0km, 9600km] for the northern hemisphere and in
* [900km, 10000km] for the southern hemisphere. (However UTM northing s * [900km, 10000km] for the southern hemisphere. (However UTM northing s
* can be continued across the equator. So the actual limits on the * can be continued across the equator. So the actual limits on the
* northings are [-9100km, 9600km] for the "northern" hemisphere and * northings are [-9100km, 9600km] for the "northern" hemisphere and
* [900km, 19600km] for the "southern" hemisphere.) * [900km, 19600km] for the "southern" hemisphere.)
* *
* UPS eastings and northings are allowed to be in the range [1200km, * UPS eastings and northings are allowed to be in the range [1200km,
* 2800km] in the northern hemisphere and in [700km, 3100km] in the * 2800km] in the northern hemisphere and in [700km, 3100km] in the
skipping to change at line 201 skipping to change at line 226
* UPS. If \e mgrslimits = true, then all the ranges are shrunk by 100 km * UPS. If \e mgrslimits = true, then all the ranges are shrunk by 100 km
* so that they agree with the stricter MGRS ranges. No checks are * so that they agree with the stricter MGRS ranges. No checks are
* performed besides these (e.g., to limit the distance outside the * performed besides these (e.g., to limit the distance outside the
* standard zone boundaries). * standard zone boundaries).
********************************************************************** / ********************************************************************** /
static void Reverse(int zone, bool northp, real x, real y, static void Reverse(int zone, bool northp, real x, real y,
real& lat, real& lon, real& gamma, real& k, real& lat, real& lon, real& gamma, real& k,
bool mgrslimits = false); bool mgrslimits = false);
/** /**
* Forward without returning convergence and scale. * UTMUPS::Forward without returning convergence and scale.
********************************************************************** / ********************************************************************** /
static void Forward(real lat, real lon, static void Forward(real lat, real lon,
int& zone, bool& northp, real& x, real& y, int& zone, bool& northp, real& x, real& y,
int setzone = STANDARD, bool mgrslimits = false) { int setzone = STANDARD, bool mgrslimits = false) {
real gamma, k; real gamma, k;
Forward(lat, lon, zone, northp, x, y, gamma, k, setzone, mgrslimits); Forward(lat, lon, zone, northp, x, y, gamma, k, setzone, mgrslimits);
} }
/** /**
* Reverse without returning convergence and scale. * UTMUPS::Reverse without returning convergence and scale.
********************************************************************** / ********************************************************************** /
static void Reverse(int zone, bool northp, real x, real y, static void Reverse(int zone, bool northp, real x, real y,
real& lat, real& lon, bool mgrslimits = false) { real& lat, real& lon, bool mgrslimits = false) {
real gamma, k; real gamma, k;
Reverse(zone, northp, x, y, lat, lon, gamma, k, mgrslimits); Reverse(zone, northp, x, y, lat, lon, gamma, k, mgrslimits);
} }
/** /**
* Decode a UTM/UPS zone string, \e zonestr, returning the resulting \e * Decode a UTM/UPS zone string.
* zone and hemisphere thru \e northp (true for northern and false for *
* southern hemispheres). For UTM, \e zonestr has the form of a zone * @param[in] zonestr string represention of zone and hemisphere.
* number in the range [UTMUPS::MINUTMZONE, UTMUPS::MAXUTMZONE] = [1, 6 * @param[out] zone the UTM zone (zero means UPS).
0] * @param[out] northp the hemisphere (true means northern, false
* followed by a hemisphere letter, N or S. For UPS, it consists just * means southern).
of *
* the hemisphere letter. The returned value of \e zone is UTMUPS::UPS * For UTM, \e zonestr has the form of a zone number in the range
= 0 * [UTMUPS::MINUTMZONE, UTMUPS::MAXUTMZONE] = [1, 60] followed by a
* for UPS. Note well that "38S" indicates the southern hemisphere of * hemisphere letter, N or S. For UPS, it consists just of the hemisph
zone ere
* 38 and not latitude band S, [32, 40]. N, 01S, 2N, 38S are legal. 0 * letter. The returned value of \e zone is UTMUPS::UPS = 0 for UPS.
N, Note
* 001S, 61N, 38P are illegal. * well that "38S" indicates the southern hemisphere of zone 38 and not
* latitude band S, [32, 40]. N, 01S, 2N, 38S are legal. 0N, 001S, 61
N,
* 38P are illegal. Throws an error is the zone string is malformed.
********************************************************************** / ********************************************************************** /
static void DecodeZone(const std::string& zonestr, int& zone, bool& nor thp); static void DecodeZone(const std::string& zonestr, int& zone, bool& nor thp);
/** /**
* Encode a UTM/UPS zone string given the \e zone and hemisphere \e * Encode a UTM/UPS zone string.
* northp. \e zone must be in the range [UTMUPS::MINZONE, *
* UTMUPS::MAXZONE] = [0, 60] with \e zone = UTMUPS::UPS, 0, indicating * @param[out] zone the UTM zone (zero means UPS).
* UPS (but the resulting string does not contain "0"). This reverses * @param[out] northp the hemisphere (true means northern, false
* DecodeZone. * means southern).
* @return string represention of zone and hemisphere.
*
* \e zone must be in the range [UTMUPS::MINZONE, UTMUPS::MAXZONE] = [0
,
* 60] with \e zone = UTMUPS::UPS, 0, indicating UPS (but the resulting
* string does not contain "0"). This reverses UTMUPS::DecodeZone.
********************************************************************** / ********************************************************************** /
static std::string EncodeZone(int zone, bool northp); static std::string EncodeZone(int zone, bool northp);
/** /**
* The shift necessary to align N and S halves of a UTM zone * @return shift (meters) necessary to align N and S halves of a UTM zo ne
* (10<sup>7</sup>). * (10<sup>7</sup>).
********************************************************************** / ********************************************************************** /
static Math::real UTMShift() throw(); static Math::real UTMShift() throw();
/** \name Inspector functions
**********************************************************************
/
///@{
/** /**
* The major radius of the ellipsoid (meters). This is the value for t * @return \e a the equatorial radius of the WGS84 ellipsoid (meters).
he *
* WGS84 ellipsoid because the UTM and UPS projections are based on thi * (The WGS84 values is returned because the UTM and UPS projections ar
s e
* ellipsoid. * based on this ellipsoid.)
********************************************************************** / ********************************************************************** /
static Math::real MajorRadius() throw() { return Constants::WGS84_a(); } static Math::real MajorRadius() throw() { return Constants::WGS84_a(); }
/** /**
* The inverse flattening of the ellipsoid. This is the value for the * @return \e r the inverse flattening of the WGS84 ellipsoid.
* WGS84 ellipsoid because the UTM and UPS projections are based on thi *
s * (The WGS84 values is returned because the UTM and UPS projections ar
* ellipsoid. e
* 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. 15 change blocks. 
68 lines changed or deleted 101 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/