diff --git a/doc/doxy/doxygen_input/groups/groups.hpp b/doc/doxy/doxygen_input/groups/groups.hpp
index 1caf7cf63..d2659b561 100644
--- a/doc/doxy/doxygen_input/groups/groups.hpp
+++ b/doc/doxy/doxygen_input/groups/groups.hpp
@@ -2,8 +2,8 @@
 //
 // Copyright (c) 2011-2012 Barend Gehrels, Amsterdam, the Netherlands.
 
-// This file was modified by Oracle on 2018.
-// Modifications copyright (c) 2018, Oracle and/or its affiliates.
+// This file was modified by Oracle on 2018-2021.
+// Modifications copyright (c) 2018-2021, Oracle and/or its affiliates.
 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
 
@@ -20,6 +20,7 @@
 \defgroup area area: calculate area of a Geometry
 \defgroup arithmetic arithmetic: arithmetic operations on points
 \defgroup assign assign: assign values to geometries
+\defgroup azimuth azimuth: calculate azimuth of a segment defined by a pair of points
 \defgroup buffer buffer: calculate buffer of a geometry
 \defgroup centroid centroid: calculate centroid (center of gravity) of a geometry
 \defgroup clear clear: clear geometries
diff --git a/doc/imports.qbk b/doc/imports.qbk
index a184f1513..a5fdacf08 100644
--- a/doc/imports.qbk
+++ b/doc/imports.qbk
@@ -5,7 +5,7 @@
   Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
   Copyright (c) 2009-2012 Bruno Lalande, Paris, France.
 
-  Copyright (c) 2018, Oracle and/or its affiliates.
+  Copyright (c) 2018-2021, Oracle and/or its affiliates.
   Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
   Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
 
@@ -24,6 +24,8 @@
 [import src/examples/algorithms/assign_3d_point.cpp]
 [import src/examples/algorithms/assign_inverse.cpp]
 [import src/examples/algorithms/assign_points.cpp]
+[import src/examples/algorithms/azimuth.cpp]
+[import src/examples/algorithms/azimuth_strategy.cpp]
 [import src/examples/algorithms/buffer_with_strategies.cpp]
 [import src/examples/algorithms/clear.cpp]
 [import src/examples/algorithms/centroid.cpp]
diff --git a/doc/make_qbk.py b/doc/make_qbk.py
index 8a5cb6100..f4d7881e3 100755
--- a/doc/make_qbk.py
+++ b/doc/make_qbk.py
@@ -6,7 +6,7 @@
 #  Copyright (c) 2009-2012 Mateusz Loskot (mateusz@loskot.net), London, UK
 #  Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland
 #
-#  Copyright (c) 2018, Oracle and/or its affiliates.
+#  Copyright (c) 2018-2021, Oracle and/or its affiliates.
 #  Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
 #  Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
 #  Use, modification and distribution is subject to the Boost Software License,
@@ -96,7 +96,7 @@ def cs_to_quickbook(section):
 call_doxygen()
 
 algorithms = ["append", "assign", "make", "clear"
-    , "area", "buffer", "centroid", "convert", "correct", "covered_by"
+    , "area", "azimuth", "buffer", "centroid", "convert", "correct", "covered_by"
     , "convex_hull", "crosses", "densify", "difference"
     , "discrete_frechet_distance", "discrete_hausdorff_distance", "disjoint"
     , "distance", "envelope", "equals", "expand", "for_each", "is_empty"
diff --git a/doc/quickref.xml b/doc/quickref.xml
index bad6b5097..092a7f2ae 100644
--- a/doc/quickref.xml
+++ b/doc/quickref.xml
@@ -364,6 +364,10 @@
         coordinate values)
      </member>
     </simplelist>
+    <bridgehead renderas="sect3">Azimuth</bridgehead>
+    <simplelist type="vert" columns="1">
+     <member><link linkend="geometry.reference.algorithms.azimuth">azimuth</link></member>    
+    </simplelist>
     <bridgehead renderas="sect3">Buffer</bridgehead>
     <simplelist type="vert" columns="1">
      <member><link linkend="geometry.reference.algorithms.buffer">buffer</link></member>
diff --git a/doc/reference.qbk b/doc/reference.qbk
index aa0be56c0..a9270a5fb 100644
--- a/doc/reference.qbk
+++ b/doc/reference.qbk
@@ -6,8 +6,8 @@
   Copyright (c) 2009-2017 Bruno Lalande, Paris, France.
   Copyright (c) 2013-2017 Adam Wulkiewicz, Lodz, Poland.
 
-  This file was modified by Oracle on 2014, 2017, 2018.
-  Modifications copyright (c) 2014-2018, Oracle and/or its affiliates.
+  This file was modified by Oracle on 2014-2021.
+  Modifications copyright (c) 2014-2021, Oracle and/or its affiliates.
 
   Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
   Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
@@ -80,6 +80,10 @@
 
 [include generated/append.qbk]
 
+[section:azimuth azimuth]
+[include generated/azimuth.qbk]
+[endsect]
+
 [section:buffer buffer]
 [include generated/buffer.qbk]
 [endsect]
diff --git a/doc/reference/algorithms/azimuth.qbk b/doc/reference/algorithms/azimuth.qbk
new file mode 100644
index 000000000..5b2606661
--- /dev/null
+++ b/doc/reference/algorithms/azimuth.qbk
@@ -0,0 +1,22 @@
+[/============================================================================
+  Boost.Geometry
+
+  Copyright (c) 2021, Oracle and/or its affiliates.
+  Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
+  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)
+=============================================================================/]
+
+[def __this_function__ azimuth]
+
+[heading_conformance_no_ogc __this_function__]
+[note PostGIS contains an algorithm ST_Azimuth with the same functionality.
+  See the [@https://postgis.net/docs/ST_Azimuth.html PostGIS documentation]. 
+]
+
+[heading Behavior]
+The algorithm calculates the azimuth of a segment defined by a pair of points.
+
+[note The result is in radians.]
diff --git a/doc/src/examples/algorithms/azimuth.cpp b/doc/src/examples/algorithms/azimuth.cpp
new file mode 100644
index 000000000..0f799920a
--- /dev/null
+++ b/doc/src/examples/algorithms/azimuth.cpp
@@ -0,0 +1,42 @@
+// Boost.Geometry
+// QuickBook Example
+
+// Copyright (c) 2021, Oracle and/or its affiliates
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
+// 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)
+
+//[azimuth
+//` Shows how to calculate azimuth
+
+#include <iostream>
+
+#include <boost/geometry.hpp>
+#include <boost/geometry/geometries/point_xy.hpp>
+
+int main()
+{
+    typedef boost::geometry::model::d2::point_xy<double> point_type;
+
+    point_type p1(0, 0);
+    point_type p2(1, 1);
+
+    auto azimuth = boost::geometry::azimuth(p1, p2);
+
+    std::cout << "azimuth: " << azimuth << std::endl;
+
+    return 0;
+}
+
+//]
+
+//[azimuth_output
+/*`
+Output:
+[pre
+azimuth: 0.785398
+]
+*/
+//]
diff --git a/doc/src/examples/algorithms/azimuth_strategy.cpp b/doc/src/examples/algorithms/azimuth_strategy.cpp
new file mode 100644
index 000000000..2a8ef5178
--- /dev/null
+++ b/doc/src/examples/algorithms/azimuth_strategy.cpp
@@ -0,0 +1,46 @@
+// Boost.Geometry
+// QuickBook Example
+
+// Copyright (c) 2021, Oracle and/or its affiliates
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
+// 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)
+
+//[azimuth_strategy
+//` Shows how to calculate azimuth in geographic coordinate system
+
+#include <iostream>
+
+#include <boost/geometry.hpp>
+#include <boost/geometry/geometries/point.hpp>
+
+int main()
+{
+    namespace bg = boost::geometry;
+    typedef bg::model::point<double, 2, bg::cs::geographic<bg::degree> > point_type;
+
+    point_type p1(0, 0);
+    point_type p2(1, 1);
+
+    bg::srs::spheroid<double> spheroid(6378137.0, 6356752.3142451793);
+    bg::strategies::azimuth::geographic<> strategy(spheroid);
+
+    auto azimuth = boost::geometry::azimuth(p1, p2, strategy);
+
+    std::cout << "azimuth: " << azimuth << std::endl;
+
+    return 0;
+}
+
+//]
+
+//[azimuth_strategy_output
+/*`
+Output:
+[pre
+azimuth: 0.788674
+]
+*/
+//]
diff --git a/include/boost/geometry/algorithms/azimuth.hpp b/include/boost/geometry/algorithms/azimuth.hpp
new file mode 100644
index 000000000..2e97eed2f
--- /dev/null
+++ b/include/boost/geometry/algorithms/azimuth.hpp
@@ -0,0 +1,226 @@
+// Boost.Geometry
+
+// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
+
+// This file was modified by Oracle on 2014-2021.
+// Modifications copyright (c) 2014-2021, Oracle and/or its affiliates.
+// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
+// 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_ALGORITHMS_AZIMUTH_HPP
+#define BOOST_GEOMETRY_ALGORITHMS_AZIMUTH_HPP
+
+
+#include <boost/geometry/algorithms/not_implemented.hpp>
+
+#include <boost/geometry/core/cs.hpp>
+#include <boost/geometry/core/access.hpp>
+#include <boost/geometry/core/radian_access.hpp>
+#include <boost/geometry/core/tag.hpp>
+#include <boost/geometry/core/tags.hpp>
+
+#include <boost/geometry/strategies/azimuth.hpp>
+#include <boost/geometry/strategies/default_strategy.hpp>
+#include <boost/geometry/strategies/azimuth/cartesian.hpp>
+#include <boost/geometry/strategies/azimuth/geographic.hpp>
+#include <boost/geometry/strategies/azimuth/spherical.hpp>
+
+#include <boost/geometry/util/math.hpp>
+
+
+namespace boost { namespace geometry
+{
+
+
+#ifndef DOXYGEN_NO_DETAIL
+namespace detail
+{
+       
+} // namespace detail
+#endif // DOXYGEN_NO_DETAIL
+
+
+#ifndef DOXYGEN_NO_DISPATCH
+namespace dispatch
+{
+
+template
+<
+    typename Geometry1, typename Geometry2,
+    typename Tag1 = typename tag<Geometry1>::type,
+    typename Tag2 = typename tag<Geometry2>::type
+>
+struct azimuth : not_implemented<Tag1, Tag2>
+{};
+
+template <typename Point1, typename Point2>
+struct azimuth<Point1, Point2, point_tag, point_tag>
+{
+    template <typename Strategy>
+    static auto apply(Point1 const& p1, Point2 const& p2, Strategy const& strategy)
+    {
+        typedef typename decltype(strategy.azimuth())::template result_type
+            <
+                typename coordinate_type<Point1>::type,
+                typename coordinate_type<Point2>::type
+            >::type calc_t;
+
+        calc_t result = 0;
+        calc_t const x1 = geometry::get_as_radian<0>(p1);
+        calc_t const y1 = geometry::get_as_radian<1>(p1);
+        calc_t const x2 = geometry::get_as_radian<0>(p2);
+        calc_t const y2 = geometry::get_as_radian<1>(p2);
+
+        strategy.azimuth().apply(x1, y1, x2, y2, result);
+
+        // NOTE: It is not clear which units we should use for the result.
+        //   For now radians are always returned but a user could expect
+        //   e.g. something like this:
+        /*
+        bool const both_degree = std::is_same
+                <
+                    typename detail::cs_angular_units<Point1>::type,
+                    geometry::degree
+                >::value
+            && std::is_same
+                <
+                    typename detail::cs_angular_units<Point2>::type,
+                    geometry::degree
+                >::value;
+        if (both_degree)
+        {
+            result *= math::r2d<calc_t>();
+        }
+        */
+
+        return result;
+    }
+};
+
+} // namespace dispatch
+#endif // DOXYGEN_NO_DISPATCH
+
+
+namespace resolve_strategy
+{
+
+template
+<
+    typename Strategy,
+    bool IsUmbrella = strategies::detail::is_umbrella_strategy<Strategy>::value
+>
+struct azimuth
+{
+    template <typename P1, typename P2>
+    static auto apply(P1 const& p1, P2 const& p2, Strategy const& strategy)
+    {
+        return dispatch::azimuth<P1, P2>::apply(p1, p2, strategy);
+    }
+};
+
+template <typename Strategy>
+struct azimuth<Strategy, false>
+{
+    template <typename P1, typename P2>
+    static auto apply(P1 const& p1, P2 const& p2, Strategy const& strategy)
+    {
+        using strategies::azimuth::services::strategy_converter;
+        return dispatch::azimuth
+            <
+                P1, P2
+            >::apply(p1, p2, strategy_converter<Strategy>::get(strategy));
+    }
+};
+
+template <>
+struct azimuth<default_strategy, false>
+{
+    template <typename P1, typename P2>
+    static auto apply(P1 const& p1, P2 const& p2, default_strategy)
+    {
+        typedef typename strategies::azimuth::services::default_strategy
+            <
+                P1, P2
+            >::type strategy_type;
+
+        return dispatch::azimuth<P1, P2>::apply(p1, p2, strategy_type());
+    }
+};
+
+
+} // namespace resolve_strategy
+
+
+namespace resolve_variant
+{
+} // namespace resolve_variant
+
+
+/*!
+\brief Calculate azimuth of a segment defined by a pair of points.
+\ingroup azimuth
+\tparam Point1 Type of the first point of a segment.
+\tparam Point2 Type of the second point of a segment.
+\param point1 First point of a segment.
+\param point2 Second point of a segment.
+\return Azimuth in radians.
+
+\qbk{[include reference/algorithms/azimuth.qbk]}
+
+\qbk{
+[heading Example]
+[azimuth]
+[azimuth_output]
+}
+*/
+template <typename Point1, typename Point2>
+inline auto azimuth(Point1 const& point1, Point2 const& point2)
+{
+    concepts::check<Point1 const>();
+    concepts::check<Point2 const>();
+    
+    return resolve_strategy::azimuth
+            <
+                default_strategy
+            >::apply(point1, point2, default_strategy());
+}
+
+
+/*!
+\brief Calculate azimuth of a segment defined by a pair of points.
+\ingroup azimuth
+\tparam Point1 Type of the first point of a segment.
+\tparam Point2 Type of the second point of a segment.
+\tparam Strategy Type of an umbrella strategy defining azimuth strategy.
+\param point1 First point of a segment.
+\param point2 Second point of a segment.
+\param strategy Umbrella strategy defining azimuth strategy.
+\return Azimuth in radians.
+
+\qbk{distinguish,with strategy}
+\qbk{[include reference/algorithms/azimuth.qbk]}
+
+\qbk{
+[heading Example]
+[azimuth_strategy]
+[azimuth_strategy_output]
+}
+*/
+template <typename Point1, typename Point2, typename Strategy>
+inline auto azimuth(Point1 const& point1, Point2 const& point2, Strategy const& strategy)
+{
+    concepts::check<Point1 const>();
+    concepts::check<Point2 const>();
+
+    return resolve_strategy::azimuth<Strategy>::apply(point1, point2, strategy);
+}
+
+
+}} // namespace boost::geometry
+
+
+#endif // BOOST_GEOMETRY_ALGORITHMS_AZIMUTH_HPP
diff --git a/include/boost/geometry/algorithms/detail/azimuth.hpp b/include/boost/geometry/algorithms/detail/azimuth.hpp
deleted file mode 100644
index 21007778b..000000000
--- a/include/boost/geometry/algorithms/detail/azimuth.hpp
+++ /dev/null
@@ -1,152 +0,0 @@
-// Boost.Geometry (aka GGL, Generic Geometry Library)
-
-// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
-
-// This file was modified by Oracle on 2014, 2016, 2017.
-// Modifications copyright (c) 2014-2017, Oracle and/or its affiliates.
-
-// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
-// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
-
-// 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_ALGORITHMS_DETAIL_AZIMUTH_HPP
-#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_AZIMUTH_HPP
-
-
-#include <boost/geometry/algorithms/not_implemented.hpp>
-
-#include <boost/geometry/core/cs.hpp>
-#include <boost/geometry/core/access.hpp>
-#include <boost/geometry/core/radian_access.hpp>
-#include <boost/geometry/core/tags.hpp>
-
-#include <boost/geometry/formulas/spherical.hpp>
-#include <boost/geometry/formulas/vincenty_inverse.hpp>
-
-#include <boost/geometry/srs/spheroid.hpp>
-
-#include <boost/geometry/util/math.hpp>
-
-
-namespace boost { namespace geometry
-{
-
-// An azimuth is an angle between a vector/segment from origin to a point of
-// interest and a reference vector. Typically north-based azimuth is used.
-// North direction is used as a reference, angle is measured clockwise
-// (North - 0deg, East - 90deg). For consistency in 2d cartesian CS
-// the reference vector is Y axis, angle is measured clockwise.
-// http://en.wikipedia.org/wiki/Azimuth
-
-#ifndef DOXYGEN_NO_DISPATCH
-namespace detail_dispatch
-{
-
-template <typename ReturnType, typename Tag>
-struct azimuth
-    : not_implemented<Tag>
-{};
-
-template <typename ReturnType>
-struct azimuth<ReturnType, geographic_tag>
-{
-    template <typename P1, typename P2, typename Spheroid>
-    static inline ReturnType apply(P1 const& p1, P2 const& p2, Spheroid const& spheroid)
-    {
-        return geometry::formula::vincenty_inverse<ReturnType, false, true>().apply
-                    ( get_as_radian<0>(p1), get_as_radian<1>(p1),
-                      get_as_radian<0>(p2), get_as_radian<1>(p2),
-                      spheroid ).azimuth;
-    }
-
-    template <typename P1, typename P2>
-    static inline ReturnType apply(P1 const& p1, P2 const& p2)
-    {
-        return apply(p1, p2, srs::spheroid<ReturnType>());
-    }
-};
-
-template <typename ReturnType>
-struct azimuth<ReturnType, spherical_equatorial_tag>
-{
-    template <typename P1, typename P2, typename Sphere>
-    static inline ReturnType apply(P1 const& p1, P2 const& p2, Sphere const& /*unused*/)
-    {
-        return geometry::formula::spherical_azimuth<ReturnType, false>
-                    ( get_as_radian<0>(p1), get_as_radian<1>(p1),
-                      get_as_radian<0>(p2), get_as_radian<1>(p2)).azimuth;
-    }
-
-    template <typename P1, typename P2>
-    static inline ReturnType apply(P1 const& p1, P2 const& p2)
-    {
-        return apply(p1, p2, 0); // dummy model
-    }
-};
-
-template <typename ReturnType>
-struct azimuth<ReturnType, spherical_polar_tag>
-    : azimuth<ReturnType, spherical_equatorial_tag>
-{};
-
-template <typename ReturnType>
-struct azimuth<ReturnType, cartesian_tag>
-{
-    template <typename P1, typename P2, typename Plane>
-    static inline ReturnType apply(P1 const& p1, P2 const& p2, Plane const& /*unused*/)
-    {
-        ReturnType x = get<0>(p2) - get<0>(p1);
-        ReturnType y = get<1>(p2) - get<1>(p1);
-
-        // NOTE: azimuth 0 is at Y axis, increasing right
-        // as in spherical/geographic where 0 is at North axis
-        return atan2(x, y);
-    }
-
-    template <typename P1, typename P2>
-    static inline ReturnType apply(P1 const& p1, P2 const& p2)
-    {
-        return apply(p1, p2, 0); // dummy model
-    }
-};
-
-} // detail_dispatch
-#endif // DOXYGEN_NO_DISPATCH
-
-#ifndef DOXYGEN_NO_DETAIL
-namespace detail
-{
-
-/// Calculate azimuth between two points.
-/// The result is in radians.
-template <typename ReturnType, typename Point1, typename Point2>
-inline ReturnType azimuth(Point1 const& p1, Point2 const& p2)
-{
-    return detail_dispatch::azimuth
-            <
-                ReturnType,
-                typename geometry::cs_tag<Point1>::type
-            >::apply(p1, p2);
-}
-
-/// Calculate azimuth between two points.
-/// The result is in radians.
-template <typename ReturnType, typename Point1, typename Point2, typename Model>
-inline ReturnType azimuth(Point1 const& p1, Point2 const& p2, Model const& model)
-{
-    return detail_dispatch::azimuth
-            <
-                ReturnType,
-                typename geometry::cs_tag<Point1>::type
-            >::apply(p1, p2, model);
-}
-
-} // namespace detail
-#endif // DOXYGEN_NO_DETAIL
-
-}} // namespace boost::geometry
-
-#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_AZIMUTH_HPP
diff --git a/include/boost/geometry/algorithms/detail/course.hpp b/include/boost/geometry/algorithms/detail/course.hpp
deleted file mode 100644
index 37424941b..000000000
--- a/include/boost/geometry/algorithms/detail/course.hpp
+++ /dev/null
@@ -1,40 +0,0 @@
-// Boost.Geometry (aka GGL, Generic Geometry Library)
-
-// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
-
-// This file was modified by Oracle on 2014.
-// Modifications copyright (c) 2014, Oracle and/or its affiliates.
-
-// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
-
-// 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_ALGORITHMS_DETAIL_COURSE_HPP
-#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_COURSE_HPP
-
-#include <boost/geometry/algorithms/detail/azimuth.hpp>
-
-namespace boost { namespace geometry
-{
-
-#ifndef DOXYGEN_NO_DETAIL
-namespace detail
-{
-
-/// Calculate course (bearing) between two points.
-///
-/// NOTE: left for convenience and temporary backward compatibility
-template <typename ReturnType, typename Point1, typename Point2>
-inline ReturnType course(Point1 const& p1, Point2 const& p2)
-{
-    return azimuth<ReturnType>(p1, p2);
-}
-
-} // namespace detail
-#endif // DOXYGEN_NO_DETAIL
-
-}} // namespace boost::geometry
-
-#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_COURSE_HPP
diff --git a/include/boost/geometry/geometry.hpp b/include/boost/geometry/geometry.hpp
index 24b12fab7..9aff956c9 100644
--- a/include/boost/geometry/geometry.hpp
+++ b/include/boost/geometry/geometry.hpp
@@ -4,8 +4,8 @@
 // Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
 // Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
 
-// This file was modified by Oracle on 2014-2020.
-// Modifications copyright (c) 2014-2020 Oracle and/or its affiliates.
+// This file was modified by Oracle on 2014-2021.
+// Modifications copyright (c) 2014-2021 Oracle and/or its affiliates.
 
 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
 // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
@@ -60,6 +60,7 @@ BOOST_PRAGMA_MESSAGE("CAUTION: Define BOOST_GEOMETRY_DISABLE_DEPRECATED_03_WARNI
 #include <boost/geometry/algorithms/append.hpp>
 #include <boost/geometry/algorithms/area.hpp>
 #include <boost/geometry/algorithms/assign.hpp>
+#include <boost/geometry/algorithms/azimuth.hpp>
 #include <boost/geometry/algorithms/buffer.hpp>
 #include <boost/geometry/algorithms/centroid.hpp>
 #include <boost/geometry/algorithms/clear.hpp>
diff --git a/include/boost/geometry/strategies/azimuth.hpp b/include/boost/geometry/strategies/azimuth.hpp
index 1573ad22a..7d1e105d9 100644
--- a/include/boost/geometry/strategies/azimuth.hpp
+++ b/include/boost/geometry/strategies/azimuth.hpp
@@ -1,6 +1,6 @@
 // Boost.Geometry (aka GGL, Generic Geometry Library)
 
-// Copyright (c) 2016-2020 Oracle and/or its affiliates.
+// Copyright (c) 2016-2021 Oracle and/or its affiliates.
 // Contributed and/or modified by Vissarion Fisikopoulos, on behalf of Oracle
 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
 
@@ -26,7 +26,7 @@ namespace strategy { namespace azimuth { namespace services
 \tparam CSTag tag of coordinate system
 \tparam CalculationType \tparam_calculation
 */
-template <typename CSTag, typename CalculationType = void>
+template <typename CSTag>
 struct default_strategy
 {
     BOOST_GEOMETRY_STATIC_ASSERT_FALSE(
diff --git a/include/boost/geometry/strategies/azimuth/cartesian.hpp b/include/boost/geometry/strategies/azimuth/cartesian.hpp
new file mode 100644
index 000000000..ac6136f89
--- /dev/null
+++ b/include/boost/geometry/strategies/azimuth/cartesian.hpp
@@ -0,0 +1,63 @@
+// Boost.Geometry
+
+// Copyright (c) 2021, Oracle and/or its affiliates.
+
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
+// Licensed under the Boost Software License version 1.0.
+// http://www.boost.org/users/license.html
+
+#ifndef BOOST_GEOMETRY_STRATEGIES_AZIMUTH_CARTESIAN_HPP
+#define BOOST_GEOMETRY_STRATEGIES_AZIMUTH_CARTESIAN_HPP
+
+
+// TODO: move this file to boost/geometry/strategy
+#include <boost/geometry/strategies/cartesian/azimuth.hpp>
+
+#include <boost/geometry/strategies/azimuth/services.hpp>
+#include <boost/geometry/strategies/detail.hpp>
+
+
+namespace boost { namespace geometry
+{
+
+namespace strategies { namespace azimuth
+{
+
+template <typename CalculationType = void>
+struct cartesian : strategies::detail::cartesian_base
+{
+    static auto azimuth()
+    {
+        return strategy::azimuth::cartesian<CalculationType>();
+    }
+};
+
+
+namespace services
+{
+
+template <typename Point1, typename Point2>
+struct default_strategy<Point1, Point2, cartesian_tag, cartesian_tag>
+{
+    using type = strategies::azimuth::cartesian<>;
+};
+
+
+template <typename CT>
+struct strategy_converter<strategy::azimuth::cartesian<CT> >
+{
+    static auto get(strategy::azimuth::cartesian<CT> const&)
+    {
+        return strategies::azimuth::cartesian<CT>();
+    }
+};
+
+
+} // namespace services
+
+}} // namespace strategies::azimuth
+
+}} // namespace boost::geometry
+
+#endif // BOOST_GEOMETRY_STRATEGIES_AZIMUTH_CARTESIAN_HPP
diff --git a/include/boost/geometry/strategies/azimuth/geographic.hpp b/include/boost/geometry/strategies/azimuth/geographic.hpp
new file mode 100644
index 000000000..a4fac5e60
--- /dev/null
+++ b/include/boost/geometry/strategies/azimuth/geographic.hpp
@@ -0,0 +1,81 @@
+// Boost.Geometry
+
+// Copyright (c) 2021, Oracle and/or its affiliates.
+
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
+// Licensed under the Boost Software License version 1.0.
+// http://www.boost.org/users/license.html
+
+#ifndef BOOST_GEOMETRY_STRATEGIES_AZIMUTH_GEOGRAPHIC_HPP
+#define BOOST_GEOMETRY_STRATEGIES_AZIMUTH_GEOGRAPHIC_HPP
+
+
+// TODO: move this file to boost/geometry/strategy
+#include <boost/geometry/strategies/geographic/azimuth.hpp>
+
+#include <boost/geometry/strategies/azimuth/services.hpp>
+#include <boost/geometry/strategies/detail.hpp>
+
+
+namespace boost { namespace geometry
+{
+
+namespace strategies { namespace azimuth
+{
+
+template
+<
+    typename FormulaPolicy = strategy::andoyer,
+    typename Spheroid = srs::spheroid<double>,
+    typename CalculationType = void
+>
+class geographic : strategies::detail::geographic_base<Spheroid>
+{
+    using base_t = strategies::detail::geographic_base<Spheroid>;
+
+public:
+    geographic()
+        : base_t()
+    {}
+
+    explicit geographic(Spheroid const& spheroid)
+        : base_t(spheroid)
+    {}
+
+    auto azimuth() const
+    {
+        return strategy::azimuth::geographic
+            <
+                FormulaPolicy, Spheroid, CalculationType
+            >(base_t::m_spheroid);
+    }
+};
+
+
+namespace services
+{
+
+template <typename Point1, typename Point2>
+struct default_strategy<Point1, Point2, geographic_tag, geographic_tag>
+{
+    using type = strategies::azimuth::geographic<>;
+};
+
+
+template <typename FP, typename S, typename CT>
+struct strategy_converter<strategy::azimuth::geographic<FP, S, CT> >
+{
+    static auto get(strategy::azimuth::geographic<FP, S, CT> const& strategy)
+    {
+        return strategies::azimuth::geographic<FP, S, CT>(strategy.model());
+    }
+};
+
+} // namespace services
+
+}} // namespace strategies::azimuth
+
+}} // namespace boost::geometry
+
+#endif // BOOST_GEOMETRY_STRATEGIES_AZIMUTH_GEOGRAPHIC_HPP
diff --git a/include/boost/geometry/strategies/azimuth/services.hpp b/include/boost/geometry/strategies/azimuth/services.hpp
new file mode 100644
index 000000000..11980e5a4
--- /dev/null
+++ b/include/boost/geometry/strategies/azimuth/services.hpp
@@ -0,0 +1,56 @@
+// Boost.Geometry
+
+// Copyright (c) 2021, Oracle and/or its affiliates.
+
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
+// Licensed under the Boost Software License version 1.0.
+// http://www.boost.org/users/license.html
+
+#ifndef BOOST_GEOMETRY_STRATEGIES_AZIMUTH_SERVICES_HPP
+#define BOOST_GEOMETRY_STRATEGIES_AZIMUTH_SERVICES_HPP
+
+
+#include <boost/geometry/core/cs.hpp>
+#include <boost/geometry/core/static_assert.hpp>
+
+
+namespace boost { namespace geometry
+{
+
+
+namespace strategies { namespace azimuth {
+
+namespace services
+{
+
+template
+<
+    typename Point1, typename Point2,
+    typename CSTag1 = typename geometry::cs_tag<Point1>::type,
+    typename CSTag2 = typename geometry::cs_tag<Point2>::type
+>
+struct default_strategy
+{
+    BOOST_GEOMETRY_STATIC_ASSERT_FALSE(
+        "Not implemented for these Points' coordinate systems.",
+        Point1, Point2, CSTag1, CSTag2);
+};
+
+template <typename Strategy>
+struct strategy_converter
+{
+    BOOST_GEOMETRY_STATIC_ASSERT_FALSE(
+        "Not implemented for this Strategy.",
+        Strategy);
+};
+
+
+} // namespace services
+
+}} // namespace strategies::azimuth
+
+
+}} // namespace boost::geometry
+
+#endif // BOOST_GEOMETRY_STRATEGIES_AZIMUTH_SERVICES_HPP
diff --git a/include/boost/geometry/strategies/azimuth/spherical.hpp b/include/boost/geometry/strategies/azimuth/spherical.hpp
new file mode 100644
index 000000000..a87cf927b
--- /dev/null
+++ b/include/boost/geometry/strategies/azimuth/spherical.hpp
@@ -0,0 +1,67 @@
+// Boost.Geometry
+
+// Copyright (c) 2021, Oracle and/or its affiliates.
+
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
+// Licensed under the Boost Software License version 1.0.
+// http://www.boost.org/users/license.html
+
+#ifndef BOOST_GEOMETRY_STRATEGIES_AZIMUTH_SPHERICAL_HPP
+#define BOOST_GEOMETRY_STRATEGIES_AZIMUTH_SPHERICAL_HPP
+
+
+// TODO: move this file to boost/geometry/strategy
+#include <boost/geometry/strategies/spherical/azimuth.hpp>
+
+#include <boost/geometry/strategies/azimuth/services.hpp>
+#include <boost/geometry/strategies/detail.hpp>
+
+
+namespace boost { namespace geometry
+{
+
+namespace strategies { namespace azimuth
+{
+
+template <typename CalculationType = void>
+class spherical : strategies::detail::spherical_base<void>
+{
+    using base_t = strategies::detail::spherical_base<void>;
+
+public:
+    
+    static auto azimuth()
+    {
+        return strategy::azimuth::spherical<CalculationType>();
+    }
+};
+
+
+namespace services
+{
+
+
+template <typename Point1, typename Point2>
+struct default_strategy<Point1, Point2, spherical_equatorial_tag, spherical_equatorial_tag>
+{
+    using type = strategies::azimuth::spherical<>;
+};
+
+template <typename CT>
+struct strategy_converter<strategy::azimuth::spherical<CT> >
+{
+    static auto get(strategy::azimuth::spherical<CT> const&)
+    {
+        return strategies::azimuth::spherical<CT>();
+    }
+};
+
+
+} // namespace services
+
+}} // namespace strategies::azimuth
+
+}} // namespace boost::geometry
+
+#endif // BOOST_GEOMETRY_STRATEGIES_AZIMUTH_SPHERICAL_HPP
diff --git a/include/boost/geometry/strategies/cartesian.hpp b/include/boost/geometry/strategies/cartesian.hpp
index 76d84caf3..9b4559ad4 100644
--- a/include/boost/geometry/strategies/cartesian.hpp
+++ b/include/boost/geometry/strategies/cartesian.hpp
@@ -1,6 +1,6 @@
 // Boost.Geometry
 
-// Copyright (c) 2020, Oracle and/or its affiliates.
+// Copyright (c) 2020-2021, Oracle and/or its affiliates.
 
 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
 
@@ -12,6 +12,7 @@
 
 
 #include <boost/geometry/strategies/area/cartesian.hpp>
+#include <boost/geometry/strategies/azimuth/cartesian.hpp>
 #include <boost/geometry/strategies/convex_hull/cartesian.hpp>
 #include <boost/geometry/strategies/envelope/cartesian.hpp>
 #include <boost/geometry/strategies/expand/cartesian.hpp>
@@ -35,7 +36,10 @@ class cartesian
 {
 public:
 
-    // point_order
+    static auto azimuth()
+    {
+        return strategy::azimuth::cartesian<CalculationType>();
+    }
 
     static auto point_order()
     {
diff --git a/include/boost/geometry/strategies/cartesian/azimuth.hpp b/include/boost/geometry/strategies/cartesian/azimuth.hpp
index 2ff16e945..74a9ea12b 100644
--- a/include/boost/geometry/strategies/cartesian/azimuth.hpp
+++ b/include/boost/geometry/strategies/cartesian/azimuth.hpp
@@ -1,6 +1,6 @@
 // Boost.Geometry (aka GGL, Generic Geometry Library)
 
-// Copyright (c) 2016-2018 Oracle and/or its affiliates.
+// Copyright (c) 2016-2021 Oracle and/or its affiliates.
 // Contributed and/or modified by Vissarion Fisikopoulos, on behalf of Oracle
 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
 
@@ -11,32 +11,81 @@
 #ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_AZIMUTH_HPP
 #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_AZIMUTH_HPP
 
+#include <cmath>
+
 #include <boost/geometry/core/tags.hpp>
 
 #include <boost/geometry/strategies/azimuth.hpp>
 
+#include <boost/geometry/util/promote_floating_point.hpp>
+#include <boost/geometry/util/select_most_precise.hpp>
+
 namespace boost { namespace geometry
 {
 
 namespace strategy { namespace azimuth
 {
 
-template
-<
-    typename CalculationType = void
->
+template <typename CalculationType = void>
 class cartesian
-{};
+{
+public:
+    template <typename T1, typename T2>
+    struct result_type
+        : geometry::select_most_precise
+              <
+                  // NOTE: this promotes any integer type to double
+                  typename geometry::promote_floating_point<T1, double>::type,
+                  typename geometry::promote_floating_point<T2, double>::type,
+                  CalculationType
+              >
+    {};
+
+    template <typename T1, typename T2, typename Result>
+    static inline void apply(T1 const& x1, T1 const& y1,
+                             T2 const& x2, T2 const& y2,
+                             Result& a1, Result& a2)
+    {
+        compute(x1, y1, x2, y2, a1, a2);
+    }
+    template <typename T1, typename T2, typename Result>
+    static inline void apply(T1 const& x1, T1 const& y1,
+                             T2 const& x2, T2 const& y2,
+                             Result& a1)
+    {
+        compute(x1, y1, x2, y2, a1, a1);
+    }
+    template <typename T1, typename T2, typename Result>
+    static inline void apply_reverse(T1 const& x1, T1 const& y1,
+                                     T2 const& x2, T2 const& y2,
+                                     Result& a2)
+    {
+        compute(x1, y1, x2, y2, a2, a2);
+    }
+
+private:
+    template <typename T1, typename T2, typename Result>
+    static inline void compute(T1 const& x1, T1 const& y1,
+                               T2 const& x2, T2 const& y2,
+                               Result& a1, Result& a2)
+    {
+        typedef typename result_type<T1, T2>::type calc_t;
+
+        // NOTE: azimuth 0 is at Y axis, increasing right
+        // as in spherical/geographic where 0 is at North axis
+        a1 = a2 = atan2(calc_t(x2) - calc_t(x1), calc_t(y2) - calc_t(y1));
+    }
+};
 
 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 
 namespace services
 {
 
-template <typename CalculationType>
-struct default_strategy<cartesian_tag, CalculationType>
+template <>
+struct default_strategy<cartesian_tag>
 {
-    typedef strategy::azimuth::cartesian<CalculationType> type;
+    typedef strategy::azimuth::cartesian<> type;
 };
 
 }
diff --git a/include/boost/geometry/strategies/geographic.hpp b/include/boost/geometry/strategies/geographic.hpp
index b20b2574a..815f1e3f3 100644
--- a/include/boost/geometry/strategies/geographic.hpp
+++ b/include/boost/geometry/strategies/geographic.hpp
@@ -1,6 +1,6 @@
 // Boost.Geometry
 
-// Copyright (c) 2020, Oracle and/or its affiliates.
+// Copyright (c) 2020-2021, Oracle and/or its affiliates.
 
 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
 
@@ -12,6 +12,7 @@
 
 
 #include <boost/geometry/strategies/area/geographic.hpp>
+#include <boost/geometry/strategies/azimuth/geographic.hpp>
 #include <boost/geometry/strategies/convex_hull/geographic.hpp>
 #include <boost/geometry/strategies/envelope/geographic.hpp>
 #include <boost/geometry/strategies/expand/geographic.hpp>
@@ -50,6 +51,14 @@ public:
         : base_t(spheroid)
     {}
 
+    auto azimuth() const
+    {
+        return strategy::azimuth::geographic
+            <
+                FormulaPolicy, Spheroid, CalculationType
+            >(base_t::m_spheroid);
+    }
+
     auto point_order() const
     {
         return strategy::point_order::geographic
diff --git a/include/boost/geometry/strategies/geographic/azimuth.hpp b/include/boost/geometry/strategies/geographic/azimuth.hpp
index ca280313d..9e6914d92 100644
--- a/include/boost/geometry/strategies/geographic/azimuth.hpp
+++ b/include/boost/geometry/strategies/geographic/azimuth.hpp
@@ -1,6 +1,6 @@
 // Boost.Geometry (aka GGL, Generic Geometry Library)
 
-// Copyright (c) 2016-2020 Oracle and/or its affiliates.
+// Copyright (c) 2016-2021 Oracle and/or its affiliates.
 // Contributed and/or modified by Vissarion Fisikopoulos, on behalf of Oracle
 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
 
@@ -19,6 +19,8 @@
 #include <boost/geometry/strategies/azimuth.hpp>
 #include <boost/geometry/strategies/geographic/parameters.hpp>
 
+#include <boost/geometry/util/select_most_precise.hpp>
+
 
 namespace boost { namespace geometry
 {
@@ -34,7 +36,14 @@ template
 >
 class geographic
 {
-public :
+public:
+    template <typename T1, typename T2>
+    struct result_type
+        : geometry::select_most_precise
+              <
+                  T1, T2, CalculationType
+              >
+    {};
 
     typedef Spheroid model_type;
 
@@ -51,28 +60,28 @@ public :
         return m_spheroid;
     }
 
-    template <typename T>
-    inline void apply(T const& lon1_rad, T const& lat1_rad,
-                      T const& lon2_rad, T const& lat2_rad,
-                      T& a1, T& a2) const
+    template <typename T1, typename T2, typename Result>
+    inline void apply(T1 const& lon1_rad, T1 const& lat1_rad,
+                      T2 const& lon2_rad, T2 const& lat2_rad,
+                      Result& a1, Result& a2) const
     {
         compute<true, true>(lon1_rad, lat1_rad,
                             lon2_rad, lat2_rad,
                             a1, a2);
     }
-    template <typename T>
-    inline void apply(T const& lon1_rad, T const& lat1_rad,
-                      T const& lon2_rad, T const& lat2_rad,
-                      T& a1) const
+    template <typename T1, typename T2, typename Result>
+    inline void apply(T1 const& lon1_rad, T1 const& lat1_rad,
+                      T2 const& lon2_rad, T2 const& lat2_rad,
+                      Result& a1) const
     {
         compute<true, false>(lon1_rad, lat1_rad,
                              lon2_rad, lat2_rad,
                              a1, a1);
     }
-    template <typename T>
-    inline void apply_reverse(T const& lon1_rad, T const& lat1_rad,
-                              T const& lon2_rad, T const& lat2_rad,
-                              T& a2) const
+    template <typename T1, typename T2, typename Result>
+    inline void apply_reverse(T1 const& lon1_rad, T1 const& lat1_rad,
+                              T2 const& lon2_rad, T2 const& lat2_rad,
+                              Result& a2) const
     {
         compute<false, true>(lon1_rad, lat1_rad,
                              lon2_rad, lat2_rad,
@@ -85,16 +94,13 @@ private :
     <
         bool EnableAzimuth,
         bool EnableReverseAzimuth,
-        typename T
+        typename T1, typename T2, typename Result
     >
-    inline void compute(T const& lon1_rad, T const& lat1_rad,
-                        T const& lon2_rad, T const& lat2_rad,
-                        T& a1, T& a2) const
+    inline void compute(T1 const& lon1_rad, T1 const& lat1_rad,
+                        T2 const& lon2_rad, T2 const& lat2_rad,
+                        Result& a1, Result& a2) const
     {
-        typedef std::conditional_t
-            <
-                std::is_void<CalculationType>::value, T, CalculationType
-            > calc_t;
+        typedef typename result_type<T1, T2>::type calc_t;
 
         typedef typename FormulaPolicy::template inverse
             <
@@ -127,14 +133,13 @@ private :
 namespace services
 {
 
-template <typename CalculationType>
-struct default_strategy<geographic_tag, CalculationType>
+template <>
+struct default_strategy<geographic_tag>
 {
     typedef strategy::azimuth::geographic
         <
             strategy::andoyer,
-            srs::spheroid<double>,
-            CalculationType
+            srs::spheroid<double>
         > type;
 };
 
diff --git a/include/boost/geometry/strategies/spherical.hpp b/include/boost/geometry/strategies/spherical.hpp
index f7d4cc8d7..07495e857 100644
--- a/include/boost/geometry/strategies/spherical.hpp
+++ b/include/boost/geometry/strategies/spherical.hpp
@@ -1,6 +1,6 @@
 // Boost.Geometry
 
-// Copyright (c) 2020, Oracle and/or its affiliates.
+// Copyright (c) 2020-2021, Oracle and/or its affiliates.
 
 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
 
@@ -12,6 +12,7 @@
 
 
 #include <boost/geometry/strategies/area/spherical.hpp>
+#include <boost/geometry/strategies/azimuth/spherical.hpp>
 #include <boost/geometry/strategies/convex_hull/spherical.hpp>
 #include <boost/geometry/strategies/envelope/spherical.hpp>
 #include <boost/geometry/strategies/expand/spherical.hpp>
@@ -53,7 +54,10 @@ public:
             >::apply(radius_or_sphere);
     }
 
-    // point_order
+    static auto azimuth()
+    {
+        return strategy::azimuth::spherical<CalculationType>();
+    }
 
     static auto point_order()
     {
diff --git a/include/boost/geometry/strategies/spherical/azimuth.hpp b/include/boost/geometry/strategies/spherical/azimuth.hpp
index 32339057b..6fcb15ac9 100644
--- a/include/boost/geometry/strategies/spherical/azimuth.hpp
+++ b/include/boost/geometry/strategies/spherical/azimuth.hpp
@@ -1,6 +1,6 @@
 // Boost.Geometry (aka GGL, Generic Geometry Library)
 
-// Copyright (c) 2016-2020 Oracle and/or its affiliates.
+// Copyright (c) 2016-2021 Oracle and/or its affiliates.
 // Contributed and/or modified by Vissarion Fisikopoulos, on behalf of Oracle
 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
 
@@ -14,9 +14,12 @@
 
 #include <type_traits>
 
-#include <boost/geometry/strategies/azimuth.hpp>
 #include <boost/geometry/formulas/spherical.hpp>
 
+#include <boost/geometry/strategies/azimuth.hpp>
+
+#include <boost/geometry/util/select_most_precise.hpp>
+
 
 namespace boost { namespace geometry
 {
@@ -24,61 +27,58 @@ namespace boost { namespace geometry
 namespace strategy { namespace azimuth
 {
 
-template
-<
-    typename CalculationType = void
->
+template <typename CalculationType = void>
 class spherical
 {
-public :
+public:
+    template <typename T1, typename T2>
+    struct result_type
+        : geometry::select_most_precise
+              <
+                  T1, T2, CalculationType
+              >
+    {};
 
-    inline spherical()
-    {}
-
-    template <typename T>
-    inline void apply(T const& lon1_rad, T const& lat1_rad,
-                      T const& lon2_rad, T const& lat2_rad,
-                      T& a1, T& a2) const
+    template <typename T1, typename T2, typename Result>
+    static inline void apply(T1 const& lon1_rad, T1 const& lat1_rad,
+                             T2 const& lon2_rad, T2 const& lat2_rad,
+                             Result& a1, Result& a2)
     {
         compute<true, true>(lon1_rad, lat1_rad,
                             lon2_rad, lat2_rad,
                             a1, a2);
     }
-    template <typename T>
-    inline void apply(T const& lon1_rad, T const& lat1_rad,
-                      T const& lon2_rad, T const& lat2_rad,
-                      T& a1) const
+    template <typename T1, typename T2, typename Result>
+    static inline void apply(T1 const& lon1_rad, T1 const& lat1_rad,
+                             T2 const& lon2_rad, T2 const& lat2_rad,
+                             Result& a1)
     {
         compute<true, false>(lon1_rad, lat1_rad,
                              lon2_rad, lat2_rad,
                              a1, a1);
     }
-    template <typename T>
-    inline void apply_reverse(T const& lon1_rad, T const& lat1_rad,
-                              T const& lon2_rad, T const& lat2_rad,
-                              T& a2) const
+    template <typename T1, typename T2, typename Result>
+    static inline void apply_reverse(T1 const& lon1_rad, T1 const& lat1_rad,
+                                     T2 const& lon2_rad, T2 const& lat2_rad,
+                                     Result& a2)
     {
         compute<false, true>(lon1_rad, lat1_rad,
                              lon2_rad, lat2_rad,
                              a2, a2);
     }
 
-private :
-
+private:
     template
     <
         bool EnableAzimuth,
         bool EnableReverseAzimuth,
-        typename T
+        typename T1, typename T2, typename Result
     >
-    inline void compute(T const& lon1_rad, T const& lat1_rad,
-                        T const& lon2_rad, T const& lat2_rad,
-                        T& a1, T& a2) const
+    static inline void compute(T1 const& lon1_rad, T1 const& lat1_rad,
+                               T2 const& lon2_rad, T2 const& lat2_rad,
+                               Result& a1, Result& a2)
     {
-        typedef std::conditional_t
-            <
-                std::is_void<CalculationType>::value, T, CalculationType
-            > calc_t;
+        typedef typename result_type<T1, T2>::type calc_t;
 
         geometry::formula::result_spherical<calc_t>
             result = geometry::formula::spherical_azimuth
@@ -104,19 +104,12 @@ private :
 namespace services
 {
 
-template <typename CalculationType>
-struct default_strategy<spherical_equatorial_tag, CalculationType>
+template <>
+struct default_strategy<spherical_equatorial_tag>
 {
-    typedef strategy::azimuth::spherical<CalculationType> type;
+    typedef strategy::azimuth::spherical<> type;
 };
 
-/*
-template <typename CalculationType>
-struct default_strategy<spherical_polar_tag, CalculationType>
-{
-    typedef strategy::azimuth::spherical<CalculationType> type;
-};
-*/
 }
 
 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
diff --git a/test/algorithms/Jamfile b/test/algorithms/Jamfile
index 3846cdb12..db2d1ce67 100644
--- a/test/algorithms/Jamfile
+++ b/test/algorithms/Jamfile
@@ -5,8 +5,8 @@
 # Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
 # Copyright (c) 2018 Adam Wulkiewicz, Lodz, Poland.
 #
-# This file was modified by Oracle on 2014, 2015, 2016, 2017.
-# Modifications copyright (c) 2014-2017, Oracle and/or its affiliates.
+# This file was modified by Oracle on 2014-2021.
+# Modifications copyright (c) 2014-2021, Oracle and/or its affiliates.
 #
 # Contributed and/or modified by Vissarion Fisikopoulos, on behalf of Oracle
 # Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
@@ -20,6 +20,7 @@ test-suite boost-geometry-algorithms
     :
     [ run append.cpp                   : : : : algorithms_append ]
     [ run assign.cpp                   : : : : algorithms_assign ]
+    [ run azimuth.cpp                  : : : : algorithms_azimuth ]
     [ run centroid.cpp                 : : : : algorithms_centroid ]
     [ run centroid_multi.cpp           : : : : algorithms_centroid_multi ]
     [ run comparable_distance.cpp      : : : : algorithms_comparable_distance ]
diff --git a/test/algorithms/azimuth.cpp b/test/algorithms/azimuth.cpp
new file mode 100644
index 000000000..541884771
--- /dev/null
+++ b/test/algorithms/azimuth.cpp
@@ -0,0 +1,88 @@
+// Boost.Geometry
+// Unit Test
+
+// Copyright (c) 2021, Oracle and/or its affiliates.
+
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
+// Licensed under the Boost Software License version 1.0.
+// http://www.boost.org/users/license.html
+
+
+#include <geometry_test_common.hpp>
+
+#include <boost/geometry/geometries/geometries.hpp>
+
+#include <boost/geometry/algorithms/azimuth.hpp>
+
+
+template <typename P, typename S>
+auto call_azimuth(P const& p1, P const& p2, S const& s)
+{
+    return bg::azimuth(p1, p2, s);
+}
+
+template <typename P>
+auto call_azimuth(P const& p1, P const& p2, bg::default_strategy const&)
+{
+    return bg::azimuth(p1, p2);
+}
+
+template <typename P, typename S = bg::default_strategy>
+void test_one(P const& p1, P const& p2, double expected, S const& s = S())
+{
+    double const result = call_azimuth(p1, p2, s) * bg::math::r2d<double>();
+    BOOST_CHECK_CLOSE(result, expected, 0.0001);
+}
+
+template <typename P>
+void test_car()
+{
+    test_one(P(0, 0), P(0, 0), 0);
+    test_one(P(0, 0), P(1, 1), 45);
+    test_one(P(0, 0), P(100, 1), 89.427061302316517);
+    test_one(P(0, 0), P(-1, 1), -45);
+    test_one(P(0, 0), P(-100, 1), -89.427061302316517);
+}
+
+template <typename P>
+void test_sph()
+{
+    test_one(P(0, 0), P(0, 0), 0);
+    test_one(P(0, 0), P(1, 1), 44.995636455344851);
+    test_one(P(0, 0), P(100, 1), 88.984576593576293);
+    test_one(P(0, 0), P(-1, 1), -44.995636455344851);
+    test_one(P(0, 0), P(-100, 1), -88.984576593576293);
+}
+
+template <typename P>
+void test_geo()
+{
+    test_one(P(0, 0), P(0, 0), 0);
+    test_one(P(0, 0), P(1, 1), 45.187718848049521);
+    test_one(P(0, 0), P(100, 1), 88.986933066023497);
+    test_one(P(0, 0), P(-1, 1), -45.187718848049521);
+    test_one(P(0, 0), P(-100, 1), -88.986933066023497);
+}
+
+template <typename P>
+void test_geo_v()
+{
+    bg::strategies::azimuth::geographic<bg::strategy::vincenty> s;
+
+    test_one(P(0, 0), P(0, 0), 0, s);
+    test_one(P(0, 0), P(1, 1), 45.188040229339755, s);
+    test_one(P(0, 0), P(100, 1), 88.986914518230208, s);
+    test_one(P(0, 0), P(-1, 1), -45.188040229339755, s);
+    test_one(P(0, 0), P(-100, 1), -88.986914518230208, s);
+}
+
+int test_main(int, char* [])
+{
+    test_car< bg::model::point<double, 2, bg::cs::cartesian> >();
+    test_sph< bg::model::point<double, 2, bg::cs::spherical_equatorial<bg::degree> > >();
+    test_geo< bg::model::point<double, 2, bg::cs::geographic<bg::degree> > >();
+    test_geo_v< bg::model::point<double, 2, bg::cs::geographic<bg::degree> > >();
+    
+    return 0;
+}