diff options
Diffstat (limited to 'boost/geometry/formulas/karney_direct.hpp')
-rw-r--r-- | boost/geometry/formulas/karney_direct.hpp | 49 |
1 files changed, 19 insertions, 30 deletions
diff --git a/boost/geometry/formulas/karney_direct.hpp b/boost/geometry/formulas/karney_direct.hpp index 463697e4bc..bfe35ddc49 100644 --- a/boost/geometry/formulas/karney_direct.hpp +++ b/boost/geometry/formulas/karney_direct.hpp @@ -1,12 +1,14 @@ // Boost.Geometry // Copyright (c) 2018 Adeel Ahmad, Islamabad, Pakistan. +// Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland. // 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, @@ -30,8 +32,6 @@ #define BOOST_GEOMETRY_FORMULAS_KARNEY_DIRECT_HPP -#include <boost/array.hpp> - #include <boost/math/constants/constants.hpp> #include <boost/math/special_functions/hypot.hpp> @@ -82,10 +82,10 @@ public: { result_type result; - CT lon1 = lo1; - CT const lat1 = la1; + CT lon1 = lo1 * math::r2d<CT>(); + CT const lat1 = la1 * math::r2d<CT>(); - Azi azi12 = azimuth12; + Azi azi12 = azimuth12 * math::r2d<CT>(); math::normalize_azimuth<degree, Azi>(azi12); CT const c0 = 0; @@ -151,10 +151,9 @@ public: // Index zero element of coeffs_C1p is unused. se::coeffs_C1p<SeriesOrder, CT> 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 +168,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<CT>(); } if (BOOST_GEOMETRY_CONDITION(CalcCoordinates)) @@ -182,9 +178,6 @@ public: result.lat2 = atan2(sin_beta2, one_minus_f * cos_beta2); - // Convert the coordinate to radians. - result.lat2 /= math::d2r<CT>(); - // Find the longitude at the second point. CT const sin_omega2 = sin_alpha0 * sin_sigma2; CT const cos_omega2 = cos_sigma2; @@ -201,15 +194,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<CT>(); + // Convert to degrees to get the longitudinal difference. + CT lon12 = lam12 * math::r2d<CT>(); // Add the longitude at first point to the longitudinal // difference and normalize the result. @@ -224,6 +213,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<degree, CT>(result.lon2); + + result.lon2 *= math::d2r<CT>(); } if (BOOST_GEOMETRY_CONDITION(CalcQuantities)) @@ -252,12 +243,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; |