CassiniSoldner.hpp | CassiniSoldner.hpp | |||
---|---|---|---|---|
/** | /** | |||
* \file CassiniSoldner.hpp | * \file CassiniSoldner.hpp | |||
* \brief Header for GeographicLib::CassiniSoldner class | * \brief Header for GeographicLib::CassiniSoldner class | |||
* | * | |||
* Copyright (c) Charles Karney (2009, 2010) <charles@karney.com> | * Copyright (c) Charles Karney (2009, 2010) <charles@karney.com> | |||
* and licensed under the LGPL. For more information, see | * and licensed under the LGPL. For more information, see | |||
* http://geographiclib.sourceforge.net/ | * http://geographiclib.sourceforge.net/ | |||
**********************************************************************/ | **********************************************************************/ | |||
#if !defined(GEOGRAPHICLIB_CASSINISOLDNER_HPP) | #if !defined(GEOGRAPHICLIB_CASSINISOLDNER_HPP) | |||
#define GEOGRAPHICLIB_CASSINISOLDNER_HPP "$Id: CassiniSoldner.hpp 6785 2010 -01-05 22:15:42Z karney $" | #define GEOGRAPHICLIB_CASSINISOLDNER_HPP "$Id: CassiniSoldner.hpp 6827 2010 -05-20 19:56:18Z karney $" | |||
#include "GeographicLib/Geodesic.hpp" | #include "GeographicLib/Geodesic.hpp" | |||
#include "GeographicLib/Constants.hpp" | #include "GeographicLib/Constants.hpp" | |||
namespace GeographicLib { | namespace GeographicLib { | |||
/** | /** | |||
* \brief Cassini-Soldner Projection. | * \brief Cassini-Soldner Projection. | |||
* | * | |||
* Cassini-Soldner projection centered at an arbitrary position, \e lat0, \e | * Cassini-Soldner projection centered at an arbitrary position, \e lat0, \e | |||
skipping to change at line 68 | skipping to change at line 68 | |||
* geodesics between that point and, respectively, (\e lat1, \e lon1) and (\e | * geodesics between that point and, respectively, (\e lat1, \e lon1) and (\e | |||
* lat, \e lon). | * lat, \e lon). | |||
**********************************************************************/ | **********************************************************************/ | |||
class CassiniSoldner { | class CassiniSoldner { | |||
private: | private: | |||
typedef Math::real real; | typedef Math::real real; | |||
const Geodesic _earth; | const Geodesic _earth; | |||
GeodesicLine _meridian; | GeodesicLine _meridian; | |||
real _sbet0, _cbet0; | real _sbet0, _cbet0; | |||
real Scale(const GeodesicLine& perp, real sig12) const throw(); | ||||
static const real eps1, eps2; | static const real eps1, eps2; | |||
static const unsigned maxit = 10; | static const unsigned maxit = 10; | |||
static inline real sq(real x) throw() { return x * x; } | static inline real sq(real x) throw() { return x * x; } | |||
// The following private helper functions are copied from Geodesic. | ||||
static inline real AngNormalize(real x) throw() { | ||||
// Place angle in [-180, 180). Assumes x is in [-540, 540). | ||||
return x >= 180 ? x - 360 : x < -180 ? x + 360 : x; | ||||
} | ||||
static inline real AngRound(real x) throw() { | ||||
// The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^ | ||||
57 | ||||
// for reals = 0.7 pm on the earth if x is an angle in degrees. (Thi | ||||
s | ||||
// is about 1000 times more resolution than we get with angles around | ||||
90 | ||||
// degrees.) We use this to avoid having to deal with near singular | ||||
// cases when x is non-zero but tiny (e.g., 1.0e-200). | ||||
const real z = real(0.0625); // 1/16 | ||||
volatile real y = std::abs(x); | ||||
// The compiler mustn't "simplify" z - (z - y) to y | ||||
y = y < z ? z - (z - y) : y; | ||||
return x < 0 ? -y : y; | ||||
} | ||||
static inline void SinCosNorm(real& sinx, real& cosx) throw() { | ||||
real r = Math::hypot(sinx, cosx); | ||||
sinx /= r; | ||||
cosx /= r; | ||||
} | ||||
public: | public: | |||
/** | /** | |||
* Constructor for CassiniSoldner setting the Geodesic object, \e earth , to | * Constructor for CassiniSoldner setting the Geodesic object, \e earth , to | |||
* use for geodesic calculations. By default this uses the WGS84 | * use for geodesic calculations. By default this uses the WGS84 | |||
* ellipsoid. This constructor makes an "uninitialized" object. Call | * ellipsoid. This constructor makes an "uninitialized" object. Call | |||
* Reset to set the central latitude and longuitude, prior to calling | * Reset to set the central latitude and longuitude, prior to calling | |||
* Forward and Reverse. | * Forward and Reverse. | |||
********************************************************************** / | ********************************************************************** / | |||
explicit CassiniSoldner(const Geodesic& earth = Geodesic::WGS84) throw( ) | explicit CassiniSoldner(const Geodesic& earth = Geodesic::WGS84) throw( ) | |||
End of changes. 3 change blocks. | ||||
2 lines changed or deleted | 26 lines changed or added | |||
Constants.hpp | Constants.hpp | |||
---|---|---|---|---|
/** | /** | |||
* \file Constants.hpp | * \file Constants.hpp | |||
* \brief Header for GeographicLib::Constants class | * \brief Header for GeographicLib::Constants class | |||
* | * | |||
* Copyright (c) Charles Karney (2008, 2009, 2010) <charles@karney.com> | * Copyright (c) Charles Karney (2008, 2009, 2010) <charles@karney.com> | |||
* and licensed under the LGPL. For more information, see | * and licensed under the LGPL. For more information, see | |||
* http://geographiclib.sourceforge.net/ | * http://geographiclib.sourceforge.net/ | |||
**********************************************************************/ | **********************************************************************/ | |||
#if !defined(GEOGRAPHICLIB_CONSTANTS_HPP) | #if !defined(GEOGRAPHICLIB_CONSTANTS_HPP) | |||
#define GEOGRAPHICLIB_CONSTANTS_HPP "$Id: Constants.hpp 6817 2010-02-08 14: 42:35Z karney $" | #define GEOGRAPHICLIB_CONSTANTS_HPP "$Id: Constants.hpp 6827 2010-05-20 19: 56:18Z karney $" | |||
/** | /** | |||
* A simple compile-time assert. This is designed to be compatible with th e | * A simple compile-time assert. This is designed to be compatible with th e | |||
* C++0X static_assert. | * C++0X static_assert. | |||
**********************************************************************/ | **********************************************************************/ | |||
#if !defined(STATIC_ASSERT) | #if !defined(STATIC_ASSERT) | |||
#define STATIC_ASSERT(cond,reason) { enum{ STATIC_ASSERT_ENUM=1/int(cond) } ; } | #define STATIC_ASSERT(cond,reason) { enum{ STATIC_ASSERT_ENUM=1/int(cond) } ; } | |||
#endif | #endif | |||
#if defined(__GNUC__) | #if defined(__GNUC__) | |||
skipping to change at line 114 | skipping to change at line 114 | |||
static inline double hypot(double x, double y) throw() | static inline double hypot(double x, double y) throw() | |||
{ return _hypot(x, y); } | { return _hypot(x, y); } | |||
static inline float hypot(float x, float y) throw() | static inline float hypot(float x, float y) throw() | |||
{ return _hypotf(x, y); } | { return _hypotf(x, y); } | |||
#else | #else | |||
// Use overloading to define generic versions | // Use overloading to define generic versions | |||
static inline double hypot(double x, double y) throw() | static inline double hypot(double x, double y) throw() | |||
{ return ::hypot(x, y); } | { return ::hypot(x, y); } | |||
static inline float hypot(float x, float y) throw() | static inline float hypot(float x, float y) throw() | |||
{ return ::hypotf(x, y); } | { return ::hypotf(x, y); } | |||
#if !defined(__NO_LONG_DOUBLE_MATH) | ||||
static inline long double hypot(long double x, long double y) throw() | static inline long double hypot(long double x, long double y) throw() | |||
{ return ::hypotl(x, y); } | { return ::hypotl(x, y); } | |||
#endif | #endif | |||
#endif | ||||
#if defined(DOXYGEN) || defined(_MSC_VER) | #if defined(DOXYGEN) || defined(_MSC_VER) | |||
/** | /** | |||
* Return exp(\e x) - 1 accurate near \e x = 0. This is taken from | * Return 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. | |||
********************************************************************** / | ********************************************************************** / | |||
static inline real expm1(real x) throw() { | static inline real expm1(real x) throw() { | |||
volatile real | volatile real | |||
y = std::exp(x), | y = std::exp(x), | |||
z = y - 1; | z = y - 1; | |||
// The reasoning here is similar to that for log1p. The expression | // The reasoning here is similar to that for log1p. The expression | |||
// mathematically reduces to exp(x) - 1, and the factor z/log(y) = (y - | // mathematically reduces to exp(x) - 1, and the factor z/log(y) = (y - | |||
// 1)/log(y) is a slowly varying quantity near y = 1 and is accuratel y | // 1)/log(y) is a slowly varying quantity near y = 1 and is accuratel y | |||
// computed. | // computed. | |||
return std::abs(x) > 1 ? z : z == 0 ? x : x * z / std::log(y); | return std::abs(x) > 1 ? z : z == 0 ? x : x * z / std::log(y); | |||
} | } | |||
#else | #else | |||
static inline double expm1(double x) throw() { return ::expm1(x); } | static inline double expm1(double x) throw() { return ::expm1(x); } | |||
static inline float expm1(float x) throw() { return ::expm1f(x); } | static inline float expm1(float x) throw() { return ::expm1f(x); } | |||
#if !defined(__NO_LONG_DOUBLE_MATH) | ||||
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 | ||||
#if defined(DOXYGEN) || defined(_MSC_VER) | #if defined(DOXYGEN) || defined(_MSC_VER) | |||
/** | /** | |||
* Return log(\e x + 1) accurate near \e x = 0. This is taken See | * Return log(\e x + 1) accurate near \e x = 0. This is taken See | |||
* D. Goldberg, | * D. Goldberg, | |||
* <a href="http://docs.sun.com/source/806-3568/ncg_goldberg.html"> Wha t | * <a href="http://docs.sun.com/source/806-3568/ncg_goldberg.html"> Wha t | |||
* every computer scientist should know about floating-point arithmetic </a> | * every computer scientist should know about floating-point arithmetic </a> | |||
* (1991), Theorem 4. See also, Higham (op. cit.), Answer to Problem 1 .5, | * (1991), Theorem 4. See also, Higham (op. cit.), Answer to Problem 1 .5, | |||
* p 528. | * p 528. | |||
********************************************************************** / | ********************************************************************** / | |||
skipping to change at line 163 | skipping to change at line 167 | |||
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; | |||
} | } | |||
#else | #else | |||
static inline double log1p(double x) throw() { return ::log1p(x); } | static inline double log1p(double x) throw() { return ::log1p(x); } | |||
static inline float log1p(float x) throw() { return ::log1pf(x); } | static inline float log1p(float x) throw() { return ::log1pf(x); } | |||
#if !defined(__NO_LONG_DOUBLE_MATH) | ||||
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 | ||||
#if defined(DOXYGEN) || defined(_MSC_VER) | #if defined(DOXYGEN) || defined(_MSC_VER) | |||
/** | /** | |||
* Return asinh(\e x). This is defined in terms of log1p(\e x) in orde | * Return asinh(\e x). This is defined in terms of Math::log1p(\e x) i | |||
r to | n | |||
* maintain accuracy near \e x = 0. In addition, the odd parity of the | * order to maintain accuracy near \e x = 0. In addition, the odd pari | |||
* function is enforced. | ty | |||
* of the function is enforced. | ||||
********************************************************************** / | ********************************************************************** / | |||
static inline real asinh(real x) throw() { | static inline real asinh(real x) throw() { | |||
real y = std::abs(x); // Enforce odd parity | real y = std::abs(x); // Enforce odd parity | |||
y = log1p(y * (1 + y/(hypot(real(1), y) + 1))); | y = log1p(y * (1 + y/(hypot(real(1), y) + 1))); | |||
return x < 0 ? -y : y; | return x < 0 ? -y : y; | |||
} | } | |||
#else | #else | |||
static inline double asinh(double x) throw() { return ::asinh(x); } | static inline double asinh(double x) throw() { return ::asinh(x); } | |||
static inline float asinh(float x) throw() { return ::asinhf(x); } | static inline float asinh(float x) throw() { return ::asinhf(x); } | |||
#if !defined(__NO_LONG_DOUBLE_MATH) | ||||
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 | ||||
#if defined(DOXYGEN) || defined(_MSC_VER) | #if defined(DOXYGEN) || defined(_MSC_VER) | |||
/** | /** | |||
* Return atanh(\e x). This is defined in terms of log1p(\e x) in orde | * Return atanh(\e x). This is defined in terms of Math::log1p(\e x) i | |||
r to | n | |||
* maintain accuracy near \e x = 0. In addition, the odd parity of the | * order to maintain accuracy near \e x = 0. In addition, the odd pari | |||
* function is enforced. | ty | |||
* of the function is enforced. | ||||
********************************************************************** / | ********************************************************************** / | |||
static inline real atanh(real x) throw() { | static inline real atanh(real x) throw() { | |||
real y = std::abs(x); // Enforce odd parity | real y = std::abs(x); // Enforce odd parity | |||
y = log1p(2 * y/(1 - y))/2; | y = log1p(2 * y/(1 - y))/2; | |||
return x < 0 ? -y : y; | return x < 0 ? -y : y; | |||
} | } | |||
#else | #else | |||
static inline double atanh(double x) throw() { return ::atanh(x); } | static inline double atanh(double x) throw() { return ::atanh(x); } | |||
static inline float atanh(float x) throw() { return ::atanhf(x); } | static inline float atanh(float x) throw() { return ::atanhf(x); } | |||
#if !defined(__NO_LONG_DOUBLE_MATH) | ||||
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 | ||||
#if defined(DOXYGEN) || defined(_MSC_VER) | #if defined(DOXYGEN) || defined(_MSC_VER) | |||
/** | /** | |||
* Return the real cube root of \e x. | * Return the real cube root of \e x. | |||
********************************************************************** / | ********************************************************************** / | |||
static inline real cbrt(real x) throw() { | static inline real cbrt(real x) throw() { | |||
real y = std::pow(std::abs(x), 1/real(3)); // Return the real cube ro ot | real y = std::pow(std::abs(x), 1/real(3)); // Return the real cube ro ot | |||
return x < 0 ? -y : y; | return x < 0 ? -y : y; | |||
} | } | |||
#else | #else | |||
static inline double cbrt(double x) throw() { return ::cbrt(x); } | static inline double cbrt(double x) throw() { return ::cbrt(x); } | |||
static inline float cbrt(float x) throw() { return ::cbrtf(x); } | static inline float cbrt(float x) throw() { return ::cbrtf(x); } | |||
#if !defined(__NO_LONG_DOUBLE_MATH) | ||||
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 | ||||
#if defined(DOXYGEN) | ||||
/** | /** | |||
* Return true if number is finite, false if NaN or infinite. | * Return true if number is finite, false if NaN or infinite. | |||
********************************************************************** / | ********************************************************************** / | |||
static inline bool isfinite(real x) throw() { | static inline bool isfinite(real x) throw() { | |||
return std::abs(x) < std::numeric_limits<real>::max(); | #if defined(DOXYGEN) | |||
} | return std::abs(x) <= std::numeric_limits<real>::max(); | |||
#elif defined(_MSC_VER) | #elif defined(_MSC_VER) | |||
static inline real isfinite(real x) throw() { return _finite(x); } | return _finite(x) != 0; | |||
#else | #else | |||
static inline real isfinite(real x) throw() { return std::isfinite(x); } | return std::isfinite(x) != 0; | |||
#endif | #endif | |||
} | ||||
}; | }; | |||
/** | /** | |||
* \brief %Constants needed by %GeographicLib | * \brief %Constants needed by %GeographicLib | |||
* | * | |||
* Define constants specifying the WGS84 ellipsoid, the UTM and UPS | * Define constants specifying the WGS84 ellipsoid, the UTM and UPS | |||
* projections, and various unit conversions. | * projections, and various unit conversions. | |||
**********************************************************************/ | **********************************************************************/ | |||
class Constants { | class Constants { | |||
private: | private: | |||
End of changes. 20 change blocks. | ||||
14 lines changed or deleted | 28 lines changed or added | |||
DMS.hpp | DMS.hpp | |||
---|---|---|---|---|
/** | /** | |||
* \file DMS.hpp | * \file DMS.hpp | |||
* \brief Header for GeographicLib::DMS class | * \brief Header for GeographicLib::DMS class | |||
* | * | |||
* Copyright (c) Charles Karney (2008, 2009, 2010) <charles@karney.com> | * Copyright (c) Charles Karney (2008, 2009, 2010) <charles@karney.com> | |||
* and licensed under the LGPL. For more information, see | * and licensed under the LGPL. For more information, see | |||
* http://geographiclib.sourceforge.net/ | * http://geographiclib.sourceforge.net/ | |||
**********************************************************************/ | **********************************************************************/ | |||
#if !defined(GEOGRAPHICLIB_DMS_HPP) | #if !defined(GEOGRAPHICLIB_DMS_HPP) | |||
#define GEOGRAPHICLIB_DMS_HPP "$Id: DMS.hpp 6785 2010-01-05 22:15:42Z karne y $" | #define GEOGRAPHICLIB_DMS_HPP "$Id: DMS.hpp 6827 2010-05-20 19:56:18Z karne y $" | |||
#include "GeographicLib/Constants.hpp" | #include "GeographicLib/Constants.hpp" | |||
#include <sstream> | #include <sstream> | |||
#include <iomanip> | ||||
namespace GeographicLib { | namespace GeographicLib { | |||
/** | /** | |||
* \brief Convert between degrees and %DMS representation. | * \brief Convert between degrees and %DMS representation. | |||
* | * | |||
* Parse a string representing degree, minutes, and seconds and return th e | * Parse a string representing degree, minutes, and seconds and return th e | |||
* angle in degrees and format an angle in degrees as degree, minutes, an d | * angle in degrees and format an angle in degrees as degree, minutes, an d | |||
* seconds. | * seconds. | |||
**********************************************************************/ | **********************************************************************/ | |||
skipping to change at line 46 | skipping to change at line 47 | |||
static const std::string signs; | static const std::string signs; | |||
static const std::string digits; | static const std::string digits; | |||
static const std::string dmsindicators; | static const std::string dmsindicators; | |||
static const std::string components[3]; | static const std::string components[3]; | |||
DMS(); // Disable constructor | DMS(); // Disable constructor | |||
public: | public: | |||
/** | /** | |||
* Indicator for presence of hemisphere indicator (N/S/E/W) on latitude s | * Indicator for presence of hemisphere indicator (N/S/E/W) on latitude s | |||
* and longitudes. AZIMUTH is used in Encode to indicate output in [00 | * and longitudes. DMS::AZIMUTH is used in Encode to indicate output i | |||
0, | n | |||
* 360) with no letter indicator. | * [000, 360) with no letter indicator. DMS::NUMBER is used in Encode | |||
to | ||||
* indicate the output of a plain number. | ||||
********************************************************************** / | ********************************************************************** / | |||
enum flag { NONE = 0, LATITUDE = 1, LONGITUDE = 2, AZIMUTH = 3 }; | enum flag { NONE = 0, LATITUDE = 1, LONGITUDE = 2, AZIMUTH = 3, | |||
NUMBER = 4, }; | ||||
/** | /** | |||
* Indicator for trailing units on an angle. | * Indicator for trailing units on an angle. | |||
********************************************************************** / | ********************************************************************** / | |||
enum component { DEGREE = 0, MINUTE = 1, SECOND = 2 }; | enum component { DEGREE = 0, MINUTE = 1, SECOND = 2 }; | |||
/** | /** | |||
* Read a string \e dms in DMS format and return the resulting angle in | * Read a string \e dms in DMS format and return the resulting angle in | |||
* degrees. Degrees, minutes, and seconds are indicated by the letters d, | * degrees. Degrees, minutes, and seconds are indicated by the letters d, | |||
* ', ", and these components may only be given in this order. Any (bu | * ', ", and these components may only be given in this order. An | |||
t | y | |||
* not all) components may be omitted. The last component indicator ma | * (but not all) components may be omitted. The last component indicat | |||
y be | or | |||
* omitted and is assumed to be tbe next smallest unit (thus 33d10 is | * may be omitted and is assumed to be tbe next smallest unit (thus 33d | |||
* interpreted as 33d10'). The final component may be a decimal fracti | 10 | |||
on | * is interpreted as 33d10'). The final component may be a decimal | |||
* but the non-final components must be integers. The integer parts of | * fraction but the non-final components must be integers. The integer | |||
the | * parts of the minutes and seconds components must be less than 60. A | |||
* minutes and seconds components must be less than 60. A single leadi | * single leading sign is permitted. A hemisphere designator (N, E, W, | |||
ng | S) | |||
* sign is permitted. A hemisphere designator (N, E, W, S) may be adde | * may be added to tbe beginning or end of the string. The result is | |||
d to | * multiplied by the implied signed of the hemisphere designator (negat | |||
* tbe beginning or end of the string. The result is multiplied by the | ive | |||
* implied signed of the hemisphere designator (negative for S and W). | * for S and W). In addition \e flag is used to indicate whether such | |||
In | a | |||
* addition \e flag is used to indicate whether such a designator was f | * designator was found and whether it implies that the angle is a lati | |||
ound | tude | |||
* and whether it implies that the angle is a latitude (N or S) or | * (N or S) or longitude (E or W). Throws an error on a malformed stri | |||
* longitude (E or W). Throws an error on a malformed string. No chec | ng. | |||
k is | * No check is performed on the range of the result. | |||
* performed on the range of the result. | ||||
********************************************************************** / | ********************************************************************** / | |||
static Math::real Decode(const std::string& dms, flag& ind); | static Math::real Decode(const std::string& dms, flag& ind); | |||
/** | /** | |||
* Convert real degrees \e d, minutes \e m, and seconds \e s, to degree s. | * Convert real degrees \e d, minutes \e m, and seconds \e s, to degree s. | |||
* This does not propagate the sign on \e d to the other components, so | * This does not propagate the sign on \e d to the other components, so | |||
* -3d20' would need to be represented as - DMS::Decode(3.0, 20.0) or | * -3d20' would need to be represented as - DMS::Decode(3.0, 20.0) or | |||
* DMS::Decode(-3.0, -20.0). | * DMS::Decode(-3.0, -20.0). | |||
********************************************************************** / | ********************************************************************** / | |||
static Math::real Decode(real d, real m = 0, real s = 0) throw() | static Math::real Decode(real d, real m = 0, real s = 0) throw() | |||
{ return d + (m + s/real(60))/real(60); } | { return d + (m + s/real(60))/real(60); } | |||
/** | /** | |||
* Convert a string \e str to a real number. | ||||
********************************************************************** | ||||
/ | ||||
static Math::real Decode(const std::string& str); | ||||
/** | ||||
* Convert two strings \e dmsa and \e dmsb to a latitude, \e lat, and | * Convert two strings \e dmsa and \e dmsb to a latitude, \e lat, and | |||
* longitude, \e lon. By default, the \e lat (resp., \e lon) is assign ed | * longitude, \e lon. By default, the \e lat (resp., \e lon) is assign ed | |||
* to the results of decoding \e dmsa (resp., \e dmsb). However this i s | * to the results of decoding \e dmsa (resp., \e dmsb). However this i s | |||
* overridden if either \e dmsa or \e dmsb contain a latitude or longit ude | * overridden if either \e dmsa or \e dmsb contain a latitude or longit ude | |||
* hemisphere designator (N, S, E, W). Throws an error if the decoded | * hemisphere designator (N, S, E, W). Throws an error if the decoded | |||
* numbers are out of the ranges [-90<sup>o</sup>, 90<sup>o</sup>] for | * numbers are out of the ranges [-90<sup>o</sup>, 90<sup>o</sup>] for | |||
* latitude and [-180<sup>o</sup>, 360<sup>o</sup>] for longitude and, in | * latitude and [-180<sup>o</sup>, 360<sup>o</sup>] for longitude and, in | |||
* which case \e lat and \e lon are unchanged. Finally the longitude i s | * which case \e lat and \e lon are unchanged. Finally the longitude i s | |||
* reduced to the range [-180<sup>o</sup>, 180<sup>o</sup>). | * reduced to the range [-180<sup>o</sup>, 180<sup>o</sup>). | |||
********************************************************************** / | ********************************************************************** / | |||
static void DecodeLatLon(const std::string& dmsa, const std::string& dm sb, | static void DecodeLatLon(const std::string& dmsa, const std::string& dm sb, | |||
real& lat, real& lon); | real& lat, real& lon); | |||
/** | /** | |||
* Convert \e degree into a DMS string. \e trailing indicates the leas | * Convert a string \e angstr to an angle in degrees. No hemisphere | |||
t | * designator is allowed and no check is done on the range of the resul | |||
* significant component of the string (and this component is given as | t. | |||
a | ********************************************************************** | |||
* decimal number if necessary). \e prec indicates the number of digit | / | |||
s | static Math::real DecodeAngle(const std::string& angstr); | |||
* after the decimal point for the trailing component. \e flag indicat | ||||
es | /** | |||
* additional formating as follows | * Convert a string \e azistr to an azimuth in degrees. A hemisphere | |||
* - flag == NONE, signed result no leading zeros on degrees except in | * designator E/W can be used; the result is multiplied by -1 if W is | |||
the | * present. Throws an error if the result is out of the range | |||
* units place, e.g., -8d03'. | * [-180<sup>o</sup>, 360<sup>o</sup>]. Finally the azimuth is reduced | |||
* - flag == LATITUDE, trailing N or S hemisphere designator, no sign, | to | |||
pad | * the range [-180<sup>o</sup>, 180<sup>o</sup>). | |||
* degress to 2 digits, e.g., 08d03'S. | ********************************************************************** | |||
* - flag == LONGITUDE, trailing E or W hemisphere designator, no sign, | / | |||
pad | static Math::real DecodeAzimuth(const std::string& azistr); | |||
* degress to 3 digits, e.g., 008d03'W. | ||||
* - flag == AZIMUTH, convert to the range [0, 360<sup>o</sup>), no sig | /** | |||
n, | * Convert \e angle (in degrees) into a DMS string. \e trailing indica | |||
* pad degrees to 3 digits, , e.g., 351d57'. | tes | |||
* the least significant component of the string (and this component is | ||||
* given as a decimal number if necessary). \e prec indicates the numb | ||||
er | ||||
* of digits after the decimal point for the trailing component. \e fl | ||||
ag | ||||
* indicates additional formating as follows | ||||
* - flag == DMS::NONE, signed result no leading zeros on degrees excep | ||||
t in | ||||
* the units place, e.g., -8d03'. | ||||
* - flag == DMS::LATITUDE, trailing N or S hemisphere designator, no s | ||||
ign, | ||||
* pad degress to 2 digits, e.g., 08d03'S. | ||||
* - flag == DMS::LONGITUDE, trailing E or W hemisphere designator, no | ||||
* sign, pad degress to 3 digits, e.g., 008d03'W. | ||||
* - flag == DMS::AZIMUTH, convert to the range [0, 360<sup>o</sup>), n | ||||
o | ||||
* sign, pad degrees to 3 digits, , e.g., 351d57'. | ||||
* . | * . | |||
* The integer parts of the minutes and seconds components are always g iven | * The integer parts of the minutes and seconds components are always g iven | |||
* with 2 digits. | * with 2 digits. | |||
********************************************************************** / | ********************************************************************** / | |||
static std::string Encode(real degree, component trailing, unsigned pre c, | static std::string Encode(real angle, component trailing, unsigned prec , | |||
flag ind = NONE); | flag ind = NONE); | |||
/** | /** | |||
* Convert \e degree into a DMS string selecting the trailing component | * Convert \e angle into a DMS string selecting the trailing component | |||
* based on \e prec. \e prec indicates the precision relative to 1 deg ree, | * based on \e prec. \e prec indicates the precision relative to 1 deg ree, | |||
* e.g., \e prec = 3 gives a result accurate to 0.1' and \e prec = 4 gi ves | * e.g., \e prec = 3 gives a result accurate to 0.1' and \e prec = 4 gi ves | |||
* a result accurate to 1". | * a result accurate to 1". If \e ind is DMS::NUMBER, then merely | |||
* format \e angle as a number in fixed format with precision \e prec. | ||||
********************************************************************** / | ********************************************************************** / | |||
static std::string Encode(real degree, unsigned prec, flag ind = NONE) | static std::string Encode(real angle, unsigned prec, flag ind = NONE) { | |||
{ | if (ind == NUMBER) { | |||
return Encode(degree, | std::ostringstream s; | |||
prec < 2 ? DEGREE : (prec < 4 ? MINUTE : SECOND), | s << std::fixed << std::setprecision(prec) << angle; | |||
prec < 2 ? prec : (prec < 4 ? prec - 2 : prec - 4), | return s.str(); | |||
ind); | } else | |||
return Encode(angle, | ||||
prec < 2 ? DEGREE : (prec < 4 ? MINUTE : SECOND), | ||||
prec < 2 ? prec : (prec < 4 ? prec - 2 : prec - 4), | ||||
ind); | ||||
} | } | |||
/** | /** | |||
* Split angle, \e ang, into degrees, \e d, and minutes \e m. | * Split angle, \e ang, into degrees, \e d, and minutes \e m. | |||
********************************************************************** / | ********************************************************************** / | |||
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); | |||
} | } | |||
/** | /** | |||
End of changes. 11 change blocks. | ||||
57 lines changed or deleted | 88 lines changed or added | |||
Geocentric.hpp | Geocentric.hpp | |||
---|---|---|---|---|
/** | /** | |||
* \file Geocentric.hpp | * \file Geocentric.hpp | |||
* \brief Header for GeographicLib::Geocentric class | * \brief Header for GeographicLib::Geocentric class | |||
* | * | |||
* Copyright (c) Charles Karney (2008, 2009, 2010) <charles@karney.com> | * Copyright (c) Charles Karney (2008, 2009, 2010) <charles@karney.com> | |||
* and licensed under the LGPL. For more information, see | * and licensed under the LGPL. For more information, see | |||
* http://geographiclib.sourceforge.net/ | * http://geographiclib.sourceforge.net/ | |||
**********************************************************************/ | **********************************************************************/ | |||
#if !defined(GEOGRAPHICLIB_GEOCENTRIC_HPP) | #if !defined(GEOGRAPHICLIB_GEOCENTRIC_HPP) | |||
#define GEOGRAPHICLIB_GEOCENTRIC_HPP "$Id: Geocentric.hpp 6785 2010-01-05 2 2:15:42Z karney $" | #define GEOGRAPHICLIB_GEOCENTRIC_HPP "$Id: Geocentric.hpp 6827 2010-05-20 1 9:56:18Z karney $" | |||
#include "GeographicLib/Constants.hpp" | #include "GeographicLib/Constants.hpp" | |||
namespace GeographicLib { | namespace GeographicLib { | |||
/** | /** | |||
* \brief %Geocentric coordinates | * \brief %Geocentric coordinates | |||
* | * | |||
* Convert between geodetic coordinates latitude = \e lat, longitude = \e | * Convert between geodetic coordinates latitude = \e lat, longitude = \e | |||
* lon, height = \e h (measured vertically from the surface of the ellips oid) | * lon, height = \e h (measured vertically from the surface of the ellips oid) | |||
skipping to change at line 48 | skipping to change at line 48 | |||
* | * | |||
* The errors in these routines are close to round-off. Specifically, fo r | * The errors in these routines are close to round-off. Specifically, fo r | |||
* points within 5000 km of the surface of the ellipsoid (either inside o r | * points within 5000 km of the surface of the ellipsoid (either inside o r | |||
* outside the ellipsoid), the error is bounded by 7 nm for the WGS84 | * outside the ellipsoid), the error is bounded by 7 nm for the WGS84 | |||
* ellipsoid. See \ref geocentric for further information on the errors. | * ellipsoid. See \ref geocentric for further information on the errors. | |||
**********************************************************************/ | **********************************************************************/ | |||
class Geocentric { | class Geocentric { | |||
private: | private: | |||
typedef Math::real real; | typedef Math::real real; | |||
const real _a, _r, _f, _e2, _e2m, _ax, _e2x, _e4x, _e2mx, _maxrad; | const real _a, _r, _f, _e2, _e2m, _e2a, _e4a, _maxrad; | |||
static inline real sq(real x) throw() { return x * x; } | static inline real sq(real x) throw() { return x * x; } | |||
public: | public: | |||
/** | /** | |||
* Constructor for a ellipsoid radius \e a (meters) and reciprocal | * Constructor for a ellipsoid radius \e a (meters) and reciprocal | |||
* flattening \e r. Setting \e r = 0 implies \e r = inf or flattening = 0 | * flattening \e r. Setting \e r = 0 implies \e r = inf or flattening = 0 | |||
* (i.e., a sphere). Negative \e r indicates a prolate spheroid. An | * (i.e., a sphere). Negative \e r indicates a prolate spheroid. An | |||
* exception is thrown if \e a is not positive. | * exception is thrown if \e a is not positive. | |||
********************************************************************** / | ********************************************************************** / | |||
Geocentric(real a, real r); | Geocentric(real a, real r); | |||
End of changes. 2 change blocks. | ||||
2 lines changed or deleted | 2 lines changed or added | |||
Geodesic.hpp | Geodesic.hpp | |||
---|---|---|---|---|
/** | /** | |||
* \file Geodesic.hpp | * \file Geodesic.hpp | |||
* \brief Header for GeographicLib::Geodesic class | * \brief Header for GeographicLib::Geodesic class | |||
* | * | |||
* Copyright (c) Charles Karney (2009, 2010) <charles@karney.com> | * Copyright (c) Charles Karney (2009, 2010) <charles@karney.com> | |||
* and licensed under the LGPL. For more information, see | * and licensed under the LGPL. For more information, see | |||
* http://geographiclib.sourceforge.net/ | * http://geographiclib.sourceforge.net/ | |||
**********************************************************************/ | **********************************************************************/ | |||
#if !defined(GEOGRAPHICLIB_GEODESIC_HPP) | #if !defined(GEOGRAPHICLIB_GEODESIC_HPP) | |||
#define GEOGRAPHICLIB_GEODESIC_HPP "$Id: Geodesic.hpp 6785 2010-01-05 22:15 :42Z karney $" | #define GEOGRAPHICLIB_GEODESIC_HPP "$Id: Geodesic.hpp 6827 2010-05-20 19:56 :18Z karney $" | |||
#if !defined(GEOD_TAU_ORD) | #include "GeographicLib/Constants.hpp" | |||
/** | ||||
* The order of the series relating \e s and \e sigma when expanded in powe | ||||
rs | ||||
* of \e k1. | ||||
**********************************************************************/ | ||||
#define GEOD_TAU_ORD \ | ||||
(sizeof(real) == sizeof(double) ? 5 : sizeof(real) == sizeof(float) ? 4 : 6 | ||||
) | ||||
#endif | ||||
#if !defined(GEOD_ETA_ORD) | #if !defined(GEOD_ORD) | |||
/** | /** | |||
* The order of the series relating \e lambda and \e eta when expanded in | * The order of the expansions used by Geodesic. | |||
* powers of \e fp. | ||||
**********************************************************************/ | **********************************************************************/ | |||
#define GEOD_ETA_ORD \ | #define GEOD_ORD (GEOGRAPHICLIB_PREC == 1 ? 6 : GEOGRAPHICLIB_PREC == 0 ? 3 | |||
(sizeof(real) == sizeof(double) ? 5 : sizeof(real) == sizeof(float) ? 4 : 6 | : 7) | |||
) | ||||
#endif | #endif | |||
#include "GeographicLib/Constants.hpp" | ||||
namespace GeographicLib { | namespace GeographicLib { | |||
class GeodesicLine; | class GeodesicLine; | |||
/** | /** | |||
* \brief %Geodesic calculations | * \brief %Geodesic calculations | |||
* | * | |||
* The shortest path between two points on a ellipsoid at (\e lat1, \e lo n1) | * The shortest path between two points on a ellipsoid at (\e lat1, \e lo n1) | |||
* and (\e lat2, \e lon2) is called the geodesic. Its length is \e s12 a nd | * and (\e lat2, \e lon2) is called the geodesic. Its length is \e s12 a nd | |||
* the geodesic from point 1 to point 2 has azimuths \e azi1 and \e azi2 at | * the geodesic from point 1 to point 2 has azimuths \e azi1 and \e azi2 at | |||
skipping to change at line 86 | skipping to change at line 75 | |||
* The calculations are accurate to better than 15 nm. (See \ref geoderr ors | * The calculations are accurate to better than 15 nm. (See \ref geoderr ors | |||
* for details.) | * for details.) | |||
* | * | |||
* For more information on geodesics see \ref geodesic. | * For more information on geodesics see \ref geodesic. | |||
**********************************************************************/ | **********************************************************************/ | |||
class Geodesic { | class Geodesic { | |||
private: | private: | |||
typedef Math::real real; | typedef Math::real real; | |||
friend class GeodesicLine; | friend class GeodesicLine; | |||
friend class CassiniSoldner; | static const int nA1 = GEOD_ORD, nC1 = GEOD_ORD, nC1p = GEOD_ORD, | |||
static const int tauord = GEOD_TAU_ORD; | nA2 = GEOD_ORD, nC2 = GEOD_ORD, nA3 = GEOD_ORD, nC3 = GEOD_ORD; | |||
static const int ntau = tauord; | ||||
static const int nsig = tauord; | ||||
static const int zetord = GEOD_TAU_ORD; | ||||
static const int nzet = zetord; | ||||
static const int etaord = GEOD_ETA_ORD; | ||||
// etaCoeff is multiplied by etaFactor which is O(f), so we reduce the | ||||
// order to which etaCoeff is computed by 1. | ||||
static const int neta = etaord > 0 ? etaord - 1 : 0; | ||||
static const unsigned maxit = 50; | static const unsigned maxit = 50; | |||
static inline real sq(real x) throw() { return x * x; } | static inline real sq(real x) throw() { return x * x; } | |||
void Lengths(real k1, real sig12, | void Lengths(real eps, real sig12, | |||
real ssig1, real csig1, real ssig2, real csig2, | real ssig1, real csig1, real ssig2, real csig2, | |||
real cbet1, real cbet2, | real cbet1, real cbet2, | |||
real& s12s, real& m12a, real& m0, | real& s12s, real& m12a, real& m0, | |||
real tc[], real zc[]) const throw(); | real tc[], real zc[]) const throw(); | |||
static void Evolute(real R, real z, real& c, real& s) throw(); | static real Astroid(real R, real z) throw(); | |||
void InverseStart(real sbet1, real cbet1, real sbet2, real cbet2, | real InverseStart(real sbet1, real cbet1, real sbet2, real cbet2, | |||
real lam12, real slam12, real clam12, | real lam12, | |||
real& salp1, real& calp1, | real& salp1, real& calp1, | |||
real tc[], real zc[]) const throw(); | real& salp2, real& calp2, | |||
real C1a[], real C2a[]) const throw(); | ||||
real Lambda12(real sbet1, real cbet1, real sbet2, real cbet2, | real Lambda12(real sbet1, real cbet1, real sbet2, real cbet2, | |||
real salp1, real calp1, | real salp1, real calp1, | |||
real& salp2, real& calp2, real& sig12, | real& salp2, real& calp2, real& sig12, | |||
real& ssig1, real& csig1, real& ssig2, real& csig2, | real& ssig1, real& csig1, real& ssig2, real& csig2, | |||
real& k1, bool diffp, real& dlam12, | real& eps, bool diffp, real& dlam12, | |||
real tc[], real zc[], real ec[]) | real C1a[], real C2a[], real C3a[]) | |||
const throw(); | const throw(); | |||
static const real eps2, tol0, tol1, tol2, xthresh; | static const real eps2, tol0, tol1, tol2, xthresh; | |||
const real _a, _r, _f, _f1, _e2, _ep2, _n, _b; | const real _a, _r, _f, _f1, _e2, _ep2, _n, _b, _etol2; | |||
static real SinSeries(real sinx, real cosx, const real c[], int n) | static real SinSeries(real sinx, real cosx, const real c[], int n) | |||
throw(); | throw(); | |||
static inline real AngNormalize(real x) throw() { | static inline real AngNormalize(real x) throw() { | |||
// Place angle in [-180, 180). Assumes x is in [-540, 540). | // Place angle in [-180, 180). Assumes x is in [-540, 540). | |||
return x >= 180 ? x - 360 : x < -180 ? x + 360 : x; | return x >= 180 ? x - 360 : x < -180 ? x + 360 : x; | |||
} | } | |||
static inline real AngRound(real x) throw() { | static inline real AngRound(real x) throw() { | |||
// The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^ 57 | // The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^ 57 | |||
// for reals = 0.7 pm on the earth if x is an angle in degrees. (Thi s | // for reals = 0.7 pm on the earth if x is an angle in degrees. (Thi s | |||
// is about 1000 times more resolution than we get with angles around 90 | // is about 1000 times more resolution than we get with angles around 90 | |||
skipping to change at line 142 | skipping to change at line 124 | |||
volatile real y = std::abs(x); | volatile real y = std::abs(x); | |||
// The compiler mustn't "simplify" z - (z - y) to y | // The compiler mustn't "simplify" z - (z - y) to y | |||
y = y < z ? z - (z - y) : y; | y = y < z ? z - (z - y) : y; | |||
return x < 0 ? -y : y; | return x < 0 ? -y : y; | |||
} | } | |||
static inline void SinCosNorm(real& sinx, real& cosx) throw() { | static inline void SinCosNorm(real& sinx, real& cosx) throw() { | |||
real r = Math::hypot(sinx, cosx); | real r = Math::hypot(sinx, cosx); | |||
sinx /= r; | sinx /= r; | |||
cosx /= r; | cosx /= r; | |||
} | } | |||
// These are Maxima generated functions to provide series approximation s to | // These are Maxima generated functions to provide series approximation s to | |||
// the integrals for the ellipsoidal geodesic. | // the integrals for the ellipsoidal geodesic. | |||
static real tauFactorm1(real k1) throw(); | static real A1m1f(real eps) throw(); | |||
static void tauCoeff(real k1, real t[]) throw(); | static void C1f(real eps, real c[]) throw(); | |||
static void sigCoeff(real k1, real tp[]) throw(); | static void C1pf(real eps, real c[]) throw(); | |||
static real zetFactorm1(real k1) throw(); | static real A2m1f(real eps) throw(); | |||
static void zetCoeff(real k1, real t[]) throw(); | static void C2f(real eps, real c[]) throw(); | |||
static real etaFactor(real f, real k1) throw(); | static real A3f(real f, real eps) throw(); | |||
static void etaCoeff(real f, real k1, real h[]) throw(); | static void C3f(real f, real eps, real c[]) throw(); | |||
public: | public: | |||
/** | /** | |||
* Constructor for a ellipsoid with equatorial radius \e a (meters) and | * Constructor for a ellipsoid with equatorial radius \e a (meters) and | |||
* reciprocal flattening \e r. Setting \e r = 0 implies \e r = inf or | * reciprocal flattening \e r. Setting \e r = 0 implies \e r = inf or | |||
* flattening = 0 (i.e., a sphere). Negative \e r indicates a prolate | * flattening = 0 (i.e., a sphere). Negative \e r indicates a prolate | |||
* ellipsoid. An exception is thrown if \e a is not positive. | * ellipsoid. An exception is thrown if \e a is not positive. | |||
********************************************************************** / | ********************************************************************** / | |||
Geodesic(real a, real r); | Geodesic(real a, real r); | |||
skipping to change at line 226 | skipping to change at line 207 | |||
* used in the constructor. A value of 0 is returned for a sphere | * used in the constructor. A value of 0 is returned for a sphere | |||
* (infinite inverse flattening). | * (infinite inverse flattening). | |||
********************************************************************** / | ********************************************************************** / | |||
Math::real InverseFlattening() const throw() { return _r; } | Math::real InverseFlattening() const throw() { return _r; } | |||
/** | /** | |||
* A global instantiation of Geodesic with the parameters for the WGS84 | * A global instantiation of Geodesic with the parameters for the WGS84 | |||
* ellipsoid. | * ellipsoid. | |||
********************************************************************** / | ********************************************************************** / | |||
const static Geodesic WGS84; | const static Geodesic WGS84; | |||
}; | }; | |||
/** | /** | |||
* \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. Geodesic.Line returns a GeodesicLine object with the | * single geodesic. Geodesic.Line returns a GeodesicLine object with the | |||
* geodesic defined by by \e lat1, \e lon1, and \e azi1. | * geodesic defined by by \e lat1, \e lon1, and \e azi1. | |||
* GeodesicLine.Position returns the \e lat2, \e lon2, \e azi2, and \e m1 2 | * GeodesicLine.Position returns the \e lat2, \e lon2, \e azi2, and \e m1 2 | |||
* given \e s12. An example of use of this class is: | * given \e s12. An example of use of this class is: | |||
skipping to change at line 268 | skipping to change at line 250 | |||
* Similarly, a vector can be used to hold GeodesicLine objects. | * Similarly, a vector can be used to hold GeodesicLine objects. | |||
* | * | |||
* The calculations are accurate to better than 12 nm. (See \ref geoderr ors | * The calculations are accurate to better than 12 nm. (See \ref geoderr ors | |||
* for details.) | * for details.) | |||
**********************************************************************/ | **********************************************************************/ | |||
class GeodesicLine { | class GeodesicLine { | |||
private: | private: | |||
typedef Math::real real; | typedef Math::real real; | |||
friend class Geodesic; | friend class Geodesic; | |||
friend class CassiniSoldner; | static const int nC1 = Geodesic::nC1, nC1p = Geodesic::nC1p, | |||
static const int ntau = Geodesic::ntau; | nC2 = Geodesic::nC2, nC3 = Geodesic::nC3; | |||
static const int nsig = Geodesic::nsig; | ||||
static const int nzet = Geodesic::nzet; | ||||
static const int neta = Geodesic::neta; | ||||
real _lat1, _lon1, _azi1; | real _lat1, _lon1, _azi1; | |||
real _a, _r, _b, _f1, _salp0, _calp0, _u2, | real _a, _r, _b, _f1, _salp0, _calp0, _k2, | |||
_ssig1, _csig1, _stau1, _ctau1, _somg1, _comg1, | _ssig1, _csig1, _stau1, _ctau1, _somg1, _comg1, | |||
_taufm1, _zetfm1, _etaf, _dtau1, _dzet1, _dlam1; | _A1m1, _A2m1, _A3c, _B11, _B21, _B31; | |||
real _tauCoeff[ntau ? ntau : 1], _sigCoeff[nsig ? nsig : 1], | // index zero elements of these arrays are unused | |||
_zetCoeff[nzet ? nzet : 1], _etaCoeff[neta ? neta : 1]; | real _C1a[nC1 + 1], _C1pa[nC1p + 1], _C2a[nC2 + 1], _C3a[nC3]; | |||
static inline real sq(real x) throw() { return x * x; } | static inline real sq(real x) throw() { return x * x; } | |||
GeodesicLine(const Geodesic& g, real lat1, real lon1, real azi1) | GeodesicLine(const Geodesic& g, real lat1, real lon1, real azi1) | |||
throw(); | throw(); | |||
public: | public: | |||
/** | /** | |||
* A default constructor. If GeodesicLine::Position is called on the | * A default constructor. If GeodesicLine::Position is called on the | |||
* resulting object, it returns immediately (without doing any | * resulting object, it returns immediately (without doing any | |||
* calculations). The object should be set with a call to Geodesic::Li ne. | * calculations). The object should be set with a call to Geodesic::Li ne. | |||
skipping to change at line 308 | skipping to change at line 287 | |||
* \e s12 can be signed. If \e arcmode (default false) is set to true, \e | * \e s12 can be signed. If \e arcmode (default false) is set to true, \e | |||
* s12 is interpreted as the arc length \e a12 (in degrees) on the | * s12 is interpreted as the arc length \e a12 (in degrees) on the | |||
* auxiliary sphere. Returned value is the arc length \e a12 (degrees) if | * auxiliary sphere. Returned value is the arc length \e a12 (degrees) if | |||
* \e arcmode is false, otherwise it is the distance \e s12 (meters). | * \e arcmode is false, otherwise it is the distance \e s12 (meters). | |||
********************************************************************** / | ********************************************************************** / | |||
Math::real Position(real s12, real& lat2, real& lon2, | Math::real Position(real s12, real& lat2, real& lon2, | |||
real& azi2, real &m12, bool arcmode = false) | real& azi2, real &m12, bool arcmode = false) | |||
const throw(); | const throw(); | |||
/** | /** | |||
* Return the scale of the geodesic line extending an arc length \e a12 | ||||
* (degrees) from point 1 to point 2. \e M12 (a number) measures the | ||||
* convergence of initially parallel geodesics. It is defined by the | ||||
* following construction: starting at point 1 proceed at azimuth \e az | ||||
i1 + | ||||
* 90<sup>o</sup> a small distance \e dt; turn -90<sup>o</sup> and proc | ||||
eed | ||||
* a distance \e s12 (\e not the arc length \e a12); the distance to po | ||||
int | ||||
* 2 is given by \e M12 \e dt. \e M21 is defined analogously. | ||||
********************************************************************** | ||||
/ | ||||
void Scale(real a12, real& M12, real& M21) const throw(); | ||||
/** | ||||
* Has this object been initialized so that Position can be called? | * Has this object been initialized so that Position can be called? | |||
********************************************************************** / | ********************************************************************** / | |||
bool Init() const throw() { return _b > 0; } | bool Init() const throw() { return _b > 0; } | |||
/** | /** | |||
* Return the latitude of point 1 (degrees). | * Return the latitude of point 1 (degrees). | |||
********************************************************************** / | ********************************************************************** / | |||
Math::real Latitude() const throw() { return Init() ? _lat1 : 0; } | Math::real Latitude() const throw() { return Init() ? _lat1 : 0; } | |||
/** | /** | |||
skipping to change at line 329 | skipping to change at line 319 | |||
********************************************************************** / | ********************************************************************** / | |||
Math::real Longitude() const throw() { return Init() ? _lon1 : 0; } | Math::real Longitude() const throw() { return Init() ? _lon1 : 0; } | |||
/** | /** | |||
* Return the azimuth (degrees) of the geodesic line as it passes throu gh | * Return the azimuth (degrees) of the geodesic line as it passes throu gh | |||
* point 1. | * point 1. | |||
********************************************************************** / | ********************************************************************** / | |||
Math::real Azimuth() const throw() { return Init() ? _azi1 : 0; } | Math::real Azimuth() const throw() { return Init() ? _azi1 : 0; } | |||
/** | /** | |||
* Return the azimuth (degrees) of the geodesic line as it crosses the | ||||
* equator in a northward direction. | ||||
********************************************************************** | ||||
/ | ||||
Math::real EquatorialAzimuth() const throw() { | ||||
return Init() ? atan2(_salp0, _calp0) / Constants::degree() : 0; | ||||
} | ||||
/** | ||||
* Return the arc length (degrees) between the northward equatorial | ||||
* crossing and point 1. | ||||
********************************************************************** | ||||
/ | ||||
Math::real EquatorialArc() const throw() { | ||||
return Init() ? atan2(_ssig1, _csig1) / Constants::degree() : 0; | ||||
} | ||||
/** | ||||
* The major radius of the ellipsoid (meters). This is that value of \ e a | * The major radius of the ellipsoid (meters). This is that value of \ e a | |||
* inherited from the Geodesic object used in the constructor. | * inherited from the Geodesic object used in the constructor. | |||
********************************************************************** / | ********************************************************************** / | |||
Math::real MajorRadius() const throw() { return Init() ? _a : 0; } | Math::real MajorRadius() const throw() { return Init() ? _a : 0; } | |||
/** | /** | |||
* The inverse flattening of the ellipsoid. This is that value of \e r | * The inverse flattening of the ellipsoid. This is that value of \e r | |||
* inherited from the Geodesic object used in the constructor. A value of | * inherited from the Geodesic object used in the constructor. A value of | |||
* 0 is returned for a sphere (infinite inverse flattening). | * 0 is returned for a sphere (infinite inverse flattening). | |||
********************************************************************** / | ********************************************************************** / | |||
End of changes. 20 change blocks. | ||||
54 lines changed or deleted | 64 lines changed or added | |||
TransverseMercator.hpp | TransverseMercator.hpp | |||
---|---|---|---|---|
/** | /** | |||
* \file TransverseMercator.hpp | * \file TransverseMercator.hpp | |||
* \brief Header for GeographicLib::TransverseMercator class | * \brief Header for GeographicLib::TransverseMercator class | |||
* | * | |||
* Copyright (c) Charles Karney (2008, 2009, 2010) <charles@karney.com> | * Copyright (c) Charles Karney (2008, 2009, 2010) <charles@karney.com> | |||
* and licensed under the LGPL. For more information, see | * and licensed under the LGPL. For more information, see | |||
* http://geographiclib.sourceforge.net/ | * http://geographiclib.sourceforge.net/ | |||
**********************************************************************/ | **********************************************************************/ | |||
#if !defined(GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP) | #if !defined(GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP) | |||
#define GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP "$Id: TransverseMercator.hpp 6 807 2010-02-01 11:26:34Z karney $" | #define GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP "$Id: TransverseMercator.hpp 6 824 2010-04-19 14:25:10Z karney $" | |||
#include "GeographicLib/Constants.hpp" | #include "GeographicLib/Constants.hpp" | |||
#if !defined(TM_TX_MAXPOW) | #if !defined(TM_TX_MAXPOW) | |||
/** | /** | |||
* The order of the series approximation used in | * The order of the series approximation used in | |||
* GeographicLib::TransverseMercator. TM_TX_MAXPOW can be set to any integ er | * GeographicLib::TransverseMercator. TM_TX_MAXPOW can be set to any integ er | |||
* in [4, 8]. | * in [4, 8]. | |||
**********************************************************************/ | **********************************************************************/ | |||
#define TM_TX_MAXPOW \ | #define TM_TX_MAXPOW \ | |||
(sizeof(real) == sizeof(double) ? 6 : sizeof(real) == sizeof(float) ? 4 : 8 ) | (GEOGRAPHICLIB_PREC == 1 ? 6 : GEOGRAPHICLIB_PREC == 0 ? 4 : 8) | |||
#endif | #endif | |||
namespace GeographicLib { | namespace GeographicLib { | |||
/** | /** | |||
* \brief Transverse Mercator Projection | * \brief Transverse Mercator Projection | |||
* | * | |||
* This uses Krüger's method which evaluates the projection and its | * This uses Krüger's method which evaluates the projection and its | |||
* inverse in terms of a series. See | * inverse in terms of a series. See | |||
* - L. Krüger, | * - L. Krüger, | |||
* <a href="http://dx.doi.org/10.2312/GFZ.b103-krueger28"> Konforme | * <a href="http://dx.doi.org/10.2312/GFZ.b103-krueger28"> Konforme | |||
* Abbildung des Erdellipsoids in der Ebene</a> (Conformal mapping of the | * Abbildung des Erdellipsoids in der Ebene</a> (Conformal mapping of the | |||
* ellipsoidal earth to the plane), Royal Prussian Geodetic Institute, New | * ellipsoidal earth to the plane), Royal Prussian Geodetic Institute, New | |||
* Series 52, 172 pp. (1912). | * Series 52, 172 pp. (1912). | |||
* | * | |||
* Krüger's method has been extended from 4th to 6th order. The max imum | * Krüger's method has been extended from 4th to 6th order. The max imum | |||
* errors is 5 nm (ground distance) for all positions within 35 degrees o f | * errors is 5 nm (ground distance) for all positions within 35 degrees o f | |||
* the central meridian. The error in the convergence is 2e-15" and the | * the central meridian. The error in the convergence is 2e-15" and the | |||
* relative error in the scale is 6e-12%%. (See \ref tmerrors for the we asel | * relative error in the scale is 6e-12%%. (See \ref tmerrors for the we asel | |||
* words.) The speed penalty in going to 6th order is only about 1%. | * words.) The speed penalty in going to 6th order is only about 1%. | |||
* GeographicLib::TransverseMercatorExact is an alternative implementatio n of | * GeographicLib::TransverseMercatorExact is an alternative implementatio n of | |||
* the projection using exact formulas which yield accurate (to 8 nm) | * the projection using exact formulas which yield accurate (to 8 nm) res | |||
* results over the entire ellipsoid. | ults | |||
* over the entire ellipsoid. | ||||
* | * | |||
* The ellipsoid parameters and the central scale are set in the construc tor. | * The ellipsoid parameters and the central scale are set in the construc tor. | |||
* The central meridian (which is a trivial shift of the longitude) is | * The central meridian (which is a trivial shift of the longitude) is | |||
* specified as the \e lon0 argument of the TransverseMercator::Forward a nd | * specified as the \e lon0 argument of the TransverseMercator::Forward a nd | |||
* TransverseMercator::Reverse functions. The latitude of origin is take n to | * TransverseMercator::Reverse functions. The latitude of origin is take n to | |||
* be the equator. There is no provision in this class for specifying a | * be the equator. There is no provision in this class for specifying a | |||
* false easting or false northing or a different latitude of origin. | * false easting or false northing or a different latitude of origin. | |||
* However these are can be simply included by the calling funtcion. For | * However these are can be simply included by the calling funtcion. For | |||
* example, the UTMUPS class applies the false easting and false northing for | * example, the UTMUPS class applies the false easting and false northing for | |||
* the UTM projections. A more complicated example is the British Nation al | * the UTM projections. A more complicated example is the British Nation al | |||
End of changes. 4 change blocks. | ||||
5 lines changed or deleted | 6 lines changed or added | |||
TransverseMercatorExact.hpp | TransverseMercatorExact.hpp | |||
---|---|---|---|---|
/** | /** | |||
* \file TransverseMercatorExact.hpp | * \file TransverseMercatorExact.hpp | |||
* \brief Header for GeographicLib::TransverseMercatorExact class | * \brief Header for GeographicLib::TransverseMercatorExact class | |||
* | * | |||
* Copyright (c) Charles Karney (2008, 2009, 2010) <charles@karney.com> | * Copyright (c) Charles Karney (2008, 2009, 2010) <charles@karney.com> | |||
* and licensed under the LGPL. For more information, see | * and licensed under the LGPL. For more information, see | |||
* http://geographiclib.sourceforge.net/ | * http://geographiclib.sourceforge.net/ | |||
**********************************************************************/ | **********************************************************************/ | |||
#if !defined(GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_HPP) | #if !defined(GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_HPP) | |||
#define GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_HPP "$Id: TransverseMercatorE xact.hpp 6807 2010-02-01 11:26:34Z karney $" | #define GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_HPP "$Id: TransverseMercatorE xact.hpp 6824 2010-04-19 14:25:10Z karney $" | |||
#include "GeographicLib/Constants.hpp" | #include "GeographicLib/Constants.hpp" | |||
#include "GeographicLib/EllipticFunction.hpp" | #include "GeographicLib/EllipticFunction.hpp" | |||
namespace GeographicLib { | namespace GeographicLib { | |||
/** | /** | |||
* \brief An exact implementation of the Transverse Mercator Projection | * \brief An exact implementation of the Transverse Mercator Projection | |||
* | * | |||
* Implementation of the Transverse Mercator Projection given in | * Implementation of the Transverse Mercator Projection given in | |||
skipping to change at line 33 | skipping to change at line 33 | |||
* <a href="http://dx.doi.org/10.3138/X687-1574-4325-WM62"> Conformal | * <a href="http://dx.doi.org/10.3138/X687-1574-4325-WM62"> Conformal | |||
* Projections Based On Jacobian Elliptic Functions</a>, Part V of | * Projections Based On Jacobian Elliptic Functions</a>, Part V of | |||
* Conformal Projections Based on Elliptic Functions, | * Conformal Projections Based on Elliptic Functions, | |||
* (B. V. Gutsell, Toronto, 1976), 128pp., | * (B. V. Gutsell, Toronto, 1976), 128pp., | |||
* ISBN: 0919870163 | * ISBN: 0919870163 | |||
* (also appeared as: | * (also appeared as: | |||
* Monograph 16, Suppl. No. 1 to Canadian Cartographer, Vol 13). | * Monograph 16, Suppl. No. 1 to Canadian Cartographer, Vol 13). | |||
* | * | |||
* This method gives the correct results for forward and reverse | * This method gives the correct results for forward and reverse | |||
* transformations subject to the branch cut rules (see the description o f | * transformations subject to the branch cut rules (see the description o f | |||
* the \e extendp argument to the constructor).. The maximum error is ab | * the \e extendp argument to the constructor). The maximum error is abo | |||
out | ut 8 | |||
* 8 nm (ground distance) for the forward and reverse transformations. T | * nm (ground distance) for the forward and reverse transformations. The | |||
he | * error in the convergence is 2e-15", the relative error in the sca | |||
* error in the convergence is 2e-15", the relative error in the scale is | le | |||
* 7e-12%%. (See \ref tmerrors for the weasel words.) The method is "ex | * is 7e-12%%. (See \ref tmerrors for the weasel words.) The method is | |||
act" | * "exact" in the sense that the errors are close to the round-off limit | |||
* in the sense that the errors are close to the round-off limit and that | and | |||
no | * that no changes are needed in the algorithms for them to be used with | |||
* changes are needed in the algorithms for them to be used with reals of | * reals of a higher precision. Thus the errors using long double (with | |||
a | a | |||
* higher precision. Thus the errors using long double (with a 64-bit | * 64-bit fraction) are about 2000 times smaller than using double (with | |||
* fraction) are about 2000 times smaller than using double (with a 53-bi | a | |||
t | * 53-bit fraction). | |||
* fraction). | ||||
* | * | |||
* This algorithm is about 4.5 times slower than the 6th-order Krüge r | * This algorithm is about 4.5 times slower than the 6th-order Krüge r | |||
* method, GeographicLib::TransverseMercator, taking about 11 us for a | * method, GeographicLib::TransverseMercator, taking about 11 us for a | |||
* combined forward and reverse projection on a 2.6 GHz Intel machine (g+ +, | * combined forward and reverse projection on a 2.6 GHz Intel machine (g+ +, | |||
* version 4.3.0, -O3). | * version 4.3.0, -O3). | |||
* | * | |||
* The ellipsoid parameters and the central scale are set in the construc tor. | * The ellipsoid parameters and the central scale are set in the construc tor. | |||
* The central meridian (which is a trivial shift of the longitude) is | * The central meridian (which is a trivial shift of the longitude) is | |||
* specified as the \e lon0 argument of the TransverseMercatorExact::Forw ard | * specified as the \e lon0 argument of the TransverseMercatorExact::Forw ard | |||
* and TransverseMercatorExact::Reverse functions. The latitude of origi n is | * and TransverseMercatorExact::Reverse functions. The latitude of origi n is | |||
End of changes. 2 change blocks. | ||||
16 lines changed or deleted | 15 lines changed or added | |||