// Boost.Geometry (aka GGL, Generic Geometry Library) // Copyright (c) 2011-2012 Barend Gehrels, Amsterdam, the Netherlands. // Use, modification and distribution is subject to the Boost Software License, // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) #ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_SSF_HPP #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_SSF_HPP #include #include #include #include #include #include #include #include //#include namespace boost { namespace geometry { namespace strategy { namespace side { /*! \brief Check at which side of a Great Circle segment a point lies left of segment (> 0), right of segment (< 0), on segment (0) \ingroup strategies \tparam CalculationType \tparam_calculation */ template class spherical_side_formula { public : template static inline int apply(P1 const& p1, P2 const& p2, P const& p) { typedef typename boost::mpl::if_c < boost::is_void::type::value, // Select at least a double... typename select_most_precise < typename select_most_precise < typename select_most_precise < typename coordinate_type::type, typename coordinate_type::type >::type, typename coordinate_type

::type >::type, double >::type, CalculationType >::type coordinate_type; // Convenient shortcuts typedef coordinate_type ct; ct const lambda1 = get_as_radian<0>(p1); ct const delta1 = get_as_radian<1>(p1); ct const lambda2 = get_as_radian<0>(p2); ct const delta2 = get_as_radian<1>(p2); ct const lambda = get_as_radian<0>(p); ct const delta = get_as_radian<1>(p); // Create temporary points (vectors) on unit a sphere ct const cos_delta1 = cos(delta1); ct const c1x = cos_delta1 * cos(lambda1); ct const c1y = cos_delta1 * sin(lambda1); ct const c1z = sin(delta1); ct const cos_delta2 = cos(delta2); ct const c2x = cos_delta2 * cos(lambda2); ct const c2y = cos_delta2 * sin(lambda2); ct const c2z = sin(delta2); // (Third point is converted directly) ct const cos_delta = cos(delta); // Apply the "Spherical Side Formula" as presented on my blog ct const dist = (c1y * c2z - c1z * c2y) * cos_delta * cos(lambda) + (c1z * c2x - c1x * c2z) * cos_delta * sin(lambda) + (c1x * c2y - c1y * c2x) * sin(delta); ct zero = ct(); return dist > zero ? 1 : dist < zero ? -1 : 0; } }; #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS namespace services { /*template struct default_strategy { typedef spherical_side_formula type; };*/ template struct default_strategy { typedef spherical_side_formula type; }; template struct default_strategy { typedef spherical_side_formula type; }; } #endif }} // namespace strategy::side }} // namespace boost::geometry #endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_SSF_HPP