Chris@102: // Boost.Geometry (aka GGL, Generic Geometry Library)
Chris@102: 
Chris@102: // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
Chris@102: 
Chris@102: // This file was modified by Oracle on 2014.
Chris@102: // Modifications copyright (c) 2014 Oracle and/or its affiliates.
Chris@102: 
Chris@102: // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
Chris@102: 
Chris@102: // Use, modification and distribution is subject to the Boost Software License,
Chris@102: // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
Chris@102: // http://www.boost.org/LICENSE_1_0.txt)
Chris@102: 
Chris@102: #ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_ANDOYER_HPP
Chris@102: #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_ANDOYER_HPP
Chris@102: 
Chris@102: 
Chris@102: #include <boost/geometry/core/coordinate_type.hpp>
Chris@102: #include <boost/geometry/core/radian_access.hpp>
Chris@102: #include <boost/geometry/core/radius.hpp>
Chris@102: #include <boost/geometry/core/srs.hpp>
Chris@102: 
Chris@102: #include <boost/geometry/algorithms/detail/flattening.hpp>
Chris@102: 
Chris@102: #include <boost/geometry/strategies/distance.hpp>
Chris@102: 
Chris@102: #include <boost/geometry/util/math.hpp>
Chris@102: #include <boost/geometry/util/promote_floating_point.hpp>
Chris@102: #include <boost/geometry/util/select_calculation_type.hpp>
Chris@102: 
Chris@102: 
Chris@102: namespace boost { namespace geometry
Chris@102: {
Chris@102: 
Chris@102: namespace strategy { namespace distance
Chris@102: {
Chris@102: 
Chris@102: 
Chris@102: /*!
Chris@102: \brief Point-point distance approximation taking flattening into account
Chris@102: \ingroup distance
Chris@102: \tparam Spheroid The reference spheroid model
Chris@102: \tparam CalculationType \tparam_calculation
Chris@102: \author After Andoyer, 19xx, republished 1950, republished by Meeus, 1999
Chris@102: \note Although not so well-known, the approximation is very good: in all cases the results
Chris@102: are about the same as Vincenty. In my (Barend's) testcases the results didn't differ more than 6 m
Chris@102: \see http://nacc.upc.es/tierra/node16.html
Chris@102: \see http://sci.tech-archive.net/Archive/sci.geo.satellite-nav/2004-12/2724.html
Chris@102: \see http://home.att.net/~srschmitt/great_circle_route.html (implementation)
Chris@102: \see http://www.codeguru.com/Cpp/Cpp/algorithms/article.php/c5115 (implementation)
Chris@102: \see http://futureboy.homeip.net/frinksamp/navigation.frink (implementation)
Chris@102: \see http://www.voidware.com/earthdist.htm (implementation)
Chris@102: */
Chris@102: template
Chris@102: <
Chris@102:     typename Spheroid,
Chris@102:     typename CalculationType = void
Chris@102: >
Chris@102: class andoyer
Chris@102: {
Chris@102: public :
Chris@102:     template <typename Point1, typename Point2>
Chris@102:     struct calculation_type
Chris@102:         : promote_floating_point
Chris@102:           <
Chris@102:               typename select_calculation_type
Chris@102:                   <
Chris@102:                       Point1,
Chris@102:                       Point2,
Chris@102:                       CalculationType
Chris@102:                   >::type
Chris@102:           >
Chris@102:     {};
Chris@102: 
Chris@102:     typedef Spheroid model_type;
Chris@102: 
Chris@102:     inline andoyer()
Chris@102:         : m_spheroid()
Chris@102:     {}
Chris@102: 
Chris@102:     explicit inline andoyer(Spheroid const& spheroid)
Chris@102:         : m_spheroid(spheroid)
Chris@102:     {}
Chris@102: 
Chris@102: 
Chris@102:     template <typename Point1, typename Point2>
Chris@102:     inline typename calculation_type<Point1, Point2>::type
Chris@102:     apply(Point1 const& point1, Point2 const& point2) const
Chris@102:     {
Chris@102:         return calc<typename calculation_type<Point1, Point2>::type>
Chris@102:             (
Chris@102:                 get_as_radian<0>(point1), get_as_radian<1>(point1),
Chris@102:                 get_as_radian<0>(point2), get_as_radian<1>(point2)
Chris@102:             );
Chris@102:     }
Chris@102: 
Chris@102:     inline Spheroid const& model() const
Chris@102:     {
Chris@102:         return m_spheroid;
Chris@102:     }
Chris@102: 
Chris@102: private :
Chris@102:     template <typename CT, typename T>
Chris@102:     inline CT calc(T const& lon1,
Chris@102:                 T const& lat1,
Chris@102:                 T const& lon2,
Chris@102:                 T const& lat2) const
Chris@102:     {
Chris@102:         CT const G = (lat1 - lat2) / 2.0;
Chris@102:         CT const lambda = (lon1 - lon2) / 2.0;
Chris@102: 
Chris@102:         if (geometry::math::equals(lambda, 0.0)
Chris@102:             && geometry::math::equals(G, 0.0))
Chris@102:         {
Chris@102:             return 0.0;
Chris@102:         }
Chris@102: 
Chris@102:         CT const F = (lat1 + lat2) / 2.0;
Chris@102: 
Chris@102:         CT const sinG2 = math::sqr(sin(G));
Chris@102:         CT const cosG2 = math::sqr(cos(G));
Chris@102:         CT const sinF2 = math::sqr(sin(F));
Chris@102:         CT const cosF2 = math::sqr(cos(F));
Chris@102:         CT const sinL2 = math::sqr(sin(lambda));
Chris@102:         CT const cosL2 = math::sqr(cos(lambda));
Chris@102: 
Chris@102:         CT const S = sinG2 * cosL2 + cosF2 * sinL2;
Chris@102:         CT const C = cosG2 * cosL2 + sinF2 * sinL2;
Chris@102: 
Chris@102:         CT const c0 = 0;
Chris@102:         CT const c1 = 1;
Chris@102:         CT const c2 = 2;
Chris@102:         CT const c3 = 3;
Chris@102: 
Chris@102:         if (geometry::math::equals(S, c0) || geometry::math::equals(C, c0))
Chris@102:         {
Chris@102:             return c0;
Chris@102:         }
Chris@102: 
Chris@102:         CT const radius_a = CT(get_radius<0>(m_spheroid));
Chris@102:         CT const flattening = geometry::detail::flattening<CT>(m_spheroid);
Chris@102: 
Chris@102:         CT const omega = atan(math::sqrt(S / C));
Chris@102:         CT const r3 = c3 * math::sqrt(S * C) / omega; // not sure if this is r or greek nu
Chris@102:         CT const D = c2 * omega * radius_a;
Chris@102:         CT const H1 = (r3 - c1) / (c2 * C);
Chris@102:         CT const H2 = (r3 + c1) / (c2 * S);
Chris@102: 
Chris@102:         return D * (c1 + flattening * (H1 * sinF2 * cosG2 - H2 * cosF2 * sinG2) );
Chris@102:     }
Chris@102: 
Chris@102:     Spheroid m_spheroid;
Chris@102: };
Chris@102: 
Chris@102: 
Chris@102: #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
Chris@102: namespace services
Chris@102: {
Chris@102: 
Chris@102: template <typename Spheroid, typename CalculationType>
Chris@102: struct tag<andoyer<Spheroid, CalculationType> >
Chris@102: {
Chris@102:     typedef strategy_tag_distance_point_point type;
Chris@102: };
Chris@102: 
Chris@102: 
Chris@102: template <typename Spheroid, typename CalculationType, typename P1, typename P2>
Chris@102: struct return_type<andoyer<Spheroid, CalculationType>, P1, P2>
Chris@102:     : andoyer<Spheroid, CalculationType>::template calculation_type<P1, P2>
Chris@102: {};
Chris@102: 
Chris@102: 
Chris@102: template <typename Spheroid, typename CalculationType>
Chris@102: struct comparable_type<andoyer<Spheroid, CalculationType> >
Chris@102: {
Chris@102:     typedef andoyer<Spheroid, CalculationType> type;
Chris@102: };
Chris@102: 
Chris@102: 
Chris@102: template <typename Spheroid, typename CalculationType>
Chris@102: struct get_comparable<andoyer<Spheroid, CalculationType> >
Chris@102: {
Chris@102:     static inline andoyer<Spheroid, CalculationType> apply(andoyer<Spheroid, CalculationType> const& input)
Chris@102:     {
Chris@102:         return input;
Chris@102:     }
Chris@102: };
Chris@102: 
Chris@102: template <typename Spheroid, typename CalculationType, typename P1, typename P2>
Chris@102: struct result_from_distance<andoyer<Spheroid, CalculationType>, P1, P2>
Chris@102: {
Chris@102:     template <typename T>
Chris@102:     static inline typename return_type<andoyer<Spheroid, CalculationType>, P1, P2>::type
Chris@102:         apply(andoyer<Spheroid, CalculationType> const& , T const& value)
Chris@102:     {
Chris@102:         return value;
Chris@102:     }
Chris@102: };
Chris@102: 
Chris@102: 
Chris@102: template <typename Point1, typename Point2>
Chris@102: struct default_strategy<point_tag, point_tag, Point1, Point2, geographic_tag, geographic_tag>
Chris@102: {
Chris@102:     typedef strategy::distance::andoyer
Chris@102:                 <
Chris@102:                     srs::spheroid
Chris@102:                         <
Chris@102:                             typename select_coordinate_type<Point1, Point2>::type
Chris@102:                         >
Chris@102:                 > type;
Chris@102: };
Chris@102: 
Chris@102: 
Chris@102: } // namespace services
Chris@102: #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
Chris@102: 
Chris@102: 
Chris@102: }} // namespace strategy::distance
Chris@102: 
Chris@102: 
Chris@102: }} // namespace boost::geometry
Chris@102: 
Chris@102: 
Chris@102: #endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_ANDOYER_HPP