Subject: [Boost-commit] svn:boost r63601 - in sandbox/geometry/boost/geometry: extensions/algorithms extensions/gis/geographic/strategies extensions/gis/latlong/detail extensions/gis/projections/impl strategies
From: barend.gehrels_at_[hidden]
Date: 2010-07-04 12:03:37
Author: barendgehrels
Date: 2010-07-04 12:03:33 EDT (Sun, 04 Jul 2010)
New Revision: 63601
URL: http://svn.boost.org/trac/boost/changeset/63601
Log:
Removed more std:: occurances
Added high precision for Vincenty
Text files modified:
sandbox/geometry/boost/geometry/extensions/algorithms/remove_spikes.hpp | 4
sandbox/geometry/boost/geometry/extensions/gis/geographic/strategies/vincenty.hpp | 243 ++++++++++++++++++++++++---------------
sandbox/geometry/boost/geometry/extensions/gis/latlong/detail/graticule.hpp | 4
sandbox/geometry/boost/geometry/extensions/gis/projections/impl/aasincos.hpp | 15 +
sandbox/geometry/boost/geometry/extensions/gis/projections/impl/adjlon.hpp | 3
sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_auth.hpp | 4
sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_ell_set.hpp | 4
sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_fwd.hpp | 10
sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_gauss.hpp | 15 +
sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_init.hpp | 5
sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_inv.hpp | 3
sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_mlfn.hpp | 6
sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_phi2.hpp | 5
sandbox/geometry/boost/geometry/extensions/gis/projections/impl/proj_mdist.hpp | 6
sandbox/geometry/boost/geometry/strategies/strategy_transform.hpp | 20 +-
15 files changed, 216 insertions(+), 131 deletions(-)
Modified: sandbox/geometry/boost/geometry/extensions/algorithms/remove_spikes.hpp
==============================================================================
--- sandbox/geometry/boost/geometry/extensions/algorithms/remove_spikes.hpp (original)
+++ sandbox/geometry/boost/geometry/extensions/algorithms/remove_spikes.hpp 2010-07-04 12:03:33 EDT (Sun, 04 Jul 2010)
@@ -343,8 +343,8 @@
// Normalize the vectors -> this results in points+direction
// and is comparible between geometries
- coordinate_type const magnitude1 = std::sqrt(dx1 * dx1 + dy1 * dy1);
- coordinate_type const magnitude2 = std::sqrt(dx2 * dx2 + dy2 * dy2);
+ coordinate_type const magnitude1 = sqrt(dx1 * dx1 + dy1 * dy1);
+ coordinate_type const magnitude2 = sqrt(dx2 * dx2 + dy2 * dy2);
if (magnitude1 > m_zero && magnitude2 > m_zero)
{
Modified: sandbox/geometry/boost/geometry/extensions/gis/geographic/strategies/vincenty.hpp
==============================================================================
--- sandbox/geometry/boost/geometry/extensions/gis/geographic/strategies/vincenty.hpp (original)
+++ sandbox/geometry/boost/geometry/extensions/gis/geographic/strategies/vincenty.hpp 2010-07-04 12:03:33 EDT (Sun, 04 Jul 2010)
@@ -15,10 +15,15 @@
#include <boost/geometry/strategies/distance.hpp>
#include <boost/geometry/core/radian_access.hpp>
#include <boost/geometry/core/coordinate_type.hpp>
+#include <boost/geometry/util/select_calculation_type.hpp>
+#include <boost/geometry/util/promote_floating_point.hpp>
+#include <boost/geometry/util/math.hpp>
#include <boost/geometry/extensions/gis/geographic/detail/ellipsoid.hpp>
+
+
namespace boost { namespace geometry
{
@@ -28,8 +33,8 @@
/*!
\brief Distance calculation formulae on latlong coordinates, after Vincenty, 1975
\ingroup distance
- \tparam P1 first point type
- \tparam P2 optional second point type
+ \tparam Point1 first point type
+ \tparam Point2 optional second point type
\author See http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
\author Adapted from various implementations to get it close to the original document
- http://www.movable-type.co.uk/scripts/LatLongVincenty.html
@@ -37,108 +42,158 @@
- http://futureboy.homeip.net/fsp/colorize.fsp?fileName=navigation.frink
*/
-template <typename P1, typename P2 = P1>
+template
+<
+ typename Point1,
+ typename Point2 = Point1,
+ typename CalculationType = void
+>
class vincenty
{
- public :
- //typedef spherical_distance return_type;
- typedef P1 first_point_type;
- typedef P2 second_point_type;
- typedef double return_type;
-
- inline vincenty()
- {}
-
- explicit inline vincenty(geometry::detail::ellipsoid const& e)
- : m_ellipsoid(e)
- {}
+public :
+ typedef typename promote_floating_point
+ <
+ typename select_most_precise
+ <
+ typename select_calculation_type
+ <
+ Point1,
+ Point2,
+ CalculationType
+ >::type,
+ double // it should be at least double otherwise Vincenty does not run
+ >::type
+ >::type return_type;
+
+ typedef Point1 first_point_type;
+ typedef Point2 second_point_type;
+
+ inline vincenty()
+ {}
+
+ explicit inline vincenty(geometry::detail::ellipsoid const& e)
+ : m_ellipsoid(e)
+ {}
- inline return_type apply(P1 const& p1, P2 const& p2) const
- {
- return calculate(get_as_radian<0>(p1), get_as_radian<1>(p1),
- get_as_radian<0>(p2), get_as_radian<1>(p2));
- }
+ inline return_type apply(Point1 const& p1, Point2 const& p2) const
+ {
+ return calculate(get_as_radian<0>(p1), get_as_radian<1>(p1),
+ get_as_radian<0>(p2), get_as_radian<1>(p2));
+ }
- inline geometry::detail::ellipsoid ellipsoid() const
- {
- return m_ellipsoid;
- }
+ inline geometry::detail::ellipsoid ellipsoid() const
+ {
+ return m_ellipsoid;
+ }
- private :
- typedef typename coordinate_type<P1>::type T1;
- typedef typename coordinate_type<P2>::type T2;
- geometry::detail::ellipsoid m_ellipsoid;
+private :
+ typedef return_type promoted_type;
+ geometry::detail::ellipsoid m_ellipsoid;
+
+ inline return_type calculate(promoted_type const& lon1,
+ promoted_type const& lat1,
+ promoted_type const& lon2,
+ promoted_type const& lat2) const
+ {
+ namespace mc = boost::math::constants;
+
+ promoted_type const c2 = 2;
+ promoted_type const two_pi = c2 * mc::pi<promoted_type>();
+
+ // lambda: difference in longitude on an auxiliary sphere
+ promoted_type L = lon2 - lon1;
+ promoted_type lambda = L;
- inline return_type calculate(T1 const& lon1, T1 const& lat1, T2 const& lon2, T2 const& lat2) const
+ if (L < -mc::pi<promoted_type>()) L += two_pi;
+ if (L > mc::pi<promoted_type>()) L -= two_pi;
+
+ if (math::equals(lat1, lat2) && math::equals(lon1, lon2))
{
- namespace mc = boost::math::constants;
+ return return_type(0);
+ }
- double const two_pi = 2.0 * mc::pi<double>();
+ // TODO: give ellipsoid a template-parameter
+ promoted_type const ellipsoid_f = m_ellipsoid.f();
+ promoted_type const ellipsoid_a = m_ellipsoid.a();
+ promoted_type const ellipsoid_b = m_ellipsoid.b();
+
+ // U: reduced latitude, defined by tan U = (1-f) tan phi
+ promoted_type const c1 = 1;
+ promoted_type const one_min_f = c1 - ellipsoid_f;
+
+ promoted_type const U1 = atan(one_min_f * tan(lat1)); // above (1)
+ promoted_type const U2 = atan(one_min_f * tan(lat2)); // above (1)
+
+ promoted_type const cos_U1 = cos(U1);
+ promoted_type const cos_U2 = cos(U2);
+ promoted_type const sin_U1 = sin(U1);
+ promoted_type const sin_U2 = sin(U2);
+
+ // alpha: azimuth of the geodesic at the equator
+ promoted_type cos2_alpha;
+ promoted_type sin_alpha;
+
+ // sigma: angular distance p1,p2 on the sphere
+ // sigma1: angular distance on the sphere from the equator to p1
+ // sigma_m: angular distance on the sphere from the equator to the midpoint of the line
+ promoted_type sigma;
+ promoted_type sin_sigma;
+ promoted_type cos2_sigma_m;
+
+ promoted_type previous_lambda;
+
+ promoted_type const c3 = 3;
+ promoted_type const c4 = 4;
+ promoted_type const c6 = 6;
+ promoted_type const c16 = 16;
- // lambda: difference in longitude on an auxiliary sphere
- double L = lon2 - lon1;
- double lambda = L;
-
- if (L < -mc::pi<double>()) L += two_pi;
- if (L > mc::pi<double>()) L -= two_pi;
-
- if (lat1 == lat2 && lon1 == lon2)
- {
- return return_type(0);
- }
-
- // U: reduced latitude, defined by tan U = (1-f) tan phi
- double U1 = atan((1-m_ellipsoid.f()) * tan(lat1)); // above (1)
- double U2 = atan((1-m_ellipsoid.f()) * tan(lat2)); // above (1)
-
- double cos_U1 = cos(U1);
- double cos_U2 = cos(U2);
- double sin_U1 = sin(U1);
- double sin_U2 = sin(U2);
-
- // alpha: azimuth of the geodesic at the equator
- double cos2_alpha;
- double sin_alpha;
-
- // sigma: angular distance p1,p2 on the sphere
- // sigma1: angular distance on the sphere from the equator to p1
- // sigma_m: angular distance on the sphere from the equator to the midpoint of the line
- double sigma;
- double sin_sigma;
- double cos2_sigma_m;
-
- double previous_lambda;
-
- do
- {
- previous_lambda = lambda; // (13)
- double sin_lambda = sin(lambda);
- double cos_lambda = cos(lambda);
- sin_sigma = sqrt(math::sqr(cos_U2 * sin_lambda) + math::sqr(cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_lambda)); // (14)
- double cos_sigma = sin_U1 * sin_U2 + cos_U1 * cos_U2 * cos_lambda; // (15)
- sin_alpha = cos_U1 * cos_U2 * sin_lambda / sin_sigma; // (17)
- cos2_alpha = 1.0 - math::sqr(sin_alpha);
- cos2_sigma_m = cos2_alpha == 0 ? 0 : cos_sigma - 2.0 * sin_U1 * sin_U2 / cos2_alpha; // (18)
-
- double C = m_ellipsoid.f()/16.0 * cos2_alpha * (4.0 + m_ellipsoid.f() * (4.0 - 3.0 * cos2_alpha)); // (10)
- sigma = atan2(sin_sigma, cos_sigma); // (16)
- lambda = L + (1.0 - C) * m_ellipsoid.f() * sin_alpha *
- (sigma + C * sin_sigma * ( cos2_sigma_m + C * cos_sigma * (-1.0 + 2.0 * math::sqr(cos2_sigma_m)))); // (11)
-
- } while (fabs(previous_lambda - lambda) > 1e-12 && fabs(lambda) < mc::pi<double>());
-
- double sqr_u = cos2_alpha * (math::sqr(m_ellipsoid.a()) - math::sqr(m_ellipsoid.b())) / math::sqr(m_ellipsoid.b()); // above (1)
-
- double A = 1.0 + sqr_u/16384.0 * (4096 + sqr_u * (-768.0 + sqr_u * (320.0 - 175.0 * sqr_u))); // (3)
- double B = sqr_u/1024.0 * (256.0 + sqr_u * ( -128.0 + sqr_u * (74.0 - 47.0 * sqr_u))); // (4)
- double delta_sigma = B * sin_sigma * ( cos2_sigma_m + (B/4.0) * (cos(sigma)* (-1.0 + 2.0 * cos2_sigma_m)
- - (B/6.0) * cos2_sigma_m * (-3.0 + 4.0 * math::sqr(sin_sigma)) * (-3.0 + 4.0 * cos2_sigma_m))); // (6)
+ promoted_type const c_e_12 = 1e-12;
- double dist = m_ellipsoid.b() * A * (sigma - delta_sigma); // (19)
+ do
+ {
+ previous_lambda = lambda; // (13)
+ promoted_type sin_lambda = sin(lambda);
+ promoted_type cos_lambda = cos(lambda);
+ sin_sigma = sqrt(math::sqr(cos_U2 * sin_lambda) + math::sqr(cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_lambda)); // (14)
+ promoted_type cos_sigma = sin_U1 * sin_U2 + cos_U1 * cos_U2 * cos_lambda; // (15)
+ sin_alpha = cos_U1 * cos_U2 * sin_lambda / sin_sigma; // (17)
+ cos2_alpha = c1 - math::sqr(sin_alpha);
+ cos2_sigma_m = cos2_alpha == 0 ? 0 : cos_sigma - c2 * sin_U1 * sin_U2 / cos2_alpha; // (18)
+
+
+ promoted_type C = ellipsoid_f/c16 * cos2_alpha * (c4 + ellipsoid_f * (c4 - c3 * cos2_alpha)); // (10)
+ sigma = atan2(sin_sigma, cos_sigma); // (16)
+ lambda = L + (c1 - C) * ellipsoid_f * sin_alpha *
+ (sigma + C * sin_sigma * ( cos2_sigma_m + C * cos_sigma * (-c1 + c2 * math::sqr(cos2_sigma_m)))); // (11)
+
+ } while (geometry::math::abs(previous_lambda - lambda) > c_e_12
+ && geometry::math::abs(lambda) < mc::pi<promoted_type>());
+
+ promoted_type sqr_u = cos2_alpha * (math::sqr(ellipsoid_a) - math::sqr(ellipsoid_b)) / math::sqr(ellipsoid_b); // above (1)
+
+ // Oops getting hard here
+ // (again, problem is that ttmath cannot divide by doubles, which is OK)
+ promoted_type const c47 = 47;
+ promoted_type const c74 = 74;
+ promoted_type const c128 = 128;
+ promoted_type const c256 = 256;
+ promoted_type const c175 = 175;
+ promoted_type const c320 = 320;
+ promoted_type const c768 = 768;
+ promoted_type const c1024 = 1024;
+ promoted_type const c4096 = 4096;
+ promoted_type const c16384 = 16384;
+
+ promoted_type A = c1 + sqr_u/c16384 * (c4096 + sqr_u * (-c768 + sqr_u * (c320 - c175 * sqr_u))); // (3)
+ promoted_type B = sqr_u/c1024 * (c256 + sqr_u * ( -c128 + sqr_u * (c74 - c47 * sqr_u))); // (4)
+ promoted_type delta_sigma = B * sin_sigma * ( cos2_sigma_m + (B/c4) * (cos(sigma)* (-c1 + c2 * cos2_sigma_m)
+ - (B/c6) * cos2_sigma_m * (-c3 + c4 * math::sqr(sin_sigma)) * (-c3 + c4 * cos2_sigma_m))); // (6)
- return return_type(dist);
- }
+ promoted_type dist = ellipsoid_b * A * (sigma - delta_sigma); // (19)
+
+ return return_type(dist);
+ }
};
#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
@@ -209,8 +264,8 @@
#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
-template <typename P1, typename P2>
-struct strategy_tag<strategy::distance::vincenty<P1, P2> >
+template <typename Point1, typename Point2>
+struct strategy_tag<strategy::distance::vincenty<Point1, Point2> >
{
typedef strategy_tag_distance_point_point type;
};
Modified: sandbox/geometry/boost/geometry/extensions/gis/latlong/detail/graticule.hpp
==============================================================================
--- sandbox/geometry/boost/geometry/extensions/gis/latlong/detail/graticule.hpp (original)
+++ sandbox/geometry/boost/geometry/extensions/gis/latlong/detail/graticule.hpp 2010-07-04 12:03:33 EDT (Sun, 04 Jul 2010)
@@ -14,6 +14,8 @@
#include <string>
#include <boost/numeric/conversion/cast.hpp>
+#include <boost/geometry/util/math.hpp>
+
namespace boost { namespace geometry
{
@@ -109,7 +111,7 @@
: (CardinalDir == south) ? 'S'
: ' ');
- value = std::fabs(value);
+ value = geometry::math::abs(value);
// Calculate the values
double fraction = 0;
Modified: sandbox/geometry/boost/geometry/extensions/gis/projections/impl/aasincos.hpp
==============================================================================
--- sandbox/geometry/boost/geometry/extensions/gis/projections/impl/aasincos.hpp (original)
+++ sandbox/geometry/boost/geometry/extensions/gis/projections/impl/aasincos.hpp 2010-07-04 12:03:33 EDT (Sun, 04 Jul 2010)
@@ -35,8 +35,12 @@
#ifndef BOOST_GEOMETRY_PROJECTIONS_IMPL_AASINCOS_HPP
#define BOOST_GEOMETRY_PROJECTIONS_IMPL_AASINCOS_HPP
+
#include <cmath>
+#include <boost/geometry/util/math.hpp>
+
+
namespace boost { namespace geometry { namespace projection
{
@@ -54,7 +58,7 @@
{
double av = 0;
- if ((av = std::fabs(v)) >= 1.0)
+ if ((av = geometry::math::abs(v)) >= 1.0)
{
if (av > aasincos::ONE_TOL)
{
@@ -63,14 +67,14 @@
return (v < 0.0 ? -HALFPI : HALFPI);
}
- return std::asin(v);
+ return asin(v);
}
inline double aacos(double v)
{
double av = 0;
- if ((av = std::fabs(v)) >= 1.0)
+ if ((av = geometry::math::abs(v)) >= 1.0)
{
if (av > aasincos::ONE_TOL)
{
@@ -84,12 +88,13 @@
inline double asqrt(double v)
{
- return ((v <= 0) ? 0. : std::sqrt(v));
+ return ((v <= 0) ? 0 : sqrt(v));
}
inline double aatan2(double n, double d)
{
- return ((std::fabs(n) < aasincos::ATOL && std::fabs(d) < aasincos::ATOL) ? 0.0 : std::atan2(n, d));
+ return ((geometry::math::abs(n) < aasincos::ATOL
+ && geometry::math::abs(d) < aasincos::ATOL) ? 0.0 : atan2(n, d));
}
Modified: sandbox/geometry/boost/geometry/extensions/gis/projections/impl/adjlon.hpp
==============================================================================
--- sandbox/geometry/boost/geometry/extensions/gis/projections/impl/adjlon.hpp (original)
+++ sandbox/geometry/boost/geometry/extensions/gis/projections/impl/adjlon.hpp 2010-07-04 12:03:33 EDT (Sun, 04 Jul 2010)
@@ -38,6 +38,7 @@
#include <cmath>
#include <boost/math/constants/constants.hpp>
+#include <boost/geometry/util/math.hpp>
#include <boost/geometry/extensions/gis/projections/impl/projects.hpp>
@@ -54,7 +55,7 @@
static const double TWOPI = 6.2831853071795864769;
static const double ONEPI = 3.14159265358979323846;
- if (std::fabs(lon) <= SPI)
+ if (geometry::math::abs(lon) <= SPI)
{
return lon;
}
Modified: sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_auth.hpp
==============================================================================
--- sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_auth.hpp (original)
+++ sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_auth.hpp 2010-07-04 12:03:33 EDT (Sun, 04 Jul 2010)
@@ -74,9 +74,9 @@
{
assert(0 != APA);
- const double t = beta + beta;
+ double const t = beta + beta;
- return(beta + APA[0] * std::sin(t) + APA[1] * std::sin(t + t) + APA[2] * std::sin(t + t + t));
+ return(beta + APA[0] * sin(t) + APA[1] * sin(t + t) + APA[2] * sin(t + t + t));
}
} // namespace detail
Modified: sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_ell_set.hpp
==============================================================================
--- sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_ell_set.hpp (original)
+++ sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_ell_set.hpp 2010-07-04 12:03:33 EDT (Sun, 04 Jul 2010)
@@ -38,6 +38,8 @@
#include <string>
#include <vector>
+#include <boost/geometry/util/math.hpp>
+
#include <boost/geometry/extensions/gis/projections/impl/pj_ellps.hpp>
#include <boost/geometry/extensions/gis/projections/impl/pj_param.hpp>
@@ -130,7 +132,7 @@
double tmp;
tmp = sin(pj_param(parameters, i ? "rR_lat_a" : "rR_lat_g").f);
- if (fabs(tmp) > HALFPI) {
+ if (geometry::math::abs(tmp) > HALFPI) {
throw proj_exception(-11);
}
tmp = 1. - es * tmp * tmp;
Modified: sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_fwd.hpp
==============================================================================
--- sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_fwd.hpp (original)
+++ sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_fwd.hpp 2010-07-04 12:03:33 EDT (Sun, 04 Jul 2010)
@@ -38,6 +38,8 @@
#include <cmath>
#include <boost/geometry/core/radian_access.hpp>
+#include <boost/geometry/util/math.hpp>
+
#include <boost/geometry/extensions/gis/projections/impl/adjlon.hpp>
/* general forward projection */
@@ -59,21 +61,21 @@
double lp_lon = geometry::get_as_radian<0>(ll);
double lp_lat = geometry::get_as_radian<1>(ll);
- const double t = std::fabs(lp_lat) - HALFPI;
+ const double t = geometry::math::abs(lp_lat) - HALFPI;
/* check for forward and latitude or longitude overange */
- if (t > forwrd::EPS || std::fabs(lp_lon) > 10.)
+ if (t > forwrd::EPS || geometry::math::abs(lp_lon) > 10.)
{
throw proj_exception();
}
- if (std::fabs(t) <= forwrd::EPS)
+ if (geometry::math::abs(t) <= forwrd::EPS)
{
lp_lat = lp_lat < 0. ? -HALFPI : HALFPI;
}
else if (par.geoc)
{
- lp_lat = std::atan(par.rone_es * std::tan(lp_lat));
+ lp_lat = atan(par.rone_es * tan(lp_lat));
}
lp_lon -= par.lam0; /* compute del lp.lam */
Modified: sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_gauss.hpp
==============================================================================
--- sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_gauss.hpp (original)
+++ sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_gauss.hpp 2010-07-04 12:03:33 EDT (Sun, 04 Jul 2010)
@@ -36,6 +36,9 @@
#define BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_GAUSS_HPP
+#include <boost/geometry/util/math.hpp>
+
+
namespace boost { namespace geometry { namespace projection {
namespace detail { namespace gauss {
@@ -55,7 +58,7 @@
inline double srat(double esinp, double exp)
{
- return (std::pow((1.0 - esinp) / (1.0 + esinp), exp));
+ return (pow((1.0 - esinp) / (1.0 + esinp), exp));
}
inline GAUSS gauss_ini(double e, double phi0, double &chi, double &rc)
@@ -90,8 +93,8 @@
template <typename T>
inline void gauss(GAUSS const& en, T& lam, T& phi)
{
- phi = 2.0 * std::atan(en.K * std::pow(std::tan(0.5 * phi + FORTPI), en.C)
- * srat(en.e * std::sin(phi), en.ratexp) ) - HALFPI;
+ phi = 2.0 * atan(en.K * pow(tan(0.5 * phi + FORTPI), en.C)
+ * srat(en.e * sin(phi), en.ratexp) ) - HALFPI;
lam *= en.C;
}
@@ -100,14 +103,14 @@
inline void inv_gauss(GAUSS const& en, T& lam, T& phi)
{
lam /= en.C;
- const double num = std::pow(std::tan(0.5 * phi + FORTPI) / en.K, 1.0 / en.C);
+ const double num = pow(tan(0.5 * phi + FORTPI) / en.K, 1.0 / en.C);
int i = 0;
for (i = MAX_ITER; i; --i)
{
- const double elp_phi = 2.0 * std::atan(num * srat(en.e * std::sin(phi), - 0.5 * en.e)) - HALFPI;
+ const double elp_phi = 2.0 * atan(num * srat(en.e * sin(phi), - 0.5 * en.e)) - HALFPI;
- if (std::fabs(elp_phi - phi) < DEL_TOL)
+ if (geometry::math::abs(elp_phi - phi) < DEL_TOL)
{
break;
}
Modified: sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_init.hpp
==============================================================================
--- sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_init.hpp (original)
+++ sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_init.hpp 2010-07-04 12:03:33 EDT (Sun, 04 Jul 2010)
@@ -41,6 +41,9 @@
#include <boost/algorithm/string.hpp>
#include <boost/range.hpp>
+#include <boost/geometry/util/math.hpp>
+
+
#include <boost/geometry/extensions/gis/projections/impl/pj_datum_set.hpp>
#include <boost/geometry/extensions/gis/projections/impl/pj_datums.hpp>
#include <boost/geometry/extensions/gis/projections/impl/pj_ell_set.hpp>
@@ -140,7 +143,7 @@
&& pin.datum_params[1] == 0.0
&& pin.datum_params[2] == 0.0
&& pin.a == 6378137.0
- && fabs(pin.es - 0.006694379990) < 0.000000000050 )/*WGS84/GRS80*/
+ && geometry::math::abs(pin.es - 0.006694379990) < 0.000000000050 )/*WGS84/GRS80*/
{
pin.datum_type = PJD_WGS84;
}
Modified: sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_inv.hpp
==============================================================================
--- sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_inv.hpp (original)
+++ sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_inv.hpp 2010-07-04 12:03:33 EDT (Sun, 04 Jul 2010)
@@ -38,6 +38,7 @@
#include <boost/geometry/extensions/gis/projections/impl/adjlon.hpp>
#include <boost/geometry/core/radian_access.hpp>
+#include <boost/geometry/util/math.hpp>
/* general inverse projection */
@@ -65,7 +66,7 @@
lon += par.lam0; /* reduce from del lp.lam */
if (!par.over)
lon = adjlon(lon); /* adjust longitude to CM */
- if (par.geoc && fabs(fabs(lat)-HALFPI) > inv::EPS)
+ if (par.geoc && geometry::math::abs(fabs(lat)-HALFPI) > inv::EPS)
lat = atan(par.one_es * tan(lat));
geometry::set_from_radian<0>(ll, lon);
Modified: sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_mlfn.hpp
==============================================================================
--- sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_mlfn.hpp (original)
+++ sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_mlfn.hpp 2010-07-04 12:03:33 EDT (Sun, 04 Jul 2010)
@@ -35,6 +35,10 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
+
+#include <boost/geometry/util/math.hpp>
+
+
namespace boost { namespace geometry { namespace projection {
namespace detail {
@@ -93,7 +97,7 @@
s = sin(phi);
t = 1. - es * s * s;
phi -= t = (pj_mlfn(phi, s, cos(phi), en) - arg) * (t * sqrt(t)) * k;
- if (fabs(t) < EPS)
+ if (geometry::math::abs(t) < EPS)
return phi;
}
throw proj_exception(-17);
Modified: sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_phi2.hpp
==============================================================================
--- sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_phi2.hpp (original)
+++ sandbox/geometry/boost/geometry/extensions/gis/projections/impl/pj_phi2.hpp 2010-07-04 12:03:33 EDT (Sun, 04 Jul 2010)
@@ -36,6 +36,9 @@
// DEALINGS IN THE SOFTWARE.
+#include <boost/geometry/util/math.hpp>
+
+
namespace boost { namespace geometry { namespace projection {
namespace detail {
@@ -58,7 +61,7 @@
dphi = HALFPI - 2. * atan (ts * pow((1. - con) /
(1. + con), eccnth)) - Phi;
Phi += dphi;
- } while ( fabs(dphi) > phi2::TOL && --i);
+ } while ( geometry::math::abs(dphi) > phi2::TOL && --i);
if (i <= 0)
throw proj_exception(-18);
return Phi;
Modified: sandbox/geometry/boost/geometry/extensions/gis/projections/impl/proj_mdist.hpp
==============================================================================
--- sandbox/geometry/boost/geometry/extensions/gis/projections/impl/proj_mdist.hpp (original)
+++ sandbox/geometry/boost/geometry/extensions/gis/projections/impl/proj_mdist.hpp 2010-07-04 12:03:33 EDT (Sun, 04 Jul 2010)
@@ -35,6 +35,10 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
+
+#include <boost/geometry/util/math.hpp>
+
+
namespace boost { namespace geometry { namespace projection
{
namespace detail
@@ -120,7 +124,7 @@
t = 1. - b.es * s * s;
phi -= t = (proj_mdist(phi, s, cos(phi), b) - dist) *
(t * sqrt(t)) * k;
- if (fabs(t) < TOL) /* that is no change */
+ if (geometry::math::abs(t) < TOL) /* that is no change */
return phi;
}
/* convergence failed */
Modified: sandbox/geometry/boost/geometry/strategies/strategy_transform.hpp
==============================================================================
--- sandbox/geometry/boost/geometry/strategies/strategy_transform.hpp (original)
+++ sandbox/geometry/boost/geometry/strategies/strategy_transform.hpp 2010-07-04 12:03:33 EDT (Sun, 04 Jul 2010)
@@ -154,10 +154,10 @@
// http://en.wikipedia.org/wiki/List_of_canonical_coordinate_transformations#From_spherical_coordinates
// Phi = first, theta is second, r is third, see documentation on cs::spherical
- double const sin_theta = std::sin(theta);
- set<0>(p, r * sin_theta * std::cos(phi));
- set<1>(p, r * sin_theta * std::sin(phi));
- set<2>(p, r * std::cos(theta));
+ double const sin_theta = sin(theta);
+ set<0>(p, r * sin_theta * cos(phi));
+ set<1>(p, r * sin_theta * sin(phi));
+ set<2>(p, r * cos(theta));
}
/// Helper function for conversion
@@ -170,7 +170,7 @@
#if defined(BOOST_GEOMETRY_TRANSFORM_CHECK_UNIT_SPHERE)
// TODO: MAYBE ONLY IF TO BE CHECKED?
- double const r = /*std::sqrt not necessary, sqrt(1)=1*/ (x * x + y * y + z * z);
+ double const r = /*sqrt not necessary, sqrt(1)=1*/ (x * x + y * y + z * z);
// Unit sphere, so r should be 1
if (geometry::math::abs(r - 1.0) > double(1e-6))
@@ -180,8 +180,8 @@
// end todo
#endif
- set_from_radian<0>(p, std::atan2(y, x));
- set_from_radian<1>(p, std::acos(z));
+ set_from_radian<0>(p, atan2(y, x));
+ set_from_radian<1>(p, acos(z));
return true;
}
@@ -191,12 +191,12 @@
assert_dimension<P, 3>();
// http://en.wikipedia.org/wiki/List_of_canonical_coordinate_transformations#From_Cartesian_coordinates
- double const r = std::sqrt(x * x + y * y + z * z);
+ double const r = sqrt(x * x + y * y + z * z);
set<2>(p, r);
- set_from_radian<0>(p, std::atan2(y, x));
+ set_from_radian<0>(p, atan2(y, x));
if (r > 0.0)
{
- set_from_radian<1>(p, std::acos(z / r));
+ set_from_radian<1>(p, acos(z / r));
return true;
}
return false;