From e314795dbb0a0d1218dc4cbdf3449772ad5eeb2e Mon Sep 17 00:00:00 2001
From: Adam Wulkiewicz <adam.wulkiewicz@gmail.com>
Date: Fri, 15 Mar 2019 17:28:48 +0100
Subject: [PATCH] [formulas] Add asserts for validity of azimuth passed to
 thomas_direct and differential_quantities.

---
 include/boost/geometry/formulas/differential_quantities.hpp | 4 ++++
 include/boost/geometry/formulas/thomas_direct.hpp           | 3 +++
 2 files changed, 7 insertions(+)

diff --git a/include/boost/geometry/formulas/differential_quantities.hpp b/include/boost/geometry/formulas/differential_quantities.hpp
index 6d4c22328..6e1658abe 100644
--- a/include/boost/geometry/formulas/differential_quantities.hpp
+++ b/include/boost/geometry/formulas/differential_quantities.hpp
@@ -11,6 +11,7 @@
 #ifndef BOOST_GEOMETRY_FORMULAS_INVERSE_DIFFERENTIAL_QUANTITIES_HPP
 #define BOOST_GEOMETRY_FORMULAS_INVERSE_DIFFERENTIAL_QUANTITIES_HPP
 
+#include <boost/geometry/core/assert.hpp>
 
 #include <boost/geometry/util/condition.hpp>
 #include <boost/geometry/util/math.hpp>
@@ -73,6 +74,9 @@ public:
             CT const sig_12 = dlon / one_minus_f;
             if (BOOST_GEOMETRY_CONDITION(EnableReducedLength))
             {
+                CT const pi = math::pi<CT>();
+                BOOST_GEOMETRY_ASSERT(-pi <= azimuth && azimuth <= pi);
+
                 int azi_sign = math::sign(azimuth) >= 0 ? 1 : -1; // for antipodal
                 CT m12 = azi_sign * sin(sig_12) * b;
                 reduced_length = m12;
diff --git a/include/boost/geometry/formulas/thomas_direct.hpp b/include/boost/geometry/formulas/thomas_direct.hpp
index 7b61616af..2218ea995 100644
--- a/include/boost/geometry/formulas/thomas_direct.hpp
+++ b/include/boost/geometry/formulas/thomas_direct.hpp
@@ -15,6 +15,7 @@
 
 #include <boost/math/constants/constants.hpp>
 
+#include <boost/geometry/core/assert.hpp>
 #include <boost/geometry/core/radius.hpp>
 
 #include <boost/geometry/util/condition.hpp>
@@ -80,6 +81,8 @@ public:
         CT const pi = math::pi<CT>();
         CT const pi_half = pi / c2;
 
+        BOOST_GEOMETRY_ASSERT(-pi <= azimuth12 && azimuth12 <= pi);
+
         // keep azimuth small - experiments show low accuracy
         // if the azimuth is closer to (+-)180 deg.
         CT azi12_alt = azimuth12;