[formulas][test] Add Karney's inverse method in inverse test cases

The compilation is successful with gcc version (7.2.0),
but not with version (5.4.1). The accepted tolerance
is set to (0.0000001). Currently, all tests are not
passing, which indicates an error in the calculation.

Additionally, some changes have been made in
karney_inverse.hpp

include/boost/geometry/formulas/karney_inverse.hpp

@@ -98,7 +98,7 @@ class karney_inverse
         CT const tol_bisection = tol0 * tol2;
 
         CT const etol2 = c0_1 * tol2 /
-            sqrt(std::max(c0_001, std::abs(f)) * std::min(c1, c1 - f / c2) / c2);
+            sqrt(std::max(CT(0.001), std::abs(f)) * std::min(CT(1), CT(1) - f / CT(2)) / c2);
 
         CT tiny = std::sqrt(std::numeric_limits<CT>::min());
 
@@ -226,8 +226,8 @@ class karney_inverse
             CT sin_sigma2 = sin_beta2;
             CT cos_sigma2 = cos_alpha2 * cos_beta2;
 
-            CT sigma12 = std::atan2(std::max(c0, cos_sigma1 * sin_sigma2 - sin_sigma1 * cos_sigma2),
-                                                 cos_sigma1 * cos_sigma2 + sin_sigma1 * sin_sigma2);
+            CT sigma12 = std::atan2(std::max(CT(0), cos_sigma1 * sin_sigma2 - sin_sigma1 * cos_sigma2),
+                                                    cos_sigma1 * cos_sigma2 + sin_sigma1 * sin_sigma2);
 
             CT dummy;
             meridian_length(n, ep2, sigma12, sin_sigma1, cos_sigma1, dn1,
@@ -670,8 +670,8 @@ class karney_inverse
                 // Strip near cut.
                 if (f >= c0)
                 {
-                    sin_alpha1 = std::min(c1, -CT(x));
-                    cos_alpha1 = - std::sqrt(c1 - math::sqr(sin_alpha1));
+                    sin_alpha1 = std::min(CT(1), -CT(x));
+                    cos_alpha1 = - std::sqrt(CT(1) - math::sqr(sin_alpha1));
                 }
                 else
                 {
@@ -780,7 +780,6 @@ class karney_inverse
             // For y small, positive root is k = abs(y)/sqrt(1-x^2).
             k = c0;
         }
-
         return k;
     }
 
@@ -844,11 +843,11 @@ class karney_inverse
         math::normalize_values<CT>(sin_sigma2, cos_sigma2);
 
         // sig12 = sig2 - sig1, limit to [0, pi].
-        sigma12 = atan2(std::max(c0, cos_sigma1 * sin_sigma2 - sin_sigma1 * cos_sigma2),
-                                     cos_sigma1 * cos_sigma2 + sin_sigma1 * sin_sigma2);
+        sigma12 = atan2(std::max(CT(0), cos_sigma1 * sin_sigma2 - sin_sigma1 * cos_sigma2),
+                                        cos_sigma1 * cos_sigma2 + sin_sigma1 * sin_sigma2);
 
         // omg12 = omg2 - omg1, limit to [0, pi].
-        sin_omega12 = std::max(c0, cos_omega1 * sin_omega2 - sin_omega1 * cos_omega2);
+        sin_omega12 = std::max(CT(0), cos_omega1 * sin_omega2 - sin_omega1 * cos_omega2);
         cos_omega12 = cos_omega1 * cos_omega2 + sin_omega1 * sin_omega2;
 
         // eta = omg12 - lam120.
@@ -900,7 +899,6 @@ class karney_inverse
                 diff_lam12 *= one_minus_f / (cos_alpha2 * cos_beta2);
             }
         }
-
         return lam12;
     }
 

test/formulas/inverse.cpp

@@ -18,6 +18,7 @@
 #include <boost/geometry/formulas/vincenty_inverse.hpp>
 #include <boost/geometry/formulas/thomas_inverse.hpp>
 #include <boost/geometry/formulas/andoyer_inverse.hpp>
+#include <boost/geometry/formulas/karney_inverse.hpp>
 
 #include <boost/geometry/srs/spheroid.hpp>
 
@@ -53,10 +54,15 @@ void test_all(expected_results const& results)
     double lon2r = results.p2.lon * d2r;
     double lat2r = results.p2.lat * d2r;
 
+    double lon1d = results.p1.lon;
+    double lat1d = results.p1.lat;
+    double lon2d = results.p2.lon;
+    double lat2d = results.p2.lat;
+
     // WGS84
     bg::srs::spheroid<double> spheroid(6378137.0, 6356752.3142451793);
 
-    bg::formula::result_inverse<double> result_v, result_t, result_a;
+    bg::formula::result_inverse<double> result_v, result_t, result_a, result_k;
 
     typedef bg::formula::vincenty_inverse<double, true, true, true, true, true> vi_t;
     result_v = vi_t::apply(lon1r, lat1r, lon2r, lat2r, spheroid);
@@ -75,6 +81,10 @@ void test_all(expected_results const& results)
     result_a.azimuth *= r2d;
     result_a.reverse_azimuth *= r2d;
     check_inverse("andoyer", results, result_a, results.andoyer, results.reference, 0.001);
+
+    typedef bg::formula::karney_inverse<double, true, true, true, true, true, 8> ka_t;
+    result_k = ka_t::apply(lon1d, lat1d, lon2d, lat2d, spheroid);
+    check_inverse("karney", results, result_k, results.vincenty, results.reference, 0.0000001);
 }
 
 int test_main(int, char*[])