Chris@16: // Boost.Geometry (aka GGL, Generic Geometry Library)
Chris@16: 
Chris@16: // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
Chris@16: 
Chris@16: // Use, modification and distribution is subject to the Boost Software License,
Chris@16: // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
Chris@16: // http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_HAVERSINE_HPP
Chris@16: #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_HAVERSINE_HPP
Chris@16: 
Chris@16: 
Chris@16: #include <boost/geometry/core/cs.hpp>
Chris@16: #include <boost/geometry/core/access.hpp>
Chris@16: #include <boost/geometry/core/radian_access.hpp>
Chris@16: 
Chris@101: #include <boost/geometry/util/math.hpp>
Chris@16: #include <boost/geometry/util/select_calculation_type.hpp>
Chris@16: #include <boost/geometry/util/promote_floating_point.hpp>
Chris@16: 
Chris@16: #include <boost/geometry/strategies/distance.hpp>
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: namespace boost { namespace geometry
Chris@16: {
Chris@16: 
Chris@16: 
Chris@16: namespace strategy { namespace distance
Chris@16: {
Chris@16: 
Chris@16: 
Chris@16: namespace comparable
Chris@16: {
Chris@16: 
Chris@16: // Comparable haversine.
Chris@16: // To compare distances, we can avoid:
Chris@16: // - multiplication with radius and 2.0
Chris@16: // - applying sqrt
Chris@16: // - applying asin (which is strictly (monotone) increasing)
Chris@16: template
Chris@16: <
Chris@16:     typename RadiusType,
Chris@16:     typename CalculationType = void
Chris@16: >
Chris@16: class haversine
Chris@16: {
Chris@16: public :
Chris@16:     template <typename Point1, typename Point2>
Chris@16:     struct calculation_type
Chris@16:         : promote_floating_point
Chris@16:           <
Chris@16:               typename select_calculation_type
Chris@16:                   <
Chris@16:                       Point1,
Chris@16:                       Point2,
Chris@16:                       CalculationType
Chris@16:                   >::type
Chris@16:           >
Chris@16:     {};
Chris@16: 
Chris@16:     typedef RadiusType radius_type;
Chris@16: 
Chris@16:     explicit inline haversine(RadiusType const& r = 1.0)
Chris@16:         : m_radius(r)
Chris@16:     {}
Chris@16: 
Chris@16:     template <typename Point1, typename Point2>
Chris@16:     static inline typename calculation_type<Point1, Point2>::type
Chris@16:     apply(Point1 const& p1, Point2 const& p2)
Chris@16:     {
Chris@16:         return calculate<typename calculation_type<Point1, Point2>::type>(
Chris@16:                    get_as_radian<0>(p1), get_as_radian<1>(p1),
Chris@16:                    get_as_radian<0>(p2), get_as_radian<1>(p2)
Chris@16:                );
Chris@16:     }
Chris@16: 
Chris@16:     inline RadiusType radius() const
Chris@16:     {
Chris@16:         return m_radius;
Chris@16:     }
Chris@16: 
Chris@16: 
Chris@16: private :
Chris@16:     template <typename R, typename T1, typename T2>
Chris@16:     static inline R calculate(T1 const& lon1, T1 const& lat1,
Chris@16:                               T2 const& lon2, T2 const& lat2)
Chris@16:     {
Chris@16:         return math::hav(lat2 - lat1)
Chris@16:                 + cos(lat1) * cos(lat2) * math::hav(lon2 - lon1);
Chris@16:     }
Chris@16: 
Chris@16:     RadiusType m_radius;
Chris@16: };
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: } // namespace comparable
Chris@16: 
Chris@16: /*!
Chris@16: \brief Distance calculation for spherical coordinates
Chris@16: on a perfect sphere using haversine
Chris@16: \ingroup strategies
Chris@16: \tparam RadiusType \tparam_radius
Chris@16: \tparam CalculationType \tparam_calculation
Chris@16: \author Adapted from: http://williams.best.vwh.net/avform.htm
Chris@16: \see http://en.wikipedia.org/wiki/Great-circle_distance
Chris@16: \note (from Wiki:) The great circle distance d between two
Chris@16: points with coordinates {lat1,lon1} and {lat2,lon2} is given by:
Chris@16:     d=acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lon1-lon2))
Chris@16: A mathematically equivalent formula, which is less subject
Chris@16:     to rounding error for short distances is:
Chris@16:     d=2*asin(sqrt((sin((lat1-lat2) / 2))^2
Chris@16:     + cos(lat1)*cos(lat2)*(sin((lon1-lon2) / 2))^2))
Chris@101: 
Chris@16: 
Chris@16: \qbk{
Chris@16: [heading See also]
Chris@16: [link geometry.reference.algorithms.distance.distance_3_with_strategy distance (with strategy)]
Chris@16: }
Chris@16: 
Chris@16: */
Chris@16: template
Chris@16: <
Chris@16:     typename RadiusType,
Chris@16:     typename CalculationType = void
Chris@16: >
Chris@16: class haversine
Chris@16: {
Chris@16:     typedef comparable::haversine<RadiusType, CalculationType> comparable_type;
Chris@16: 
Chris@16: public :
Chris@16:     template <typename Point1, typename Point2>
Chris@16:     struct calculation_type
Chris@16:         : services::return_type<comparable_type, Point1, Point2>
Chris@16:     {};
Chris@16: 
Chris@16:     typedef RadiusType radius_type;
Chris@16: 
Chris@16:     /*!
Chris@16:     \brief Constructor
Chris@16:     \param radius radius of the sphere, defaults to 1.0 for the unit sphere
Chris@16:     */
Chris@16:     inline haversine(RadiusType const& radius = 1.0)
Chris@16:         : m_radius(radius)
Chris@16:     {}
Chris@16: 
Chris@16:     /*!
Chris@16:     \brief applies the distance calculation
Chris@16:     \return the calculated distance (including multiplying with radius)
Chris@16:     \param p1 first point
Chris@16:     \param p2 second point
Chris@16:     */
Chris@16:     template <typename Point1, typename Point2>
Chris@16:     inline typename calculation_type<Point1, Point2>::type
Chris@16:     apply(Point1 const& p1, Point2 const& p2) const
Chris@16:     {
Chris@16:         typedef typename calculation_type<Point1, Point2>::type calculation_type;
Chris@16:         calculation_type const a = comparable_type::apply(p1, p2);
Chris@101:         calculation_type const c = calculation_type(2.0) * asin(math::sqrt(a));
Chris@101:         return calculation_type(m_radius) * c;
Chris@16:     }
Chris@16: 
Chris@16:     /*!
Chris@16:     \brief access to radius value
Chris@16:     \return the radius
Chris@16:     */
Chris@16:     inline RadiusType radius() const
Chris@16:     {
Chris@16:         return m_radius;
Chris@16:     }
Chris@16: 
Chris@16: private :
Chris@16:     RadiusType m_radius;
Chris@16: };
Chris@16: 
Chris@16: 
Chris@16: #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
Chris@16: namespace services
Chris@16: {
Chris@16: 
Chris@16: template <typename RadiusType, typename CalculationType>
Chris@16: struct tag<haversine<RadiusType, CalculationType> >
Chris@16: {
Chris@16:     typedef strategy_tag_distance_point_point type;
Chris@16: };
Chris@16: 
Chris@16: 
Chris@16: template <typename RadiusType, typename CalculationType, typename P1, typename P2>
Chris@16: struct return_type<haversine<RadiusType, CalculationType>, P1, P2>
Chris@16:     : haversine<RadiusType, CalculationType>::template calculation_type<P1, P2>
Chris@16: {};
Chris@16: 
Chris@16: 
Chris@16: template <typename RadiusType, typename CalculationType>
Chris@16: struct comparable_type<haversine<RadiusType, CalculationType> >
Chris@16: {
Chris@16:     typedef comparable::haversine<RadiusType, CalculationType> type;
Chris@16: };
Chris@16: 
Chris@16: 
Chris@16: template <typename RadiusType, typename CalculationType>
Chris@16: struct get_comparable<haversine<RadiusType, CalculationType> >
Chris@16: {
Chris@16: private :
Chris@16:     typedef haversine<RadiusType, CalculationType> this_type;
Chris@16:     typedef comparable::haversine<RadiusType, CalculationType> comparable_type;
Chris@16: public :
Chris@16:     static inline comparable_type apply(this_type const& input)
Chris@16:     {
Chris@16:         return comparable_type(input.radius());
Chris@16:     }
Chris@16: };
Chris@16: 
Chris@16: template <typename RadiusType, typename CalculationType, typename P1, typename P2>
Chris@16: struct result_from_distance<haversine<RadiusType, CalculationType>, P1, P2>
Chris@16: {
Chris@16: private :
Chris@16:     typedef haversine<RadiusType, CalculationType> this_type;
Chris@16:     typedef typename return_type<this_type, P1, P2>::type return_type;
Chris@16: public :
Chris@16:     template <typename T>
Chris@16:     static inline return_type apply(this_type const& , T const& value)
Chris@16:     {
Chris@16:         return return_type(value);
Chris@16:     }
Chris@16: };
Chris@16: 
Chris@16: 
Chris@16: // Specializations for comparable::haversine
Chris@16: template <typename RadiusType, typename CalculationType>
Chris@16: struct tag<comparable::haversine<RadiusType, CalculationType> >
Chris@16: {
Chris@16:     typedef strategy_tag_distance_point_point type;
Chris@16: };
Chris@16: 
Chris@16: 
Chris@16: template <typename RadiusType, typename CalculationType, typename P1, typename P2>
Chris@16: struct return_type<comparable::haversine<RadiusType, CalculationType>, P1, P2>
Chris@16:     : comparable::haversine<RadiusType, CalculationType>::template calculation_type<P1, P2>
Chris@16: {};
Chris@16: 
Chris@16: 
Chris@16: template <typename RadiusType, typename CalculationType>
Chris@16: struct comparable_type<comparable::haversine<RadiusType, CalculationType> >
Chris@16: {
Chris@16:     typedef comparable::haversine<RadiusType, CalculationType> type;
Chris@16: };
Chris@16: 
Chris@16: 
Chris@16: template <typename RadiusType, typename CalculationType>
Chris@16: struct get_comparable<comparable::haversine<RadiusType, CalculationType> >
Chris@16: {
Chris@16: private :
Chris@16:     typedef comparable::haversine<RadiusType, CalculationType> this_type;
Chris@16: public :
Chris@16:     static inline this_type apply(this_type const& input)
Chris@16:     {
Chris@16:         return input;
Chris@16:     }
Chris@16: };
Chris@16: 
Chris@16: 
Chris@16: template <typename RadiusType, typename CalculationType, typename P1, typename P2>
Chris@16: struct result_from_distance<comparable::haversine<RadiusType, CalculationType>, P1, P2>
Chris@16: {
Chris@16: private :
Chris@16:     typedef comparable::haversine<RadiusType, CalculationType> strategy_type;
Chris@16:     typedef typename return_type<strategy_type, P1, P2>::type return_type;
Chris@16: public :
Chris@16:     template <typename T>
Chris@16:     static inline return_type apply(strategy_type const& strategy, T const& distance)
Chris@16:     {
Chris@16:         return_type const s = sin((distance / strategy.radius()) / return_type(2));
Chris@16:         return s * s;
Chris@16:     }
Chris@16: };
Chris@16: 
Chris@16: 
Chris@101: // Register it as the default for point-types
Chris@16: // in a spherical equatorial coordinate system
Chris@16: template <typename Point1, typename Point2>
Chris@101: struct default_strategy
Chris@101:     <
Chris@101:         point_tag, point_tag, Point1, Point2,
Chris@101:         spherical_equatorial_tag, spherical_equatorial_tag
Chris@101:     >
Chris@16: {
Chris@16:     typedef strategy::distance::haversine<typename select_coordinate_type<Point1, Point2>::type> type;
Chris@16: };
Chris@16: 
Chris@16: // Note: spherical polar coordinate system requires "get_as_radian_equatorial"
Chris@16: 
Chris@16: 
Chris@16: } // namespace services
Chris@16: #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
Chris@16: 
Chris@16: 
Chris@16: }} // namespace strategy::distance
Chris@16: 
Chris@16: 
Chris@16: }} // namespace boost::geometry
Chris@16: 
Chris@16: 
Chris@16: #endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_HAVERSINE_HPP