Accumulator.hpp   Accumulator.hpp 
/** /**
* \file Accumulator.hpp * \file Accumulator.hpp
* \brief Header for GeographicLib::Accumulator class * \brief Header for GeographicLib::Accumulator class
* *
* Copyright (c) Charles Karney (2010, 2011) <charles@karney.com> and licen sed * Copyright (c) Charles Karney (2010, 2011) <charles@karney.com> and licen sed
* under the MIT/X11 License. For more information, see * under the MIT/X11 License. For more information, see
* http://geographiclib.sourceforge.net/ * http://geographiclib.sourceforge.net/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_ACCUMULATOR_HPP) #if !defined(GEOGRAPHICLIB_ACCUMULATOR_HPP)
#define GEOGRAPHICLIB_ACCUMULATOR_HPP "$Id: d84e4e2bc5c07f7956205da62656e9b #define GEOGRAPHICLIB_ACCUMULATOR_HPP \
bb3765968 $" "$Id: 39476b3fc8d7c5abcc980a7f434bb132d24f53ef $"
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief An accumulator for sums. * \brief An accumulator for sums.
* *
* This allow many numbers of floating point type \e T to be added togeth er * This allow many numbers of floating point type \e T to be added togeth er
* with twice the normal precision. Thus if \e T is double, the effectiv e * with twice the normal precision. Thus if \e T is double, the effectiv e
 End of changes. 1 change blocks. 
2 lines changed or deleted 2 lines changed or added


 AlbersEqualArea.hpp   AlbersEqualArea.hpp 
/** /**
* \file AlbersEqualArea.hpp * \file AlbersEqualArea.hpp
* \brief Header for GeographicLib::AlbersEqualArea class * \brief Header for GeographicLib::AlbersEqualArea class
* *
* Copyright (c) Charles Karney (2010, 2011) <charles@karney.com> and licen sed * Copyright (c) Charles Karney (2010, 2011) <charles@karney.com> and licen sed
* under the MIT/X11 License. For more information, see * under the MIT/X11 License. For more information, see
* http://geographiclib.sourceforge.net/ * http://geographiclib.sourceforge.net/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_ALBERSEQUALAREA_HPP) #if !defined(GEOGRAPHICLIB_ALBERSEQUALAREA_HPP)
#define GEOGRAPHICLIB_ALBERSEQUALAREA_HPP "$Id: eee2b79bf8f65888c9a3b75a6b5 #define GEOGRAPHICLIB_ALBERSEQUALAREA_HPP \
be04ada04e985 $" "$Id: 554a5a83b2feb1f029cc7a7e3c94bdc86ccc36d1 $"
#include <algorithm> #include <algorithm>
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief Albers Equal Area Conic Projection * \brief Albers Equal Area Conic Projection
* *
* Implementation taken from the report, * Implementation taken from the report,
skipping to change at line 114 skipping to change at line 115
// DDatanhee(x,y) = (Datanhee(1,y) - Datanhee(1,x))/(y-x) // DDatanhee(x,y) = (Datanhee(1,y) - Datanhee(1,x))/(y-x)
real DDatanhee(real x, real y) const throw(); real DDatanhee(real x, real y) const throw();
void Init(real sphi1, real cphi1, real sphi2, real cphi2, real k1) thro w(); void Init(real sphi1, real cphi1, real sphi2, real cphi2, real k1) thro w();
real txif(real tphi) const throw(); real txif(real tphi) const throw();
real tphif(real txi) const throw(); real tphif(real txi) const throw();
public: public:
/** /**
* Constructor with a single standard parallel. * Constructor with a single standard parallel.
* *
* @param[in] a equatorial radius of ellipsoid (meters) * @param[in] a equatorial radius of ellipsoid (meters).
* @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphe re. * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphe re.
* Negative \e f gives a prolate ellipsoid. If \e f > 1, set flatten ing * Negative \e f gives a prolate ellipsoid. If \e f > 1, set flatten ing
* to 1/\e f. * to 1/\e f.
* @param[in] stdlat standard parallel (degrees), the circle of tangenc y. * @param[in] stdlat standard parallel (degrees), the circle of tangenc y.
* @param[in] k0 azimuthal scale on the standard parallel. * @param[in] k0 azimuthal scale on the standard parallel.
* *
* An exception is thrown if \e a or \e k0 is not positive or if \e std lat * 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]. * is not in the range [-90, 90].
********************************************************************** / ********************************************************************** /
AlbersEqualArea(real a, real f, real stdlat, real k0); AlbersEqualArea(real a, real f, real stdlat, real k0);
/** /**
* Constructor with two standard parallels. * Constructor with two standard parallels.
* *
* @param[in] a equatorial radius of ellipsoid (meters) * @param[in] a equatorial radius of ellipsoid (meters).
* @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphe re. * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphe re.
* Negative \e f gives a prolate ellipsoid. If \e f > 1, set flatten ing * Negative \e f gives a prolate ellipsoid. If \e f > 1, set flatten ing
* to 1/\e f. * to 1/\e f.
* @param[in] stdlat1 first standard parallel (degrees). * @param[in] stdlat1 first standard parallel (degrees).
* @param[in] stdlat2 second standard parallel (degrees). * @param[in] stdlat2 second standard parallel (degrees).
* @param[in] k1 azimuthal scale on the standard parallels. * @param[in] k1 azimuthal scale on the standard parallels.
* *
* An exception is thrown if \e a or \e k0 is not positive or if \e std lat1 * An exception is thrown if \e a or \e k0 is not positive or if \e std lat1
* or \e stdlat2 is not in the range [-90, 90]. In addition, an except ion * or \e stdlat2 is not in the range [-90, 90]. In addition, an except ion
* is thrown if \e stdlat1 and \e stdlat2 are opposite poles. * is thrown if \e stdlat1 and \e stdlat2 are opposite poles.
********************************************************************** / ********************************************************************** /
AlbersEqualArea(real a, real f, real stdlat1, real stdlat2, real k1); AlbersEqualArea(real a, real f, real stdlat1, real stdlat2, real k1);
/** /**
* Constructor with two standard parallels specified by sines and cosin es. * Constructor with two standard parallels specified by sines and cosin es.
* *
* @param[in] a equatorial radius of ellipsoid (meters) * @param[in] a equatorial radius of ellipsoid (meters).
* @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphe re. * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphe re.
* Negative \e f gives a prolate ellipsoid. If \e f > 1, set flatten ing * Negative \e f gives a prolate ellipsoid. If \e f > 1, set flatten ing
* to 1/\e f. * to 1/\e f.
* @param[in] sinlat1 sine of first standard parallel. * @param[in] sinlat1 sine of first standard parallel.
* @param[in] coslat1 cosine of first standard parallel. * @param[in] coslat1 cosine of first standard parallel.
* @param[in] sinlat2 sine of second standard parallel. * @param[in] sinlat2 sine of second standard parallel.
* @param[in] coslat2 cosine of second standard parallel. * @param[in] coslat2 cosine of second standard parallel.
* @param[in] k1 azimuthal scale on the standard parallels. * @param[in] k1 azimuthal scale on the standard parallels.
* *
* This allows parallels close to the poles to be specified accurately. * This allows parallels close to the poles to be specified accurately.
 End of changes. 4 change blocks. 
5 lines changed or deleted 5 lines changed or added


 AzimuthalEquidistant.hpp   AzimuthalEquidistant.hpp 
/** /**
* \file AzimuthalEquidistant.hpp * \file AzimuthalEquidistant.hpp
* \brief Header for GeographicLib::AzimuthalEquidistant class * \brief Header for GeographicLib::AzimuthalEquidistant class
* *
* Copyright (c) Charles Karney (2009, 2010, 2011) <charles@karney.com> and * Copyright (c) Charles Karney (2009, 2010, 2011) <charles@karney.com> and
* licensed under the MIT/X11 License. For more information, see * licensed under the MIT/X11 License. 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: 48f5a84645ea5925a2e782 #define GEOGRAPHICLIB_AZIMUTHALEQUIDISTANT_HPP \
6adb07a42fe091d09f $" "$Id: 457b3566d14500cdc67b82e6422e0e621f136b3e $"
#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
 End of changes. 1 change blocks. 
2 lines changed or deleted 2 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, 2011) <charles@karney.com> and * Copyright (c) Charles Karney (2009, 2010, 2011) <charles@karney.com> and
* licensed under the MIT/X11 License. For more information, see * licensed under the MIT/X11 License. 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: 2527a0e24486c29c443a79173a18 #define GEOGRAPHICLIB_CASSINISOLDNER_HPP \
def243c091c3 $" "$Id: 5714da5714d03e293c707371c888faefb2bf5fd2 $"
#include <GeographicLib/Geodesic.hpp> #include <GeographicLib/Geodesic.hpp>
#include <GeographicLib/GeodesicLine.hpp> #include <GeographicLib/GeodesicLine.hpp>
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief Cassini-Soldner Projection. * \brief Cassini-Soldner Projection.
* *
skipping to change at line 57 skipping to change at line 58
* 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 Geodesic object (by default * The conversions all take place using a Geodesic object (by default
* Geodesic::WGS84). For more information on geodesics see \ref geodesic . * Geodesic::WGS84). For more information on geodesics see \ref geodesic .
* The determination of (\e lat1, \e lon1) in the forward projection is b y * The determination of (\e lat1, \e lon1) in the forward projection is b y
* solving the inverse geodesic problem for (\e lat, \e lon) and its twin * solving the inverse geodesic problem for (\e lat, \e lon) and its twin
* obtained by reflection in the meridional plane. The scale is found by * obtained by reflection in the meridional plane. The scale is found by
* determining where two neighboring geodesics intersecting the central * determining where two neighboring geodesics intersecting the central
* meridan at \e lat1 and \e lat1 + \e dlat1 intersect and taking the rat io * meridian at \e lat1 and \e lat1 + \e dlat1 intersect and taking the ra tio
* of the reduced lengths for the two geodesics between that point and, * of the reduced lengths for the two geodesics between that point and,
* respectively, (\e lat1, \e lon1) and (\e lat, \e lon). * respectively, (\e lat1, \e lon1) and (\e lat, \e lon).
**********************************************************************/ **********************************************************************/
class GEOGRAPHIC_EXPORT CassiniSoldner { class GEOGRAPHIC_EXPORT CassiniSoldner {
private: private:
typedef Math::real real; typedef Math::real real;
Geodesic _earth; Geodesic _earth;
GeodesicLine _meridian; GeodesicLine _meridian;
real _sbet0, _cbet0; real _sbet0, _cbet0;
skipping to change at line 103 skipping to change at line 104
} }
public: public:
/** /**
* Constructor for CassiniSoldner. * Constructor for CassiniSoldner.
* *
* @param[in] earth the Geodesic object to use for geodesic calculation s. * @param[in] earth the Geodesic object to use for geodesic calculation s.
* By default this uses the WGS84 ellipsoid. * By default this uses the WGS84 ellipsoid.
* *
* This constructor makes an "uninitialized" object. Call Reset to set the * This constructor makes an "uninitialized" object. Call Reset to set the
* central latitude and longuitude, prior to calling Forward and Revers e. * central latitude and longitude, prior to calling Forward and Reverse .
********************************************************************** / ********************************************************************** /
explicit CassiniSoldner(const Geodesic& earth = Geodesic::WGS84) throw( ) explicit CassiniSoldner(const Geodesic& earth = Geodesic::WGS84) throw( )
: _earth(earth) {} : _earth(earth) {}
/** /**
* Constructor for CassiniSoldner specifying a center point. * Constructor for CassiniSoldner specifying a center point.
* *
* @param[in] lat0 latitude of center point of projection (degrees). * @param[in] lat0 latitude of center point of projection (degrees).
* @param[in] lon0 longitude of center point of projection (degrees). * @param[in] lon0 longitude of center point of projection (degrees).
* @param[in] earth the Geodesic object to use for geodesic calculation s. * @param[in] earth the Geodesic object to use for geodesic calculation s.
 End of changes. 3 change blocks. 
4 lines changed or deleted 4 lines changed or added


 CircularEngine.hpp   CircularEngine.hpp 
/** /**
* \file CircularEngine.hpp * \file CircularEngine.hpp
* \brief Header for GeographicLib::CircularEngine class * \brief Header for GeographicLib::CircularEngine class
* *
* Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed un der * Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed un der
* the MIT/X11 License. For more information, see * the MIT/X11 License. For more information, see
* http://geographiclib.sourceforge.net/ * http://geographiclib.sourceforge.net/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_CIRCULARENGINE_HPP) #if !defined(GEOGRAPHICLIB_CIRCULARENGINE_HPP)
#define GEOGRAPHICLIB_CIRCULARENGINE_HPP "$Id: bc8b718e6ab4981abe34da1b8d11 #define GEOGRAPHICLIB_CIRCULARENGINE_HPP \
2dc5f5037465 $" "$Id: 6694dcc59d2fe64a0971fdb0d0d78ff084b0b484 $"
#include <vector> #include <vector>
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
#include <GeographicLib/SphericalEngine.hpp> #include <GeographicLib/SphericalEngine.hpp>
#if defined(_MSC_VER) #if defined(_MSC_VER)
// Squelch warnings about dll vs vector // Squelch warnings about dll vs vector
#pragma warning (push) #pragma warning (push)
#pragma warning (disable: 4251) #pragma warning (disable: 4251)
#endif #endif
skipping to change at line 53 skipping to change at line 54
* CircularEngine stores the coefficients needed to allow the summation o ver * CircularEngine stores the coefficients needed to allow the summation o ver
* order to be performed in 2 or 6 vectors of length \e M + 1 (depending on * order to be performed in 2 or 6 vectors of length \e M + 1 (depending on
* whether gradients are to be calculated). For this reason the construc tor * whether gradients are to be calculated). For this reason the construc tor
* may throw a bad_alloc exception. * may throw a bad_alloc exception.
**********************************************************************/ **********************************************************************/
class GEOGRAPHIC_EXPORT CircularEngine { class GEOGRAPHIC_EXPORT CircularEngine {
private: private:
typedef Math::real real; typedef Math::real real;
enum normalization { enum normalization {
full = SphericalEngine::full, FULL = SphericalEngine::FULL,
schmidt = SphericalEngine::schmidt, SCHMIDT = SphericalEngine::SCHMIDT,
}; };
int _M; int _M;
bool _gradp; bool _gradp;
normalization _norm; unsigned _norm;
real _scale, _a, _r, _u, _t; real _a, _r, _u, _t;
std::vector<real> _wc, _ws, _wrc, _wrs, _wtc, _wts; std::vector<real> _wc, _ws, _wrc, _wrs, _wtc, _wts;
real _q, _uq, _uq2; real _q, _uq, _uq2;
Math::real Value(bool gradp, real cl, real sl, Math::real Value(bool gradp, real cl, real sl,
real& gradx, real& grady, real& gradz) const throw(); real& gradx, real& grady, real& gradz) const throw();
static inline void cossin(real x, real& cosx, real& sinx) { static inline void cossin(real x, real& cosx, real& sinx) {
x = x >= 180 ? x - 360 : (x < -180 ? x + 360 : x); x = x >= 180 ? x - 360 : (x < -180 ? x + 360 : x);
real xi = x * Math::degree<real>(); real xi = x * Math::degree<real>();
cosx = std::abs(x) == 90 ? 0 : cos(xi); cosx = std::abs(x) == 90 ? 0 : cos(xi);
sinx = x == -180 ? 0 : sin(xi); sinx = x == -180 ? 0 : sin(xi);
} }
friend class SphericalEngine; friend class SphericalEngine;
CircularEngine(int M, bool gradp, SphericalEngine::normalization norm, friend class GravityCircle; // Access to cossin
real scale, real a, real r, real u, real t) friend class MagneticCircle; // Access to cossin
CircularEngine(int M, bool gradp, unsigned norm,
real a, real r, real u, real t)
: _M(M) : _M(M)
, _gradp(gradp) , _gradp(gradp)
, _norm(normalization(norm)) , _norm(norm)
, _scale(scale)
, _a(a) , _a(a)
, _r(r) , _r(r)
, _u(u) , _u(u)
, _t(t) , _t(t)
, _wc(std::vector<real>(_M + 1, 0)) , _wc(std::vector<real>(_M + 1, 0))
, _ws(std::vector<real>(_M + 1, 0)) , _ws(std::vector<real>(_M + 1, 0))
, _wrc(std::vector<real>(_gradp ? _M + 1 : 0, 0)) , _wrc(std::vector<real>(_gradp ? _M + 1 : 0, 0))
, _wrs(std::vector<real>(_gradp ? _M + 1 : 0, 0)) , _wrs(std::vector<real>(_gradp ? _M + 1 : 0, 0))
, _wtc(std::vector<real>(_gradp ? _M + 1 : 0, 0)) , _wtc(std::vector<real>(_gradp ? _M + 1 : 0, 0))
, _wts(std::vector<real>(_gradp ? _M + 1 : 0, 0)) , _wts(std::vector<real>(_gradp ? _M + 1 : 0, 0))
skipping to change at line 109 skipping to change at line 111
void SetCoeff(int m, real wc, real ws, void SetCoeff(int m, real wc, real ws,
real wrc, real wrs, real wtc, real wts) { real wrc, real wrs, real wtc, real wts) {
_wc[m] = wc; _ws[m] = ws; _wc[m] = wc; _ws[m] = ws;
if (_gradp) { if (_gradp) {
_wrc[m] = wrc; _wrs[m] = wrs; _wrc[m] = wrc; _wrs[m] = wrs;
_wtc[m] = wtc; _wts[m] = wts; _wtc[m] = wtc; _wts[m] = wts;
} }
} }
public: public:
/**
* A default constructor. CircularEngine::operator()() on the resultin
g
* object returns zero. The resulting object can be assigned to the re
sult
* of SphericalHarmonic::Circle.
**********************************************************************
/
CircularEngine()
: _M(-1)
, _gradp(true)
, _u(0)
, _t(1)
{}
/** /**
* Evaluate the sum for a particular longitude given in terms of its * Evaluate the sum for a particular longitude given in terms of its
* cosine and sine. * cosine and sine.
* *
* @param[in] coslon the cosine of the longitude. * @param[in] coslon the cosine of the longitude.
* @param[in] sinlon the sine of the longitude. * @param[in] sinlon the sine of the longitude.
* @return[in] \e V the value of the sum. * @return \e V the value of the sum.
* *
* The arguments must satisfy \e coslon<sup>2</sup> + \e sinlon<sup>2</ * The arguments must satisfy <i>coslon</i><sup>2</sup> +
sup> * <i>sinlon</i><sup>2</sup> = 1.
* = 1.
********************************************************************** / ********************************************************************** /
Math::real operator()(real coslon, real sinlon) const throw() { Math::real operator()(real coslon, real sinlon) const throw() {
real dummy; real dummy;
return Value(false, coslon, sinlon, dummy, dummy, dummy); return Value(false, coslon, sinlon, dummy, dummy, dummy);
} }
/** /**
* Evaluate the sum for a particular longitude. * Evaluate the sum for a particular longitude.
* *
* @param[in] lon the longitude (degrees). * @param[in] lon the longitude (degrees).
* @return[in] \e V the value of the sum. * @return \e V the value of the sum.
********************************************************************** / ********************************************************************** /
Math::real operator()(real lon) const throw() { Math::real operator()(real lon) const throw() {
real coslon, sinlon; real coslon, sinlon;
cossin(lon, coslon, sinlon); cossin(lon, coslon, sinlon);
return (*this)(coslon, sinlon); return (*this)(coslon, sinlon);
} }
/** /**
* Evaluate the sum and its gradient for a particular longitude given i n * Evaluate the sum and its gradient for a particular longitude given i n
* terms of its cosine and sine. * terms of its cosine and sine.
* *
* @param[in] coslon the cosine of the longitude. * @param[in] coslon the cosine of the longitude.
* @param[in] sinlon the sine of the longitude. * @param[in] sinlon the sine of the longitude.
* @param[out] gradx \e x component of the gradient. * @param[out] gradx \e x component of the gradient.
* @param[out] grady \e y component of the gradient. * @param[out] grady \e y component of the gradient.
* @param[out] gradz \e z component of the gradient. * @param[out] gradz \e z component of the gradient.
* @return[in] \e V the value of the sum. * @return \e V the value of the sum.
* *
* The gradients will only be computed if the CircularEngine object was * The gradients will only be computed if the CircularEngine object was
* created with this capability (e.g., via \e gradp = true in * created with this capability (e.g., via \e gradp = true in
* SphericalHarmonic::Circle). If not, \e gradx, etc., will not be * SphericalHarmonic::Circle). If not, \e gradx, etc., will not be
* touched. The arguments must satisfy \e coslon<sup>2</sup> + \e * touched. The arguments must satisfy <i>coslon</i><sup>2</sup> +
* sinlon<sup>2</sup> = 1. * <i>sinlon</i><sup>2</sup> = 1.
********************************************************************** / ********************************************************************** /
Math::real operator()(real coslon, real sinlon, Math::real operator()(real coslon, real sinlon,
real& gradx, real& grady, real& gradz) const thro w() { real& gradx, real& grady, real& gradz) const thro w() {
return Value(true, coslon, sinlon, gradx, grady, gradz); return Value(true, coslon, sinlon, gradx, grady, gradz);
} }
/** /**
* Evaluate the sum and its gradient for a particular longitude. * Evaluate the sum and its gradient for a particular longitude.
* *
* @param[in] lon the longitude (degrees). * @param[in] lon the longitude (degrees).
* @param[out] gradx \e x component of the gradient. * @param[out] gradx \e x component of the gradient.
* @param[out] grady \e y component of the gradient. * @param[out] grady \e y component of the gradient.
* @param[out] gradz \e z component of the gradient. * @param[out] gradz \e z component of the gradient.
* @return[in] \e V the value of the sum. * @return \e V the value of the sum.
* *
* The gradients will only be computed if the CircularEngine object was * The gradients will only be computed if the CircularEngine object was
* created with this capability (e.g., via \e gradp = true in * created with this capability (e.g., via \e gradp = true in
* SphericalHarmonic::Circle). If not, \e gradx, etc., will not be * SphericalHarmonic::Circle). If not, \e gradx, etc., will not be
* touched. * touched.
********************************************************************** / ********************************************************************** /
Math::real operator()(real lon, Math::real operator()(real lon,
real& gradx, real& grady, real& gradz) const thro w() { real& gradx, real& grady, real& gradz) const thro w() {
real coslon, sinlon; real coslon, sinlon;
cossin(lon, coslon, sinlon); cossin(lon, coslon, sinlon);
 End of changes. 12 change blocks. 
19 lines changed or deleted 35 lines changed or added


 Config.h   Config.h 
#define HAVE_LONG_DOUBLE 1 #define HAVE_LONG_DOUBLE 1
#define GEOGRAPHICLIB_VERSION_STRING "1.15" #define GEOGRAPHICLIB_VERSION_STRING "1.16"
/* # undef WORDS_BIGENDIAN */ /* # undef WORDS_BIGENDIAN */
 End of changes. 1 change blocks. 
1 lines changed or deleted 1 lines changed or added


 Constants.hpp   Constants.hpp 
/** /**
* \file Constants.hpp * \file Constants.hpp
* \brief Header for GeographicLib::Constants class * \brief Header for GeographicLib::Constants class
* *
* Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co m> * Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co m>
* and licensed under the MIT/X11 License. For more information, see * and licensed under the MIT/X11 License. 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: 6201f1cbf92bd8bf866f879fdd95a3cbf #define GEOGRAPHICLIB_CONSTANTS_HPP \
5e3780f $" "$Id: b5139b6cb944bbfb59ea89fef8ab7f60306057c6 $"
#include <GeographicLib/Config.h> #include <GeographicLib/Config.h>
/** /**
* A compile-time assert. Use C++0X static_assert, if available. * A compile-time assert. Use C++0X static_assert, if available.
**********************************************************************/ **********************************************************************/
#if !defined(STATIC_ASSERT) #if !defined(STATIC_ASSERT)
# if defined(__GXX_EXPERIMENTAL_CXX0X__) # if defined(__GXX_EXPERIMENTAL_CXX0X__)
# define STATIC_ASSERT static_assert # define STATIC_ASSERT static_assert
# elif defined(_MSC_VER) && _MSC_VER >= 1600 # elif defined(_MSC_VER) && _MSC_VER >= 1600
skipping to change at line 57 skipping to change at line 58
# define GEOGRAPHIC_EXPORT # define GEOGRAPHIC_EXPORT
#endif #endif
#include <stdexcept> #include <stdexcept>
#include <GeographicLib/Math.hpp> #include <GeographicLib/Math.hpp>
/** /**
* \brief Namespace for %GeographicLib * \brief Namespace for %GeographicLib
* *
* All of %GeographicLib is defined within the GeographicLib namespace. In * All of %GeographicLib is defined within the GeographicLib namespace. In
* addtion all the header files are included via %GeographicLib/filename. * addition all the header files are included via %GeographicLib/filename.
This * This minimizes the likelihood of conflicts with other packages.
* minimizes the likelihood of conflicts with other packages.
**********************************************************************/ **********************************************************************/
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief %Constants needed by %GeographicLib * \brief %Constants needed by %GeographicLib
* *
* Define constants specifying the WGS84 ellipsoid, the UTM and UPS * Define constants specifying the WGS84 ellipsoid, the UTM and UPS
* projections, and various unit conversions. * projections, and various unit conversions.
**********************************************************************/ **********************************************************************/
class GEOGRAPHIC_EXPORT Constants { class GEOGRAPHIC_EXPORT Constants {
skipping to change at line 96 skipping to change at line 97
static inline Math::real arcsecond() throw() static inline Math::real arcsecond() throw()
{ return Math::degree<real>() / 3600; } { return Math::degree<real>() / 3600; }
/** \name Ellipsoid parameters /** \name Ellipsoid parameters
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* @tparam T the type of the returned value. * @tparam T the type of the returned value.
* @return the equatorial radius of WGS84 ellipsoid (6378137 m). * @return the equatorial radius of WGS84 ellipsoid (6378137 m).
********************************************************************** / ********************************************************************** /
template<typename T> template<typename T> static inline T WGS84_a() throw()
static inline T WGS84_a() throw() { return T(6378137) * meter<T>(); } { return T(6378137) * meter<T>(); }
/** /**
* A synonym for WGS84_a<real>(). * A synonym for WGS84_a<real>().
********************************************************************** / ********************************************************************** /
static inline Math::real WGS84_a() throw() { return WGS84_a<real>(); } static inline Math::real WGS84_a() throw() { return WGS84_a<real>(); }
/** /**
* @tparam T the type of the returned value. * @tparam T the type of the returned value.
* @return the flattening of WGS84 ellipsoid (1/298.257223563). * @return the flattening of WGS84 ellipsoid (1/298.257223563).
********************************************************************** / ********************************************************************** /
template<typename T> template<typename T> static inline T WGS84_f() throw()
static inline T WGS84_f() throw() { { return T(1) / ( T(298) + T(257223563) / T(1000000000) ); }
return T(1) / ( T(298) + T(257223563) / T(1000000000) );
}
/** /**
* A synonym for WGS84_f<real>(). * A synonym for WGS84_f<real>().
********************************************************************** / ********************************************************************** /
static inline Math::real WGS84_f() throw() { return WGS84_f<real>(); } static inline Math::real WGS84_f() throw() { return WGS84_f<real>(); }
/** /**
* @tparam T the type of the returned value.
* @return the gravitational constant of the WGS84 ellipsoid, \e GM, in
* m<sup>3</sup> s<sup>-2</sup>.
**********************************************************************
/
template<typename T> static inline T WGS84_GM() throw()
{ return T(3986004) * T(100000000) + T(41800000); }
/**
* @tparam T the type of the returned value.
* @return the angular velocity of the the WGS84 ellipsoid, \e omega, i
n
* rad s<sup>-1</sup>.
**********************************************************************
/
template<typename T> static inline T WGS84_omega() throw()
{ return T(7292115) / (T(1000000) * T(100000)); }
/**
* <b>DEPRECATED</b> * <b>DEPRECATED</b>
* @return the reciprocal flattening of WGS84 ellipsoid. * @return the reciprocal flattening of WGS84 ellipsoid.
********************************************************************** / ********************************************************************** /
template<typename T> template<typename T> static inline T WGS84_r() throw()
static inline T WGS84_r() throw() { return 1/WGS84_f<T>(); } { return 1/WGS84_f<T>(); }
/** /**
* <b>DEPRECATED</b> * <b>DEPRECATED</b>
* A synonym for WGS84_r<real>(). * A synonym for WGS84_r<real>().
********************************************************************** / ********************************************************************** /
static inline Math::real WGS84_r() throw() { return WGS84_r<real>(); } static inline Math::real WGS84_r() throw() { return WGS84_r<real>(); }
/** /**
* @tparam T the type of the returned value. * @tparam T the type of the returned value.
* @return the equatorial radius of GRS80 ellipsoid, \e a, in m.
**********************************************************************
/
template<typename T> static inline T GRS80_a() throw()
{ return T(6378137); }
/**
* @tparam T the type of the returned value.
* @return the gravitational constant of the GRS80 ellipsoid, \e GM, in
* m<sup>3</sup> s<sup>-2</sup>.
**********************************************************************
/
template<typename T> static inline T GRS80_GM() throw()
{ return T(3986005) * T(100000000); }
/**
* @tparam T the type of the returned value.
* @return the angular velocity of the the GRS80 ellipsoid, \e omega, i
n
* rad s<sup>-1</sup>.
*
* This is about 2*pi*366.25 / (365.25*24*3600) rad s<sup>-1</sup>. 36
5.25
* is the number of days in a Julian year and 365.35/366.25 converts fr
om
* solar days to sidereal days. Using the number of days in a Gregoria
n
* year (365.2425) results in a worse approximation (because the Gregor
ian
* year includes the precession of the earth's axis).
**********************************************************************
/
template<typename T> static inline T GRS80_omega() throw()
{ return T(7292115) / (T(1000000) * T(100000)); }
/**
* @tparam T the type of the returned value.
* @return the dynamical form factor of the GRS80 ellipsoid,
* <i>J</i><sub>2</sub>.
**********************************************************************
/
template<typename T> static inline T GRS80_J2() throw()
{ return T(108263) / T(100000000); }
/**
* @tparam T the type of the returned value.
* @return the central scale factor for UTM (0.9996). * @return the central scale factor for UTM (0.9996).
********************************************************************** / ********************************************************************** /
template<typename T> template<typename T> static inline T UTM_k0() throw()
static inline T UTM_k0() throw() {return T(9996) / T(10000); } {return T(9996) / T(10000); }
/** /**
* A synonym for UTM_k0<real>(). * A synonym for UTM_k0<real>().
********************************************************************** / ********************************************************************** /
static inline Math::real UTM_k0() throw() { return UTM_k0<real>(); } static inline Math::real UTM_k0() throw() { return UTM_k0<real>(); }
/** /**
* @tparam T the type of the returned value. * @tparam T the type of the returned value.
* @return the central scale factor for UPS (0.994). * @return the central scale factor for UPS (0.994).
********************************************************************** / ********************************************************************** /
template<typename T> template<typename T> static inline T UPS_k0() throw()
static inline T UPS_k0() throw() { return T(994) / T(1000); } { return T(994) / T(1000); }
/** /**
* A synonym for UPS_k0<real>(). * A synonym for UPS_k0<real>().
********************************************************************** / ********************************************************************** /
static inline Math::real UPS_k0() throw() { return UPS_k0<real>(); } static inline Math::real UPS_k0() throw() { return UPS_k0<real>(); }
///@} ///@}
/** \name SI units /** \name SI units
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* @tparam T the type of the returned value. * @tparam T the type of the returned value.
* @return the number of meters in a meter. * @return the number of meters in a meter.
* *
* This is unity, but this lets the internal system of units be changed if * This is unity, but this lets the internal system of units be changed if
* necessary. * necessary.
********************************************************************** / ********************************************************************** /
template<typename T> template<typename T> static inline T meter() throw() { return T(1); }
static inline T meter() throw() { return T(1); }
/** /**
* A synonym for meter<real>(). * A synonym for meter<real>().
********************************************************************** / ********************************************************************** /
static inline Math::real meter() throw() { return meter<real>(); } static inline Math::real meter() throw() { return meter<real>(); }
/** /**
* @return the number of meters in a kilometer. * @return the number of meters in a kilometer.
********************************************************************** / ********************************************************************** /
static inline Math::real kilometer() throw() static inline Math::real kilometer() throw()
{ return 1000 * meter<real>(); } { return 1000 * meter<real>(); }
/** /**
skipping to change at line 182 skipping to change at line 227
static inline Math::real nauticalmile() throw() static inline Math::real nauticalmile() throw()
{ return 1852 * meter<real>(); } { return 1852 * meter<real>(); }
/** /**
* @tparam T the type of the returned value. * @tparam T the type of the returned value.
* @return the number of square meters in a square meter. * @return the number of square meters in a square meter.
* *
* This is unity, but this lets the internal system of units be changed if * This is unity, but this lets the internal system of units be changed if
* necessary. * necessary.
********************************************************************** / ********************************************************************** /
template<typename T> template<typename T> static inline T square_meter() throw()
static inline T square_meter() throw()
{ return meter<real>() * meter<real>(); } { return meter<real>() * meter<real>(); }
/** /**
* A synonym for square_meter<real>(). * A synonym for square_meter<real>().
********************************************************************** / ********************************************************************** /
static inline Math::real square_meter() throw() static inline Math::real square_meter() throw()
{ return square_meter<real>(); } { return square_meter<real>(); }
/** /**
* @return the number of square meters in a hectare. * @return the number of square meters in a hectare.
********************************************************************** / ********************************************************************** /
static inline Math::real hectare() throw() static inline Math::real hectare() throw()
 End of changes. 11 change blocks. 
21 lines changed or deleted 75 lines changed or added


 DMS.hpp   DMS.hpp 
/** /**
* \file DMS.hpp * \file DMS.hpp
* \brief Header for GeographicLib::DMS class * \brief Header for GeographicLib::DMS class
* *
* Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co m> * Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co m>
* and licensed under the MIT/X11 License. For more information, see * and licensed under the MIT/X11 License. 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: bc4a072cb427ec5f1648616fdceea58c6fccae7 a $" #define GEOGRAPHICLIB_DMS_HPP "$Id: 83292a5889c09c7f444a587d501a4145bcdf993 c $"
#include <sstream> #include <sstream>
#include <iomanip> #include <iomanip>
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
#include <GeographicLib/Utility.hpp> #include <GeographicLib/Utility.hpp>
#if defined(_MSC_VER) #if defined(_MSC_VER)
// Squelch warnings about dll vs string // Squelch warnings about dll vs string
#pragma warning (push) #pragma warning (push)
#pragma warning (disable: 4251) #pragma warning (disable: 4251)
skipping to change at line 114 skipping to change at line 114
* Trailing unit is arc seconds. * Trailing unit is arc seconds.
* @hideinitializer * @hideinitializer
******************************************************************** **/ ******************************************************************** **/
SECOND = 2, SECOND = 2,
}; };
/** /**
* Convert a string in DMS to an angle. * Convert a string in DMS to an angle.
* *
* @param[in] dms string input. * @param[in] dms string input.
* @param[out] ind a DMS::flag value signalling the presence of a * @param[out] ind a DMS::flag value signaling the presence of a
* hemisphere indicator. * hemisphere indicator.
* @return angle (degrees). * @return angle (degrees).
* *
* Degrees, minutes, and seconds are indicated by the characters d, ' * Degrees, minutes, and seconds are indicated by the characters d, '
* (single quote), &quot; (double quote), and these components may only be * (single quote), &quot; (double quote), and these components may only be
* given in this order. Any (but not all) components may be omitted an d * given in this order. Any (but not all) components may be omitted an d
* other symbols (e.g., the <sup>o</sup> symbol for degrees and the uni code * other symbols (e.g., the <sup>o</sup> symbol for degrees and the uni code
* prime and double prime symbols for minutes and seconds) may be * prime and double prime symbols for minutes and seconds) may be
* substituted. The last component indicator may be omitted and is ass umed * substituted. The last component indicator may be omitted and is ass umed
* to be the next smallest unit (thus 33d10 is interpreted as 33d10'). The * to be the next smallest unit (thus 33d10 is interpreted as 33d10'). The
skipping to change at line 311 skipping to change at line 311
********************************************************************** / ********************************************************************** /
static std::string Encode(real angle, unsigned prec, flag ind = NONE) { static std::string Encode(real angle, unsigned prec, flag ind = NONE) {
return ind == NUMBER ? Utility::str<real>(angle, int(prec)) : return ind == NUMBER ? Utility::str<real>(angle, int(prec)) :
Encode(angle, 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 into degrees and mintues * Split angle into degrees and minutes
* *
* @param[in] ang angle (degrees) * @param[in] ang angle (degrees)
* @param[out] d degrees (an integer returned as a real) * @param[out] d degrees (an integer returned as a real)
* @param[out] m arc minutes. * @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 into degrees and mintues and seconds. * Split angle into degrees and minutes and seconds.
* *
* @param[in] ang angle (degrees) * @param[in] ang angle (degrees)
* @param[out] d degrees (an integer returned as a real) * @param[out] d degrees (an integer returned as a real)
* @param[out] m arc minutes (an integer returned as a real) * @param[out] m arc minutes (an integer returned as a real)
* @param[out] s arc seconds. * @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);
} }
 End of changes. 4 change blocks. 
4 lines changed or deleted 4 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, 2011) <charles@karney.com> * Copyright (c) Charles Karney (2008, 2009, 2011) <charles@karney.com>
* and licensed under the MIT/X11 License. For more information, see * and licensed under the MIT/X11 License. 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: 1ed2e0020530f29db210a5816b #define GEOGRAPHICLIB_ELLIPTICFUNCTION_HPP \
d9990ce5cd05f5 $" "$Id: fcee9dee564eaf87b75e742f16d1954091a77dbf $"
#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<sup>2</sup> and 1 - \e e<sup>2</sup> are needed. This class taken t * <i>e</i><sup>2</sup> and 1 - <i>e</i><sup>2</sup> are needed. This cl
he ass
* parameter as a constructor parameters and caches the values of the * taken the parameter as a constructor parameters and caches the values
* required complete integrals. A method is provided for Jacobi elliptic of
* the required complete integrals. A method is provided for Jacobi elli
ptic
* functions and for the incomplete elliptic integral of the second kind in * functions and for the incomplete elliptic integral of the second kind in
* terms of the amplitude. * terms of the amplitude.
* *
* The computation of the elliptic integrals uses the algorithms given in * The computation of the elliptic integrals uses the algorithms given in
* - B. C. Carlson, * - B. C. Carlson,
* <a href="http://dx.doi.org/10.1007/BF02198293"> Computation of ellip tic * <a href="http://dx.doi.org/10.1007/BF02198293"> Computation of ellip tic
* integrals</a>, Numerical Algorithms 10, 13&ndash;26 (1995). * integrals</a>, Numerical Algorithms 10, 13&ndash;26 (1995).
* . * .
* The computation of the Jacobi elliptic functions uses the algorithm gi ven * The computation of the Jacobi elliptic functions uses the algorithm gi ven
* in * in
skipping to change at line 122 skipping to change at line 123
/** /**
* The incomplete integral of the second kind in terms of Jacobi ellipt ic * The incomplete integral of the second kind in terms of Jacobi ellipt ic
* functions * functions
* *
* @param[in] sn * @param[in] sn
* @param[in] cn * @param[in] cn
* @param[in] dn * @param[in] dn
* @return int dn(\e w)<sup>2</sup> \e dw (A+S 17.2.10). * @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 * Instead of specifying the amplitude \e phi, we provide \e sn = sin(\ e
* phi), \e cn = cos(\e phi), \e dn = sqrt(1 - \e m sin<sup>2</sup>(\e * phi), \e cn = cos(\e phi), \e dn = sqrt(1 - \e m sin<sup>2</sup>(\e
* phi)). * phi)).
********************************************************************** / ********************************************************************** /
Math::real E(real sn, real cn, real dn) const throw(); Math::real E(real sn, real cn, real dn) const throw();
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif // GEOGRAPHICLIB_ELLIPTICFUNCTION_HPP #endif // GEOGRAPHICLIB_ELLIPTICFUNCTION_HPP
 End of changes. 3 change blocks. 
8 lines changed or deleted 10 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, 2011) <charles@karney.co m> * Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co m>
* and licensed under the MIT/X11 License. For more information, see * and licensed under the MIT/X11 License. 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: 0498603ef805f63465c670a5ded816492 #define GEOGRAPHICLIB_GEOCOORDS_HPP \
d3bc631 $" "$Id: 8f389fbd5bd0231892039ef8ffae9ded5736f013 $"
#include <GeographicLib/UTMUPS.hpp> #include <GeographicLib/UTMUPS.hpp>
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief Conversion between geographic coordinates * \brief Conversion between geographic coordinates
* *
* This class stores a geographic position which may be set via the * This class stores a geographic position which may be set via the
* constructors or Reset via * constructors or Reset via
* - latitude and longitude * - latitude and longitude
* - UTM or UPS coordinates * - UTM or UPS coordinates
* - a string represention of these or an MGRS coordinate string * - a string representation of these or an MGRS coordinate string
* *
* The state consists of the latitude and longitude and the supplied UTM or * The state consists of the latitude and longitude and the supplied UTM or
* UPS coordinates (possibly derived from the MGRS coordinates). If lati tude * UPS coordinates (possibly derived from the MGRS coordinates). If lati tude
* and longitude were given then the UTM/UPS coordinates follows the stan dard * and longitude were given then the UTM/UPS coordinates follows the stan dard
* conventions. * conventions.
* *
* The mutable state consists of the UTM or UPS coordinates for a alterna te * The mutable state consists of the UTM or UPS coordinates for a alterna te
* zone. A method SetAltZone is provided to set the alternate UPS/UTM zo ne. * zone. A method SetAltZone is provided to set the alternate UPS/UTM zo ne.
* *
* Methods are provided to return the geographic coordinates, the input U TM * Methods are provided to return the geographic coordinates, the input U TM
skipping to change at line 68 skipping to change at line 69
} }
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 /** \name Initializing the GeoCoords object
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* The default contructor is equivalent to \e latitude = 90<sup>o</sup> * The default constructor is equivalent to \e latitude = 90<sup>o</sup
, \e >,
* longitude = 0<sup>o</sup>. * \e 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(); }
skipping to change at line 123 skipping to change at line 124
* N, S, E, W hemisphere designator is used on one or both coordinates. If * N, S, E, W hemisphere designator is used on one or both coordinates. If
* \e swaplatlong = true (default is false), then longitude precedes * \e swaplatlong = true (default is false), then longitude precedes
* latitude in the absence of a hemisphere designator. Thus (with \e * latitude in the absence of a hemisphere designator. Thus (with \e
* swaplatlong = false) * swaplatlong = false)
* - 40 -75 * - 40 -75
* - N40 W75 * - N40 W75
* - -75 N40 * - -75 N40
* - 75W 40N * - 75W 40N
* - E-75 -40S * - E-75 -40S
* . * .
* are all the same position. The coodinates may be given in decimal * are all the same position. The coordinates may be given in decimal
* degrees, degrees and decimal minutes, degrees, minutes, seconds, etc . * degrees, degrees and decimal minutes, degrees, minutes, seconds, etc .
* Use d, ', and " to make off the degrees, minutes and seconds. Thus * Use d, ', and " to make off the degrees, minutes and seconds. Thus
* - 40d30'30" * - 40d30'30"
* - 40d30'30 * - 40d30'30
* - 40d30.5' * - 40d30.5'
* - 40d30.5 * - 40d30.5
* - 40.508333333 * - 40.508333333
* . * .
* all specify the same angle. The leading sign applies to all compone nts * all specify the same angle. The leading sign applies to all compone nts
* so -1d30 is -(1+30/60) = -1.5. Latitudes must be in the range [-90, 90] * so -1d30 is -(1+30/60) = -1.5. Latitudes must be in the range [-90, 90]
skipping to change at line 318 skipping to change at line 319
* @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 alternate 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 alternate zone.
********************************************************************** / ********************************************************************** /
Math::real AltScale() const throw() { return _alt_k; } Math::real AltScale() const throw() { return _alt_k; }
///@} ///@}
/** \name String representations of the GeoCoords object /** \name String representations of the GeoCoords object
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* String representation with latitude and longitude as signed decimal * String representation with latitude and longitude as signed decimal
* degrees. * degrees.
skipping to change at line 376 skipping to change at line 377
/** /**
* MGRS string. * MGRS string.
* *
* @param[in] prec precision (relative to about 1m). * @param[in] prec precision (relative to about 1m).
* @return MGRS string. * @return MGRS string.
* *
* This gives the coordinates of the enclosing grid square with size gi ven * This gives the coordinates of the enclosing grid square with size gi ven
* by the precision. Thus 38N 444180 3684790 converted to a MGRS * by the precision. Thus 38N 444180 3684790 converted to a MGRS
* coordinate at precision -2 (100m) is 38SMB441847 and not 38SMB442848 . * coordinate at precision -2 (100m) is 38SMB441847 and not 38SMB442848 .
* \e prec specifies the precision of the MSGRS string as follows: * \e prec specifies the precision of the MGRS 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;
 End of changes. 7 change blocks. 
10 lines changed or deleted 10 lines changed or added


 Geocentric.hpp   Geocentric.hpp 
/** /**
* \file Geocentric.hpp * \file Geocentric.hpp
* \brief Header for GeographicLib::Geocentric class * \brief Header for GeographicLib::Geocentric class
* *
* Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co m> * Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co m>
* and licensed under the MIT/X11 License. For more information, see * and licensed under the MIT/X11 License. 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: 5e13a145df558737b41cb38be92d4e01 #define GEOGRAPHICLIB_GEOCENTRIC_HPP \
03120e6d $" "$Id: 869e8cba7a96857eef318b2b02e4135d88bfbc75 $"
#include <vector> #include <vector>
#include <algorithm> #include <algorithm>
#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)
* to geocentric coordinates (\e x, \e y, \e z). The origin of geocentri * to geocentric coordinates (\e X, \e Y, \e Z). The origin of geocentri
c c
* coordinates is at the center of the earth. The \e z axis goes thru th * coordinates is at the center of the earth. The \e Z axis goes thru th
e e
* north pole, \e lat = 90<sup>o</sup>. The \e x axis goes thru \e lat = * north pole, \e lat = 90<sup>o</sup>. The \e X axis goes thru \e lat =
0, 0,
* \e lon = 0. %Geocentric coordinates are also known as earth centered, * \e lon = 0. %Geocentric coordinates are also known as earth centered,
* earth fixed (ECEF) coordinates. * earth fixed (ECEF) coordinates.
* *
* The conversion from geographic to geocentric coordinates is * The conversion from geographic to geocentric coordinates is
* straightforward. For the reverse transformation we use * straightforward. For the reverse transformation we use
* - H. Vermeille, * - H. Vermeille,
* <a href="http://dx.doi.org/10.1007/s00190-002-0273-6"> Direct * <a href="http://dx.doi.org/10.1007/s00190-002-0273-6"> Direct
* transformation from geocentric coordinates to geodetic coordinates</ a>, * transformation from geocentric coordinates to geodetic coordinates</ a>,
* J. Geodesy 76, 451&ndash;454 (2002). * J. Geodesy 76, 451&ndash;454 (2002).
* . * .
skipping to change at line 61 skipping to change at line 62
* outside the ellipsoid), the error is bounded by 7 nm for the WGS84 * outside the ellipsoid), the error is bounded by 7 nm for the WGS84
* ellipsoid. See \ref geocentric for further information on the errors. * ellipsoid. See \ref geocentric for further information on the errors.
**********************************************************************/ **********************************************************************/
class GEOGRAPHIC_EXPORT Geocentric { class GEOGRAPHIC_EXPORT Geocentric {
private: private:
typedef Math::real real; typedef Math::real real;
friend class LocalCartesian; friend class LocalCartesian;
friend class MagneticCircle; // MagneticCircle uses Rotation friend class MagneticCircle; // MagneticCircle uses Rotation
friend class MagneticModel; // MagneticModel uses IntForward friend class MagneticModel; // MagneticModel uses IntForward
friend class GravityCircle; // GravityCircle uses Rotation
friend class GravityModel; // GravityModel uses IntForward
friend class NormalGravity; // NormalGravity uses IntForward
friend class SphericalHarmonic;
friend class SphericalHarmonic1;
friend class SphericalHarmonic2;
static const size_t dim_ = 3; static const size_t dim_ = 3;
static const size_t dim2_ = dim_ * dim_; static const size_t dim2_ = dim_ * dim_;
real _a, _f, _e2, _e2m, _e2a, _e4a, _maxrad; real _a, _f, _e2, _e2m, _e2a, _e4a, _maxrad;
static void Rotation(real sphi, real cphi, real slam, real clam, static void Rotation(real sphi, real cphi, real slam, real clam,
real M[dim2_]) throw(); real M[dim2_]) throw();
void IntForward(real lat, real lon, real h, real& x, real& y, real& z, static void Rotate(real M[dim2_], real x, real y, real z,
real& X, real& Y, real& Z) throw() {
// Perform [X,Y,Z]^t = M.[x,y,z]^t
// (typically local cartesian to geocentric)
X = M[0] * x + M[1] * y + M[2] * z;
Y = M[3] * x + M[4] * y + M[5] * z;
Z = M[6] * x + M[7] * y + M[8] * z;
}
static void Unrotate(real M[dim2_], real X, real Y, real Z,
real& x, real& y, real& z) throw() {
// Perform [x,y,z]^t = M^t.[X,Y,Z]^t
// (typically geocentric to local cartesian)
x = M[0] * X + M[3] * Y + M[6] * Z;
y = M[1] * X + M[4] * Y + M[7] * Z;
z = M[2] * X + M[5] * Y + M[8] * Z;
}
void IntForward(real lat, real lon, real h, real& X, real& Y, real& Z,
real M[dim2_]) const throw(); real M[dim2_]) const throw();
void IntReverse(real x, real y, real z, real& lat, real& lon, real& h, void IntReverse(real X, real Y, real Z, real& lat, real& lon, real& h,
real M[dim2_]) const throw(); real M[dim2_]) const throw();
public: public:
/** /**
* Constructor for a ellipsoid with * Constructor for a ellipsoid with
* *
* @param[in] a equatorial radius (meters) * @param[in] a equatorial radius (meters).
* @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphe re. * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphe re.
* Negative \e f gives a prolate ellipsoid. If \e f > 1, set flatten ing * Negative \e f gives a prolate ellipsoid. If \e f > 1, set flatten ing
* to 1/\e f. * to 1/\e f.
* *
* An exception is thrown if either of the axes of the ellipsoid is * An exception is thrown if either of the axes of the ellipsoid is
* non-positive. * non-positive.
********************************************************************** / ********************************************************************** /
Geocentric(real a, real f); Geocentric(real a, real f);
/** /**
* A default constructor (for use by NormalGravity).
**********************************************************************
/
Geocentric() : _a(-1) {}
/**
* Convert from geodetic to geocentric coordinates. * Convert from geodetic to geocentric coordinates.
* *
* @param[in] lat latitude of point (degrees). * @param[in] lat latitude of point (degrees).
* @param[in] lon longitude of point (degrees). * @param[in] lon longitude of point (degrees).
* @param[in] h height of point above the ellipsoid (meters). * @param[in] h height of point above the ellipsoid (meters).
* @param[out] x geocentric coordinate (meters). * @param[out] X geocentric coordinate (meters).
* @param[out] y geocentric coordinate (meters). * @param[out] Y geocentric coordinate (meters).
* @param[out] z 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 * \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() {
IntForward(lat, lon, h, x, y, z, NULL); if (Init())
IntForward(lat, lon, h, X, Y, Z, NULL);
} }
/** /**
* Convert from geodetic to geocentric coordinates and return rotation * Convert from geodetic to geocentric coordinates and return rotation
* matrix. * matrix.
* *
* @param[in] lat latitude of point (degrees). * @param[in] lat latitude of point (degrees).
* @param[in] lon longitude of point (degrees). * @param[in] lon longitude of point (degrees).
* @param[in] h height of point above the ellipsoid (meters). * @param[in] h height of point above the ellipsoid (meters).
* @param[out] x geocentric coordinate (meters). * @param[out] X geocentric coordinate (meters).
* @param[out] y geocentric coordinate (meters). * @param[out] Y geocentric coordinate (meters).
* @param[out] z geocentric coordinate (meters). * @param[out] Z geocentric coordinate (meters).
* @param[out] M if the length of the vector is 9, fill with the rotati on * @param[out] M if the length of the vector is 9, fill with the rotati on
* matrix in row-major order. * matrix in row-major order.
* *
* Let \e v be a unit vector located at (\e lat, \e lon, \e h). We can * Let \e v be a unit vector located at (\e lat, \e lon, \e h). We can
* express \e v as \e column vectors in one of two ways * express \e v as \e column vectors in one of two ways
* - in east, north, up coordinates (where the components are relative to a * - in east, north, up coordinates (where the components are relative to a
* local coordinate system at (\e lat, \e lon, \e h)); call this * local coordinate system at (\e lat, \e lon, \e h)); call this
* representation \e v1. * representation \e v1.
* - in geocentric \e x, \e y, \e z coordinates; call this representati on * - in geocentric \e X, \e Y, \e Z coordinates; call this representati on
* \e v0. * \e v0.
* . * .
* Then we have \e v0 = \e M . \e v1. * Then we have \e v0 = \e M . \e v1.
********************************************************************** / ********************************************************************** /
void Forward(real lat, real lon, real h, real& x, real& y, real& z, void Forward(real lat, real lon, real h, real& X, real& Y, real& Z,
std::vector<real>& M) std::vector<real>& M)
const throw() { const throw() {
if (!Init())
return;
if (M.end() == M.begin() + dim2_) { if (M.end() == M.begin() + dim2_) {
real t[dim2_]; real t[dim2_];
IntForward(lat, lon, h, x, y, z, t); IntForward(lat, lon, h, X, Y, Z, t);
copy(t, t + dim2_, M.begin()); copy(t, t + dim2_, M.begin());
} else } else
IntForward(lat, lon, h, x, y, z, NULL); IntForward(lat, lon, h, X, Y, Z, NULL);
} }
/** /**
* Convert from geocentric to geodetic to coordinates. * Convert from geocentric to geodetic to coordinates.
* *
* @param[in] x geocentric coordinate (meters). * @param[in] X geocentric coordinate (meters).
* @param[in] y geocentric coordinate (meters). * @param[in] Y geocentric coordinate (meters).
* @param[in] z geocentric coordinate (meters). * @param[in] Z geocentric coordinate (meters).
* @param[out] lat latitude of point (degrees). * @param[out] lat latitude of point (degrees).
* @param[out] lon longitude of point (degrees). * @param[out] lon longitude of point (degrees).
* @param[out] h height of point above the ellipsoid (meters). * @param[out] h height of point above the ellipsoid (meters).
* *
* In general there are multiple solutions and the result which maximiz es * 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 * latitudes (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 - * returned. The value of \e h returned satisfies \e h >= - \e a (1 -
\e * <i>e</i><sup>2</sup>) / sqrt(1 - <i>e</i><sup>2</sup> sin<sup>2</sup
* e<sup>2</sup>) / sqrt(1 - \e e<sup>2</sup> sin<sup>2</sup>\e lat). >\e
The * 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() {
IntReverse(x, y, z, lat, lon, h, NULL); if (Init())
IntReverse(X, Y, Z, lat, lon, h, NULL);
} }
/** /**
* Convert from geocentric to geodetic to coordinates. * Convert from geocentric to geodetic to coordinates.
* *
* @param[in] x geocentric coordinate (meters). * @param[in] X geocentric coordinate (meters).
* @param[in] y geocentric coordinate (meters). * @param[in] Y geocentric coordinate (meters).
* @param[in] z geocentric coordinate (meters). * @param[in] Z geocentric coordinate (meters).
* @param[out] lat latitude of point (degrees). * @param[out] lat latitude of point (degrees).
* @param[out] lon longitude of point (degrees). * @param[out] lon longitude of point (degrees).
* @param[out] h height of point above the ellipsoid (meters). * @param[out] h height of point above the ellipsoid (meters).
* @param[out] M if the length of the vector is 9, fill with the rotati on * @param[out] M if the length of the vector is 9, fill with the rotati on
* matrix in row-major order. * matrix in row-major order.
* *
* Let \e v be a unit vector located at (\e lat, \e lon, \e h). We can * Let \e v be a unit vector located at (\e lat, \e lon, \e h). We can
* express \e v as \e column vectors in one of two ways * express \e v as \e column vectors in one of two ways
* - in east, north, up coordinates (where the components are relative to a * - in east, north, up coordinates (where the components are relative to a
* local coordinate system at (\e lat, \e lon, \e h)); call this * local coordinate system at (\e lat, \e lon, \e h)); call this
* representation \e v1. * representation \e v1.
* - in geocentric \e x, \e y, \e z coordinates; call this representati on * - in geocentric \e X, \e Y, \e Z coordinates; call this representati on
* \e v0. * \e v0.
* . * .
* Then we have \e v1 = \e M^T . \e v0, where \e M^T is the transpose o f \e * Then we have \e v1 = \e M^T . \e v0, where \e M^T is the transpose o f \e
* M. * M.
********************************************************************** / ********************************************************************** /
void Reverse(real x, real y, real z, real& lat, real& lon, real& h, void Reverse(real X, real Y, real Z, real& lat, real& lon, real& h,
std::vector<real>& M) std::vector<real>& M)
const throw() { const throw() {
if (!Init())
return;
if (M.end() == M.begin() + dim2_) { if (M.end() == M.begin() + dim2_) {
real t[dim2_]; real t[dim2_];
IntReverse(x, y, z, lat, lon, h, t); IntReverse(X, Y, Z, lat, lon, h, t);
copy(t, t + dim2_, M.begin()); copy(t, t + dim2_, M.begin());
} else } else
IntReverse(x, y, z, lat, lon, h, NULL); IntReverse(X, Y, Z, lat, lon, h, NULL);
} }
/** \name Inspector functions /** \name Inspector functions
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* @return true if the object has been initialized.
**********************************************************************
/
bool Init() const throw() { return _a > 0; }
/**
* @return \e a the equatorial radius of the ellipsoid (meters). This is * @return \e a the equatorial radius of the ellipsoid (meters). This is
* the value used in the constructor. * the value used in the constructor.
********************************************************************** / ********************************************************************** /
Math::real MajorRadius() const throw() { return _a; } Math::real MajorRadius() const throw()
{ return Init() ? _a : Math::NaN<real>(); }
/** /**
* @return \e f the flattening of the ellipsoid. This is the * @return \e f the flattening of the ellipsoid. This is the
* value used in the constructor. * value used in the constructor.
********************************************************************** / ********************************************************************** /
Math::real Flattening() const throw() { return _f; } Math::real Flattening() const throw()
{ return Init() ? _f : Math::NaN<real>(); }
/** /**
* <b>DEPRECATED</b> * <b>DEPRECATED</b>
* @return \e r the inverse flattening of the ellipsoid. * @return \e r the inverse flattening of the ellipsoid.
********************************************************************** / ********************************************************************** /
Math::real InverseFlattening() const throw() { return 1/_f; } Math::real InverseFlattening() const throw()
{ return Init() ? 1/_f : Math::NaN<real>(); }
///@} ///@}
/** /**
* A global instantiation of Geocentric with the parameters for the WGS 84 * A global instantiation of Geocentric with the parameters for the WGS 84
* ellipsoid. * ellipsoid.
********************************************************************** / ********************************************************************** /
static const Geocentric WGS84; static const Geocentric WGS84;
}; };
} // namespace GeographicLib } // namespace GeographicLib
 End of changes. 31 change blocks. 
45 lines changed or deleted 86 lines changed or added


 Geodesic.hpp   Geodesic.hpp 
/** /**
* \file Geodesic.hpp * \file Geodesic.hpp
* \brief Header for GeographicLib::Geodesic class * \brief Header for GeographicLib::Geodesic class
* *
* Copyright (c) Charles Karney (2009, 2010, 2011) <charles@karney.com> * Copyright (c) Charles Karney (2009, 2010, 2011) <charles@karney.com>
* and licensed under the MIT/X11 License. For more information, see * and licensed under the MIT/X11 License. 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: e8ec844ee6192e0fd4e0e4105664cc5362 #define GEOGRAPHICLIB_GEODESIC_HPP \
a1660b $" "$Id: 81785396c1ad962e10595c8b81d9673658669591 $"
#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 \ #define GEOD_ORD \
(GEOGRAPHICLIB_PREC == 1 ? 6 : (GEOGRAPHICLIB_PREC == 0 ? 3 : 7)) (GEOGRAPHICLIB_PREC == 1 ? 6 : (GEOGRAPHICLIB_PREC == 0 ? 3 : 7))
#endif #endif
skipping to change at line 48 skipping to change at line 49
* 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, and \e azi2. This is the \e direct geodesic problem and its * lon2, and \e azi2. This is the \e direct geodesic problem and its
* solution is given by the function Geodesic::Direct. (If \e s12 is * solution is given by the function Geodesic::Direct. (If \e s12 is
* sufficiently large that the geodesic wraps more than halfway around th e * sufficiently large that the geodesic wraps more than halfway around th e
* earth, there will be another geodesic between the points with a smalle r \e * earth, there will be another geodesic between the points with a smalle r \e
* s12.) * 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, and \e s12. This is the \e inverse geodesic problem, whose solu tion * azi2, and \e s12. This is the \e inverse geodesic problem, whose solu tion
* is given by Geodesic::Inverse. Usually, the solution to the inverse * is given by Geodesic::Inverse. Usually, the solution to the inverse
* problem is unique. In cases where there are muliple solutions (all wi th * problem is unique. In cases where there are multiple solutions (all w ith
* the same \e s12, of course), all the solutions can be easily generated * the same \e s12, of course), all the solutions 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 * The standard way of specifying the direct problem is the specify the
* distance \e s12 to the second point. However it is sometimes useful * 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 * 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 * sphere. This is a mathematical construct used in solving the geodesic
* problems. The solution of the direct problem in this form is provide by * 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 * 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 * that the geodesic is not a shortest path. In addition, the arc length
skipping to change at line 70 skipping to change at line 71
* geodesic is 90<sup>o</sup>. * geodesic is 90<sup>o</sup>.
* *
* This class can also calculate several other quantities related to * This class can also calculate several other quantities related to
* geodesics. These are: * geodesics. These are:
* - <i>reduced length</i>. If we fix the first point and increase \e az i1 * - <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 * 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 * 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 * called the "reduced length" and is symmetric under interchange of th e
* two points. On a curved surface the reduced length obeys a symmetry * two points. On a curved surface the reduced length obeys a symmetry
* relation, \e m12 + \e m21 = 0. On a flat surface, we have \e m12 = \e * relation, \e m12 + \e m21 = 0. On a flat surface, we have \e m12 = \e
* s12. The ratio \e s12/\e m12 gives the azimuthal scale for an azimu * s12. The ratio <i>s12</i>/\e m12 gives the azimuthal scale for an
thal * azimuthal equidistant projection.
* equidistant projection.
* - <i>geodesic scale</i>. Consider a reference geodesic and a second * - <i>geodesic scale</i>. Consider a reference geodesic and a second
* geodesic parallel to this one at point 1 and separated by a small * 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 * 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 * 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 * 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 * flat surface, we have \e M12 = \e M21 = 1. The quantity 1/\e M12 gi ves
* the scale of the Cassini-Soldner projection. * the scale of the Cassini-Soldner projection.
* - <i>area</i>. Consider the quadrilateral bounded by the following li nes: * - <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 * 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, the equator from \e lon2 to \e lon1, the meridian from the
skipping to change at line 291 skipping to change at line 292
******************************************************************** **/ ******************************************************************** **/
ALL = OUT_ALL| CAP_ALL, ALL = OUT_ALL| CAP_ALL,
}; };
/** \name Constructor /** \name Constructor
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* Constructor for a ellipsoid with * Constructor for a ellipsoid with
* *
* @param[in] a equatorial radius (meters) * @param[in] a equatorial radius (meters).
* @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphe re. * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphe re.
* Negative \e f gives a prolate ellipsoid. If \e f > 1, set flatten ing * Negative \e f gives a prolate ellipsoid. If \e f > 1, set flatten ing
* to 1/\e f. * to 1/\e f.
* *
* An exception is thrown if either of the axes of the ellipsoid is * An exception is thrown if either of the axes of the ellipsoid is
* non-positive. * non-positive.
********************************************************************** / ********************************************************************** /
Geodesic(real a, real f); Geodesic(real a, real f);
///@} ///@}
 End of changes. 4 change blocks. 
7 lines changed or deleted 6 lines changed or added


 GeodesicLine.hpp   GeodesicLine.hpp 
/** /**
* \file GeodesicLine.hpp * \file GeodesicLine.hpp
* \brief Header for GeographicLib::GeodesicLine class * \brief Header for GeographicLib::GeodesicLine class
* *
* Copyright (c) Charles Karney (2009, 2010, 2011) <charles@karney.com> * Copyright (c) Charles Karney (2009, 2010, 2011) <charles@karney.com>
* and licensed under the MIT/X11 License. For more information, see * and licensed under the MIT/X11 License. For more information, see
* http://geographiclib.sourceforge.net/ * http://geographiclib.sourceforge.net/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_GEODESICLINE_HPP) #if !defined(GEOGRAPHICLIB_GEODESICLINE_HPP)
#define GEOGRAPHICLIB_GEODESICLINE_HPP "$Id: 24c298d5d343f35eaeeee27e3c65f2 #define GEOGRAPHICLIB_GEODESICLINE_HPP \
2cae729faa $" "$Id: b18e66a0299f857bf2d71c6debd98b7f14da8a64 $"
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
#include <GeographicLib/Geodesic.hpp> #include <GeographicLib/Geodesic.hpp>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief A geodesic line. * \brief A geodesic line.
* *
* GeodesicLine facilitates the determination of a series of points on a * GeodesicLine facilitates the determination of a series of points on a
* single geodesic. The starting point (\e lat1, \e lon1) and the azimut h \e * single geodesic. The starting point (\e lat1, \e lon1) and the azimut h \e
* azi1 are specified in the constructor. GeodesicLine.Position returns the * azi1 are specified in the constructor. GeodesicLine.Position returns the
* location of point 2 a distance \e s12 along the geodesic. Alternative ly * location of point 2 a distance \e s12 along the geodesic. Alternative ly
* GeodesicLine.ArcPosition gives the position of point 2 an arc length \ e * GeodesicLine.ArcPosition gives the position of point 2 an arc length \ e
* a12 along the geodesic. An example of use of this class is: * a12 along the geodesic. An example of use of this class is:
\code \code
#include <iostream>
#include <iomanip>
#include <cmath>
#include <GeographicLib/Geodesic.hpp> #include <GeographicLib/Geodesic.hpp>
#include <GeographicLib/GeodesicLine.hpp> #include <GeographicLib/GeodesicLine.hpp>
int main() { int main() {
// Print waypoints between JFK and SIN at approximately // Print waypoints between JFK and SIN at approximately
// equal 100km intervals // equal 100km intervals
double double
lat1 = 40.640, lon1 = -73.779, // JFK lat1 = 40.640, lon1 = -73.779, // JFK
lat2 = 1.359, lon2 = 177.486; // SIN lat2 = 1.359, lon2 = 177.486; // SIN
const GeographicLib::Geodesic& const GeographicLib::Geodesic&
skipping to change at line 116 skipping to change at line 114
CAP_C1p = Geodesic::CAP_C1p, CAP_C1p = Geodesic::CAP_C1p,
CAP_C2 = Geodesic::CAP_C2, CAP_C2 = Geodesic::CAP_C2,
CAP_C3 = Geodesic::CAP_C3, CAP_C3 = Geodesic::CAP_C3,
CAP_C4 = Geodesic::CAP_C4, CAP_C4 = Geodesic::CAP_C4,
CAP_ALL = Geodesic::CAP_ALL, CAP_ALL = Geodesic::CAP_ALL,
OUT_ALL = Geodesic::OUT_ALL, OUT_ALL = Geodesic::OUT_ALL,
}; };
public: public:
/** /**
* Bit masks for what calculations to do. These masks do double duty. * Bit masks for what calculations to do. They signify to the
* They signify to the GeodesicLine::GeodesicLine constructor and to * GeodesicLine::GeodesicLine constructor and to Geodesic::Line what
* Geodesic::Line what capabilities should be included in the GeodesicL * capabilities should be included in the GeodesicLine object. This is
ine
* object. They also specify which results to return in the general
* routines Geodesic::GenDirect and Geodesic::GenInverse routines. Thi
s is
* merely a duplication of Geodesic::mask. * merely a duplication of Geodesic::mask.
********************************************************************** / ********************************************************************** /
enum mask { enum mask {
/** /**
* No capabilities, no output. * No capabilities, no output.
* @hideinitializer * @hideinitializer
******************************************************************** **/ ******************************************************************** **/
NONE = Geodesic::NONE, NONE = Geodesic::NONE,
/** /**
* Calculate latitude \e lat2. (It's not necessary to include this a s a * Calculate latitude \e lat2. (It's not necessary to include this a s a
skipping to change at line 186 skipping to change at line 182
******************************************************************** **/ ******************************************************************** **/
ALL = Geodesic::ALL, ALL = Geodesic::ALL,
}; };
/** \name Constructors /** \name Constructors
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* Constructor for a geodesic line staring at latitude \e lat1, longitu de * Constructor for a geodesic line staring at latitude \e lat1, longitu de
* \e lon1, and aziumuth \e azi1 (all in degrees). * \e lon1, and azimuth \e azi1 (all in degrees).
* *
* @param[in] g A Geodesic object used to compute the necessary informa tion * @param[in] g A Geodesic object used to compute the necessary informa tion
* about the GeodesicLine. * about the GeodesicLine.
* *
* @param[in] lat1 latitude of point 1 (degrees). * @param[in] lat1 latitude of point 1 (degrees).
* @param[in] lon1 longitude of point 1 (degrees). * @param[in] lon1 longitude of point 1 (degrees).
* @param[in] azi1 azimuth at point 1 (degrees). * @param[in] azi1 azimuth at point 1 (degrees).
* @param[in] caps bitor'ed combination of GeodesicLine::mask values * @param[in] caps bitor'ed combination of GeodesicLine::mask values
* specifying the capabilities the GeodesicLine object should possess , * specifying the capabilities the GeodesicLine object should possess ,
* i.e., which quantities can be returned in calls to * i.e., which quantities can be returned in calls to
 End of changes. 4 change blocks. 
13 lines changed or deleted 6 lines changed or added


 Geoid.hpp   Geoid.hpp 
/** /**
* \file Geoid.hpp * \file Geoid.hpp
* \brief Header for GeographicLib::Geoid class * \brief Header for GeographicLib::Geoid class
* *
* Copyright (c) Charles Karney (2009, 2010, 2011) <charles@karney.com> * Copyright (c) Charles Karney (2009, 2010, 2011) <charles@karney.com>
* and licensed under the MIT/X11 License. For more information, see * and licensed under the MIT/X11 License. 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: ea21a7ef2bf421621967b0c442da427106c8d #define GEOGRAPHICLIB_GEOID_HPP \
6c2 $" "$Id: 73a380e4dcfed4a6ec72c818622ec45f8daf3352 $"
#include <string> #include <string>
#include <vector> #include <vector>
#include <fstream> #include <fstream>
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
#if defined(_MSC_VER) #if defined(_MSC_VER)
// Squelch warnings about dll vs vector // Squelch warnings about dll vs vector
#pragma warning (push) #pragma warning (push)
#pragma warning (disable: 4251) #pragma warning (disable: 4251)
skipping to change at line 41 skipping to change at line 42
* grid of data. These geoid models are documented in * grid of data. These geoid models are documented in
* - EGM84: * - EGM84:
* http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180 .html * http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180 .html
* - EGM96: * - EGM96:
* http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html * http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html
* - EGM2008: * - EGM2008:
* http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008 * http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008
* *
* The geoids are defined in terms of spherical harmonics. However in or der * The geoids are defined in terms of spherical harmonics. However in or der
* to provide a quick and flexible method of evaluating the geoid heights , * to provide a quick and flexible method of evaluating the geoid heights ,
* this class evaluates the height by interpolation inot a grid of * this class evaluates the height by interpolation into a grid of
* precomputed values. * precomputed values.
* *
* See \ref geoid for details of how to install the data sets, the data * See \ref geoid for details of how to install the data sets, the data
* format, estimates of the interpolation errors, and how to use caching. * format, estimates of the interpolation errors, and how to use caching.
* *
* In addition to returning the geoid height, the gradient of the geoid c an * In addition to returning the geoid height, the gradient of the geoid c an
* be calculated. The gradient is defined as the rate of change of the g eoid * be calculated. The gradient is defined as the rate of change of the g eoid
* as a function of position on the ellipsoid. This uses the parameters for * as a function of position on the ellipsoid. This uses the parameters for
* the WGS84 ellipsoid. The gradient defined in terms of the interpolate d * the WGS84 ellipsoid. The gradient defined in terms of the interpolate d
* heights. As a result of the way that the geoid data is stored, the * heights. As a result of the way that the geoid data is stored, the
* calculation of gradients can result in large quantization errors. Thi s is * calculation of gradients can result in large quantization errors. Thi s is
* particularly acute for fine grids, at high latitudes, and for the east erly * particularly acute for fine grids, at high latitudes, and for the east erly
* gradient. * gradient.
* *
* This class is typically \e not thread safe in that a single instantiat ion * This class is typically \e not thread safe in that a single instantiat ion
* cannot be safely used by multiple threads because of the way the objec t * cannot be safely used by multiple threads because of the way the objec t
* reads the data set and because it maintains a single-cell cache. If * reads the data set and because it maintains a single-cell cache. If
* multiple threads need to calculate geoid heights they should all const ruct * multiple threads need to calculate geoid heights they should all const ruct
* thread-local instantiations. Alternatively, set the optional \e * thread-local instantiations. Alternatively, set the optional \e
* threadsafe parameter to true in the constuctor. This causes the * threadsafe parameter to true in the constructor. This causes the
* constructor to read all the data into memory and to turn off the * constructor to read all the data into memory and to turn off the
* single-cell caching which results in a Geoid object which \e is thread * single-cell caching which results in a Geoid object which \e is thread
* safe. * safe.
**********************************************************************/ **********************************************************************/
class GEOGRAPHIC_EXPORT Geoid { class GEOGRAPHIC_EXPORT Geoid {
private: private:
typedef Math::real real; typedef Math::real real;
static const unsigned stencilsize_ = 12; static const unsigned stencilsize_ = 12;
static const unsigned nterms_ = ((3 + 1) * (3 + 2))/2; // for a cubic f it static const unsigned nterms_ = ((3 + 1) * (3 + 2))/2; // for a cubic f it
skipping to change at line 183 skipping to change at line 184
* *
* @param[in] name the name of the geoid. * @param[in] name the name of the geoid.
* @param[in] path (optional) directory for data file. * @param[in] path (optional) directory for data file.
* @param[in] cubic (optional) interpolation method; false means biline ar, * @param[in] cubic (optional) interpolation method; false means biline ar,
* true (the default) means cubic. * true (the default) means cubic.
* @param[in] threadsafe (optional), if true, construct a thread safe * @param[in] threadsafe (optional), if true, construct a thread safe
* object. The default is false * object. The default is false
* *
* The data file is formed by appending ".pgm" to the name. If \e path is * The data file is formed by appending ".pgm" to the name. If \e path is
* specified (and is non-empty), then the file is loaded from directory , \e * specified (and is non-empty), then the file is loaded from directory , \e
* path. Otherwise the path is given by the GEOID_PATH environment * path. Otherwise the path is given by DefaultGeoidPath(). This may
* variable. If that is undefined, a compile-time default path is used * throw an exception because the file does not exist, is unreadable, o
* (/usr/local/share/GeographicLib/geoids on non-Windows systems and r is
* C:/Documents and Settings/All Users/Application * corrupt. If the \e threadsafe parameter is true, the data set is re
* Data/GeographicLib/geoids on Windows systems). This may throw an ad
* exception because the file does not exist, is unreadable, or is corr * into memory (which this may also cause an exception to be thrown), t
upt. he
* If the \e threadsafe parameter is true, the data set is read into me * data file is closed, and single-cell caching is turned off; this res
mory ults
* (which this may also cause an exception to be thrown), the data file * in a Geoid object which \e is thread safe.
is
* closed, and single-cell caching is turned off; this results in a Geo
id
* object which \e is thread safe.
********************************************************************** / ********************************************************************** /
explicit Geoid(const std::string& name, const std::string& path = "", explicit Geoid(const std::string& name, const std::string& path = "",
bool cubic = true, bool threadsafe = false); bool cubic = true, bool threadsafe = false);
/** /**
* Set up a cache. * Set up a cache.
* *
* @param[in] south latitude (degrees) of the south edge of the cached area. * @param[in] south latitude (degrees) of the south edge of the cached area.
* @param[in] west longitude (degrees) of the west edge of the cached a rea. * @param[in] west longitude (degrees) of the west edge of the cached a rea.
* @param[in] north latitude (degrees) of the north edge of the cached area. * @param[in] north latitude (degrees) of the north edge of the cached area.
* @param[in] east longitude (degrees) of the east edge of the cached a rea. * @param[in] east longitude (degrees) of the east edge of the cached a rea.
* *
* Cache the data for the specified "rectangular" area bounded by the * Cache the data for the specified "rectangular" area bounded by the
* parallels \e south and \e north and the meridians \e west and \e eas t. * parallels \e south and \e north and the meridians \e west and \e eas t.
* \e east is always interpreted as being east of \e west, if necessary by * \e east is always interpreted as being east of \e west, if necessary by
* adding 360<sup>o</sup> to its value. This may throw an error becaus e of * adding 360<sup>o</sup> to its value. This may throw an error becaus e of
* insufficent memory or because of an error reading the data from the * insufficient 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]. An exception is thrown if this * should be in the range [-180, 360]. An exception is thrown if this
* routine is called on a thread safe Geoid. * routine is called on a thread safe Geoid.
********************************************************************** / ********************************************************************** /
void CacheArea(real south, real west, real north, real east) const; void CacheArea(real south, real west, real north, real east) const;
/** /**
* Cache all the data. On most computers, this is fast for data sets w ith * Cache all the data. On most computers, this is fast for data sets w ith
* grid resolution of 5' or coarser. For a 1' grid, the required RAM i s * grid resolution of 5' or coarser. For a 1' grid, the required RAM i s
* 450MB; a 2.5' grid needs 72MB; and a 5' grid needs 18MB. This may t hrow * 450MB; a 2.5' grid needs 72MB; and a 5' grid needs 18MB. This may t hrow
* an error because of insufficent memory or because of an error readin g * an error because of insufficient memory or because of an error readi ng
* the data from the file. In this case, you can catch the error and * the data from the file. In this case, you can catch the error and
* either do nothing (you will have no cache in this case) or try using * either do nothing (you will have no cache in this case) or try using
* Geoid::CacheArea on a specific area. An exception is thrown if this * Geoid::CacheArea on a specific area. An exception is thrown if this
* routine is called on a thread safe Geoid. * routine is called on a thread safe Geoid.
********************************************************************** / ********************************************************************** /
void CacheAll() const { CacheArea(real(-90), real(0), void CacheAll() const { CacheArea(real(-90), real(0),
real(90), real(360)); } real(90), real(360)); }
/** /**
* Clear the cache. This never throws an error. (This does nothing wi th a * Clear the cache. This never throws an error. (This does nothing wi th a
skipping to change at line 249 skipping to change at line 246
/** \name Compute geoid heights /** \name Compute geoid heights
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* Compute the geoid height at a point * Compute the geoid height at a point
* *
* @param[in] lat latitude of the point (degrees). * @param[in] lat latitude of the point (degrees).
* @param[in] lon longitude of the point (degrees). * @param[in] lon longitude of the point (degrees).
* @return geoid height (meters). * @return geoid height (meters).
* *
* The latitude should be in [-90, 90] and longitude should bein * The latitude should be in [-90, 90] and longitude should be in
* [-180,360]. This may throw an error because of an error reading dat a * [-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 * from disk. However, it will not throw if (\e lat, \e lon) is within a
* successfully cached area. * 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);
} }
/** /**
skipping to change at line 440 skipping to change at line 437
* <b>DEPRECATED</b> * <b>DEPRECATED</b>
* @return \e r the inverse flattening of the WGS84 ellipsoid. * @return \e r the inverse flattening of the WGS84 ellipsoid.
********************************************************************** / ********************************************************************** /
Math::real InverseFlattening() const throw() Math::real InverseFlattening() const throw()
{ return 1/Constants::WGS84_f<real>(); } { return 1/Constants::WGS84_f<real>(); }
///@} ///@}
/** /**
* @return the default path for geoid data files. * @return the default path for geoid data files.
* *
* This is the value of the environment variable GEOID_PATH, if set, * This is the value of the environment variable GEOID_PATH, if set;
* otherwise, it is a compile-time default. * otherwise, it is $GEOGRAPHICLIB_DATA/geoids if the environment varia
ble
* GEOGRAPHICLIB_DATA is set; otherwise, it is a compile-time default
* (/usr/local/share/GeographicLib/geoids on non-Windows systems and
* C:/Documents and Settings/All Users/Application
* Data/GeographicLib/geoids on Windows systems).
********************************************************************** / ********************************************************************** /
static std::string DefaultGeoidPath(); static std::string DefaultGeoidPath();
/** /**
* @return the default name for the geoid. * @return the default name for the geoid.
* *
* This is the value of the environment variable GEOID_NAME, if set, * This is the value of the environment variable GEOID_NAME, if set,
* otherwise, it is "egm96-5". The Geoid class does not use this funct ion; * otherwise, it is "egm96-5". The Geoid class does not use this funct ion;
* it is just provided as a convenience for a calling program when * it is just provided as a convenience for a calling program when
* constructing a Geoid object. * constructing a Geoid object.
 End of changes. 8 change blocks. 
23 lines changed or deleted 24 lines changed or added


 Gnomonic.hpp   Gnomonic.hpp 
/** /**
* \file Gnomonic.hpp * \file Gnomonic.hpp
* \brief Header for GeographicLib::Gnomonic class * \brief Header for GeographicLib::Gnomonic class
* *
* Copyright (c) Charles Karney (2010, 2011) <charles@karney.com> and licen sed * Copyright (c) Charles Karney (2010, 2011) <charles@karney.com> and licen sed
* under the MIT/X11 License. For more information, see * under the MIT/X11 License. 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: 35168bbc58f320e1e117301fd7653b8835 #define GEOGRAPHICLIB_GNOMONIC_HPP \
407792 $" "$Id: 3ae0cfd77aad984813d955a03726248ac6b427a9 $"
#include <GeographicLib/Geodesic.hpp> #include <GeographicLib/Geodesic.hpp>
#include <GeographicLib/GeodesicLine.hpp> #include <GeographicLib/GeodesicLine.hpp>
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief %Gnomonic Projection. * \brief %Gnomonic Projection.
* *
skipping to change at line 40 skipping to change at line 41
* . * .
* See also Section 8 of * See also Section 8 of
* - C. F. F. Karney, * - C. F. F. Karney,
* <a href="http://arxiv.org/abs/1109.4448">Algorithms for geodesics</a >, * <a href="http://arxiv.org/abs/1109.4448">Algorithms for geodesics</a >,
* Sept. 2011; * Sept. 2011;
* preprint * preprint
* <a href="http://arxiv.org/abs/1109.4448">arxiv:1109.4448</a>. * <a href="http://arxiv.org/abs/1109.4448">arxiv:1109.4448</a>.
* . * .
* The projection of \e P is defined as follows: compute the * The projection of \e P is defined as follows: compute the
* geodesic line from \e C to \e P; compute the reduced length \e m12, * geodesic line from \e C to \e P; compute the reduced length \e m12,
* geodesic scale \e M12, and \e rho = \e m12/\e M12; finally \e x = \e r * geodesic scale \e M12, and \e rho = <i>m12</i>/\e M12; finally \e x =
ho \e
* sin \e azi1; \e y = \e rho cos \e azi1, where \e azi1 is the azimuth o * rho sin \e azi1; \e y = \e rho cos \e azi1, where \e azi1 is the azimu
f th
* the geodesic at \e C. The Gnomonic::Forward and Gnomonic::Reverse met * of the geodesic at \e C. The Gnomonic::Forward and Gnomonic::Reverse
hods * methods also return the azimuth \e azi of the geodesic at \e P and
* also return the azimuth \e azi of the geodesic at \e P and reciprocal * reciprocal scale \e rk in the azimuthal direction. The scale in the
* scale \e rk in the azimuthal direction. The scale in the radial direc * radial direction if 1/<i>rk</i><sup>2</sup>.
tion *
* if 1/\e rk<sup>2</sup>. * For a sphere, \e rho is reduces to \e a tan(<i>s12</i>/<i>a</i>), wher
* e \e
* For a sphere, \e rho is reduces to \e a tan(\e s12/\e a), where \e s12 * s12 is the length of the geodesic from \e C to \e P, and the gnomonic
is * projection has the property that all geodesics appear as straight line
* the length of the geodesic from \e C to \e P, and the gnomonic project s.
ion * For an ellipsoid, this property holds only for geodesics interesting t
* has the property that all geodesics appear as straight lines. For an he
* ellipsoid, this property holds only for geodesics interesting the cent * centers. However geodesic segments close to the center are approximat
ers. ely
* However geodesic segments close to the center are approximately straig * straight.
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 (extended if necessary) closest to \e C the center of the * geodesic (extended if necessary) closest to \e C the center of the
* projection and \e t be the distance \e CT. To lowest order, the maxim um * projection and \e t be the distance \e CT. To lowest order, the maxim um
* deviation (as a true distance) of the corresponding gnomonic line segm ent * deviation (as a true distance) of the corresponding gnomonic line segm ent
* (i.e., with the same end points) from the geodesic is<br> * (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> * (<i>K</i>(<i>T</i>) - <i>K</i>(<i>C</i>))
* <i>l</i><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 revolution, * This result applies for any surface. For an ellipsoid 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. * (<i>r</i>/2<i>a</i>)<sup>3</sup> \e r.
* *
* The conversions all take place using a Geodesic object (by default * The conversions all take place using a Geodesic object (by default
* Geodesic::WGS84). For more information on geodesics see \ref geodesic . * Geodesic::WGS84). For more information on geodesics see \ref geodesic .
* *
* <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;
skipping to change at line 134 skipping to change at line 137
* @param[in] lon0 longitude of center point of projection (degrees). * @param[in] lon0 longitude of center point of projection (degrees).
* @param[in] lat latitude of point (degrees). * @param[in] lat latitude of point (degrees).
* @param[in] lon longitude of point (degrees). * @param[in] lon longitude of point (degrees).
* @param[out] x easting of point (meters). * @param[out] x easting of point (meters).
* @param[out] y northing of point (meters). * @param[out] y northing of point (meters).
* @param[out] azi azimuth of geodesic at point (degrees). * @param[out] azi azimuth of geodesic at point (degrees).
* @param[out] rk reciprocal of azimuthal scale at point. * @param[out] rk reciprocal of azimuthal scale at point.
* *
* \e lat0 and \e lat should be in the range [-90, 90] and \e lon0 and \e * \e lat0 and \e lat should be in the range [-90, 90] and \e lon0 and \e
* lon should be in the range [-180, 360]. The scale of the projection is * lon should be in the range [-180, 360]. The scale of the projection is
* 1/\e rk<sup>2</sup> in the "radial" direction, \e azi clockwise from * 1/<i>rk</i><sup>2</sup> in the "radial" direction, \e azi clockwise from
* true north, and is 1/\e rk in the direction perpendicular to this. If * true north, and is 1/\e rk in the direction perpendicular to this. If
* the point lies "over the horizon", i.e., if \e rk <= 0, then NaNs ar e * the point lies "over the horizon", i.e., if \e rk <= 0, then NaNs ar e
* returned for \e x and \e y (the correct values are returned for \e a zi * 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 * 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 * the original (\e lat, \e lon) (to within roundoff) provided the poin t in
* not over the horizon. * 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();
 End of changes. 6 change blocks. 
26 lines changed or deleted 26 lines changed or added


 LambertConformalConic.hpp   LambertConformalConic.hpp 
/** /**
* \file LambertConformalConic.hpp * \file LambertConformalConic.hpp
* \brief Header for GeographicLib::LambertConformalConic class * \brief Header for GeographicLib::LambertConformalConic class
* *
* Copyright (c) Charles Karney (2010, 2011) <charles@karney.com> and licen sed * Copyright (c) Charles Karney (2010, 2011) <charles@karney.com> and licen sed
* under the MIT/X11 License. For more information, see * under the MIT/X11 License. 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: 22e444fc38e3e8b49a56c #define GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_HPP \
7a22121ddeb78c230b9 $" "$Id: 7975cff875e88fd42f81f9e1d7117c1dd36667d2 $"
#include <algorithm> #include <algorithm>
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief Lambert Conformal Conic Projection * \brief Lambert Conformal Conic Projection
* *
* Implementation taken from the report, * Implementation taken from the report,
* - J. P. Snyder, * - J. P. Snyder,
* <a href="http://pubs.er.usgs.gov/usgspubs/pp/pp1395"> Map Projection s: A * <a href="http://pubs.er.usgs.gov/usgspubs/pp/pp1395"> Map Projection s: A
* Working Manual</a>, USGS Professional Paper 1395 (1987), * Working Manual</a>, USGS Professional Paper 1395 (1987),
* pp. 107&ndash;109. * pp. 107&ndash;109.
* *
* This is a implementation of the equations in Snyder except that divide d * This is a implementation of the equations in Snyder except that divide d
* differences have been used to transform the expressions into ones whic h * differences have been used to transform the expressions into ones whic h
* may be evaluated accurately and that Newton's method is used to invert the * may be evaluated accurately and that Newton's method is used to invert the
* projection. In this implementation, the projection correctly becomes the * projection. In this implementation, the projection correctly becomes the
* Mercator projection or the polar sterographic projection when the stan * Mercator projection or the polar stereographic projection when the
dard * standard latitude is the equator or a pole. The accuracy of the
* latitude is the equator or a pole. The accuracy of the projections is * projections is about 10 nm.
* about 10 nm.
* *
* The ellipsoid parameters, the standard parallels, and the scale on the * The ellipsoid parameters, the standard parallels, and the scale on the
* standard parallels are set in the constructor. Internally, the case w ith * standard parallels are set in the constructor. Internally, the case w ith
* two standard parallels is converted into a single standard parallel, t he * two standard parallels is converted into a single standard parallel, t he
* latitude of tangency (also the latitude of minimum scale), with a scal e * latitude of tangency (also the latitude of minimum scale), with a scal e
* specified on this parallel. This latitude is also used as the latitud e of * specified on this parallel. This latitude is also used as the latitud e of
* origin which is returned by LambertConformalConic::OriginLatitude. Th e * origin which is returned by LambertConformalConic::OriginLatitude. Th e
* scale on the latitude of origin is given by * scale on the latitude of origin is given by
* LambertConformalConic::CentralScale. The case with two distinct stand ard * LambertConformalConic::CentralScale. The case with two distinct stand ard
* parallels where one is a pole is singular and is disallowed. The cent ral * parallels where one is a pole is singular and is disallowed. The cent ral
skipping to change at line 166 skipping to change at line 167
inline real Deatanhe(real x, real y) const throw() { inline real Deatanhe(real x, real y) const throw() {
real t = x - y, d = 1 - _e2 * x * y; real t = x - y, d = 1 - _e2 * x * y;
return t != 0 ? eatanhe(t / d) / t : _e2 / d; return t != 0 ? eatanhe(t / d) / t : _e2 / d;
} }
void Init(real sphi1, real cphi1, real sphi2, real cphi2, real k1) thro w(); void Init(real sphi1, real cphi1, real sphi2, real cphi2, real k1) thro w();
public: public:
/** /**
* Constructor with a single standard parallel. * Constructor with a single standard parallel.
* *
* @param[in] a equatorial radius of ellipsoid (meters) * @param[in] a equatorial radius of ellipsoid (meters).
* @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphe re. * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphe re.
* Negative \e f gives a prolate ellipsoid. If \e f > 1, set flatten ing * Negative \e f gives a prolate ellipsoid. If \e f > 1, set flatten ing
* to 1/\e f. * to 1/\e f.
* @param[in] stdlat standard parallel (degrees), the circle of tangenc y. * @param[in] stdlat standard parallel (degrees), the circle of tangenc y.
* @param[in] k0 scale on the standard parallel. * @param[in] k0 scale on the standard parallel.
* *
* An exception is thrown if \e a or \e k0 is not positive or if \e std lat * 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]. * is not in the range [-90, 90].
********************************************************************** / ********************************************************************** /
LambertConformalConic(real a, real f, real stdlat, real k0); LambertConformalConic(real a, real f, real stdlat, real k0);
/** /**
* Constructor with two standard parallels. * Constructor with two standard parallels.
* *
* @param[in] a equatorial radius of ellipsoid (meters) * @param[in] a equatorial radius of ellipsoid (meters).
* @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphe re. * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphe re.
* Negative \e f gives a prolate ellipsoid. If \e f > 1, set flatten ing * Negative \e f gives a prolate ellipsoid. If \e f > 1, set flatten ing
* to 1/\e f. * to 1/\e f.
* @param[in] stdlat1 first standard parallel (degrees). * @param[in] stdlat1 first standard parallel (degrees).
* @param[in] stdlat2 second standard parallel (degrees). * @param[in] stdlat2 second standard parallel (degrees).
* @param[in] k1 scale on the standard parallels. * @param[in] k1 scale on the standard parallels.
* *
* An exception is thrown if \e a or \e k0 is not positive or if \e std lat1 * An exception is thrown if \e a or \e k0 is not positive or if \e std lat1
* or \e stdlat2 is not in the range [-90, 90]. In addition, if either \e * 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 or \e stdlat2 is a pole, then an exception is thrown if \e
* stdlat1 is not equal \e stdlat2. * stdlat1 is not equal \e stdlat2.
********************************************************************** / ********************************************************************** /
LambertConformalConic(real a, real f, real stdlat1, real stdlat2, real k1); LambertConformalConic(real a, real f, real stdlat1, real stdlat2, real k1);
/** /**
* Constructor with two standard parallels specified by sines and cosin es. * Constructor with two standard parallels specified by sines and cosin es.
* *
* @param[in] a equatorial radius of ellipsoid (meters) * @param[in] a equatorial radius of ellipsoid (meters).
* @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphe re. * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphe re.
* Negative \e f gives a prolate ellipsoid. If \e f > 1, set flatten ing * Negative \e f gives a prolate ellipsoid. If \e f > 1, set flatten ing
* to 1/\e f. * to 1/\e f.
* @param[in] sinlat1 sine of first standard parallel. * @param[in] sinlat1 sine of first standard parallel.
* @param[in] coslat1 cosine of first standard parallel. * @param[in] coslat1 cosine of first standard parallel.
* @param[in] sinlat2 sine of second standard parallel. * @param[in] sinlat2 sine of second standard parallel.
* @param[in] coslat2 cosine of second standard parallel. * @param[in] coslat2 cosine of second standard parallel.
* @param[in] k1 scale on the standard parallels. * @param[in] k1 scale on the standard parallels.
* *
* This allows parallels close to the poles to be specified accurately. * This allows parallels close to the poles to be specified accurately.
* This routine computes the latitude of origin and the scale at this * This routine computes the latitude of origin and the scale at this
* latitude. In the case where \e lat1 and \e lat2 are different, the * latitude. In the case where \e lat1 and \e lat2 are different, the
* errors in this routines are as follows: if \e dlat = abs(\e lat2 - \ e * errors in this routines are as follows: if \e dlat = abs(\e lat2 - \ e
* lat1) <= 160<sup>o</sup> and max(abs(\e lat1), abs(\e lat2)) <= 90 - * lat1) <= 160<sup>o</sup> and max(abs(\e lat1), abs(\e lat2)) <= 90 -
* min(0.0002, 2.2e-6(180 - \e dlat), 6e-8\e dlat<sup>2</sup>) (in * min(0.0002, 2.2e-6(180 - \e dlat), 6e-8 <i>dlat</i><sup>2</sup>) (in
* degrees), then the error in the latitude of origin is less than * degrees), then the error in the latitude of origin is less than
* 4.5e-14<sup>o</sup> and the relative error in the scale is less than * 4.5e-14<sup>o</sup> and the relative error in the scale is less than
* 7e-15. * 7e-15.
********************************************************************** / ********************************************************************** /
LambertConformalConic(real a, real f, LambertConformalConic(real a, real f,
real sinlat1, real coslat1, real sinlat1, real coslat1,
real sinlat2, real coslat2, real sinlat2, real coslat2,
real k1); real k1);
/** /**
 End of changes. 6 change blocks. 
10 lines changed or deleted 9 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, 2011) <charles@karney.co m> * Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co m>
* and licensed under the MIT/X11 License. For more information, see * and licensed under the MIT/X11 License. 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: 66ad3a5333b92c5ff6f373080e62 #define GEOGRAPHICLIB_LOCALCARTESIAN_HPP \
c5392841afd3 $" "$Id: 0200a2653b8473d3fbf1f1382291fbff5294821c $"
#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
 End of changes. 1 change blocks. 
2 lines changed or deleted 2 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, 2011) <charles@karney.co m> * Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co m>
* and licensed under the MIT/X11 License. For more information, see * and licensed under the MIT/X11 License. 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: b15aa92f13463a8432dc6ff095cbb2650eb4a9 a2 $" #define GEOGRAPHICLIB_MGRS_HPP "$Id: 9f5fb4813d4ed429fee1570bf9378ad7b8f495 9b $"
#include <sstream> #include <sstream>
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
#include <GeographicLib/UTMUPS.hpp> #include <GeographicLib/UTMUPS.hpp>
#if defined(_MSC_VER) #if defined(_MSC_VER)
// Squelch warnings about dll vs string // Squelch warnings about dll vs string
#pragma warning (push) #pragma warning (push)
#pragma warning (disable: 4251) #pragma warning (disable: 4251)
#endif #endif
skipping to change at line 133 skipping to change at line 133
/** /**
* Convert UTM or UPS coordinate to an MGRS coordinate. * Convert UTM or UPS coordinate to an MGRS coordinate.
* *
* @param[in] zone UTM zone (zero means UPS). * @param[in] zone UTM zone (zero means UPS).
* @param[in] northp hemisphere (true means north, false means south). * @param[in] northp hemisphere (true means north, false means south).
* @param[in] x easting of point (meters). * @param[in] x easting of point (meters).
* @param[in] y northing of point (meters). * @param[in] y northing of point (meters).
* @param[in] prec precision relative to 100 km. * @param[in] prec precision relative to 100 km.
* @param[out] mgrs MGRS string. * @param[out] mgrs MGRS string.
* *
* \e prec specifies the precision of the MSGRS string as follows: * \e prec specifies the precision of the MGRS string as follows:
* - prec = 0 (min), 100 km * - prec = 0 (min), 100 km
* - prec = 1, 10 km * - prec = 1, 10 km
* - prec = 2, 1 km * - prec = 2, 1 km
* - prec = 3, 100 m * - prec = 3, 100 m
* - prec = 4, 10 m * - prec = 4, 10 m
* - prec = 5, 1 m * - prec = 5, 1 m
* - prec = 6, 0.1 m * - prec = 6, 0.1 m
* - prec = 11 (max), 1 um * - prec = 11 (max), 1 um
* *
* 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
skipping to change at line 161 skipping to change at line 161
* 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
* hemisphere. * hemisphere.
* *
* The ranges are 100 km more restrictive that for the conversion betwe en * The ranges are 100 km more restrictive that for the conversion betwe en
* geographic coordinates and UTM and UPS given by UTMUPS. These * geographic coordinates and UTM and UPS given by UTMUPS. These
* restrictions are dictated by the allowed letters in MGRS coordinates . * restrictions are dictated by the allowed letters in MGRS coordinates .
* The choice of 9500 km for the maximum northing for northern hemisphe re * The choice of 9500 km for the maximum northing for northern hemisphe re
* and of 1000 km as the minimum northing for southern hemisphere provi de * and of 1000 km as the minimum northing for southern hemisphere provi de
* at least 0.5 degree extension into standard UPS zones. The upper en ds * at least 0.5 degree extension into standard UPS zones. The upper en ds
* of the ranges for the UPS coordinates is dictated by requiring symme try * of the ranges for the UPS coordinates is dictated by requiring symme try
* about the meridans 0E and 90E. * about the meridians 0E and 90E.
* *
* All allowed UTM and UPS coordinates may now be converted to legal MG RS * All allowed UTM and UPS coordinates may now be converted to legal MG RS
* coordinates with the proviso that eastings and northings on the uppe r * coordinates with the proviso that eastings and northings on the uppe r
* boundaries are silently reduced by about 4 nm to place them \e withi n the * boundaries are silently reduced by about 4 nm to place them \e withi n the
* allowed range. (This includes reducing a southern hemisphere northi ng * allowed range. (This includes reducing a southern hemisphere northi ng
* of 10000 km by 4 nm so that it is placed in latitude band M.) The U TM or * of 10000 km by 4 nm so that it is placed in latitude band M.) The U TM or
* UPS coordinates are truncated to requested precision to determine th e * UPS coordinates are truncated to requested precision to determine th e
* MGRS coordinate. Thus in UTM zone 38N, the square area with easting in * MGRS coordinate. Thus in UTM zone 38N, the square area with easting in
* [444 km, 445 km) and northing in [3688 km, 3689 km) maps to MGRS * [444 km, 445 km) and northing in [3688 km, 3689 km) maps to MGRS
* coordinate 38SMB4488 (at \e prec = 2, 1 km), Khulani Sq., Baghdad. * coordinate 38SMB4488 (at \e prec = 2, 1 km), Khulani Sq., Baghdad.
 End of changes. 3 change blocks. 
3 lines changed or deleted 3 lines changed or added


 MagneticCircle.hpp   MagneticCircle.hpp 
/** /**
* \file MagneticCircle.hpp * \file MagneticCircle.hpp
* \brief Header for GeographicLib::MagneticCircle class * \brief Header for GeographicLib::MagneticCircle class
* *
* Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed un der * Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed un der
* the MIT/X11 License. For more information, see * the MIT/X11 License. For more information, see
* http://geographiclib.sourceforge.net/ * http://geographiclib.sourceforge.net/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_MAGNETICCIRCLE_HPP) #if !defined(GEOGRAPHICLIB_MAGNETICCIRCLE_HPP)
#define GEOGRAPHICLIB_MAGNETICCIRCLE_HPP "$Id: de919c519a73ddbaa944ebdd377e #define GEOGRAPHICLIB_MAGNETICCIRCLE_HPP \
514d0e7899cb $" "$Id: 070a0f02ac2dae39473f225ebf9b827ac75c826d $"
#include <string> #include <string>
#include <vector> #include <vector>
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
#include <GeographicLib/CircularEngine.hpp> #include <GeographicLib/CircularEngine.hpp>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief Geomagnetic field on a circle of latitude * \brief Geomagnetic field on a circle of latitude
* *
* Evaluate the earth's magnetic field on a circle of constant height and * Evaluate the earth's magnetic field on a circle of constant height and
* latitude. This uses a CircleEngine to pre-evaluate the inner sum of t he * latitude. This uses a CircleEngine to pre-evaluate the inner sum of t he
* spherical harmonic sum, allowing the values of the field at different * spherical harmonic sum, allowing the values of the field at several
* latitude to be evaluated rapidly. * different longitudes to be evaluated rapidly.
* *
* Use MagneticModel::Circle to create a MagneticCircle object. (The * Use MagneticModel::Circle to create a MagneticCircle object. (The
* constructor for this class is private.) * constructor for this class is private.)
**********************************************************************/ **********************************************************************/
class GEOGRAPHIC_EXPORT MagneticCircle { class GEOGRAPHIC_EXPORT MagneticCircle {
private: private:
typedef Math::real real; typedef Math::real real;
real _a, _cphi, _sphi, _t, _dt0; real _a, _f, _lat, _h, _t, _cphi, _sphi, _t1, _dt0;
bool _interpolate; bool _interpolate;
CircularEngine _circ0, _circ1; CircularEngine _circ0, _circ1;
MagneticCircle(real a, real cphi, real sphi, real t, real dt0, MagneticCircle(real a, real f, real lat, real h, real t,
real cphi, real sphi, real t1, real dt0,
bool interpolate, bool interpolate,
const CircularEngine& circ0, const CircularEngine& circ1 ) const CircularEngine& circ0, const CircularEngine& circ1 )
: _a(a) : _a(a)
, _f(f)
, _lat(lat)
, _h(h)
, _t(t)
, _cphi(cphi) , _cphi(cphi)
, _sphi(sphi) , _sphi(sphi)
, _t(t) , _t1(t1)
, _dt0(dt0) , _dt0(dt0)
, _interpolate(interpolate) , _interpolate(interpolate)
, _circ0(circ0) , _circ0(circ0)
, _circ1(circ1) , _circ1(circ1)
{} {}
void Field(real lon, bool diffp, void Field(real lon, bool diffp,
real& Bx, real& By, real& Bz, real& Bx, real& By, real& Bz,
real& Bxt, real& Byt, real& Bzt) const; real& Bxt, real& Byt, real& Bzt) const throw();
friend class MagneticModel; // MagneticModel calls the private construc tor friend class MagneticModel; // MagneticModel calls the private construc tor
public: public:
/** /**
* A default constructor for the normal gravity. This sets up an
* uninitialized object which can be later replaced by the
* MagneticModel::Circle.
**********************************************************************
/
MagneticCircle() : _a(-1) {}
/** \name Compute the magnetic field
**********************************************************************
/
///@{
/**
* Evaluate the components of the geomagnetic field at a particular * Evaluate the components of the geomagnetic field at a particular
* longitude. * longitude.
* *
* @param[in] lon longitude of the point (degrees). * @param[in] lon longitude of the point (degrees).
* @param[out] Bx the easterly component of the magnetic field (nanotes la). * @param[out] Bx the easterly component of the magnetic field (nanotes la).
* @param[out] By the northerly component of the magnetic field (nanote sla). * @param[out] By the northerly component of the magnetic field (nanote sla).
* @param[out] Bz the vertical (up) component of the magnetic field * @param[out] Bz the vertical (up) component of the magnetic field
* (nanotesla). * (nanotesla).
********************************************************************** / ********************************************************************** /
void operator()(real lon, real& Bx, real& By, real& Bz) const { void operator()(real lon, real& Bx, real& By, real& Bz) const throw() {
real dummy; real dummy;
Field(lon, false, Bx, By, Bz, dummy, dummy, dummy); Field(lon, false, Bx, By, Bz, dummy, dummy, dummy);
} }
/** /**
* Evaluate the components of the geomagnetic field and their time * Evaluate the components of the geomagnetic field and their time
* derivatives at a particular longitude. * derivatives at a particular longitude.
* *
* @param[in] lon longitude of the point (degrees). * @param[in] lon longitude of the point (degrees).
* @param[out] Bx the easterly component of the magnetic field (nanotes la). * @param[out] Bx the easterly component of the magnetic field (nanotes la).
* @param[out] By the northerly component of the magnetic field (nanote sla). * @param[out] By the northerly component of the magnetic field (nanote sla).
* @param[out] Bz the vertical (up) component of the magnetic field * @param[out] Bz the vertical (up) component of the magnetic field
* (nanotesla). * (nanotesla).
* @param[out] Bxt the rate of change of \e Bx (nT/yr). * @param[out] Bxt the rate of change of \e Bx (nT/yr).
* @param[out] Byt the rate of change of \e By (nT/yr). * @param[out] Byt the rate of change of \e By (nT/yr).
* @param[out] Bzt the rate of change of \e Bz (nT/yr). * @param[out] Bzt the rate of change of \e Bz (nT/yr).
********************************************************************** / ********************************************************************** /
void operator()(real lon, real& Bx, real& By, real& Bz, void operator()(real lon, real& Bx, real& By, real& Bz,
real& Bxt, real& Byt, real& Bzt) const { real& Bxt, real& Byt, real& Bzt) const throw() {
Field(lon, true, Bx, By, Bz, Bxt, Byt, Bzt); Field(lon, true, Bx, By, Bz, Bxt, Byt, Bzt);
} }
///@}
/** \name Inspector functions
**********************************************************************
/
///@{
/**
* @return true if the object has been initialized.
**********************************************************************
/
bool Init() const throw() { return _a > 0; }
/**
* @return \e a the equatorial radius of the ellipsoid (meters). This
is
* the value inherited from the MagneticModel object used in the
* constructor.
**********************************************************************
/
Math::real MajorRadius() const throw()
{ return Init() ? _a : Math::NaN<real>(); }
/**
* @return \e f the flattening of the ellipsoid. This is the value
* inherited from the MagneticModel object used in the constructor.
**********************************************************************
/
Math::real Flattening() const throw()
{ return Init() ? _f : Math::NaN<real>(); }
/**
* @return the latitude of the circle (degrees).
**********************************************************************
/
Math::real Latitude() const throw()
{ return Init() ? _lat : Math::NaN<real>(); }
/**
* @return the height of the circle (meters).
**********************************************************************
/
Math::real Height() const throw()
{ return Init() ? _h : Math::NaN<real>(); }
/**
* @return the time (fractional years).
**********************************************************************
/
Math::real Time() const throw()
{ return Init() ? _t : Math::NaN<real>(); }
///@}
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif // GEOGRAPHICLIB_MAGNETICCIRCLE_HPP #endif // GEOGRAPHICLIB_MAGNETICCIRCLE_HPP
 End of changes. 11 change blocks. 
10 lines changed or deleted 74 lines changed or added


 MagneticModel.hpp   MagneticModel.hpp 
/** /**
* \file MagneticModel.hpp * \file MagneticModel.hpp
* \brief Header for GeographicLib::MagneticModel class * \brief Header for GeographicLib::MagneticModel class
* *
* Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed un der * Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed un der
* the MIT/X11 License. For more information, see * the MIT/X11 License. For more information, see
* http://geographiclib.sourceforge.net/ * http://geographiclib.sourceforge.net/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_MAGNETICMODEL_HPP) #if !defined(GEOGRAPHICLIB_MAGNETICMODEL_HPP)
#define GEOGRAPHICLIB_MAGNETICMODEL_HPP "$Id: fc3e51862699ad9f6cf56a3aa2c35 #define GEOGRAPHICLIB_MAGNETICMODEL_HPP \
883f0103820 $" "$Id: 4da5b2c9d2a2ced31182528eb47b129520ae484b $"
#include <string> #include <string>
#include <sstream> #include <sstream>
#include <vector> #include <vector>
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
#include <GeographicLib/Geocentric.hpp> #include <GeographicLib/Geocentric.hpp>
#include <GeographicLib/SphericalHarmonic.hpp> #include <GeographicLib/SphericalHarmonic.hpp>
#include <GeographicLib/SphericalHarmonic1.hpp>
#if defined(_MSC_VER) #if defined(_MSC_VER)
// Squelch warnings about dll vs vector // Squelch warnings about dll vs vector
#pragma warning (push) #pragma warning (push)
#pragma warning (disable: 4251) #pragma warning (disable: 4251)
#endif #endif
namespace GeographicLib { namespace GeographicLib {
class MagneticCircle; class MagneticCircle;
skipping to change at line 74 skipping to change at line 74
int _Nmodels; int _Nmodels;
SphericalHarmonic::normalization _norm; SphericalHarmonic::normalization _norm;
Geocentric _earth; Geocentric _earth;
std::vector< std::vector<real> > _G; std::vector< std::vector<real> > _G;
std::vector< std::vector<real> > _H; std::vector< std::vector<real> > _H;
std::vector<SphericalHarmonic> _harm; std::vector<SphericalHarmonic> _harm;
void Field(real t, real lat, real lon, real h, bool diffp, void Field(real t, real lat, real lon, real h, bool diffp,
real& Bx, real& By, real& Bz, real& Bx, real& By, real& Bz,
real& Bxt, real& Byt, real& Bzt) const throw(); real& Bxt, real& Byt, real& Bzt) const throw();
void ReadMetadata(const std::string& name); void ReadMetadata(const std::string& name);
static bool ParseLine(const std::string& line, MagneticModel(const MagneticModel&); // copy constructor not allowed
std::string& key, std::string& val); MagneticModel& operator=(const MagneticModel&); // nor copy assignment
public: public:
/** \name Setting up the magnetic model /** \name Setting up the magnetic model
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* Construct a magnetic model. * Construct a magnetic model.
* *
* @param[in] name the name of the model. * @param[in] name the name of the model.
* @param[in] path (optional) directory for data file. * @param[in] path (optional) directory for data file.
* @param[in] earth (optional) Geocentric object for converting * @param[in] earth (optional) Geocentric object for converting
* coordinates; default Geocentric::WGS84. * coordinates; default Geocentric::WGS84.
* *
* A filename is formed by appending ".wwm" (World Magnetic Model) to * A filename is formed by appending ".wmm" (World Magnetic Model) to t
* the name. If \e path is specified (and is non-empty), then the file he
is * 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 * loaded from directory, \e path. Otherwise the path is given by the
* MAGNETIC_PATH environment variable. If that is undefined, a * DefaultMagneticPath(). This may throw an exception because the file
* compile-time default path is used * does not exist, is unreadable, or is corrupt.
* (/usr/local/share/GeographicLib/magnetic on non-Windows systems and
* C:/Documents and Settings/All Users/Application
* Data/GeographicLib/magnetic on Windows systems). This may throw an
* exception because the file does not exist, is unreadable, or is corr
upt.
* *
* This file contains the metadata which specifies the properties of th e * This file contains the metadata which specifies the properties of th e
* model. The coefficients for the spherical harmonic sums are obtaine d * model. The coefficients for the spherical harmonic sums are obtaine d
* from a file obtained by appending ".cof" to metadata file (so the * from a file obtained by appending ".cof" to metadata file (so the
* filename ends in ".wwm.cof"). * filename ends in ".wwm.cof").
* *
* The model is not tied to a particular ellipsoidal model of the earth . * The model is not tied to a particular ellipsoidal model of the earth .
* The final earth argument to the constructor specify an ellipsoid to * The final earth argument to the constructor specify an ellipsoid to
* allow geodetic coordinates to the transformed into the spherical * allow geodetic coordinates to the transformed into the spherical
* coordinates used in the spherical harmonic sum. * coordinates used in the spherical harmonic sum.
********************************************************************** / ********************************************************************** /
MagneticModel(const std::string& name, const std::string& path = "", explicit MagneticModel(const std::string& name,
const Geocentric& earth = Geocentric::WGS84); const std::string& path = "",
const Geocentric& earth = Geocentric::WGS84);
///@} ///@}
/** \name Compute the magnetic field /** \name Compute the magnetic field
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* Evaluate the components of the geomagnetic field. * Evaluate the components of the geomagnetic field.
* *
* @param[in] t the time (years). * @param[in] t the time (years).
* @param[in] lat latitude of the point (degrees). * @param[in] lat latitude of the point (degrees).
skipping to change at line 159 skipping to change at line 156
********************************************************************** / ********************************************************************** /
void operator()(real t, real lat, real lon, real h, void operator()(real t, real lat, real lon, real h,
real& Bx, real& By, real& Bz, real& Bx, real& By, real& Bz,
real& Bxt, real& Byt, real& Bzt) const throw() { real& Bxt, real& Byt, real& Bzt) const throw() {
Field(t, lat, lon, h, true, Bx, By, Bz, Bxt, Byt, Bzt); Field(t, lat, lon, h, true, Bx, By, Bz, Bxt, Byt, Bzt);
} }
/** /**
* Create a MagneticCircle object to allow the geomagnetic field at man y * Create a MagneticCircle object to allow the geomagnetic field at man y
* points with constant \e lat, \e h, and \e t and varying \e lon to be * points with constant \e lat, \e h, and \e t and varying \e lon to be
* computed efficienty. * computed efficiently.
* *
* @param[in] t the time (years).
* @param[in] lat latitude of the point (degrees). * @param[in] lat latitude of the point (degrees).
* @param[in] h the height of the point above the ellipsoid (meters). * @param[in] h the height of the point above the ellipsoid (meters).
* @param[in] t the time (years).
* @return a MagneticCircle object whose MagneticCircle::operator()(rea l * @return a MagneticCircle object whose MagneticCircle::operator()(rea l
* lon) member function computes the field at a particular \e lon. * lon) member function computes the field at particular values of \e
* lon.
* *
* If the field at several points on a circle of latitude need to be * If the field at several points on a circle of latitude need to be
* calculated then instead of * calculated then instead of
\code \code
SphericalModel m(...); // Create a magnetic model MagneticModel m(...); // Create a magnetic model
double lat = 33, lon0 = 44, dlon = 0.01, h = 1000, t = 2012; double lat = 33, lon0 = 44, dlon = 0.01, h = 1000, t = 2012;
for (int i = 0; i <= 100; ++i) { for (int i = 0; i <= 100; ++i) {
real real
lon = lon0 + i * dlon, Bx, By, Bz; lon = lon0 + i * dlon,
m(lat, lon, h, t, Bx, By, Bz); Bx, By, Bz;
m(t, lat, lon, h, Bx, By, Bz);
std::cout << lon << " " << Bx << " " << By << " " << Bz << "\n"; std::cout << lon << " " << Bx << " " << By << " " << Bz << "\n";
} }
\endcode \endcode
* use a MagneticCircle as in * use a MagneticCircle as in
\code \code
SphericalModel m(...); // Create a magnetic model MagneticModel m(...); // Create a magnetic model
double lat = 33, lon0 = 44, dlon = 0.01, h = 1000, t = 2012; double lat = 33, lon0 = 44, dlon = 0.01, h = 1000, t = 2012;
MagneticCircle c(m.Circle(lat, h, t)); // the MagneticCircle object MagneticCircle c(m.Circle(t, lat, h)); // the MagneticCircle object
for (int i = 0; i <= 100; ++i) { for (int i = 0; i <= 100; ++i) {
real real
lon = lon0 + i * dlon, Bx, By, Bz; lon = lon0 + i * dlon, Bx, By, Bz;
c(lon, Bx, By, Bz); c(lon, Bx, By, Bz);
std::cout << lon << " " << Bx << " " << By << " " << Bz << "\n"; std::cout << lon << " " << Bx << " " << By << " " << Bz << "\n";
} }
\endcode \endcode
* For high-degree models, this will be substantially faster. * For high-degree models, this will be substantially faster.
********************************************************************** / ********************************************************************** /
MagneticCircle Circle(real lat, real h, real t) const; MagneticCircle Circle(real t, real lat, real h) const;
/** /**
* Compute various quantities dependent on the magnetic field. * Compute various quantities dependent on the magnetic field.
* *
* @param[in] Bx the \e x (easterly) component of the magnetic field (n T). * @param[in] Bx the \e x (easterly) component of the magnetic field (n T).
* @param[in] By the \e y (northerly) component of the magnetic field ( nT). * @param[in] By the \e y (northerly) component of the magnetic field ( nT).
* @param[in] Bz the \e z (vertical, up positive) component of the magn etic * @param[in] Bz the \e z (vertical, up positive) component of the magn etic
* field (nT). * field (nT).
* @param[out] H the horizontal magnetic field (nT). * @param[out] H the horizontal magnetic field (nT).
* @param[out] F the total magnetic field (nT). * @param[out] F the total magnetic field (nT).
skipping to change at line 246 skipping to change at line 245
static void FieldComponents(real Bx, real By, real Bz, static void FieldComponents(real Bx, real By, real Bz,
real Bxt, real Byt, real Bzt, real Bxt, real Byt, real Bzt,
real& H, real& F, real& D, real& I, real& H, real& F, real& D, real& I,
real& Ht, real& Ft, real& Dt, real& It) thr ow(); real& Ht, real& Ft, real& Dt, real& It) thr ow();
///@} ///@}
/** \name Inspector functions /** \name Inspector functions
********************************************************************** / ********************************************************************** /
///@{ ///@{
/** /**
* @return the description of the magnetic model, if available, in the * @return the description of the magnetic model, if available, from th
data e
* file; if absent, return "NONE". * Description file 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 date of the model; if absent, return "UNKNOWN". * @return date of the model, if available, from the ReleaseDate field
in
* the data file; if absent, return "UNKNOWN".
********************************************************************** / ********************************************************************** /
const std::string& DateTime() const throw() { return _date; } const std::string& DateTime() const throw() { return _date; }
/** /**
* @return full file name used to load the magnetic model. * @return full file name used to load the magnetic model.
********************************************************************** / ********************************************************************** /
const std::string& MagneticFile() const throw() { return _filename; } const std::string& MagneticFile() const throw() { return _filename; }
/** /**
* @return "name" used to load the magnetic model (from the first argum ent * @return "name" used to load the magnetic model (from the first argum ent
skipping to change at line 279 skipping to change at line 279
********************************************************************** / ********************************************************************** /
const std::string& MagneticModelDirectory() const throw() { return _dir ; } const std::string& MagneticModelDirectory() const throw() { return _dir ; }
/** /**
* @return the minimum height above the ellipsoid (in meters) for which * @return the minimum height above the ellipsoid (in meters) for which
* this MagneticModel should be used. * this MagneticModel should be used.
* *
* Because the model will typically provide useful results * Because the model will typically provide useful results
* slightly outside the range of allowed heights, no check of \e t * slightly outside the range of allowed heights, no check of \e t
* argument is made by MagneticModel::operator()() or * argument is made by MagneticModel::operator()() or
* MagnetgicModel::Circle. * MagneticModel::Circle.
********************************************************************** / ********************************************************************** /
Math::real MinHeight() const throw() { return _hmin; } Math::real MinHeight() const throw() { return _hmin; }
/** /**
* @return the maximum height above the ellipsoid (in meters) for which * @return the maximum height above the ellipsoid (in meters) for which
* this MagneticModel should be used. * this MagneticModel should be used.
* *
* Because the model will typically provide useful results * Because the model will typically provide useful results
* slightly outside the range of allowed heights, no check of \e t * slightly outside the range of allowed heights, no check of \e t
* argument is made by MagneticModel::operator()() or * argument is made by MagneticModel::operator()() or
* MagnetgicModel::Circle. * MagneticModel::Circle.
********************************************************************** / ********************************************************************** /
Math::real MaxHeight() const throw() { return _hmax; } Math::real MaxHeight() const throw() { return _hmax; }
/** /**
* @return the minimum time (in years) for which this MagneticModel sho uld * @return the minimum time (in years) for which this MagneticModel sho uld
* be used. * be used.
* *
* Because the model will typically provide useful results * Because the model will typically provide useful results
* slightly outside the range of allowed times, no check of \e t * slightly outside the range of allowed times, no check of \e t
* argument is made by MagneticModel::operator()() or * argument is made by MagneticModel::operator()() or
* MagnetgicModel::Circle. * MagneticModel::Circle.
********************************************************************** / ********************************************************************** /
Math::real MinTime() const throw() { return _tmin; } Math::real MinTime() const throw() { return _tmin; }
/** /**
* @return the maximum time (in years) for which this MagneticModel sho uld * @return the maximum time (in years) for which this MagneticModel sho uld
* be used. * be used.
* *
* Because the model will typically provide useful results * Because the model will typically provide useful results
* slightly outside the range of allowed times, no check of \e t * slightly outside the range of allowed times, no check of \e t
* argument is made by MagneticModel::operator()() or * argument is made by MagneticModel::operator()() or
* MagnetgicModel::Circle. * MagneticModel::Circle.
********************************************************************** / ********************************************************************** /
Math::real MaxTime() const throw() { return _tmax; } Math::real MaxTime() const throw() { return _tmax; }
/** /**
* @return \e a the equatorial radius of the ellipsoid (meters). This is * @return \e a the equatorial radius of the ellipsoid (meters). This is
* the value of \e a inherited from the Geocentric object used in the * the value of \e a inherited from the Geocentric object used in the
* constructor. * constructor.
********************************************************************** / ********************************************************************** /
Math::real MajorRadius() const throw() { return _earth.MajorRadius(); } Math::real MajorRadius() const throw() { return _earth.MajorRadius(); }
/** /**
* @return \e f the flattening of the ellipsoid. This is the value * @return \e f the flattening of the ellipsoid. This is the value
* inherited from the Geocentric object used in the constructor. * inherited from the Geocentric object used in the constructor.
********************************************************************** / ********************************************************************** /
Math::real Flattening() const throw() { return _earth.Flattening(); } Math::real Flattening() const throw() { return _earth.Flattening(); }
///@} ///@}
/** /**
* @return the default path for magnetic model data files. * @return the default path for magnetic model data files.
* *
* This is the value of the environment variable MAGNETIC_PATH, if set, * This is the value of the environment variable MAGNETIC_PATH, if set;
* otherwise, it is a compile-time default. * otherwise, it is $GEOGRAPHICLIB_DATA/magnetic if the environment
* variable GEOGRAPHICLIB_DATA is set; otherwise, it is a compile-time
* default (/usr/local/share/GeographicLib/magnetic on non-Windows syst
ems
* and C:/Documents and Settings/All Users/Application
* Data/GeographicLib/magnetic on Windows systems).
********************************************************************** / ********************************************************************** /
static std::string DefaultMagneticPath(); static std::string DefaultMagneticPath();
/** /**
* @return the default name for the magnetic model. * @return the default name for the magnetic model.
* *
* This is the value of the environment variable MAGNETIC_NAME, if set, * This is the value of the environment variable MAGNETIC_NAME, if set,
* otherwise, it is "wmm2010-12". The MagneticModel class does not use * otherwise, it is "wmm2010". The MagneticModel class does not use th
* this function; it is just provided as a convenience for a calling is
* program when constructing a MagneticModel object. * function; it is just provided as a convenience for a calling program
* when constructing a MagneticModel object.
********************************************************************** / ********************************************************************** /
static std::string DefaultMagneticName(); static std::string DefaultMagneticName();
}; };
} // namespace GeographicLib } // namespace GeographicLib
#if defined(_MSC_VER) #if defined(_MSC_VER)
#pragma warning (pop) #pragma warning (pop)
#endif #endif
 End of changes. 24 change blocks. 
40 lines changed or deleted 44 lines changed or added


 Math.hpp   Math.hpp 
skipping to change at line 15 skipping to change at line 15
* Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co m> * Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co m>
* and licensed under the MIT/X11 License. For more information, see * and licensed under the MIT/X11 License. For more information, see
* http://geographiclib.sourceforge.net/ * http://geographiclib.sourceforge.net/
**********************************************************************/ **********************************************************************/
// Constants.hpp includes Math.hpp. Place this include outside Math.hpp's // Constants.hpp includes Math.hpp. Place this include outside Math.hpp's
// include guard to enforce this ordering. // include guard to enforce this ordering.
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
#if !defined(GEOGRAPHICLIB_MATH_HPP) #if !defined(GEOGRAPHICLIB_MATH_HPP)
#define GEOGRAPHICLIB_MATH_HPP "$Id: dd6afb3ee17e4c12342f7691e714c5f6dc78f0 12 $" #define GEOGRAPHICLIB_MATH_HPP "$Id: 0bc5689770853cdf3c6dad20eb38980624ffcd 61 $"
/** /**
* Are C++0X math functions available? * Are C++0X math functions available?
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_CPLUSPLUS0X_MATH) #if !defined(GEOGRAPHICLIB_CPLUSPLUS0X_MATH)
# if defined(__GXX_EXPERIMENTAL_CXX0X__) # if defined(__GXX_EXPERIMENTAL_CXX0X__)
# define GEOGRAPHICLIB_CPLUSPLUS0X_MATH 1 # define GEOGRAPHICLIB_CPLUSPLUS0X_MATH 1
# else # else
# define GEOGRAPHICLIB_CPLUSPLUS0X_MATH 0 # define GEOGRAPHICLIB_CPLUSPLUS0X_MATH 0
# endif # endif
skipping to change at line 53 skipping to change at line 53
#include <cmath> #include <cmath>
#include <limits> #include <limits>
#include <algorithm> #include <algorithm>
#include <vector> #include <vector>
/** /**
* \brief Namespace for %GeographicLib * \brief Namespace for %GeographicLib
* *
* All of %GeographicLib is defined within the GeographicLib namespace. In * All of %GeographicLib is defined within the GeographicLib namespace. In
* addtion all the header files are included via %GeographicLib/filename. * addtiion all the header files are included via %GeographicLib/filename.
This * This minimizes the likelihood of conflicts with other packages.
* minimizes the likelihood of conflicts with other packages.
**********************************************************************/ **********************************************************************/
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief Mathematical functions needed by %GeographicLib * \brief Mathematical functions needed by %GeographicLib
* *
* Define mathematical functions in order to localize system dependencies and * Define mathematical functions in order to localize system dependencies and
* to provide generic versions of the functions. In addition define a re al * to provide generic versions of the functions. In addition define a re al
* type to be used by %GeographicLib. * type to be used by %GeographicLib.
**********************************************************************/ **********************************************************************/
skipping to change at line 109 skipping to change at line 109
/** /**
* true if the machine is big-endian * true if the machine is big-endian
********************************************************************** / ********************************************************************** /
static const bool bigendian = WORDS_BIGENDIAN; static const bool bigendian = WORDS_BIGENDIAN;
/** /**
* @tparam T the type of the returned value. * @tparam T the type of the returned value.
* @return \e pi. * @return \e pi.
********************************************************************** / ********************************************************************** /
template<typename T> template<typename T> static inline T pi() throw()
static inline T pi() throw() { return std::atan2(T(0), -T(1)); } { return std::atan2(T(0), -T(1)); }
/** /**
* A synonym for pi<real>(). * A synonym for pi<real>().
********************************************************************** / ********************************************************************** /
static inline real pi() throw() { return pi<real>(); } static inline real pi() throw() { return pi<real>(); }
/** /**
* @tparam T the type of the returned value. * @tparam T the type of the returned value.
* @return the number of radians in a degree. * @return the number of radians in a degree.
********************************************************************** / ********************************************************************** /
template<typename T> template<typename T> static inline T degree() throw()
static inline T degree() throw() { return pi<T>() / T(180); } { return pi<T>() / T(180); }
/** /**
* A synonym for degree<real>(). * A synonym for degree<real>().
********************************************************************** / ********************************************************************** /
static inline real degree() throw() { return degree<real>(); } static inline real degree() throw() { return degree<real>(); }
/** /**
* Square a number. * Square a number.
* @tparam T the type of the argument and the returned value. * @tparam T the type of the argument and the returned value.
* @param[in] x * @param[in] x
* @return \e x<sup>2</sup>. * @return <i>x</i><sup>2</sup>.
********************************************************************** / ********************************************************************** /
template<typename T> template<typename T> static inline T sq(T x) throw()
static inline T sq(T x) throw() { return x * x; } { return x * x; }
#if defined(DOXYGEN) #if defined(DOXYGEN)
/** /**
* The hypotenuse function avoiding underflow and overflow. * The hypotenuse function avoiding underflow and overflow.
* *
* @tparam T the type of the arguments and the returned value. * @tparam T the type of the arguments and the returned value.
* @param[in] x * @param[in] x
* @param[in] y * @param[in] y
* @return sqrt(\e x<sup>2</sup> + \e y<sup>2</sup>). * @return sqrt(<i>x</i><sup>2</sup> + <i>y</i><sup>2</sup>).
********************************************************************** / ********************************************************************** /
template<typename T> template<typename T> static inline T hypot(T x, T y) throw() {
static inline T hypot(T x, T y) throw() {
x = std::abs(x); x = std::abs(x);
y = std::abs(y); y = std::abs(y);
T a = (std::max)(x, y), T a = (std::max)(x, y),
b = (std::min)(x, y) / (a ? a : 1); b = (std::min)(x, y) / (a ? a : 1);
return a * std::sqrt(1 + b * b); return a * std::sqrt(1 + b * b);
} }
#elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH
template<typename T> template<typename T> static inline T hypot(T x, T y) throw()
static inline T hypot(T x, T y) throw() { return std::hypot(x, y); } { return std::hypot(x, y); }
#elif defined(_MSC_VER) #elif defined(_MSC_VER)
static inline double hypot(double x, double y) throw() static inline double hypot(double x, double y) throw()
{ return _hypot(x, y); } { return _hypot(x, y); }
#if (_MSC_VER < 1400) #if (_MSC_VER < 1400)
// Visual C++ 7.1/VS .NET 2003 does not have _hypotf() // Visual C++ 7.1/VS .NET 2003 does not have _hypotf()
static inline float hypot(float x, float y) throw() static inline float hypot(float x, float y) throw()
{ return float(_hypot(x, y)); } { return float(_hypot(x, y)); }
#else #else
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); }
skipping to change at line 194 skipping to change at line 193
#if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MA TH) #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MA TH)
/** /**
* exp(\e x) - 1 accurate near \e x = 0. This is taken from * exp(\e x) - 1 accurate near \e x = 0. This is taken from
* N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd * N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd
* Edition (SIAM, 2002), Sec 1.14.1, p 19. * Edition (SIAM, 2002), Sec 1.14.1, p 19.
* *
* @tparam T the type of the argument and the returned value. * @tparam T the type of the argument and the returned value.
* @param[in] x * @param[in] x
* @return exp(\e x) - 1. * @return exp(\e x) - 1.
********************************************************************** / ********************************************************************** /
template<typename T> template<typename T> static inline T expm1(T x) throw() {
static inline T expm1(T x) throw() {
volatile T volatile T
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));
} }
#elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH
template<typename T> template<typename T> static inline T expm1(T x) throw()
static inline T expm1(T x) throw() { return std::expm1(x); } { return std::expm1(x); }
#else #else
static inline double expm1(double x) throw() { return ::expm1(x); } static inline double expm1(double x) throw() { return ::expm1(x); }
static inline float expm1(float x) throw() { return ::expm1f(x); } static inline float expm1(float x) throw() { return ::expm1f(x); }
#if defined(HAVE_LONG_DOUBLE) #if defined(HAVE_LONG_DOUBLE)
static inline long double expm1(long double x) throw() static inline long double expm1(long double x) throw()
{ return ::expm1l(x); } { return ::expm1l(x); }
#endif #endif
#endif #endif
#if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MA TH) #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MA TH)
skipping to change at line 230 skipping to change at line 228
* *
* This is taken from D. Goldberg, * This is taken from D. Goldberg,
* <a href="http://dx.doi.org/10.1145/103162.103163">What every compute r * <a href="http://dx.doi.org/10.1145/103162.103163">What every compute r
* scientist should know about floating-point arithmetic</a> (1991), * scientist should know about floating-point arithmetic</a> (1991),
* Theorem 4. See also, Higham (op. cit.), Answer to Problem 1.5, p 52 8. * Theorem 4. See also, Higham (op. cit.), Answer to Problem 1.5, p 52 8.
* *
* @tparam T the type of the argument and the returned value. * @tparam T the type of the argument and the returned value.
* @param[in] x * @param[in] x
* @return log(1 + \e x). * @return log(1 + \e x).
********************************************************************** / ********************************************************************** /
template<typename T> template<typename T> static inline T log1p(T x) throw() {
static inline T log1p(T x) throw() {
volatile T volatile T
y = 1 + x, y = 1 + x,
z = y - 1; z = y - 1;
// Here's the explanation for this magic: y = 1 + z, exactly, and z // Here's the explanation for this magic: y = 1 + z, exactly, and z
// 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;
} }
#elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH
template<typename T> template<typename T> static inline T log1p(T x) throw()
static inline T log1p(T x) throw() { return std::log1p(x); } { return std::log1p(x); }
#else #else
static inline double log1p(double x) throw() { return ::log1p(x); } static inline double log1p(double x) throw() { return ::log1p(x); }
static inline float log1p(float x) throw() { return ::log1pf(x); } static inline float log1p(float x) throw() { return ::log1pf(x); }
#if defined(HAVE_LONG_DOUBLE) #if defined(HAVE_LONG_DOUBLE)
static inline long double log1p(long double x) throw() static inline long double log1p(long double x) throw()
{ return ::log1pl(x); } { return ::log1pl(x); }
#endif #endif
#endif #endif
#if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MA TH) #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MA TH)
/** /**
* The inverse hyperbolic sine function. This is defined in terms of * The inverse hyperbolic sine function. This is defined in terms of
* Math::log1p(\e x) in order to maintain accuracy near \e x = 0. In * Math::log1p(\e x) in order to maintain accuracy near \e x = 0. In
* addition, the odd parity of the function is enforced. * addition, the odd parity of the function is enforced.
* *
* @tparam T the type of the argument and the returned value. * @tparam T the type of the argument and the returned value.
* @param[in] x * @param[in] x
* @return asinh(\e x). * @return asinh(\e x).
********************************************************************** / ********************************************************************** /
template<typename T> template<typename T> static inline T asinh(T x) throw() {
static inline T asinh(T x) throw() {
T y = std::abs(x); // Enforce odd parity T y = std::abs(x); // Enforce odd parity
y = log1p(y * (1 + y/(hypot(T(1), y) + 1))); y = log1p(y * (1 + y/(hypot(T(1), y) + 1)));
return x < 0 ? -y : y; return x < 0 ? -y : y;
} }
#elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH
template<typename T> template<typename T> static inline T asinh(T x) throw()
static inline T asinh(T x) throw() { return std::asinh(x); } { return std::asinh(x); }
#else #else
static inline double asinh(double x) throw() { return ::asinh(x); } static inline double asinh(double x) throw() { return ::asinh(x); }
static inline float asinh(float x) throw() { return ::asinhf(x); } static inline float asinh(float x) throw() { return ::asinhf(x); }
#if defined(HAVE_LONG_DOUBLE) #if defined(HAVE_LONG_DOUBLE)
static inline long double asinh(long double x) throw() static inline long double asinh(long double x) throw()
{ return ::asinhl(x); } { return ::asinhl(x); }
#endif #endif
#endif #endif
#if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MA TH) #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MA TH)
/** /**
* The inverse hyperbolic tangent function. This is defined in terms o f * The inverse hyperbolic tangent function. This is defined in terms o f
* Math::log1p(\e x) in order to maintain accuracy near \e x = 0. In * Math::log1p(\e x) in order to maintain accuracy near \e x = 0. In
* addition, the odd parity of the function is enforced. * addition, the odd parity of the function is enforced.
* *
* @tparam T the type of the argument and the returned value. * @tparam T the type of the argument and the returned value.
* @param[in] x * @param[in] x
* @return atanh(\e x). * @return atanh(\e x).
********************************************************************** / ********************************************************************** /
template<typename T> template<typename T> static inline T atanh(T x) throw() {
static inline T atanh(T x) throw() {
T y = std::abs(x); // Enforce odd parity T y = std::abs(x); // Enforce odd parity
y = log1p(2 * y/(1 - y))/2; y = log1p(2 * y/(1 - y))/2;
return x < 0 ? -y : y; return x < 0 ? -y : y;
} }
#elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH
template<typename T> template<typename T> static inline T atanh(T x) throw()
static inline T atanh(T x) throw() { return std::atanh(x); } { return std::atanh(x); }
#else #else
static inline double atanh(double x) throw() { return ::atanh(x); } static inline double atanh(double x) throw() { return ::atanh(x); }
static inline float atanh(float x) throw() { return ::atanhf(x); } static inline float atanh(float x) throw() { return ::atanhf(x); }
#if defined(HAVE_LONG_DOUBLE) #if defined(HAVE_LONG_DOUBLE)
static inline long double atanh(long double x) throw() static inline long double atanh(long double x) throw()
{ return ::atanhl(x); } { return ::atanhl(x); }
#endif #endif
#endif #endif
#if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MA TH) #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MA TH)
/** /**
* The cube root function. * The cube root function.
* *
* @tparam T the type of the argument and the returned value. * @tparam T the type of the argument and the returned value.
* @param[in] x * @param[in] x
* @return the real cube root of \e x. * @return the real cube root of \e x.
********************************************************************** / ********************************************************************** /
template<typename T> template<typename T> static inline T cbrt(T x) throw() {
static inline T cbrt(T x) throw() {
T y = std::pow(std::abs(x), 1/T(3)); // Return the real cube root T y = std::pow(std::abs(x), 1/T(3)); // Return the real cube root
return x < 0 ? -y : y; return x < 0 ? -y : y;
} }
#elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH #elif GEOGRAPHICLIB_CPLUSPLUS0X_MATH
template<typename T> template<typename T> static inline T cbrt(T x) throw()
static inline T cbrt(T x) throw() { return std::cbrt(x); } { return std::cbrt(x); }
#else #else
static inline double cbrt(double x) throw() { return ::cbrt(x); } static inline double cbrt(double x) throw() { return ::cbrt(x); }
static inline float cbrt(float x) throw() { return ::cbrtf(x); } static inline float cbrt(float x) throw() { return ::cbrtf(x); }
#if defined(HAVE_LONG_DOUBLE) #if defined(HAVE_LONG_DOUBLE)
static inline long double cbrt(long double x) throw() { return ::cbrtl( x); } static inline long double cbrt(long double x) throw() { return ::cbrtl( x); }
#endif #endif
#endif #endif
/** /**
* Test for finiteness. * Test for finiteness.
* *
* @tparam T the type of the argument. * @tparam T the type of the argument.
* @param[in] x * @param[in] x
* @return true if number is finite, false if NaN or infinite. * @return true if number is finite, false if NaN or infinite.
********************************************************************** / ********************************************************************** /
template<typename T> template<typename T> static inline bool isfinite(T x) throw() {
static inline bool isfinite(T x) throw() {
#if defined(DOXYGEN) #if defined(DOXYGEN)
return std::abs(x) <= (std::numeric_limits<T>::max)(); return std::abs(x) <= (std::numeric_limits<T>::max)();
#elif (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MATH) #elif (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MATH)
return _finite(x) != 0; return _finite(x) != 0;
#else #else
return std::isfinite(x); return std::isfinite(x);
#endif #endif
} }
/** /**
* The NaN (not a number) * The NaN (not a number)
* *
* @tparam T the type of the returned value. * @tparam T the type of the returned value.
* @return NaN if available, otherwise return the max real. * @return NaN if available, otherwise return the max real.
********************************************************************** / ********************************************************************** /
template<typename T> template<typename T> static inline T NaN() throw() {
static inline T NaN() throw() {
return std::numeric_limits<T>::has_quiet_NaN ? return std::numeric_limits<T>::has_quiet_NaN ?
std::numeric_limits<T>::quiet_NaN() : std::numeric_limits<T>::quiet_NaN() :
(std::numeric_limits<T>::max)(); (std::numeric_limits<T>::max)();
} }
/** /**
* A synonym for NaN<real>(). * A synonym for NaN<real>().
********************************************************************** / ********************************************************************** /
static inline real NaN() throw() { return NaN<real>(); } static inline real NaN() throw() { return NaN<real>(); }
/** /**
* Test for NaN. * Test for NaN.
* *
* @tparam T the type of the argument. * @tparam T the type of the argument.
* @param[in] x * @param[in] x
* @return true if argument is a NaN. * @return true if argument is a NaN.
********************************************************************** / ********************************************************************** /
template<typename T> template<typename T> static inline bool isnan(T x) throw() {
static inline bool isnan(T x) throw() {
#if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MA TH) #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS0X_MA TH)
return x != x; return x != x;
#else #else
return std::isnan(x); return std::isnan(x);
#endif #endif
} }
/** /**
* Infinity * Infinity
* *
* @tparam T the type of the returned value. * @tparam T the type of the returned value.
* @return infinity if available, otherwise return the max real. * @return infinity if available, otherwise return the max real.
********************************************************************** / ********************************************************************** /
template<typename T> template<typename T> static inline T infinity() throw() {
static inline T infinity() throw() {
return std::numeric_limits<T>::has_infinity ? return std::numeric_limits<T>::has_infinity ?
std::numeric_limits<T>::infinity() : std::numeric_limits<T>::infinity() :
(std::numeric_limits<T>::max)(); (std::numeric_limits<T>::max)();
} }
/** /**
* A synonym for infinity<real>(). * A synonym for infinity<real>().
********************************************************************** / ********************************************************************** /
static inline real infinity() throw() { return infinity<real>(); } static inline real infinity() throw() { return infinity<real>(); }
/** /**
 End of changes. 23 change blocks. 
44 lines changed or deleted 33 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, 2011) <charles@karney.co m> * Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co m>
* and licensed under the MIT/X11 License. For more information, see * and licensed under the MIT/X11 License. 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: babda3abc04257e799b84de8 #define GEOGRAPHICLIB_POLARSTEREOGRAPHIC_HPP \
aaaec20082b960e0 $" "$Id: 53979440f6d23cf95aefe41a2b9ccb8eedf5b811 $"
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief Polar Stereographic Projection * \brief Polar Stereographic Projection
* *
* Implementation taken from the report, * Implementation taken from the report,
* - J. P. Snyder, * - J. P. Snyder,
skipping to change at line 38 skipping to change at line 39
**********************************************************************/ **********************************************************************/
class GEOGRAPHIC_EXPORT PolarStereographic { class GEOGRAPHIC_EXPORT PolarStereographic {
private: private:
typedef Math::real real; typedef Math::real real;
// _Cx used to be _C but g++ 3.4 has a macro of that name // _Cx used to be _C but g++ 3.4 has a macro of that name
real _a, _f, _e2, _e, _e2m, _Cx, _c; real _a, _f, _e2, _e, _e2m, _Cx, _c;
real _k0; real _k0;
static const real tol_; static const real tol_;
static const real overflow_; static const real overflow_;
static const int numit_ = 5; static const int numit_ = 5;
// tan(x) for x in [-pi/2, pi/2] ensuring that the sign is right
static inline real tanx(real x) throw() {
real t = std::tan(x);
// Write the tests this way to ensure that tanx(NaN()) is NaN()
return x >= 0 ? (!(t < 0) ? t : overflow_) : (!(t >= 0) ? t : -overfl
ow_);
}
// Return e * atanh(e * x) for f >= 0, else return // Return e * atanh(e * x) for f >= 0, else return
// - sqrt(-e2) * atan( sqrt(-e2) * x) for f < 0 // - sqrt(-e2) * atan( sqrt(-e2) * x) for f < 0
inline real eatanhe(real x) const throw() { inline real eatanhe(real x) const throw() {
return _f >= 0 ? _e * Math::atanh(_e * x) : - _e * std::atan(_e * x); return _f >= 0 ? _e * Math::atanh(_e * x) : - _e * std::atan(_e * x);
} }
public: public:
/** /**
* Constructor for a ellipsoid with * Constructor for a ellipsoid with
* *
* @param[in] a equatorial radius (meters) * @param[in] a equatorial radius (meters).
* @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphe re. * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphe re.
* Negative \e f gives a prolate ellipsoid. If \e f > 1, set flatten ing * Negative \e f gives a prolate ellipsoid. If \e f > 1, set flatten ing
* to 1/\e f. * to 1/\e f.
* @param[in] k0 central scale factor. * @param[in] k0 central scale factor.
* *
* An exception is thrown if either of the axes of the ellipsoid is * An exception is thrown if either of the axes of the ellipsoid is
* not positive \e a or if \e k0 is not positive. * not positive \e a or if \e k0 is not positive.
********************************************************************** / ********************************************************************** /
PolarStereographic(real a, real f, real k0); PolarStereographic(real a, real f, real k0);
 End of changes. 3 change blocks. 
3 lines changed or deleted 10 lines changed or added


 PolygonArea.hpp   PolygonArea.hpp 
/** /**
* \file PolygonArea.hpp * \file PolygonArea.hpp
* \brief Header for GeographicLib::PolygonArea class * \brief Header for GeographicLib::PolygonArea class
* *
* Copyright (c) Charles Karney (2010, 2011) <charles@karney.com> and licen sed * Copyright (c) Charles Karney (2010, 2011) <charles@karney.com> and licen sed
* under the MIT/X11 License. For more information, see * under the MIT/X11 License. For more information, see
* http://geographiclib.sourceforge.net/ * http://geographiclib.sourceforge.net/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_POLYGONAREA_HPP) #if !defined(GEOGRAPHICLIB_POLYGONAREA_HPP)
#define GEOGRAPHICLIB_POLYGONAREA_HPP "$Id: 705d561a35fc61675a77d0d54d938c7 #define GEOGRAPHICLIB_POLYGONAREA_HPP \
690763bcf $" "$Id: c6b922a70d4f5ea75580f4e2746a2148b288a03e $"
#include <GeographicLib/Geodesic.hpp> #include <GeographicLib/Geodesic.hpp>
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
#include <GeographicLib/Accumulator.hpp> #include <GeographicLib/Accumulator.hpp>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief Polygon Areas. * \brief Polygon Areas.
* *
 End of changes. 1 change blocks. 
2 lines changed or deleted 2 lines changed or added


 SphericalEngine.hpp   SphericalEngine.hpp 
/** /**
* \file SphericalEngine.hpp * \file SphericalEngine.hpp
* \brief Header for GeographicLib::SphericalEngine class * \brief Header for GeographicLib::SphericalEngine class
* *
* Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed un der * Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed un der
* the MIT/X11 License. For more information, see * the MIT/X11 License. For more information, see
* http://geographiclib.sourceforge.net/ * http://geographiclib.sourceforge.net/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_SPHERICALENGINE_HPP) #if !defined(GEOGRAPHICLIB_SPHERICALENGINE_HPP)
#define GEOGRAPHICLIB_SPHERICALENGINE_HPP "$Id: 08b65902abf634bdf1a9d6685b0 #define GEOGRAPHICLIB_SPHERICALENGINE_HPP \
065e762190214 $" "$Id: 3410576dbc0276f23077ba3662817e14981ab919 $"
#include <vector> #include <vector>
#include <istream>
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
#if defined(_MSC_VER) #if defined(_MSC_VER)
// Squelch warnings about dll vs vector // Squelch warnings about dll vs vector
#pragma warning (push) #pragma warning (push)
#pragma warning (disable: 4251) #pragma warning (disable: 4251)
#endif #endif
namespace GeographicLib { namespace GeographicLib {
class CircularEngine; class CircularEngine;
/** /**
* \brief The evaluation engine for SphericalHarmonic * \brief The evaluation engine for SphericalHarmonic
* *
* This serves as the backend to SphericalHarmonic, SphericalHarmonic1, a nd * This serves as the backend to SphericalHarmonic, SphericalHarmonic1, a nd
* SphericalHarmonic2. Typically end-users will not have to access this * SphericalHarmonic2. Typically end-users will not have to access this
* class directly. * class directly.
*
* See SphericalEngine.cpp for more information on the implementation.
**********************************************************************/ **********************************************************************/
class GEOGRAPHIC_EXPORT SphericalEngine { class GEOGRAPHIC_EXPORT SphericalEngine {
private: private:
typedef Math::real real; typedef Math::real real;
// A table of the square roots of integers
static std::vector<real> root_;
friend class CircularEngine; // CircularEngine needs access to root_, s
cale_
// An internal scaling of the coefficients to avoid overflow in // An internal scaling of the coefficients to avoid overflow in
// intermediate calculations. // intermediate calculations.
static const real scale_; static const real scale_;
// Move latitudes near the pole off the axis by this amount. // Move latitudes near the pole off the axis by this amount.
static const real eps_; static const real eps_;
static const std::vector<real> Z_; static const std::vector<real> Z_;
SphericalEngine(); // Disable constructor SphericalEngine(); // Disable constructor
public: public:
/** /**
* Supported normalizations for associate Legendre polynomials. * Supported normalizations for associated Legendre polynomials.
********************************************************************** / ********************************************************************** /
enum normalization { enum normalization {
/** /**
* Fully normalized associated Legendre polynomials. See * Fully normalized associated Legendre polynomials. See
* SphericalHarmonic::full for documentation. * SphericalHarmonic::FULL for documentation.
* *
* @hideinitializer * @hideinitializer
******************************************************************** **/ ******************************************************************** **/
full = 0, FULL = 0,
/** /**
* Schmidt semi-normalized associated Legendre polynomials. See * Schmidt semi-normalized associated Legendre polynomials. See
* SphericalHarmonic::schmidt for documentation. * SphericalHarmonic::SCHMIDT for documentation.
* *
* @hideinitializer * @hideinitializer
******************************************************************** **/ ******************************************************************** **/
schmidt = 1, SCHMIDT = 1,
/// \cond SKIP
// These are deprecated...
full = FULL,
schmidt = SCHMIDT,
/// \endcond
}; };
/** /**
* \brief Package up coefficients for SphericalEngine * \brief Package up coefficients for SphericalEngine
* *
* This packages up the \e C, \e S coefficients and information about h ow * This packages up the \e C, \e S coefficients and information about h ow
* the coefficients are stored into a single structure. This allows a * the coefficients are stored into a single structure. This allows a
* vector of type SphericalEngine::coeff to be passed to * vector of type SphericalEngine::coeff to be passed to
* SphericalEngine::Value. This class also includes functions to aid * SphericalEngine::Value. This class also includes functions to aid
* indexing into \e C and \e S. * indexing into \e C and \e S.
skipping to change at line 120 skipping to change at line 131
, _nmx(nmx) , _nmx(nmx)
, _mmx(mmx) , _mmx(mmx)
, _Cnm(C.begin()) , _Cnm(C.begin())
, _Snm(S.begin() - (_N + 1)) , _Snm(S.begin() - (_N + 1))
{ {
if (!(_N >= _nmx && _nmx >= _mmx && _mmx >= -1)) if (!(_N >= _nmx && _nmx >= _mmx && _mmx >= -1))
throw GeographicErr("Bad indices for coeff"); throw GeographicErr("Bad indices for coeff");
if (!(index(_nmx, _mmx) < int(C.size()) && if (!(index(_nmx, _mmx) < int(C.size()) &&
index(_nmx, _mmx) < int(S.size()) + (_N + 1))) index(_nmx, _mmx) < int(S.size()) + (_N + 1)))
throw GeographicErr("Arrays too small in coeff"); throw GeographicErr("Arrays too small in coeff");
SphericalEngine::RootTable(_nmx);
} }
/** /**
* The constructor for full coefficient vectors. * The constructor for full coefficient vectors.
* *
* @param[in] C a vector of coefficients for the cosine terms. * @param[in] C a vector of coefficients for the cosine terms.
* @param[in] S a vector of coefficients for the sine terms. * @param[in] S a vector of coefficients for the sine terms.
* @param[in] N the maximum degree and order. * @param[in] N the maximum degree and order.
* *
* This requires \e N >= -1. \e C and \e S must also be large enough to * This requires \e N >= -1. \e C and \e S must also be large enough to
* hold the coefficients. Otherwise an exception is thrown. * hold the coefficients. Otherwise an exception is thrown.
skipping to change at line 145 skipping to change at line 157
, _nmx(N) , _nmx(N)
, _mmx(N) , _mmx(N)
, _Cnm(C.begin()) , _Cnm(C.begin())
, _Snm(S.begin() - (_N + 1)) , _Snm(S.begin() - (_N + 1))
{ {
if (!(_N >= -1)) if (!(_N >= -1))
throw GeographicErr("Bad indices for coeff"); throw GeographicErr("Bad indices for coeff");
if (!(index(_nmx, _mmx) < int(C.size()) && if (!(index(_nmx, _mmx) < int(C.size()) &&
index(_nmx, _mmx) < int(S.size()) + (_N + 1))) index(_nmx, _mmx) < int(S.size()) + (_N + 1)))
throw GeographicErr("Arrays too small in coeff"); throw GeographicErr("Arrays too small in coeff");
SphericalEngine::RootTable(_nmx);
} }
/** /**
* @return \e N the degree giving storage layout for \e C and \e S. * @return \e N the degree giving storage layout for \e C and \e S.
******************************************************************** **/ ******************************************************************** **/
inline int N() const throw() { return _N; } inline int N() const throw() { return _N; }
/** /**
* @return \e nmx the maximum degree to be used. * @return \e nmx the maximum degree to be used.
******************************************************************** **/ ******************************************************************** **/
inline int nmx() const throw() { return _nmx; } inline int nmx() const throw() { return _nmx; }
/** /**
skipping to change at line 227 skipping to change at line 240
/** /**
* The size of the coefficient vector for the sine terms. * The size of the coefficient vector for the sine terms.
* *
* @param[in] N the maximum degree. * @param[in] N the maximum degree.
* @param[in] M the maximum order. * @param[in] M the maximum order.
* @return the size of the vector of cosine terms as stored in column * @return the size of the vector of cosine terms as stored in column
* major order. * major order.
******************************************************************** **/ ******************************************************************** **/
static inline int Ssize(int N, int M) static inline int Ssize(int N, int M)
{ return Csize(N, M) - (N + 1); } { return Csize(N, M) - (N + 1); }
/**
* Load coefficients from a binary stream.
*
* @param[in] stream the input stream.
* @param[out] N The maximum degree of the coefficients.
* @param[out] M The maximum order of the coefficients.
* @param[out] C The vector of cosine coefficients.
* @param[out] S The vector of sine coefficients.
*
* \e N and \e M are read as 4-byte ints. \e C and \e S are resized
to
* accommodate all the coefficients (with the \e m = 0 coefficients f
or
* \e S excluded) and the data for these coefficients read as 8-byte
* doubles. The coefficients are stored in column major order. The
* bytes in the stream should use little-endian ordering. IEEE float
ing
* point is assumed for the coefficients.
********************************************************************
**/
static void readcoeffs(std::istream& stream, int& N, int& M,
std::vector<real>& C, std::vector<real>& S);
}; };
/** /**
* Evaluate a spherical harmonic sum and its gradient. * Evaluate a spherical harmonic sum and its gradient.
* *
* @tparam gradp should the gradient be calculated. * @tparam gradp should the gradient be calculated.
* @tparam norm the normalization for the associated Legendre polynomia ls. * @tparam norm the normalization for the associated Legendre polynomia ls.
* @tparam L the number of terms in the coefficents. * @tparam L the number of terms in the coefficients.
* @param[in] c an array of coeff objects. * @param[in] c an array of coeff objects.
* @param[in] f array of coefficient multipliers. f[0] should be 1. * @param[in] f array of coefficient multipliers. f[0] should be 1.
* @param[in] x the \e x component of the cartesian position. * @param[in] x the \e x component of the cartesian position.
* @param[in] y the \e y component of the cartesian position. * @param[in] y the \e y component of the cartesian position.
* @param[in] z the \e z component of the cartesian position. * @param[in] z the \e z component of the cartesian position.
* @param[in] a the normalizing radius. * @param[in] a the normalizing radius.
* @param[out] gradx the \e x component of the gradient. * @param[out] gradx the \e x component of the gradient.
* @param[out] grady the \e y component of the gradient. * @param[out] grady the \e y component of the gradient.
* @param[out] gradz the \e z component of the gradient. * @param[out] gradz the \e z component of the gradient.
* @result the spherical harmonic sum. * @result the spherical harmonic sum.
skipping to change at line 257 skipping to change at line 289
* The coefficients used by this function are, for example, * The coefficients used by this function are, for example,
* c[0].Cv + f[1] * c[1].Cv + ... + f[L-1] * c[L-1].Cv. (Note * c[0].Cv + f[1] * c[1].Cv + ... + f[L-1] * c[L-1].Cv. (Note
* that f[0] is \e not used.) The upper limits on the sum are * that f[0] is \e not used.) The upper limits on the sum are
* determined by c[0].nmx() and c[0].mmx(); these limits apply to * determined by c[0].nmx() and c[0].mmx(); these limits apply to
* \e all the components of the coefficients. The parameters \e * \e all the components of the coefficients. The parameters \e
* gradp, \e norm, and \e L are template parameters, to allow more * gradp, \e norm, and \e L are template parameters, to allow more
* optimization to be done at compile time. * optimization to be done at compile time.
* *
* Clenshaw summation is used which permits the evaluation of the sum * Clenshaw summation is used which permits the evaluation of the sum
* without the need to allocate temporary arrays. Thus this function n ever * without the need to allocate temporary arrays. Thus this function n ever
* throws an excpetion. * throws an exception.
********************************************************************** / ********************************************************************** /
template<bool gradp, normalization norm, int L> template<bool gradp, normalization norm, int L>
static Math::real Value(const coeff c[], const real f[], static Math::real Value(const coeff c[], const real f[],
real x, real y, real z, real a, real x, real y, real z, real a,
real& gradx, real& grady, real& gradz) throw( ); real& gradx, real& grady, real& gradz) throw( );
/** /**
* Create a CircularEngine object * Create a CircularEngine object
* *
* @tparam gradp should the gradient be calculated. * @tparam gradp should the gradient be calculated.
* @tparam norm the normalization for the associated Legendre polynomia ls. * @tparam norm the normalization for the associated Legendre polynomia ls.
* @tparam L the number of terms in the coefficents. * @tparam L the number of terms in the coefficients.
* @param[in] c an array of coeff objects. * @param[in] c an array of coeff objects.
* @param[in] f array of coefficient multipliers. f[0] should be 1. * @param[in] f array of coefficient multipliers. f[0] should be 1.
* @param[in] p the radius of the circle = sqrt(\e x<sup>2</sup> + \e * @param[in] p the radius of the circle = sqrt(<i>x</i><sup>2</sup> +
* y<sup>2</sup>). * <i>y</i><sup>2</sup>).
* @param[in] z the height of the circle. * @param[in] z the height of the circle.
* @param[in] a the normalizing radius. * @param[in] a the normalizing radius.
* @result the CircularEngine object. * @result the CircularEngine object.
* *
* If you need to evaluate the spherical harmonic sum for several point s * If you need to evaluate the spherical harmonic sum for several point s
* with constant \e f, \e p = sqrt(\e x<sup>2</sup> + \e y<sup>2</sup>) * with constant \e f, \e p = sqrt(<i>x</i><sup>2</sup> +
, \e * <i>y</i><sup>2</sup>), \e z, and \e a, it is more efficient to const
* z, and \e a, it is more efficient to construct call ruct
* SphericalEngine::Circle to give a CircularEngine object and then cal * call SphericalEngine::Circle to give a CircularEngine object and the
l n
* CircularEngine::operator()() with arguments \e x/\e p and \e y/\e p. * call CircularEngine::operator()() with arguments <i>x</i>/\e p and
* <i>y</i>/\e p.
********************************************************************** / ********************************************************************** /
template<bool gradp, SphericalEngine::normalization norm, int L> template<bool gradp, SphericalEngine::normalization norm, int L>
static CircularEngine Circle(const coeff c[], const real f[], static CircularEngine Circle(const coeff c[], const real f[],
real p, real z, real a); real p, real z, real a);
/**
* Check that the static table of square roots is big enough and enlarg
e it
* if necessary.
*
* @param[in] N the maximum degree to be used in SphericalEngine.
*
* Typically, there's no need for an end-user to call this routine, bec
ause
* the constructors for SphericalEngine::coeff do so. However, since t
his
* updates a static table, there's a possible race condition in a
* multi-threaded environment. Because this routine does nothing if th
e
* table is already large enough, one way to avoid race conditions is t
o
* call this routine at program start up (when it's still single thread
ed),
* supplying the largest degree that your program will use. E.g.,
\code
GeographicLib::SphericalEngine::RootTable(2190);
\endcode
* suffices to accommodate extant magnetic and gravity models.
**********************************************************************
/
static void RootTable(int N);
/**
* Clear the static table of square roots and release the memory. Call
* this only when you are sure you no longer will be using SphericalEng
ine.
* Your program will crash if you call SphericalEngine after calling th
is
* routine. <b>It's safest not to call this routine at all.</b> (The s
pace
* used by the table is modest.)
**********************************************************************
/
static void ClearRootTable() {
std::vector<real> temp(0);
root_.swap(temp);
}
}; };
} // namespace GeographicLib } // namespace GeographicLib
#if defined(_MSC_VER) #if defined(_MSC_VER)
#pragma warning (pop) #pragma warning (pop)
#endif #endif
#endif // GEOGRAPHICLIB_SPHERICALENGINE_HPP #endif // GEOGRAPHICLIB_SPHERICALENGINE_HPP
 End of changes. 19 change blocks. 
20 lines changed or deleted 99 lines changed or added


 SphericalHarmonic.hpp   SphericalHarmonic.hpp 
/** /**
* \file SphericalHarmonic.hpp * \file SphericalHarmonic.hpp
* \brief Header for GeographicLib::SphericalHarmonic class * \brief Header for GeographicLib::SphericalHarmonic class
* *
* Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed un der * Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed un der
* the MIT/X11 License. For more information, see * the MIT/X11 License. For more information, see
* http://geographiclib.sourceforge.net/ * http://geographiclib.sourceforge.net/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_SPHERICALHARMONIC_HPP) #if !defined(GEOGRAPHICLIB_SPHERICALHARMONIC_HPP)
#define GEOGRAPHICLIB_SPHERICALHARMONIC_HPP "$Id: 772f453ae5c0c8f69164279b9 #define GEOGRAPHICLIB_SPHERICALHARMONIC_HPP \
f5eeae42b45b42c $" "$Id: 12b23c05ab8bae57213d6762e45189a7cb46b397 $"
#include <vector> #include <vector>
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
#include <GeographicLib/SphericalEngine.hpp> #include <GeographicLib/SphericalEngine.hpp>
#include <GeographicLib/CircularEngine.hpp> #include <GeographicLib/CircularEngine.hpp>
#include <GeographicLib/Geocentric.hpp>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief Spherical Harmonic series * \brief Spherical Harmonic series
* *
* This class evaluates the spherical harmonic sum \verbatim * This class evaluates the spherical harmonic sum \verbatim
V(x, y, z) = sum(n = 0..N)[ q^(n+1) * sum(m = 0..n)[ V(x, y, z) = sum(n = 0..N)[ q^(n+1) * sum(m = 0..n)[
(C[n,m] * cos(m*lambda) + S[n,m] * sin(m*lambda)) * (C[n,m] * cos(m*lambda) + S[n,m] * sin(m*lambda)) *
P[n,m](cos(theta)) ] ] P[n,m](cos(theta)) ] ]
\endverbatim \endverbatim
* where * where
* - <i>p</i><sup>2</sup> = <i>x</i><sup>2</sup> + <i>y</i><sup>2</sup>, * - <i>p</i><sup>2</sup> = <i>x</i><sup>2</sup> + <i>y</i><sup>2</sup>,
* - <i>r</i><sup>2</sup> = <i>p</i><sup>2</sup> + <i>z</i><sup>2</sup>, * - <i>r</i><sup>2</sup> = <i>p</i><sup>2</sup> + <i>z</i><sup>2</sup>,
* - \e q = <i>a</i>/<i>r</i>, * - \e q = <i>a</i>/<i>r</i>,
* - \e theta = atan2(\e p, \e z) = the spherical \e colatitude, * - \e theta = atan2(\e p, \e z) = the spherical \e colatitude,
* - \e lambda = atan2(\e y, \e x) = the longitude. * - \e lambda = atan2(\e y, \e x) = the longitude.
* - P<sub>\e nm</sub>(\e t) is the associated Legendre polynomial of deg ree * - P<sub>\e nm</sub>(\e t) is the associated Legendre polynomial of deg ree
* \e n and order \e m. * \e n and order \e m.
* *
* One of two normalizations are supported for P<sub>\e nm</sub> * Two normalizations are supported for P<sub>\e nm</sub>
* - fully normalized denoted by SphericalHarmonic::full. * - fully normalized denoted by SphericalHarmonic::FULL.
* - Schmidt semi-normalized denoted by SphericalHarmonic::schmidt. * - Schmidt semi-normalized denoted by SphericalHarmonic::SCHMIDT.
* *
* Clenshaw summation is used for the sums over both \e n and \e m. This * Clenshaw summation is used for the sums over both \e n and \e m. This
* allows the computation to be carried out without the need for any * allows the computation to be carried out without the need for any
* temporary arrays. References: * temporary arrays. See SphericalEngine.cpp for more information on the
* implementation.
*
* References:
* - C. W. Clenshaw, A note on the summation of Chebyshev series, * - C. W. Clenshaw, A note on the summation of Chebyshev series,
* %Math. Tables Aids Comput. 9(51), 118-120 (1955). * %Math. Tables Aids Comput. 9(51), 118-120 (1955).
* - R. E. Deakin, Derivatives of the earth's potentials, Geomatics * - R. E. Deakin, Derivatives of the earth's potentials, Geomatics
* Research Australasia 68, 31-60, (June 1998). * Research Australasia 68, 31-60, (June 1998).
* - W. A. Heiskanen and H. Moritz, Physical Geodesy, (Freeman, San * - W. A. Heiskanen and H. Moritz, Physical Geodesy, (Freeman, San
* Fransisco, 1967). (See Sec. 1-14, for a definition of Pbar.) * Francisco, 1967). (See Sec. 1-14, for a definition of Pbar.)
* - S. A. Holmes and W. E. Featherstone, A unified approach to the * - S. A. Holmes and W. E. Featherstone, A unified approach to the
* Clenshaw summation and the recursive computation of very high degree * Clenshaw summation and the recursive computation of very high degree
* and order normalised associated Legendre functions, J. Geod. 76(5), * and order normalised associated Legendre functions, J. Geod. 76(5),
* 279-299 (2002). * 279-299 (2002).
* - C. C. Tscherning and K. Poder, Some geodetic applications of Clensha w * - C. C. Tscherning and K. Poder, Some geodetic applications of Clensha w
* summation, Boll. Geod. Sci. Aff. 41(4), 349-375 (1982). * summation, Boll. Geod. Sci. Aff. 41(4), 349-375 (1982).
**********************************************************************/ **********************************************************************/
class GEOGRAPHIC_EXPORT SphericalHarmonic { class GEOGRAPHIC_EXPORT SphericalHarmonic {
public: public:
skipping to change at line 82 skipping to change at line 87
* associated Legendre polynomial) http://dlmf.nist.gov/14.7.E10 and \e k * associated Legendre polynomial) http://dlmf.nist.gov/14.7.E10 and \e k
* = 1 for \e m = 0 and \e k = 2 otherwise. * = 1 for \e m = 0 and \e k = 2 otherwise.
* *
* The mean squared value of * The mean squared value of
* <i>P</i><sub><i>nm</i></sub><sup>full</sup>(cos \e theta) cos(\e m \e * <i>P</i><sub><i>nm</i></sub><sup>full</sup>(cos \e theta) cos(\e m \e
* lambda) and <i>P</i><sub><i>nm</i></sub><sup>full</sup>(cos \e the ta) * lambda) and <i>P</i><sub><i>nm</i></sub><sup>full</sup>(cos \e the ta)
* sin(\e m \e lambda) over the sphere is 1. * sin(\e m \e lambda) over the sphere is 1.
* *
* @hideinitializer * @hideinitializer
******************************************************************** **/ ******************************************************************** **/
full = SphericalEngine::full, FULL = SphericalEngine::FULL,
/** /**
* Schmidt semi-normalized associated Legendre polynomials. * Schmidt semi-normalized associated Legendre polynomials.
* *
* These are defined by <i>P</i><sub><i>nm</i></sub><sup>schmidt</sup >(\e * These are defined by <i>P</i><sub><i>nm</i></sub><sup>schmidt</sup >(\e
* z) = (-1)<sup><i>m</i></sup> sqrt(\e k (\e n - \e m)! / (\e n + \e * z) = (-1)<sup><i>m</i></sup> sqrt(\e k (\e n - \e m)! / (\e n + \e
* m)!) <b>P</b><sub><i>n</i></sub><sup><i>m</i></sup>(\e z), where * m)!) <b>P</b><sub><i>n</i></sub><sup><i>m</i></sup>(\e z), where
* <b>P</b><sub><i>n</i></sub><sup><i>m</i></sup>(\e z) is Ferrers * <b>P</b><sub><i>n</i></sub><sup><i>m</i></sup>(\e z) is Ferrers
* function (also known as the Legendre function on the cut or the * function (also known as the Legendre function on the cut or the
* associated Legendre polynomial) http://dlmf.nist.gov/14.7.E10 and \e k * associated Legendre polynomial) http://dlmf.nist.gov/14.7.E10 and \e k
* = 1 for \e m = 0 and \e k = 2 otherwise. * = 1 for \e m = 0 and \e k = 2 otherwise.
* *
* The mean squared value of * The mean squared value of
* <i>P</i><sub><i>nm</i></sub><sup>schmidt</sup>(cos \e theta) cos(\ e m * <i>P</i><sub><i>nm</i></sub><sup>schmidt</sup>(cos \e theta) cos(\ e m
* \e lambda) and <i>P</i><sub><i>nm</i></sub><sup>schmidt</sup>(cos \e * \e lambda) and <i>P</i><sub><i>nm</i></sub><sup>schmidt</sup>(cos \e
* theta) sin(\e m \e lambda) over the sphere is 1/(2\e n + 1). * theta) sin(\e m \e lambda) over the sphere is 1/(2\e n + 1).
* *
* @hideinitializer * @hideinitializer
******************************************************************** **/ ******************************************************************** **/
schmidt = SphericalEngine::schmidt, SCHMIDT = SphericalEngine::SCHMIDT,
/// \cond SKIP
// These are deprecated...
full = FULL,
schmidt = SCHMIDT,
/// \endcond
}; };
private: private:
typedef Math::real real; typedef Math::real real;
SphericalEngine::coeff _c[1]; SphericalEngine::coeff _c[1];
real _a; real _a;
normalization _norm; unsigned _norm;
public: public:
/** /**
* Constructor with a full set of coefficients specified. * Constructor with a full set of coefficients specified.
* *
* @param[in] C the coefficients \e C<sub>\e nm</sub>. * @param[in] C the coefficients \e C<sub>\e nm</sub>.
* @param[in] S the coefficients \e S<sub>\e nm</sub>. * @param[in] S the coefficients \e S<sub>\e nm</sub>.
* @param[in] N the maximum degree and order of the sum * @param[in] N the maximum degree and order of the sum
* @param[in] a the reference radius appearing in the definition of the * @param[in] a the reference radius appearing in the definition of the
* sum. * sum.
skipping to change at line 146 skipping to change at line 156
* <i>C</i><sub>33</sub>. * <i>C</i><sub>33</sub>.
* In general the (\e n,\e m) element is at index \e m*\e N - \e m*(\e m - * In general the (\e n,\e m) element is at index \e m*\e N - \e m*(\e m -
* 1)/2 + \e n. The layout of \e S is the same except that the first * 1)/2 + \e n. The layout of \e S is the same except that the first
* column is omitted (since the \e m = 0 terms never contribute to the sum) * column is omitted (since the \e m = 0 terms never contribute to the sum)
* and the 0th element is <i>S</i><sub>11</sub> * and the 0th element is <i>S</i><sub>11</sub>
* *
* The class stores <i>pointers</i> to the first elements of \e C and \ e S. * The class stores <i>pointers</i> to the first elements of \e C and \ e S.
* These arrays should not be altered or destroyed during the lifetime of a * These arrays should not be altered or destroyed during the lifetime of a
* SphericalHarmonic object. * SphericalHarmonic object.
********************************************************************** / ********************************************************************** /
SphericalHarmonic(const std::vector<double>& C, SphericalHarmonic(const std::vector<real>& C,
const std::vector<double>& S, const std::vector<real>& S,
int N, real a, normalization norm = full) int N, real a, unsigned norm = FULL)
: _a(a) : _a(a)
, _norm(norm) , _norm(norm)
{ _c[0] = SphericalEngine::coeff(C, S, N); } { _c[0] = SphericalEngine::coeff(C, S, N); }
/** /**
* Constructor with a subset of coefficients specified. * Constructor with a subset of coefficients specified.
* *
* @param[in] C the coefficients \e C<sub>\e nm</sub>. * @param[in] C the coefficients \e C<sub>\e nm</sub>.
* @param[in] S the coefficients \e S<sub>\e nm</sub>. * @param[in] S the coefficients \e S<sub>\e nm</sub>.
* @param[in] N the degree used to determine the layout of \e C and \e S. * @param[in] N the degree used to determine the layout of \e C and \e S.
* @param[in] nmx the maximum degree used in the sum. The sum over \e n is * @param[in] nmx the maximum degree used in the sum. The sum over \e n is
* from 0 thru \e nmx. * from 0 thru \e nmx.
* @param[in] mmx the maximum order used in the sum. The sum over \e m is * @param[in] mmx the maximum order used in the sum. The sum over \e m is
* from 0 thru min(\e n, \e mmx). * from 0 thru min(\e n, \e mmx).
* @param[in] a the reference radius appearing in the definition of the * @param[in] a the reference radius appearing in the definition of the
* sum. * sum.
* @param[in] norm the normalization for the associated Legendre * @param[in] norm the normalization for the associated Legendre
* polynomials, either SphericalHarmonic::full (the default) or * polynomials, either SphericalHarmonic::FULL (the default) or
* SphericalHarmonic::schmidt. * SphericalHarmonic::SCHMIDT.
* *
* The class stores <i>pointers</i> to the first elements of \e C and \ e S. * The class stores <i>pointers</i> to the first elements of \e C and \ e S.
* These arrays should not be altered or destroyed during the lifetime of a * These arrays should not be altered or destroyed during the lifetime of a
* SphericalHarmonic object. * SphericalHarmonic object.
********************************************************************** / ********************************************************************** /
SphericalHarmonic(const std::vector<double>& C, SphericalHarmonic(const std::vector<real>& C,
const std::vector<double>& S, const std::vector<real>& S,
int N, int nmx, int mmx, int N, int nmx, int mmx,
real a, normalization norm = full) real a, unsigned norm = FULL)
: _a(a) : _a(a)
, _norm(norm) , _norm(norm)
{ _c[0] = SphericalEngine::coeff(C, S, N, nmx, mmx); } { _c[0] = SphericalEngine::coeff(C, S, N, nmx, mmx); }
/** /**
* A default constructor so that the object can be created when the * A default constructor so that the object can be created when the
* constructor for another object is initialized. This default object can * constructor for another object is initialized. This default object can
* then be reset with the default copy assignment operator. * then be reset with the default copy assignment operator.
********************************************************************** / ********************************************************************** /
SphericalHarmonic() {} SphericalHarmonic() {}
skipping to change at line 204 skipping to change at line 214
* @return \e V the spherical harmonic sum. * @return \e V the spherical harmonic sum.
* *
* This routine requires constant memory and thus never throws an * This routine requires constant memory and thus never throws an
* exception. * exception.
********************************************************************** / ********************************************************************** /
Math::real operator()(real x, real y, real z) const throw() { Math::real operator()(real x, real y, real z) const throw() {
real f[] = {1}; real f[] = {1};
real v = 0; real v = 0;
real dummy; real dummy;
switch (_norm) { switch (_norm) {
case full: case FULL:
v = SphericalEngine::Value<false, SphericalEngine::full, 1> v = SphericalEngine::Value<false, SphericalEngine::FULL, 1>
(_c, f, x, y, z, _a, dummy, dummy, dummy); (_c, f, x, y, z, _a, dummy, dummy, dummy);
break; break;
case schmidt: case SCHMIDT:
v = SphericalEngine::Value<false, SphericalEngine::schmidt, 1> v = SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 1>
(_c, f, x, y, z, _a, dummy, dummy, dummy); (_c, f, x, y, z, _a, dummy, dummy, dummy);
break; break;
} }
return v; return v;
} }
/** /**
* Compute a spherical harmonic sum and its gradient. * Compute a spherical harmonic sum and its gradient.
* *
* @param[in] x cartesian coordinate. * @param[in] x cartesian coordinate.
skipping to change at line 237 skipping to change at line 247
* This is the same as the previous function, except that the component s of * This is the same as the previous function, except that the component s of
* the gradients of the sum in the \e x, \e y, and \e z directions are * the gradients of the sum in the \e x, \e y, and \e z directions are
* computed. This routine requires constant memory and thus never thro ws * computed. This routine requires constant memory and thus never thro ws
* an exception. * an exception.
********************************************************************** / ********************************************************************** /
Math::real operator()(real x, real y, real z, Math::real operator()(real x, real y, real z,
real& gradx, real& grady, real& gradz) const thro w() { real& gradx, real& grady, real& gradz) const thro w() {
real f[] = {1}; real f[] = {1};
real v = 0; real v = 0;
switch (_norm) { switch (_norm) {
case full: case FULL:
v = SphericalEngine::Value<true, SphericalEngine::full, 1> v = SphericalEngine::Value<true, SphericalEngine::FULL, 1>
(_c, f, x, y, z, _a, gradx, grady, gradz); (_c, f, x, y, z, _a, gradx, grady, gradz);
break; break;
case schmidt: case SCHMIDT:
v = SphericalEngine::Value<true, SphericalEngine::schmidt, 1> v = SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 1>
(_c, f, x, y, z, _a, gradx, grady, gradz); (_c, f, x, y, z, _a, gradx, grady, gradz);
break; break;
} }
return v; return v;
} }
/** /**
* Create a CircularEngine to allow the efficient evaluation of several * Create a CircularEngine to allow the efficient evaluation of several
* points on a circle of latitude. * points on a circle of latitude.
* *
* @param[in] p the radius of the circle. * @param[in] p the radius of the circle.
* @param[in] z the height of the circle above the equatorial plane. * @param[in] z the height of the circle above the equatorial plane.
* @param[in] gradp if true the returned object will be able to compute the * @param[in] gradp if true the returned object will be able to compute the
* gradient of the sum. * gradient of the sum.
* @return the CircularEngine object. * @return the CircularEngine object.
* *
* SphericalHarmonic::operator()() exchanges the order of the sums in t he * SphericalHarmonic::operator()() exchanges the order of the sums in t he
* definition, i.e., sum(n = 0..N)[sum(m = 0..n)[...]] becomes sum(m = * definition, i.e., sum(n = 0..N)[sum(m = 0..n)[...]] becomes sum(m =
* 0..N)[sum(n = m..N)[...]]. SphericalHarmonic::Circle performs the i nner * 0..N)[sum(n = m..N)[...]]. SphericalHarmonic::Circle performs the i nner
* sum over degree \e n (which entails about \e N<sup>2</sup> operation * sum over degree \e n (which entails about <i>N</i><sup>2</sup>
s). * operations). Calling CircularEngine::operator()() on the returned
* Calling CircularEngine::operator()() on the returned object performs * object performs the outer sum over the order \e m (about \e N
the * operations). This routine may throw a bad_alloc exception in the
* outer sum over the order \e m (about \e N operations). This routine * CircularEngine constructor.
may
* throw a bad_alloc exception in the GeographicLib::CircularEngine
* constructor.
* *
* Here's an example of computing the spherical sum at a sequence of * Here's an example of computing the spherical sum at a sequence of
* longitudes without using a CircularEngine object * longitudes without using a CircularEngine object
\code \code
SphericalHarmonic h(...); // Create the SphericalHarmonic object SphericalHarmonic h(...); // Create the SphericalHarmonic object
double r = 2, lat = 33, lon0 = 44, dlon = 0.01; double r = 2, lat = 33, lon0 = 44, dlon = 0.01;
double double
phi = lat * Math::degree<double>(), phi = lat * Math::degree<double>(),
z = r * sin(phi), p = r * cos(phi); z = r * sin(phi), p = r * cos(phi);
for (int i = 0; i <= 100; ++i) { for (int i = 0; i <= 100; ++i) {
real real
lon = lon0 + i * dlon, lon = lon0 + i * dlon,
lam = lon * Math::degree<double>(); lam = lon * Math::degree<double>();
std::cout << lon << " " << h(p * cos(lam), p * sin(lam), z) << "\n"; std::cout << lon << " " << h(p * cos(lam), p * sin(lam), z) << "\n";
} }
\endcode \endcode
* Here is the same calculation done using a CircularEngine object. Th is * Here is the same calculation done using a CircularEngine object. Th is
* will be about \e N/2 times faster. * will be about <i>N</i>/2 times faster.
\code \code
SphericalHarmonic h(...); // Create the SphericalHarmonic object SphericalHarmonic h(...); // Create the SphericalHarmonic object
double r = 2, lat = 33, lon0 = 44, dlon = 0.01; double r = 2, lat = 33, lon0 = 44, dlon = 0.01;
double double
phi = lat * Math::degree<double>(), phi = lat * Math::degree<double>(),
z = r * sin(phi), p = r * cos(phi); z = r * sin(phi), p = r * cos(phi);
CircularEngine c(h(p, z, false)); // Create the CircularEngine object CircularEngine c(h(p, z, false)); // Create the CircularEngine object
for (int i = 0; i <= 100; ++i) { for (int i = 0; i <= 100; ++i) {
real real
lon = lon0 + i * dlon; lon = lon0 + i * dlon;
std::cout << lon << " " << c(lon) << "\n"; std::cout << lon << " " << c(lon) << "\n";
} }
\endcode \endcode
********************************************************************** / ********************************************************************** /
CircularEngine Circle(real p, real z, bool gradp) const { CircularEngine Circle(real p, real z, bool gradp) const {
real f[] = {1}; real f[] = {1};
switch (_norm) { switch (_norm) {
case full: case FULL:
return gradp ? return gradp ?
SphericalEngine::Circle<true, SphericalEngine::full, 1> SphericalEngine::Circle<true, SphericalEngine::FULL, 1>
(_c, f, p, z, _a) : (_c, f, p, z, _a) :
SphericalEngine::Circle<false, SphericalEngine::full, 1> SphericalEngine::Circle<false, SphericalEngine::FULL, 1>
(_c, f, p, z, _a); (_c, f, p, z, _a);
break; break;
case schmidt: case SCHMIDT:
default: // To avoid compiler warnings default: // To avoid compiler warnings
return gradp ? return gradp ?
SphericalEngine::Circle<true, SphericalEngine::schmidt, 1> SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 1>
(_c, f, p, z, _a) : (_c, f, p, z, _a) :
SphericalEngine::Circle<false, SphericalEngine::schmidt, 1> SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 1>
(_c, f, p, z, _a); (_c, f, p, z, _a);
break; break;
} }
} }
/**
* @return the zeroth SphericalEngine::coeff object.
**********************************************************************
/
const SphericalEngine::coeff& Coefficients() const throw()
{ return _c[0]; }
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif // GEOGRAPHICLIB_SPHERICALHARMONIC_HPP #endif // GEOGRAPHICLIB_SPHERICALHARMONIC_HPP
 End of changes. 25 change blocks. 
41 lines changed or deleted 54 lines changed or added


 SphericalHarmonic1.hpp   SphericalHarmonic1.hpp 
/** /**
* \file SphericalHarmonic1.hpp * \file SphericalHarmonic1.hpp
* \brief Header for GeographicLib::SphericalHarmonic1 class * \brief Header for GeographicLib::SphericalHarmonic1 class
* *
* Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed un der * Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed un der
* the MIT/X11 License. For more information, see * the MIT/X11 License. For more information, see
* http://geographiclib.sourceforge.net/ * http://geographiclib.sourceforge.net/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_SPHERICALHARMONIC1_HPP) #if !defined(GEOGRAPHICLIB_SPHERICALHARMONIC1_HPP)
#define GEOGRAPHICLIB_SPHERICALHARMONIC1_HPP "$Id: 65b2e8b9b73bae6b1b575da2 #define GEOGRAPHICLIB_SPHERICALHARMONIC1_HPP \
2ba8bcf634d3a00b $" "$Id: f873848926021872a17b732171bb0bbbba1bb520 $"
#include <vector> #include <vector>
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
#include <GeographicLib/SphericalEngine.hpp> #include <GeographicLib/SphericalEngine.hpp>
#include <GeographicLib/CircularEngine.hpp> #include <GeographicLib/CircularEngine.hpp>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief Spherical Harmonic series with a correction to the coefficients . * \brief Spherical Harmonic series with a correction to the coefficients .
skipping to change at line 36 skipping to change at line 37
**********************************************************************/ **********************************************************************/
class GEOGRAPHIC_EXPORT SphericalHarmonic1 { class GEOGRAPHIC_EXPORT SphericalHarmonic1 {
public: public:
/** /**
* Supported normalizations for associate Legendre polynomials. * Supported normalizations for associate Legendre polynomials.
********************************************************************** / ********************************************************************** /
enum normalization { enum normalization {
/** /**
* Fully normalized associated Legendre polynomials. See * Fully normalized associated Legendre polynomials. See
* SphericalHarmonic::full for documentation. * SphericalHarmonic::FULL for documentation.
* *
* @hideinitializer * @hideinitializer
******************************************************************** **/ ******************************************************************** **/
full = SphericalEngine::full, FULL = SphericalEngine::FULL,
/** /**
* Schmidt semi-normalized associated Legendre polynomials. See * Schmidt semi-normalized associated Legendre polynomials. See
* SphericalHarmonic::schmidt for documentation. * SphericalHarmonic::SCHMIDT for documentation.
* *
* @hideinitializer * @hideinitializer
******************************************************************** **/ ******************************************************************** **/
schmidt = SphericalEngine::schmidt, SCHMIDT = SphericalEngine::SCHMIDT,
/// \cond SKIP
// These are deprecated...
full = FULL,
schmidt = SCHMIDT,
/// \endcond
}; };
private: private:
typedef Math::real real; typedef Math::real real;
SphericalEngine::coeff _c[2]; SphericalEngine::coeff _c[2];
real _a; real _a;
normalization _norm; unsigned _norm;
public: public:
/** /**
* Constructor with a full set of coefficients specified. * Constructor with a full set of coefficients specified.
* *
* @param[in] C the coefficients \e C<sub>\e nm</sub>. * @param[in] C the coefficients \e C<sub>\e nm</sub>.
* @param[in] S the coefficients \e S<sub>\e nm</sub>. * @param[in] S the coefficients \e S<sub>\e nm</sub>.
* @param[in] N the maximum degree and order of the sum * @param[in] N the maximum degree and order of the sum
* @param[in] C1 the coefficients \e C'<sub>\e nm</sub>. * @param[in] C1 the coefficients \e C'<sub>\e nm</sub>.
* @param[in] S1 the coefficients \e S'<sub>\e nm</sub>. * @param[in] S1 the coefficients \e S'<sub>\e nm</sub>.
* @param[in] N1 the maximum degree and order of the correction * @param[in] N1 the maximum degree and order of the correction
* coefficients \e C'<sub>\e nm</sub> and \e S'<sub>\e nm</sub>. * coefficients \e C'<sub>\e nm</sub> and \e S'<sub>\e nm</sub>.
* @param[in] a the reference radius appearing in the definition of the * @param[in] a the reference radius appearing in the definition of the
* sum. * sum.
* @param[in] norm the normalization for the associated Legendre * @param[in] norm the normalization for the associated Legendre
* polynomials, either SphericalHarmonic1::full (the default) or * polynomials, either SphericalHarmonic1::FULL (the default) or
* SphericalHarmonic1::schmidt. * SphericalHarmonic1::SCHMIDT.
* *
* See SphericalHarmonic for the way the coefficients should be stored. \e * See SphericalHarmonic for the way the coefficients should be stored. \e
* N1 should satisfy \e N1 <= \e N. * N1 should satisfy \e N1 <= \e N.
* *
* The class stores <i>pointers</i> to the first elements of \e C, \e S , \e * The class stores <i>pointers</i> to the first elements of \e C, \e S , \e
* C', and \e S'. These arrays should not be altered or destroyed duri ng * C', and \e S'. These arrays should not be altered or destroyed duri ng
* the lifetime of a SphericalHarmonic object. * the lifetime of a SphericalHarmonic object.
********************************************************************** / ********************************************************************** /
SphericalHarmonic1(const std::vector<double>& C, SphericalHarmonic1(const std::vector<real>& C,
const std::vector<double>& S, const std::vector<real>& S,
int N, int N,
const std::vector<double>& C1, const std::vector<real>& C1,
const std::vector<double>& S1, const std::vector<real>& S1,
int N1, int N1,
real a, normalization norm = full) real a, unsigned norm = FULL)
: _a(a) : _a(a)
, _norm(norm) { , _norm(norm) {
if (!(N1 <= N)) if (!(N1 <= N))
throw GeographicErr("N1 cannot be larger that N"); throw GeographicErr("N1 cannot be larger that N");
_c[0] = SphericalEngine::coeff(C, S, N); _c[0] = SphericalEngine::coeff(C, S, N);
_c[1] = SphericalEngine::coeff(C1, S1, N1); _c[1] = SphericalEngine::coeff(C1, S1, N1);
} }
/** /**
* Constructor with a subset of coefficients specified. * Constructor with a subset of coefficients specified.
skipping to change at line 114 skipping to change at line 120
* from 0 thru min(\e n, \e mmx). * from 0 thru min(\e n, \e mmx).
* @param[in] C1 the coefficients \e C'<sub>\e nm</sub>. * @param[in] C1 the coefficients \e C'<sub>\e nm</sub>.
* @param[in] S1 the coefficients \e S'<sub>\e nm</sub>. * @param[in] S1 the coefficients \e S'<sub>\e nm</sub>.
* @param[in] N1 the degree used to determine the layout of \e C' and \ e * @param[in] N1 the degree used to determine the layout of \e C' and \ e
* S'. * S'.
* @param[in] nmx1 the maximum degree used for \e C' and \e S'. * @param[in] nmx1 the maximum degree used for \e C' and \e S'.
* @param[in] mmx1 the maximum order used for \e C' and \e S'. * @param[in] mmx1 the maximum order used for \e C' and \e S'.
* @param[in] a the reference radius appearing in the definition of the * @param[in] a the reference radius appearing in the definition of the
* sum. * sum.
* @param[in] norm the normalization for the associated Legendre * @param[in] norm the normalization for the associated Legendre
* polynomials, either SphericalHarmonic1::full (the default) or * polynomials, either SphericalHarmonic1::FULL (the default) or
* SphericalHarmonic1::schmidt. * SphericalHarmonic1::SCHMIDT.
* *
* The class stores <i>pointers</i> to the first elements of \e C, \e S , \e * The class stores <i>pointers</i> to the first elements of \e C, \e S , \e
* C', and \e S'. These arrays should not be altered or destroyed duri ng * C', and \e S'. These arrays should not be altered or destroyed duri ng
* the lifetime of a SphericalHarmonic object. * the lifetime of a SphericalHarmonic object.
********************************************************************** / ********************************************************************** /
SphericalHarmonic1(const std::vector<double>& C, SphericalHarmonic1(const std::vector<real>& C,
const std::vector<double>& S, const std::vector<real>& S,
int N, int nmx, int mmx, int N, int nmx, int mmx,
const std::vector<double>& C1, const std::vector<real>& C1,
const std::vector<double>& S1, const std::vector<real>& S1,
int N1, int nmx1, int mmx1, int N1, int nmx1, int mmx1,
real a, normalization norm = full) real a, unsigned norm = FULL)
: _a(a) : _a(a)
, _norm(norm) { , _norm(norm) {
if (!(nmx1 <= nmx)) if (!(nmx1 <= nmx))
throw GeographicErr("nmx1 cannot be larger that nmx"); throw GeographicErr("nmx1 cannot be larger that nmx");
if (!(mmx1 <= mmx)) if (!(mmx1 <= mmx))
throw GeographicErr("mmx1 cannot be larger that mmx"); throw GeographicErr("mmx1 cannot be larger that mmx");
_c[0] = SphericalEngine::coeff(C, S, N, nmx, mmx); _c[0] = SphericalEngine::coeff(C, S, N, nmx, mmx);
_c[1] = SphericalEngine::coeff(C1, S1, N1, nmx1, mmx1); _c[1] = SphericalEngine::coeff(C1, S1, N1, nmx1, mmx1);
} }
skipping to change at line 162 skipping to change at line 168
* @return \e V the spherical harmonic sum. * @return \e V the spherical harmonic sum.
* *
* This routine requires constant memory and thus never throws * This routine requires constant memory and thus never throws
* an exception. * an exception.
********************************************************************** / ********************************************************************** /
Math::real operator()(real tau, real x, real y, real z) const throw() { Math::real operator()(real tau, real x, real y, real z) const throw() {
real f[] = {1, tau}; real f[] = {1, tau};
real v = 0; real v = 0;
real dummy; real dummy;
switch (_norm) { switch (_norm) {
case full: case FULL:
v = SphericalEngine::Value<false, SphericalEngine::full, 2> v = SphericalEngine::Value<false, SphericalEngine::FULL, 2>
(_c, f, x, y, z, _a, dummy, dummy, dummy); (_c, f, x, y, z, _a, dummy, dummy, dummy);
break; break;
case schmidt: case SCHMIDT:
v = SphericalEngine::Value<false, SphericalEngine::schmidt, 2> v = SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 2>
(_c, f, x, y, z, _a, dummy, dummy, dummy); (_c, f, x, y, z, _a, dummy, dummy, dummy);
break; break;
} }
return v; return v;
} }
/** /**
* Compute a spherical harmonic sum with a correction term and its * Compute a spherical harmonic sum with a correction term and its
* gradient. * gradient.
* *
skipping to change at line 197 skipping to change at line 203
* This is the same as the previous function, except that the component s of * This is the same as the previous function, except that the component s of
* the gradients of the sum in the \e x, \e y, and \e z directions are * the gradients of the sum in the \e x, \e y, and \e z directions are
* computed. This routine requires constant memory and thus never thro ws * computed. This routine requires constant memory and thus never thro ws
* an exception. * an exception.
********************************************************************** / ********************************************************************** /
Math::real operator()(real tau, real x, real y, real z, Math::real operator()(real tau, real x, real y, real z,
real& gradx, real& grady, real& gradz) const thro w() { real& gradx, real& grady, real& gradz) const thro w() {
real f[] = {1, tau}; real f[] = {1, tau};
real v = 0; real v = 0;
switch (_norm) { switch (_norm) {
case full: case FULL:
v = SphericalEngine::Value<true, SphericalEngine::full, 2> v = SphericalEngine::Value<true, SphericalEngine::FULL, 2>
(_c, f, x, y, z, _a, gradx, grady, gradz); (_c, f, x, y, z, _a, gradx, grady, gradz);
break; break;
case schmidt: case SCHMIDT:
v = SphericalEngine::Value<true, SphericalEngine::schmidt, 2> v = SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 2>
(_c, f, x, y, z, _a, gradx, grady, gradz); (_c, f, x, y, z, _a, gradx, grady, gradz);
break; break;
} }
return v; return v;
} }
/** /**
* Create a CircularEngine to allow the efficient evaluation of several * Create a CircularEngine to allow the efficient evaluation of several
* points on a circle of latitude at a fixed value of \e tau. * points on a circle of latitude at a fixed value of \e tau.
* *
* @param[in] tau the multiplier for the correction coefficients. * @param[in] tau the multiplier for the correction coefficients.
* @param[in] p the radius of the circle. * @param[in] p the radius of the circle.
* @param[in] z the height of the circle above the equatorial plane. * @param[in] z the height of the circle above the equatorial plane.
* @param[in] gradp if true the returned object will be able to compute the * @param[in] gradp if true the returned object will be able to compute the
* gradient of the sum. * gradient of the sum.
* @return the CircularEngine object. * @return the CircularEngine object.
* *
* SphericalHarmonic1::operator()() exchanges the order of the sums in the * SphericalHarmonic1::operator()() exchanges the order of the sums in the
* definition, i.e., sum(n = 0..N)[sum(m = 0..n)[...]] becomes sum(m = * definition, i.e., sum(n = 0..N)[sum(m = 0..n)[...]] becomes sum(m =
* 0..N)[sum(n = m..N)[...]]. SphericalHarmonic1::Circle performs the * 0..N)[sum(n = m..N)[...]]. SphericalHarmonic1::Circle performs the
* inner sum over degree \e n (which entails about \e N<sup>2</sup> * inner sum over degree \e n (which entails about <i>N</i><sup>2</sup>
* operations). Calling CircularEngine::operator()() on the returned * operations). Calling CircularEngine::operator()() on the returned
* object performs the outer sum over the order \e m (about \e N * object performs the outer sum over the order \e m (about \e N
* operations). This routine may throw a bad_alloc exception in the * operations). This routine may throw a bad_alloc exception in the
* GeographicLib::CircularEngine constructor. * CircularEngine constructor.
* *
* See SphericalHarmonic::Circle for an example of its use. * See SphericalHarmonic::Circle for an example of its use.
********************************************************************** / ********************************************************************** /
CircularEngine Circle(real tau, real p, real z, bool gradp) const { CircularEngine Circle(real tau, real p, real z, bool gradp) const {
real f[] = {1, tau}; real f[] = {1, tau};
switch (_norm) { switch (_norm) {
case full: case FULL:
return gradp ? return gradp ?
SphericalEngine::Circle<true, SphericalEngine::full, 2> SphericalEngine::Circle<true, SphericalEngine::FULL, 2>
(_c, f, p, z, _a) : (_c, f, p, z, _a) :
SphericalEngine::Circle<false, SphericalEngine::full, 2> SphericalEngine::Circle<false, SphericalEngine::FULL, 2>
(_c, f, p, z, _a); (_c, f, p, z, _a);
break; break;
case schmidt: case SCHMIDT:
default: // To avoid compiler warnings default: // To avoid compiler warnings
return gradp ? return gradp ?
SphericalEngine::Circle<true, SphericalEngine::schmidt, 2> SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 2>
(_c, f, p, z, _a) : (_c, f, p, z, _a) :
SphericalEngine::Circle<false, SphericalEngine::schmidt, 2> SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 2>
(_c, f, p, z, _a); (_c, f, p, z, _a);
break; break;
} }
} }
/**
* @return the zeroth SphericalEngine::coeff object.
**********************************************************************
/
const SphericalEngine::coeff& Coefficients() const throw()
{ return _c[0]; }
/**
* @return the first SphericalEngine::coeff object.
**********************************************************************
/
const SphericalEngine::coeff& Coefficients1() const throw()
{ return _c[1]; }
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif // GEOGRAPHICLIB_SPHERICALHARMONIC1_HPP #endif // GEOGRAPHICLIB_SPHERICALHARMONIC1_HPP
 End of changes. 27 change blocks. 
37 lines changed or deleted 55 lines changed or added


 SphericalHarmonic2.hpp   SphericalHarmonic2.hpp 
/** /**
* \file SphericalHarmonic2.hpp * \file SphericalHarmonic2.hpp
* \brief Header for GeographicLib::SphericalHarmonic2 class * \brief Header for GeographicLib::SphericalHarmonic2 class
* *
* Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed un der * Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed un der
* the MIT/X11 License. For more information, see * the MIT/X11 License. For more information, see
* http://geographiclib.sourceforge.net/ * http://geographiclib.sourceforge.net/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_SPHERICALHARMONIC2_HPP) #if !defined(GEOGRAPHICLIB_SPHERICALHARMONIC2_HPP)
#define GEOGRAPHICLIB_SPHERICALHARMONIC2_HPP "$Id: dc815a8205bbae87d3644671 #define GEOGRAPHICLIB_SPHERICALHARMONIC2_HPP \
11811b952c48af57 $" "$Id: eb59f9ca6566eedfb06c24652873a69e4784f302 $"
#include <vector> #include <vector>
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
#include <GeographicLib/SphericalEngine.hpp> #include <GeographicLib/SphericalEngine.hpp>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief Spherical Harmonic series with two corrections to the coefficie nts. * \brief Spherical Harmonic series with two corrections to the coefficie nts.
* *
skipping to change at line 36 skipping to change at line 37
**********************************************************************/ **********************************************************************/
class GEOGRAPHIC_EXPORT SphericalHarmonic2 { class GEOGRAPHIC_EXPORT SphericalHarmonic2 {
public: public:
/** /**
* Supported normalizations for associate Legendre polynomials. * Supported normalizations for associate Legendre polynomials.
********************************************************************** / ********************************************************************** /
enum normalization { enum normalization {
/** /**
* Fully normalized associated Legendre polynomials. See * Fully normalized associated Legendre polynomials. See
* SphericalHarmonic::full for documentation. * SphericalHarmonic::FULL for documentation.
* *
* @hideinitializer * @hideinitializer
******************************************************************** **/ ******************************************************************** **/
full = SphericalEngine::full, FULL = SphericalEngine::FULL,
/** /**
* Schmidt semi-normalized associated Legendre polynomials. See * Schmidt semi-normalized associated Legendre polynomials. See
* SphericalHarmonic::schmidt for documentation. * SphericalHarmonic::SCHMIDT for documentation.
* *
* @hideinitializer * @hideinitializer
******************************************************************** **/ ******************************************************************** **/
schmidt = SphericalEngine::schmidt, SCHMIDT = SphericalEngine::SCHMIDT,
/// \cond SKIP
// These are deprecated...
full = FULL,
schmidt = SCHMIDT,
/// \endcond
}; };
private: private:
typedef Math::real real; typedef Math::real real;
SphericalEngine::coeff _c[3]; SphericalEngine::coeff _c[3];
real _a; real _a;
normalization _norm; normalization _norm;
public: public:
/** /**
skipping to change at line 74 skipping to change at line 80
* @param[in] S1 the coefficients \e S'<sub>\e nm</sub>. * @param[in] S1 the coefficients \e S'<sub>\e nm</sub>.
* @param[in] N1 the maximum degree and order of the first correction * @param[in] N1 the maximum degree and order of the first correction
* coefficients \e C'<sub>\e nm</sub> and \e S'<sub>\e nm</sub>. * coefficients \e C'<sub>\e nm</sub> and \e S'<sub>\e nm</sub>.
* @param[in] C2 the coefficients \e C''<sub>\e nm</sub>. * @param[in] C2 the coefficients \e C''<sub>\e nm</sub>.
* @param[in] S2 the coefficients \e S''<sub>\e nm</sub>. * @param[in] S2 the coefficients \e S''<sub>\e nm</sub>.
* @param[in] N2 the maximum degree and order of the second correction * @param[in] N2 the maximum degree and order of the second correction
* coefficients \e C'<sub>\e nm</sub> and \e S'<sub>\e nm</sub>. * coefficients \e C'<sub>\e nm</sub> and \e S'<sub>\e nm</sub>.
* @param[in] a the reference radius appearing in the definition of the * @param[in] a the reference radius appearing in the definition of the
* sum. * sum.
* @param[in] norm the normalization for the associated Legendre * @param[in] norm the normalization for the associated Legendre
* polynomials, either SphericalHarmonic2::full (the default) or * polynomials, either SphericalHarmonic2::FULL (the default) or
* SphericalHarmonic2::schmidt. * SphericalHarmonic2::SCHMIDT.
* *
* See SphericalHarmonic for the way the coefficients should be stored. \e * See SphericalHarmonic for the way the coefficients should be stored. \e
* N1 and \e N2 should satisfy \e N1 <= \e N and \e N2 <= \e N. * N1 and \e N2 should satisfy \e N1 <= \e N and \e N2 <= \e N.
* *
* The class stores <i>pointers</i> to the first elements of \e C, \e S , \e * The class stores <i>pointers</i> to the first elements of \e C, \e S , \e
* C', \e S', \e C'', and \e S''. These arrays should not be altered o r * C', \e S', \e C'', and \e S''. These arrays should not be altered o r
* destroyed during the lifetime of a SphericalHarmonic object. * destroyed during the lifetime of a SphericalHarmonic object.
********************************************************************** / ********************************************************************** /
SphericalHarmonic2(const std::vector<double>& C, SphericalHarmonic2(const std::vector<real>& C,
const std::vector<double>& S, const std::vector<real>& S,
int N, int N,
const std::vector<double>& C1, const std::vector<real>& C1,
const std::vector<double>& S1, const std::vector<real>& S1,
int N1, int N1,
const std::vector<double>& C2, const std::vector<real>& C2,
const std::vector<double>& S2, const std::vector<real>& S2,
int N2, int N2,
real a, normalization norm = full) real a, unsigned norm = FULL)
: _a(a) : _a(a)
, _norm(norm) { , _norm(norm) {
if (!(N1 <= N && N2 <= N)) if (!(N1 <= N && N2 <= N))
throw GeographicErr("N1 and N2 cannot be larger that N"); throw GeographicErr("N1 and N2 cannot be larger that N");
_c[0] = SphericalEngine::coeff(C, S, N); _c[0] = SphericalEngine::coeff(C, S, N);
_c[1] = SphericalEngine::coeff(C1, S1, N1); _c[1] = SphericalEngine::coeff(C1, S1, N1);
_c[2] = SphericalEngine::coeff(C2, S2, N2); _c[2] = SphericalEngine::coeff(C2, S2, N2);
} }
/** /**
skipping to change at line 128 skipping to change at line 134
* @param[in] mmx1 the maximum order used for \e C' and \e S'. * @param[in] mmx1 the maximum order used for \e C' and \e S'.
* @param[in] C2 the coefficients \e C''<sub>\e nm</sub>. * @param[in] C2 the coefficients \e C''<sub>\e nm</sub>.
* @param[in] S2 the coefficients \e S''<sub>\e nm</sub>. * @param[in] S2 the coefficients \e S''<sub>\e nm</sub>.
* @param[in] N2 the degree used to determine the layout of \e C'' and \e * @param[in] N2 the degree used to determine the layout of \e C'' and \e
* S''. * S''.
* @param[in] nmx2 the maximum degree used for \e C'' and \e S''. * @param[in] nmx2 the maximum degree used for \e C'' and \e S''.
* @param[in] mmx2 the maximum order used for \e C'' and \e S''. * @param[in] mmx2 the maximum order used for \e C'' and \e S''.
* @param[in] a the reference radius appearing in the definition of the * @param[in] a the reference radius appearing in the definition of the
* sum. * sum.
* @param[in] norm the normalization for the associated Legendre * @param[in] norm the normalization for the associated Legendre
* polynomials, either SphericalHarmonic2::full (the default) or * polynomials, either SphericalHarmonic2::FULL (the default) or
* SphericalHarmonic2::schmidt. * SphericalHarmonic2::SCHMIDT.
* *
* The class stores <i>pointers</i> to the first elements of \e C, \e S , \e * The class stores <i>pointers</i> to the first elements of \e C, \e S , \e
* C', \e S', \e C'', and \e S''. These arrays should not be altered o r * C', \e S', \e C'', and \e S''. These arrays should not be altered o r
* destroyed during the lifetime of a SphericalHarmonic object. * destroyed during the lifetime of a SphericalHarmonic object.
********************************************************************** / ********************************************************************** /
SphericalHarmonic2(const std::vector<double>& C, SphericalHarmonic2(const std::vector<real>& C,
const std::vector<double>& S, const std::vector<real>& S,
int N, int nmx, int mmx, int N, int nmx, int mmx,
const std::vector<double>& C1, const std::vector<real>& C1,
const std::vector<double>& S1, const std::vector<real>& S1,
int N1, int nmx1, int mmx1, int N1, int nmx1, int mmx1,
const std::vector<double>& C2, const std::vector<real>& C2,
const std::vector<double>& S2, const std::vector<real>& S2,
int N2, int nmx2, int mmx2, int N2, int nmx2, int mmx2,
real a, normalization norm = full) real a, unsigned norm = FULL)
: _a(a) : _a(a)
, _norm(norm) { , _norm(norm) {
if (!(nmx1 <= nmx && nmx2 <= nmx)) if (!(nmx1 <= nmx && nmx2 <= nmx))
throw GeographicErr("nmx1 and nmx2 cannot be larger that nmx"); throw GeographicErr("nmx1 and nmx2 cannot be larger that nmx");
if (!(mmx1 <= mmx && mmx2 <= mmx)) if (!(mmx1 <= mmx && mmx2 <= mmx))
throw GeographicErr("mmx1 and mmx2cannot be larger that mmx"); throw GeographicErr("mmx1 and mmx2cannot be larger that mmx");
_c[0] = SphericalEngine::coeff(C, S, N, nmx, mmx); _c[0] = SphericalEngine::coeff(C, S, N, nmx, mmx);
_c[1] = SphericalEngine::coeff(C1, S1, N1, nmx1, mmx1); _c[1] = SphericalEngine::coeff(C1, S1, N1, nmx1, mmx1);
_c[2] = SphericalEngine::coeff(C2, S2, N2, nmx2, mmx2); _c[2] = SphericalEngine::coeff(C2, S2, N2, nmx2, mmx2);
} }
skipping to change at line 182 skipping to change at line 188
* *
* This routine requires constant memory and thus never throws an * This routine requires constant memory and thus never throws an
* exception. * exception.
********************************************************************** / ********************************************************************** /
Math::real operator()(real tau1, real tau2, real x, real y, real z) Math::real operator()(real tau1, real tau2, real x, real y, real z)
const throw() { const throw() {
real f[] = {1, tau1, tau2}; real f[] = {1, tau1, tau2};
real v = 0; real v = 0;
real dummy; real dummy;
switch (_norm) { switch (_norm) {
case full: case FULL:
v = SphericalEngine::Value<false, SphericalEngine::full, 3> v = SphericalEngine::Value<false, SphericalEngine::FULL, 3>
(_c, f, x, y, z, _a, dummy, dummy, dummy); (_c, f, x, y, z, _a, dummy, dummy, dummy);
break; break;
case schmidt: case SCHMIDT:
v = SphericalEngine::Value<false, SphericalEngine::schmidt, 3> v = SphericalEngine::Value<false, SphericalEngine::SCHMIDT, 3>
(_c, f, x, y, z, _a, dummy, dummy, dummy); (_c, f, x, y, z, _a, dummy, dummy, dummy);
break; break;
} }
return v; return v;
} }
/** /**
* Compute a spherical harmonic sum with two correction terms and its * Compute a spherical harmonic sum with two correction terms and its
* gradient. * gradient.
* *
skipping to change at line 218 skipping to change at line 224
* This is the same as the previous function, except that the component s of * This is the same as the previous function, except that the component s of
* the gradients of the sum in the \e x, \e y, and \e z directions are * the gradients of the sum in the \e x, \e y, and \e z directions are
* computed. This routine requires constant memory and thus never thro ws * computed. This routine requires constant memory and thus never thro ws
* an exception. * an exception.
********************************************************************** / ********************************************************************** /
Math::real operator()(real tau1, real tau2, real x, real y, real z, Math::real operator()(real tau1, real tau2, real x, real y, real z,
real& gradx, real& grady, real& gradz) const thro w() { real& gradx, real& grady, real& gradz) const thro w() {
real f[] = {1, tau1, tau2}; real f[] = {1, tau1, tau2};
real v = 0; real v = 0;
switch (_norm) { switch (_norm) {
case full: case FULL:
v = SphericalEngine::Value<true, SphericalEngine::full, 3> v = SphericalEngine::Value<true, SphericalEngine::FULL, 3>
(_c, f, x, y, z, _a, gradx, grady, gradz); (_c, f, x, y, z, _a, gradx, grady, gradz);
break; break;
case schmidt: case SCHMIDT:
v = SphericalEngine::Value<true, SphericalEngine::schmidt, 3> v = SphericalEngine::Value<true, SphericalEngine::SCHMIDT, 3>
(_c, f, x, y, z, _a, gradx, grady, gradz); (_c, f, x, y, z, _a, gradx, grady, gradz);
break; break;
} }
return v; return v;
} }
/** /**
* Create a CircularEngine to allow the efficient evaluation of several * Create a CircularEngine to allow the efficient evaluation of several
* points on a circle of latitude at fixed values of \e tau1 and \e tau 2. * points on a circle of latitude at fixed values of \e tau1 and \e tau 2.
* *
skipping to change at line 245 skipping to change at line 251
* @param[in] tau2 multiplier for correction coefficients \e C'' and \e S''. * @param[in] tau2 multiplier for correction coefficients \e C'' and \e S''.
* @param[in] p the radius of the circle. * @param[in] p the radius of the circle.
* @param[in] z the height of the circle above the equatorial plane. * @param[in] z the height of the circle above the equatorial plane.
* @param[in] gradp if true the returned object will be able to compute the * @param[in] gradp if true the returned object will be able to compute the
* gradient of the sum. * gradient of the sum.
* @return the CircularEngine object. * @return the CircularEngine object.
* *
* SphericalHarmonic2::operator()() exchanges the order of the sums in the * SphericalHarmonic2::operator()() exchanges the order of the sums in the
* definition, i.e., sum(n = 0..N)[sum(m = 0..n)[...]] becomes sum(m = * definition, i.e., sum(n = 0..N)[sum(m = 0..n)[...]] becomes sum(m =
* 0..N)[sum(n = m..N)[...]]. SphericalHarmonic2::Circle performs the * 0..N)[sum(n = m..N)[...]]. SphericalHarmonic2::Circle performs the
* inner sum over degree \e n (which entails about \e N<sup>2</sup> * inner sum over degree \e n (which entails about <i>N</i><sup>2</sup>
* operations). Calling CircularEngine::operator()() on the returned * operations). Calling CircularEngine::operator()() on the returned
* object performs the outer sum over the order \e m (about \e N * object performs the outer sum over the order \e m (about \e N
* operations). This routine may throw a bad_alloc exception in the * operations). This routine may throw a bad_alloc exception in the
* GeographicLib::CircularEngine constructor. * CircularEngine constructor.
* *
* See SphericalHarmonic::Circle for an example of its use. * See SphericalHarmonic::Circle for an example of its use.
********************************************************************** / ********************************************************************** /
CircularEngine Circle(real tau1, real tau2, real p, real z, bool gradp) CircularEngine Circle(real tau1, real tau2, real p, real z, bool gradp)
const { const {
real f[] = {1, tau1, tau2}; real f[] = {1, tau1, tau2};
switch (_norm) { switch (_norm) {
case full: case FULL:
return gradp ? return gradp ?
SphericalEngine::Circle<true, SphericalEngine::full, 3> SphericalEngine::Circle<true, SphericalEngine::FULL, 3>
(_c, f, p, z, _a) : (_c, f, p, z, _a) :
SphericalEngine::Circle<false, SphericalEngine::full, 3> SphericalEngine::Circle<false, SphericalEngine::FULL, 3>
(_c, f, p, z, _a); (_c, f, p, z, _a);
break; break;
case schmidt: case SCHMIDT:
default: // To avoid compiler warnings default: // To avoid compiler warnings
return gradp ? return gradp ?
SphericalEngine::Circle<true, SphericalEngine::schmidt, 3> SphericalEngine::Circle<true, SphericalEngine::SCHMIDT, 3>
(_c, f, p, z, _a) : (_c, f, p, z, _a) :
SphericalEngine::Circle<false, SphericalEngine::schmidt, 3> SphericalEngine::Circle<false, SphericalEngine::SCHMIDT, 3>
(_c, f, p, z, _a); (_c, f, p, z, _a);
break; break;
} }
} }; }
/**
* @return the zeroth SphericalEngine::coeff object.
**********************************************************************
/
const SphericalEngine::coeff& Coefficients() const throw()
{ return _c[0]; }
/**
* @return the first SphericalEngine::coeff object.
**********************************************************************
/
const SphericalEngine::coeff& Coefficients1() const throw()
{ return _c[1]; }
/**
* @return the second SphericalEngine::coeff object.
**********************************************************************
/
const SphericalEngine::coeff& Coefficients2() const throw()
{ return _c[2]; }
};
} // namespace GeographicLib } // namespace GeographicLib
#endif // GEOGRAPHICLIB_SPHERICALHARMONIC2_HPP #endif // GEOGRAPHICLIB_SPHERICALHARMONIC2_HPP
 End of changes. 28 change blocks. 
41 lines changed or deleted 66 lines changed or added


 TransverseMercator.hpp   TransverseMercator.hpp 
/** /**
* \file TransverseMercator.hpp * \file TransverseMercator.hpp
* \brief Header for GeographicLib::TransverseMercator class * \brief Header for GeographicLib::TransverseMercator class
* *
* Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co m> * Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co m>
* and licensed under the MIT/X11 License. For more information, see * and licensed under the MIT/X11 License. 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: c59b9988116b53613878bbe2 #define GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP \
65ae9fb3afa1d62f $" "$Id: 532c4a657f0355a74397cee3b87ab8f686205526 $"
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
#if !defined(TM_TX_MAXPOW) #if !defined(TM_TX_MAXPOW)
/** /**
* The order of the series approximation used in TransverseMercator. * The order of the series approximation used in TransverseMercator.
* TM_TX_MAXPOW can be set to any integer in [4, 8]. * TM_TX_MAXPOW can be set to any integer in [4, 8].
**********************************************************************/ **********************************************************************/
#define TM_TX_MAXPOW \ #define TM_TX_MAXPOW \
(GEOGRAPHICLIB_PREC == 1 ? 6 : (GEOGRAPHICLIB_PREC == 0 ? 4 : 8)) (GEOGRAPHICLIB_PREC == 1 ? 6 : (GEOGRAPHICLIB_PREC == 0 ? 4 : 8))
skipping to change at line 59 skipping to change at line 60
* TransverseMercatorExact is an alternative implementation of the projec tion * TransverseMercatorExact is an alternative implementation of the projec tion
* using exact formulas which yield accurate (to 8 nm) results over the * using exact formulas which yield accurate (to 8 nm) results over the
* entire ellipsoid. * 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 function. For
* example, the UTMUPS class applies the false easting and false northing for * example, the UTMUPS class applies the false easting and false northing for
* the UTM projections. A more complicated example is the British Nation al * the UTM projections. A more complicated example is the British Nation al
* Grid (<a href="http://www.spatialreference.org/ref/epsg/7405/"> * Grid (<a href="http://www.spatialreference.org/ref/epsg/7405/">
* EPSG:7405</a>) which requires the use of a latitude of origin. This i s * EPSG:7405</a>) which requires the use of a latitude of origin. This i s
* implemented by the GeographicLib::OSGB class. * implemented by the GeographicLib::OSGB class.
* *
* See TransverseMercator.cpp for more information on the implementation. * See TransverseMercator.cpp for more information on the implementation.
* *
* See \ref transversemercator for a discussion of this projection. * See \ref transversemercator for a discussion of this projection.
**********************************************************************/ **********************************************************************/
skipping to change at line 97 skipping to change at line 98
// Return e * atanh(e * x) for f >= 0, else return // Return e * atanh(e * x) for f >= 0, else return
// - sqrt(-e2) * atan( sqrt(-e2) * x) for f < 0 // - sqrt(-e2) * atan( sqrt(-e2) * x) for f < 0
inline real eatanhe(real x) const throw() { inline real eatanhe(real x) const throw() {
return _f >= 0 ? _e * Math::atanh(_e * x) : - _e * std::atan(_e * x); return _f >= 0 ? _e * Math::atanh(_e * x) : - _e * std::atan(_e * x);
} }
public: public:
/** /**
* Constructor for a ellipsoid with * Constructor for a ellipsoid with
* *
* @param[in] a equatorial radius (meters) * @param[in] a equatorial radius (meters).
* @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphe re. * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphe re.
* Negative \e f gives a prolate ellipsoid. If \e f > 1, set flatten ing * Negative \e f gives a prolate ellipsoid. If \e f > 1, set flatten ing
* to 1/\e f. * to 1/\e f.
* @param[in] k0 central scale factor. * @param[in] k0 central scale factor.
* *
* An exception is thrown if either of the axes of the ellipsoid or \e k0 * An exception is thrown if either of the axes of the ellipsoid or \e k0
* is not positive. * is not positive.
********************************************************************** / ********************************************************************** /
TransverseMercator(real a, real f, real k0); TransverseMercator(real a, real f, real k0);
 End of changes. 3 change blocks. 
4 lines changed or deleted 4 lines changed or added


 TransverseMercatorExact.hpp   TransverseMercatorExact.hpp 
/** /**
* \file TransverseMercatorExact.hpp * \file TransverseMercatorExact.hpp
* \brief Header for GeographicLib::TransverseMercatorExact class * \brief Header for GeographicLib::TransverseMercatorExact class
* *
* Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co m> * Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co m>
* and licensed under the MIT/X11 License. For more information, see * and licensed under the MIT/X11 License. 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: 779cb64ba88ff77e42e #define GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_HPP \
606e9b186649f78b73064 $" "$Id: 85d1d4b68c9f3d093838ba3801adbf038c4a1391 $"
#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 125 skipping to change at line 126
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 with * Constructor for a ellipsoid with
* *
* @param[in] a equatorial radius (meters) * @param[in] a equatorial radius (meters).
* @param[in] f flattening of ellipsoid. If \e f > 1, set flattening * @param[in] f flattening of ellipsoid. If \e f > 1, set flattening
* to 1/\e f. * to 1/\e f.
* @param[in] k0 central scale factor. * @param[in] k0 central scale factor.
* @param[in] extendp use extended domain. * @param[in] extendp use extended domain.
* *
* The transverse Mercator projection has a branch point singularity at \e * The transverse Mercator projection has a branch point singularity at \e
* lat = 0 and \e lon - \e lon0 = 90 (1 - \e e) or (for * lat = 0 and \e lon - \e lon0 = 90 (1 - \e e) or (for
* TransverseMercatorExact::UTM) x = 18381 km, y = 0m. The \e extendp * TransverseMercatorExact::UTM) x = 18381 km, y = 0m. The \e extendp
* argument governs where the branch cut is placed. With \e extendp = * argument governs where the branch cut is placed. With \e extendp =
* false, the "standard" convention is followed, namely the cut is plac ed * false, the "standard" convention is followed, namely the cut is plac ed
skipping to change at line 152 skipping to change at line 153
* m, \e y = 4235043.607933 m both map to \e lat = -2 deg, \e lon = 88 deg. * 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 * - <i>x</i>/(\e k0 \e a) in [0, inf) and
* \e y/(\e k0 \e a) in [0, E(\e e^2)] * <i>y</i>/(\e k0 \e a) in [0, E(<i>e</i><sup>2</sup>)]
* - \e x/(\e k0 \e a) in [K(1 - \e e^2) - E(1 - \e e^2), inf) and * - <i>x</i>/(\e k0 \e a) in [K(1 - <i>e</i><sup>2</sup>) - E(1 -
* \e y/(\e k0 \e a) in (-inf, 0] * <i>e</i><sup>2</sup>), inf) and <i>y</i>/(\e k0 \e a) in (-inf,
0]
* . * .
* See Sec. 5 of * See Sec. 5 of
* <a href="http://arxiv.org/abs/1002.1417">arXiv:1002.1417</a> for a f ull * <a href="http://arxiv.org/abs/1002.1417">arXiv:1002.1417</a> for a f ull
* discussion of the treatment of the branch cut. * discussion of the treatment of the branch cut.
* *
* The method will work for all ellipsoids used in terrestial geodesy. * The method will work for all ellipsoids used in terrestrial geodesy.
The * The method cannot be applied directly to the case of a sphere (\e f
* method cannot be applied directly to the case of a sphere (\e f = 0) = 0)
* 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 f should be larger than about numeric_lim its< * limit, and in practice, \e f should be larger than about numeric_lim its<
* real >::%epsilon(). However, TransverseMercator treats the sphere * real >::%epsilon(). However, TransverseMercator treats the sphere
* exactly. An exception is thrown if either axis of the ellipsoid or \e * exactly. An exception is thrown if either axis of the ellipsoid or \e
* k0 is not positive or if \e f <= 0. * k0 is not positive or if \e f <= 0.
********************************************************************** / ********************************************************************** /
TransverseMercatorExact(real a, real f, real k0, bool extendp = false); TransverseMercatorExact(real a, real f, real k0, bool extendp = false);
/** /**
* Forward projection, from geographic to transverse Mercator. * Forward projection, from geographic to transverse Mercator.
* *
* @param[in] lon0 central meridian of the projection (degrees). * @param[in] lon0 central meridian of the projection (degrees).
 End of changes. 5 change blocks. 
11 lines changed or deleted 12 lines changed or added


 UTMUPS.hpp   UTMUPS.hpp 
/** /**
* \file UTMUPS.hpp * \file UTMUPS.hpp
* \brief Header for GeographicLib::UTMUPS class * \brief Header for GeographicLib::UTMUPS class
* *
* Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co m> * Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.co m>
* and licensed under the MIT/X11 License. For more information, see * and licensed under the MIT/X11 License. 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: 7973c2fa261a964b7ccdccd2057a08612c1d #define GEOGRAPHICLIB_UTMUPS_HPP \
4cd5 $" "$Id: bfc2e2a500d4adf0cae87cbeaf4e3753e20b3019 $"
#include <sstream> #include <sstream>
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
namespace GeographicLib { namespace GeographicLib {
/** /**
* \brief Convert between Geographic coordinates and UTM/UPS * \brief Convert between Geographic coordinates and UTM/UPS
* *
* UTM and UPS are defined * UTM and UPS are defined
skipping to change at line 172 skipping to change at line 173
* @param[out] zone the UTM zone (zero means UPS). * @param[out] zone the UTM zone (zero means UPS).
* @param[out] northp hemisphere (true means north, false means south). * @param[out] northp hemisphere (true means north, false means south).
* @param[out] x easting of point (meters). * @param[out] x easting of point (meters).
* @param[out] y northing of point (meters). * @param[out] y northing of point (meters).
* @param[out] gamma meridian convergence at point (degrees). * @param[out] gamma meridian convergence at point (degrees).
* @param[out] k scale of projection at point. * @param[out] k scale of projection at point.
* @param[in] setzone zone override. * @param[in] setzone zone override.
* @param[in] mgrslimits if true enforce the stricter MGRS limits on th e * @param[in] mgrslimits if true enforce the stricter MGRS limits on th e
* coordinates (default = false). * coordinates (default = false).
* *
* The prefered zone for the result can be specified with \e setzone, s ee * The preferred zone for the result can be specified with \e setzone, see
* 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. This also returns meridian convergence \e gamma * are unchanged. This also returns meridian convergence \e gamma
* (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);
skipping to change at line 194 skipping to change at line 195
* Reverse projection, from UTM/UPS to geographic. * Reverse projection, from UTM/UPS to geographic.
* *
* @param[in] zone the UTM zone (zero means UPS). * @param[in] zone the UTM zone (zero means UPS).
* @param[in] northp hemisphere (true means north, false means south). * @param[in] northp hemisphere (true means north, false means south).
* @param[in] x easting of point (meters). * @param[in] x easting of point (meters).
* @param[in] y northing of point (meters). * @param[in] y northing of point (meters).
* @param[out] lat latitude of point (degrees). * @param[out] lat latitude of point (degrees).
* @param[out] lon longitude of point (degrees). * @param[out] lon longitude of point (degrees).
* @param[out] gamma meridian convergence at point (degrees). * @param[out] gamma meridian convergence at point (degrees).
* @param[out] k scale of projection at point. * @param[out] k scale of projection at point.
* @param[in] mgrslimits if true enforce the stricted MGRS limits on th e * @param[in] mgrslimits if true enforce the stricter MGRS limits on th e
* coordinates (default = false). * coordinates (default = false).
* *
* Throw error if easting or northing is outside the allowed range (see * Throw error if easting or northing is outside the allowed range (see
* below), in which case the arguments are unchanged. The accuracy of the * below), in which case the arguments are unchanged. The accuracy of the
* conversion is about 5nm. * 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
skipping to change at line 246 skipping to change at line 247
********************************************************************** / ********************************************************************** /
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. * Decode a UTM/UPS zone string.
* *
* @param[in] zonestr string represention of zone and hemisphere. * @param[in] zonestr string representation of zone and hemisphere.
* @param[out] zone the UTM zone (zero means UPS). * @param[out] zone the UTM zone (zero means UPS).
* @param[out] northp hemisphere (true means north, false means south). * @param[out] northp hemisphere (true means north, false means south).
* *
* For UTM, \e zonestr has the form of a zone number in the range * For UTM, \e zonestr has the form of a zone number in the range
* [UTMUPS::MINUTMZONE, UTMUPS::MAXUTMZONE] = [1, 60] followed by a * [UTMUPS::MINUTMZONE, UTMUPS::MAXUTMZONE] = [1, 60] followed by a
* hemisphere letter, N or S. For UPS, it consists just of the hemisph ere * hemisphere letter, N or S. For UPS, it consists just of the hemisph ere
* letter. The returned value of \e zone is UTMUPS::UPS = 0 for UPS. Note * letter. The returned value of \e zone is UTMUPS::UPS = 0 for UPS. Note
* well that "38S" indicates the southern hemisphere of zone 38 and not * 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, * latitude band S, [32, 40]. N, 01S, 2N, 38S are legal. 0N, 001S, 61 N,
* 38P are illegal. INV is a special value for which the returned valu e of * 38P are illegal. INV is a special value for which the returned valu e of
* \e is UTMUPS::INVALID. Throws an error is the zone string is malfor med. * \e is UTMUPS::INVALID. Throws an error is the zone string is malfor med.
********************************************************************** / ********************************************************************** /
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. * Encode a UTM/UPS zone string.
* *
* @param[out] zone the UTM zone (zero means UPS). * @param[out] zone the UTM zone (zero means UPS).
* @param[out] northp hemisphere (true means north, false means south). * @param[out] northp hemisphere (true means north, false means south).
* @return string represention of zone and hemisphere. * @return string representation of zone and hemisphere.
* *
* \e zone must be in the range [UTMUPS::MINZONE, UTMUPS::MAXZONE] = [0 , * \e zone must be in the range [UTMUPS::MINZONE, UTMUPS::MAXZONE] = [0 ,
* 60] with \e zone = UTMUPS::UPS, 0, indicating UPS (but the resulting * 60] with \e zone = UTMUPS::UPS, 0, indicating UPS (but the resulting
* string does not contain "0"). \e zone may also be UTMUPS::INVALID, in * string does not contain "0"). \e zone may also be UTMUPS::INVALID, in
* which case the returned string is "INV". This reverses * which case the returned string is "INV". This reverses
* UTMUPS::DecodeZone. * UTMUPS::DecodeZone.
********************************************************************** / ********************************************************************** /
static std::string EncodeZone(int zone, bool northp); static std::string EncodeZone(int zone, bool northp);
/** /**
 End of changes. 5 change blocks. 
6 lines changed or deleted 6 lines changed or added


 Utility.hpp   Utility.hpp 
/** /**
* \file Utility.hpp * \file Utility.hpp
* \brief Header for GeographicLib::Utility class * \brief Header for GeographicLib::Utility class
* *
* Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed un der * Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed un der
* the MIT/X11 License. For more information, see * the MIT/X11 License. For more information, see
* http://geographiclib.sourceforge.net/ * http://geographiclib.sourceforge.net/
**********************************************************************/ **********************************************************************/
#if !defined(GEOGRAPHICLIB_UTILITY_HPP) #if !defined(GEOGRAPHICLIB_UTILITY_HPP)
#define GEOGRAPHICLIB_UTILITY_HPP "$Id: 5930e0a4f55b8a9478bfca81bdc962c7f8e #define GEOGRAPHICLIB_UTILITY_HPP \
5d26f $" "$Id: a4a9a65ef1b4dc72579b2310257cde0d626480a5 $"
#include <GeographicLib/Constants.hpp> #include <GeographicLib/Constants.hpp>
#include <iomanip> #include <iomanip>
#include <vector> #include <vector>
#include <string> #include <string>
#include <sstream> #include <sstream>
#include <algorithm> #include <algorithm>
#include <cctype> #include <cctype>
namespace GeographicLib { namespace GeographicLib {
skipping to change at line 374 skipping to change at line 375
* @param[out] array the output array of type IntT (internal). * @param[out] array the output array of type IntT (internal).
* @param[in] num the size of the array. * @param[in] num the size of the array.
********************************************************************** / ********************************************************************** /
template<typename ExtT, typename IntT, bool bigendp> template<typename ExtT, typename IntT, bool bigendp>
static inline void readarray(std::istream& str, static inline void readarray(std::istream& str,
IntT array[], size_t num) { IntT array[], size_t num) {
if (sizeof(IntT) == sizeof(ExtT) && if (sizeof(IntT) == sizeof(ExtT) &&
std::numeric_limits<IntT>::is_integer == std::numeric_limits<IntT>::is_integer ==
std::numeric_limits<ExtT>::is_integer) { std::numeric_limits<ExtT>::is_integer) {
// Data is compatible (aside from the issue of endian-ness). // Data is compatible (aside from the issue of endian-ness).
str.read(reinterpret_cast<char *>(array), num * sizeof(IntT)); str.read(reinterpret_cast<char *>(array), num * sizeof(ExtT));
if (!str.good()) if (!str.good())
throw GeographicErr("Failure reading data"); throw GeographicErr("Failure reading data");
if (bigendp != Math::bigendian) { // endian mismatch -> swap bytes if (bigendp != Math::bigendian) { // endian mismatch -> swap bytes
for (int i = num; i--;) for (int i = num; i--;)
array[i] = Math::swab<IntT>(array[i]); array[i] = Math::swab<IntT>(array[i]);
} }
} else { } else {
const int bufsize = 1024; // read this many values at a time const int bufsize = 1024; // read this many values at a time
ExtT buffer[bufsize]; // temporary buffer ExtT buffer[bufsize]; // temporary buffer
int k = int(num); // data values left to read int k = int(num); // data values left to read
int i = 0; // index into output array int i = 0; // index into output array
while (k) { while (k) {
int num = (std::min)(k, bufsize); int num = (std::min)(k, bufsize);
str.read(reinterpret_cast<char *>(buffer), num * sizeof(IntT)); str.read(reinterpret_cast<char *>(buffer), num * sizeof(ExtT));
if (!str.good()) if (!str.good())
throw GeographicErr("Failure reading data"); throw GeographicErr("Failure reading data");
for (int j = 0; j < num; ++j) for (int j = 0; j < num; ++j)
// fix endian-ness and cast to IntT // fix endian-ness and cast to IntT
array[i++] = IntT(bigendp == Math::bigendian ? buffer[j] : array[i++] = IntT(bigendp == Math::bigendian ? buffer[j] :
Math::swab<IntT>(buffer[j])); Math::swab<IntT>(buffer[j]));
k -= num; k -= num;
} }
} }
return; return;
skipping to change at line 418 skipping to change at line 419
* @tparam bigendp true if the external storage format is big-endian. * @tparam bigendp true if the external storage format is big-endian.
* @param[in] str the input stream containing the data of type ExtT * @param[in] str the input stream containing the data of type ExtT
* (external). * (external).
* @param[out] array the output vector of type IntT (internal). * @param[out] array the output vector of type IntT (internal).
********************************************************************** / ********************************************************************** /
template<typename ExtT, typename IntT, bool bigendp> template<typename ExtT, typename IntT, bool bigendp>
static inline void readarray(std::istream& str, static inline void readarray(std::istream& str,
std::vector<IntT>& array) { std::vector<IntT>& array) {
readarray<ExtT, IntT, bigendp>(str, &array[0], array.size()); readarray<ExtT, IntT, bigendp>(str, &array[0], array.size());
} }
/**
* Parse a KEY VALUE line.
*
* @param[in] line the input line.
* @param[out] key the key.
* @param[out] val the value.
* @return whether a key was found.
*
* A # character and everything after it are discarded. If the results
is
* just white space, the routine returns false (and \e key and \e val a
re
* not set). Otherwise the first token is taken to be the key and the
rest
* of the line (trimmed of leading and trailing white space) is the val
ue.
**********************************************************************
/
static bool ParseLine(const std::string& line,
std::string& key, std::string& val);
}; };
} // namespace GeographicLib } // namespace GeographicLib
#endif // GEOGRAPHICLIB_UTILITY_HPP #endif // GEOGRAPHICLIB_UTILITY_HPP
 End of changes. 4 change blocks. 
4 lines changed or deleted 26 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/