diff --git a/include/boost/geometry/formulas/karney_direct.hpp b/include/boost/geometry/formulas/karney_direct.hpp index 463697e4b..20b0120d7 100644 --- a/include/boost/geometry/formulas/karney_direct.hpp +++ b/include/boost/geometry/formulas/karney_direct.hpp @@ -5,8 +5,9 @@ // Contributed and/or modified by Adeel Ahmad, // as part of Google Summer of Code 2018 program. -// This file was modified by Oracle on 2018-2020. -// Modifications copyright (c) 2018-2020 Oracle and/or its affiliates. +// This file was modified by Oracle on 2018-2022. +// Modifications copyright (c) 2018-2022 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, @@ -82,10 +83,10 @@ public: { result_type result; - CT lon1 = lo1; - CT const lat1 = la1; + CT lon1 = lo1 * math::r2d(); + CT const lat1 = la1 * math::r2d(); - Azi azi12 = azimuth12; + Azi azi12 = azimuth12 * math::r2d(); math::normalize_azimuth(azi12); CT const c0 = 0; @@ -151,10 +152,9 @@ public: // Index zero element of coeffs_C1p is unused. se::coeffs_C1p const coeffs_C1p(epsilon); - CT const B12 = - se::sin_cos_series - (sin_tau1 * cos_tau12 + cos_tau1 * sin_tau12, - cos_tau1 * cos_tau12 - sin_tau1 * sin_tau12, - coeffs_C1p); + CT const B12 = - se::sin_cos_series(sin_tau1 * cos_tau12 + cos_tau1 * sin_tau12, + cos_tau1 * cos_tau12 - sin_tau1 * sin_tau12, + coeffs_C1p); CT const sigma12 = tau12 - (B12 - B11); CT const sin_sigma12 = sin(sigma12); @@ -169,9 +169,6 @@ public: CT const cos_alpha2 = cos_alpha0 * cos_sigma2; result.reverse_azimuth = atan2(sin_alpha2, cos_alpha2); - - // Convert the angle to radians. - result.reverse_azimuth /= math::d2r(); } if (BOOST_GEOMETRY_CONDITION(CalcCoordinates)) @@ -182,9 +179,6 @@ public: result.lat2 = atan2(sin_beta2, one_minus_f * cos_beta2); - // Convert the coordinate to radians. - result.lat2 /= math::d2r(); - // Find the longitude at the second point. CT const sin_omega2 = sin_alpha0 * sin_sigma2; CT const cos_omega2 = cos_sigma2; @@ -201,15 +195,11 @@ public: CT const B31 = se::sin_cos_series(sin_sigma1, cos_sigma1, coeffs_C3); - CT const lam12 = omega12 + A3c * - (sigma12 + (se::sin_cos_series - (sin_sigma2, - cos_sigma2, - coeffs_C3) - B31)); + CT const sin_cos_res = se::sin_cos_series(sin_sigma2, cos_sigma2, coeffs_C3); + CT const lam12 = omega12 + A3c * (sigma12 + (sin_cos_res - B31)); - // Convert to radians to get the - // longitudinal difference. - CT lon12 = lam12 / math::d2r(); + // Convert to degrees to get the longitudinal difference. + CT lon12 = lam12 * math::r2d(); // Add the longitude at first point to the longitudinal // difference and normalize the result. @@ -224,6 +214,8 @@ public: // otherwise differential quantities are calculated incorrectly. // But here it's ok since result.lon2 is not used after this point. math::normalize_longitude(result.lon2); + + result.lon2 *= math::d2r(); } if (BOOST_GEOMETRY_CONDITION(CalcQuantities)) @@ -252,12 +244,10 @@ public: cos_sigma1 * cos_sigma2 * J12); // Find the geodesic scale. - CT const t = k2 * (sin_sigma2 - sin_sigma1) * - (sin_sigma2 + sin_sigma1) / (dn1 + dn2); + CT const t = k2 * (sin_sigma2 - sin_sigma1) * (sin_sigma2 + sin_sigma1) / (dn1 + dn2); - result.geodesic_scale = cos_sigma12 + - (t * sin_sigma2 - cos_sigma2 * J12) * - sin_sigma1 / dn1; + result.geodesic_scale = cos_sigma12 + (t * sin_sigma2 - cos_sigma2 * J12) * + sin_sigma1 / dn1; } return result; diff --git a/test/formulas/direct.cpp b/test/formulas/direct.cpp index 97fb170f4..40b9d5291 100644 --- a/test/formulas/direct.cpp +++ b/test/formulas/direct.cpp @@ -48,7 +48,7 @@ inline expected_results symmetric_wrt_origin(expected_results r) } template -void check_direct(Result const& result, expected_result const& expected, expected_result const& reference, +void check_direct(Result& result, expected_result const& expected, expected_result const& reference, double reference_error, bool check_reference_only = false) { check_direct_sph(result, expected, reference, reference_error, check_reference_only); @@ -57,9 +57,13 @@ void check_direct(Result const& result, expected_result const& expected, expecte } template -void check_direct_sph(Result const& result, expected_result const& expected, expected_result const& reference, +void check_direct_sph(Result& result, expected_result const& expected, expected_result const& reference, double reference_error, bool check_reference_only = false) { + double const r2d = bg::math::r2d(); + result.lon2 *= r2d; + result.lat2 *= r2d; + result.reverse_azimuth *= r2d; check_one(result.lon2, expected.lon2, reference.lon2, reference_error, true, check_reference_only); check_one(result.lat2, expected.lat2, reference.lat2, reference_error, true, check_reference_only); check_one(result.reverse_azimuth, expected.reverse_azimuth, reference.reverse_azimuth, reference_error, true, check_reference_only); @@ -68,17 +72,12 @@ void check_direct_sph(Result const& result, expected_result const& expected, exp void test_all(expected_results const& results) { double const d2r = bg::math::d2r(); - double const r2d = bg::math::r2d(); double lon1r = results.p1.lon * d2r; double lat1r = results.p1.lat * d2r; double distance = results.distance; double azi12r = results.azimuth12 * d2r; - double lon1d = results.p1.lon; - double lat1d = results.p1.lat; - double azi12d = results.azimuth12; - // WGS84 bg::srs::spheroid spheroid(6378137.0, 6356752.3142451793); bg::srs::sphere const sphere; @@ -87,23 +86,14 @@ void test_all(expected_results const& results) typedef bg::formula::vincenty_direct vi_t; result = vi_t::apply(lon1r, lat1r, distance, azi12r, spheroid); - result.lon2 *= r2d; - result.lat2 *= r2d; - result.reverse_azimuth *= r2d; check_direct(result, results.vincenty, results.karney, 0.00000001); typedef bg::formula::thomas_direct th_t; result = th_t::apply(lon1r, lat1r, distance, azi12r, spheroid); - result.lon2 *= r2d; - result.lat2 *= r2d; - result.reverse_azimuth *= r2d; check_direct(result, results.thomas, results.karney, 0.0000001); typedef bg::formula::thomas_direct th_t1st; result = th_t1st::apply(lon1r, lat1r, distance, azi12r, spheroid); - result.lon2 *= r2d; - result.lat2 *= r2d; - result.reverse_azimuth *= r2d; check_direct(result, results.thomas1st, results.karney, 0.0000001); /* typedef bg::formula::series_expansion_direct series; @@ -115,13 +105,10 @@ void test_all(expected_results const& results) */ result = bg::formula::spherical_direct(lon1r, lat1r, distance, azi12r, sphere); - result.lon2 *= r2d; - result.lat2 *= r2d; - result.reverse_azimuth *= r2d; check_direct_sph(result, results.spherical, results.karney, 0.1); typedef bg::formula::karney_direct ka_t; - result = ka_t::apply(lon1d, lat1d, distance, azi12d, spheroid); + result = ka_t::apply(lon1r, lat1r, distance, azi12r, spheroid); check_direct(result, results.karney, results.karney, 0.0000001); #ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB @@ -140,10 +127,10 @@ void test_all(expected_results const& results) void test_karney_antipodal(expected_results_antipodal const& results) { - double lon1d = results.p1.lon; - double lat1d = results.p1.lat; + double lon1r = results.p1.lon * bg::math::d2r(); + double lat1r = results.p1.lat * bg::math::d2r(); double distance = results.distance; - double azi12d = results.azimuth12; + double azi12r = results.azimuth12 * bg::math::d2r(); // WGS84 bg::srs::spheroid spheroid(6378137.0, 6356752.3142451793); @@ -151,7 +138,7 @@ void test_karney_antipodal(expected_results_antipodal const& results) bg::formula::result_direct result; typedef bg::formula::karney_direct ka_t; - result = ka_t::apply(lon1d, lat1d, distance, azi12d, spheroid); + result = ka_t::apply(lon1r, lat1r, distance, azi12r, spheroid); check_direct(result, results.karney, results.karney, 0.0000001, true); }