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_ALGORITHMS_DETAIL_AZIMUTH_HPP
Chris@102: #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_AZIMUTH_HPP
Chris@102: 
Chris@102: #include <boost/geometry/core/cs.hpp>
Chris@102: #include <boost/geometry/core/access.hpp>
Chris@102: #include <boost/geometry/core/radian_access.hpp>
Chris@102: #include <boost/geometry/core/tags.hpp>
Chris@102: 
Chris@102: #include <boost/geometry/util/math.hpp>
Chris@102: 
Chris@102: #include <boost/geometry/algorithms/not_implemented.hpp>
Chris@102: #include <boost/geometry/algorithms/detail/vincenty_inverse.hpp>
Chris@102: 
Chris@102: namespace boost { namespace geometry
Chris@102: {
Chris@102: 
Chris@102: // An azimuth is an angle between a vector/segment from origin to a point of
Chris@102: // interest and a reference vector. Typically north-based azimuth is used.
Chris@102: // North direction is used as a reference, angle is measured clockwise
Chris@102: // (North - 0deg, East - 90deg). For consistency in 2d cartesian CS
Chris@102: // the reference vector is Y axis, angle is measured clockwise.
Chris@102: // http://en.wikipedia.org/wiki/Azimuth
Chris@102: 
Chris@102: #ifndef DOXYGEN_NO_DISPATCH
Chris@102: namespace detail_dispatch
Chris@102: {
Chris@102: 
Chris@102: template <typename ReturnType, typename Tag>
Chris@102: struct azimuth
Chris@102:     : not_implemented<Tag>
Chris@102: {};
Chris@102: 
Chris@102: template <typename ReturnType>
Chris@102: struct azimuth<ReturnType, geographic_tag>
Chris@102: {
Chris@102:     template <typename P1, typename P2, typename Spheroid>
Chris@102:     static inline ReturnType apply(P1 const& p1, P2 const& p2, Spheroid const& spheroid)
Chris@102:     {
Chris@102:         return geometry::detail::vincenty_inverse<ReturnType>
Chris@102:                     ( get_as_radian<0>(p1), get_as_radian<1>(p1),
Chris@102:                       get_as_radian<0>(p2), get_as_radian<1>(p2),
Chris@102:                       spheroid ).azimuth12();
Chris@102:     }
Chris@102: 
Chris@102:     template <typename P1, typename P2>
Chris@102:     static inline ReturnType apply(P1 const& p1, P2 const& p2)
Chris@102:     {
Chris@102:         return apply(p1, p2, srs::spheroid<ReturnType>());
Chris@102:     }
Chris@102: };
Chris@102: 
Chris@102: template <typename ReturnType>
Chris@102: struct azimuth<ReturnType, spherical_equatorial_tag>
Chris@102: {
Chris@102:     template <typename P1, typename P2, typename Sphere>
Chris@102:     static inline ReturnType apply(P1 const& p1, P2 const& p2, Sphere const& /*unused*/)
Chris@102:     {
Chris@102:         // http://williams.best.vwh.net/avform.htm#Crs
Chris@102:         ReturnType dlon = get_as_radian<0>(p2) - get_as_radian<0>(p1);
Chris@102:         ReturnType cos_p2lat = cos(get_as_radian<1>(p2));
Chris@102: 
Chris@102:         // An optimization which should kick in often for Boxes
Chris@102:         //if ( math::equals(dlon, ReturnType(0)) )
Chris@102:         //if ( get<0>(p1) == get<0>(p2) )
Chris@102:         //{
Chris@102:         //    return - sin(get_as_radian<1>(p1)) * cos_p2lat);
Chris@102:         //}
Chris@102: 
Chris@102:         // "An alternative formula, not requiring the pre-computation of d"
Chris@102:         // In the formula below dlon is used as "d"
Chris@102:         return atan2(sin(dlon) * cos_p2lat,
Chris@102:             cos(get_as_radian<1>(p1)) * sin(get_as_radian<1>(p2))
Chris@102:             - sin(get_as_radian<1>(p1)) * cos_p2lat * cos(dlon));
Chris@102:     }
Chris@102: 
Chris@102:     template <typename P1, typename P2>
Chris@102:     static inline ReturnType apply(P1 const& p1, P2 const& p2)
Chris@102:     {
Chris@102:         return apply(p1, p2, 0); // dummy model
Chris@102:     }
Chris@102: };
Chris@102: 
Chris@102: template <typename ReturnType>
Chris@102: struct azimuth<ReturnType, spherical_polar_tag>
Chris@102:     : azimuth<ReturnType, spherical_equatorial_tag>
Chris@102: {};
Chris@102: 
Chris@102: template <typename ReturnType>
Chris@102: struct azimuth<ReturnType, cartesian_tag>
Chris@102: {
Chris@102:     template <typename P1, typename P2, typename Plane>
Chris@102:     static inline ReturnType apply(P1 const& p1, P2 const& p2, Plane const& /*unused*/)
Chris@102:     {
Chris@102:         ReturnType x = get<0>(p2) - get<0>(p1);
Chris@102:         ReturnType y = get<1>(p2) - get<1>(p1);
Chris@102: 
Chris@102:         // NOTE: azimuth 0 is at Y axis, increasing right
Chris@102:         // as in spherical/geographic where 0 is at North axis
Chris@102:         return atan2(x, y);
Chris@102:     }
Chris@102: 
Chris@102:     template <typename P1, typename P2>
Chris@102:     static inline ReturnType apply(P1 const& p1, P2 const& p2)
Chris@102:     {
Chris@102:         return apply(p1, p2, 0); // dummy model
Chris@102:     }
Chris@102: };
Chris@102: 
Chris@102: } // detail_dispatch
Chris@102: #endif // DOXYGEN_NO_DISPATCH
Chris@102: 
Chris@102: #ifndef DOXYGEN_NO_DETAIL
Chris@102: namespace detail
Chris@102: {
Chris@102: 
Chris@102: /// Calculate azimuth between two points.
Chris@102: /// The result is in radians.
Chris@102: template <typename ReturnType, typename Point1, typename Point2>
Chris@102: inline ReturnType azimuth(Point1 const& p1, Point2 const& p2)
Chris@102: {
Chris@102:     return detail_dispatch::azimuth
Chris@102:             <
Chris@102:                 ReturnType,
Chris@102:                 typename geometry::cs_tag<Point1>::type
Chris@102:             >::apply(p1, p2);
Chris@102: }
Chris@102: 
Chris@102: /// Calculate azimuth between two points.
Chris@102: /// The result is in radians.
Chris@102: template <typename ReturnType, typename Point1, typename Point2, typename Model>
Chris@102: inline ReturnType azimuth(Point1 const& p1, Point2 const& p2, Model const& model)
Chris@102: {
Chris@102:     return detail_dispatch::azimuth
Chris@102:             <
Chris@102:                 ReturnType,
Chris@102:                 typename geometry::cs_tag<Point1>::type
Chris@102:             >::apply(p1, p2, model);
Chris@102: }
Chris@102: 
Chris@102: } // namespace detail
Chris@102: #endif // DOXYGEN_NO_DETAIL
Chris@102: 
Chris@102: }} // namespace boost::geometry
Chris@102: 
Chris@102: #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_AZIMUTH_HPP