From 4520fbc373681745fac040c356e4f781e8a3b65f Mon Sep 17 00:00:00 2001 From: Adam Wulkiewicz Date: Sun, 25 Jan 2015 19:49:44 +0100 Subject: [PATCH] [algorithms] Improve the robustness of thomas_inverse formula. For lat equal to pi/2 or -pi/2 avoid calculating "infinite" tangent result. --- .../algorithms/detail/thomas_inverse.hpp | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/include/boost/geometry/algorithms/detail/thomas_inverse.hpp b/include/boost/geometry/algorithms/detail/thomas_inverse.hpp index d427e7697..f798dd204 100644 --- a/include/boost/geometry/algorithms/detail/thomas_inverse.hpp +++ b/include/boost/geometry/algorithms/detail/thomas_inverse.hpp @@ -59,10 +59,18 @@ public: CT const one_minus_f = CT(1) - f; - CT const tan_theta1 = one_minus_f * tan(lat1); - CT const tan_theta2 = one_minus_f * tan(lat2); - CT const theta1 = atan(tan_theta1); - CT const theta2 = atan(tan_theta2); +// CT const tan_theta1 = one_minus_f * tan(lat1); +// CT const tan_theta2 = one_minus_f * tan(lat2); +// CT const theta1 = atan(tan_theta1); +// CT const theta2 = atan(tan_theta2); + + CT const pi_half = math::pi() / CT(2); + CT const theta1 = math::equals(lat1, pi_half) ? lat1 : + math::equals(lat1, -pi_half) ? lat1 : + atan(one_minus_f * tan(lat1)); + CT const theta2 = math::equals(lat2, pi_half) ? lat2 : + math::equals(lat2, -pi_half) ? lat2 : + atan(one_minus_f * tan(lat2)); CT const theta_m = (theta1 + theta2) / CT(2); CT const d_theta_m = (theta2 - theta1) / CT(2);