[Boost-commit] svn:boost r72286 - in trunk/boost/geometry/strategies: cartesian spherical

Subject: [Boost-commit] svn:boost r72286 - in trunk/boost/geometry/strategies: cartesian spherical
From: barend.gehrels_at_[hidden]
Date: 2011-05-30 11:07:13

Author: barendgehrels
Date: 2011-05-30 11:07:12 EDT (Mon, 30 May 2011)
New Revision: 72286
URL: http://svn.boost.org/trac/boost/changeset/72286

Log:
Huiller: changed to calculation types (now supporting ttmath)
distance_projected_point.hpp: minor changes
Text files modified: 
   trunk/boost/geometry/strategies/cartesian/distance_projected_point.hpp |    25 ++++++++++-----------                   
   trunk/boost/geometry/strategies/spherical/area_huiller.hpp             |    47 +++++++++++++++++++++++---------------- 
   2 files changed, 40 insertions(+), 32 deletions(-)

Modified: trunk/boost/geometry/strategies/cartesian/distance_projected_point.hpp
==============================================================================
--- trunk/boost/geometry/strategies/cartesian/distance_projected_point.hpp	(original)
+++ trunk/boost/geometry/strategies/cartesian/distance_projected_point.hpp	2011-05-30 11:07:12 EDT (Mon, 30 May 2011)
@@ -116,13 +116,14 @@
     {
         assert_dimension_equal<Point, PointOfSegment>();
 
-        /* Algorithm
-        POINT v(x2 - x1, y2 - y1);
-        POINT w(px - x1, py - y1);
-        c1 = w . v
-        c2 = v . v
-        b = c1 / c2
-        RETURN POINT(x1 + b * vx, y1 + b * vy);
+        /* 
+            Algorithm [p1: (x1,y1), p2: (x2,y2), p: (px,py)]
+            VECTOR v(x2 - x1, y2 - y1)
+            VECTOR w(px - x1, py - y1)
+            c1 = w . v
+            c2 = v . v
+            b = c1 / c2
+            RETURN POINT(x1 + b * vx, y1 + b * vy)
         */
 
         // v is multiplied below with a (possibly) FP-value, so should be in FP
@@ -137,21 +138,20 @@
         Strategy strategy;
         boost::ignore_unused_variable_warning(strategy);
 
-        calculation_type zero = calculation_type();
-        fp_type c1 = dot_product(w, v);
+        calculation_type const zero = calculation_type();
+        fp_type const c1 = dot_product(w, v);
         if (c1 <= zero)
         {
             return strategy.apply(p, p1);
         }
-        fp_type c2 = dot_product(v, v);
+        fp_type const c2 = dot_product(v, v);
         if (c2 <= c1)
         {
             return strategy.apply(p, p2);
         }
 
         // See above, c1 > 0 AND c2 > c1 so: c2 != 0
-        fp_type b = fp_type(c1) / fp_type(c2);
-
+        fp_type const b = c1 / c2;
 
         fp_strategy_type fp_strategy
             = strategy::distance::services::get_similar
@@ -168,7 +168,6 @@
 
         return fp_strategy.apply(p, projected);
     }
-
 };
 
 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS

Modified: trunk/boost/geometry/strategies/spherical/area_huiller.hpp
==============================================================================
--- trunk/boost/geometry/strategies/spherical/area_huiller.hpp	(original)
+++ trunk/boost/geometry/strategies/spherical/area_huiller.hpp	2011-05-30 11:07:12 EDT (Mon, 30 May 2011)
@@ -10,7 +10,6 @@
 #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AREA_HUILLER_HPP
 
 
-#include <boost/math/constants/constants.hpp>
 
 #include <boost/geometry/strategies/spherical/distance_haversine.hpp>
 
@@ -66,10 +65,21 @@
 >
 class huiller
 {
+typedef typename boost::mpl::if_c
+    <
+        boost::is_void<CalculationType>::type::value,
+        typename select_most_precise
+            <
+                typename coordinate_type<PointOfSegment>::type,
+                double
+            >::type,
+        CalculationType
+    >::type calculation_type;
+
 protected :
     struct excess_sum
     {
-        double sum;
+        calculation_type sum;
 
         // Distances are calculated on unit sphere here
         strategy::distance::haversine<PointOfSegment, PointOfSegment>
@@ -80,18 +90,18 @@
             : sum(0)
             , distance_over_unit_sphere(1)
         {}
-        inline double area(double radius) const
+        inline calculation_type area(calculation_type radius) const
         {
             return - sum * radius * radius;
         }
     };
 
 public :
-    typedef double return_type;
+    typedef calculation_type return_type;
     typedef PointOfSegment segment_point_type;
     typedef excess_sum state_type;
 
-    inline huiller(double radius = 1.0)
+    inline huiller(calculation_type radius = 1.0)
         : m_radius(radius)
     {}
 
@@ -101,26 +111,25 @@
     {
         if (! geometry::math::equals(get<0>(p1), get<0>(p2)))
         {
-            namespace mc = boost::math::constants;
-
-            double const two_pi = 2.0 * mc::pi<double>();
-            double const half = 0.5;
-            double const two = 2.0;
-            double const four = 4.0;
+            calculation_type const half = 0.5;
+            calculation_type const two = 2.0;
+            calculation_type const four = 4.0;
+            calculation_type const two_pi = two * geometry::math::pi<calculation_type>();
+            calculation_type const half_pi = half * geometry::math::pi<calculation_type>();
 
             // Distance p1 p2
-            double a = state.distance_over_unit_sphere.apply(p1, p2);
+            calculation_type a = state.distance_over_unit_sphere.apply(p1, p2);
 
             // Sides on unit sphere to south pole
-            double b = half * mc::pi<double>() - geometry::get_as_radian<1>(p2);
-            double c = half * mc::pi<double>() - geometry::get_as_radian<1>(p1);
+            calculation_type b = half_pi - geometry::get_as_radian<1>(p2);
+            calculation_type c = half_pi - geometry::get_as_radian<1>(p1);
 
             // Semi parameter
-            double s = half * (a + b + c);
+            calculation_type s = half * (a + b + c);
 
             // E: spherical excess, using l'Huiller's formula
             // [tg(e / 4)]2   =   tg[s / 2]  tg[(s-a) / 2]  tg[(s-b) / 2]  tg[(s-c) / 2]
-            double E = four * atan(sqrt(geometry::math::abs(tan(s / two)
+            calculation_type E = four * atan(sqrt(geometry::math::abs(tan(s / two)
                     * tan((s - a) / two)
                     * tan((s - b) / two)
                     * tan((s - c) / two))));
@@ -132,11 +141,11 @@
             // we have to take the dateline into account.
             // TODO: check this / enhance this, should be more robust. See also the "grow" for ll
             // TODO: use minmax or "smaller"/"compare" strategy for this
-            double lon1 = geometry::get_as_radian<0>(p1) < 0
+            calculation_type lon1 = geometry::get_as_radian<0>(p1) < 0
                 ? geometry::get_as_radian<0>(p1) + two_pi
                 : geometry::get_as_radian<0>(p1);
 
-            double lon2 = geometry::get_as_radian<0>(p2) < 0
+            calculation_type lon2 = geometry::get_as_radian<0>(p2) < 0
                 ? geometry::get_as_radian<0>(p2) + two_pi
                 : geometry::get_as_radian<0>(p2);
 
@@ -156,7 +165,7 @@
 
 private :
     /// Radius of the sphere
-    double m_radius;
+    calculation_type m_radius;
 };
 
 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS