From 07c112ce085806507adbc93859bfa4b7f155083a Mon Sep 17 00:00:00 2001 From: Marcel Jacobse Date: Mon, 27 Apr 2026 20:56:44 +0200 Subject: [PATCH] Improve short-distance accuracy of Andoyer inverse Use haversine formula instead of law of cosines to avoid numerical loss of precision for close points as suggested in #1217. This had the side-effect of returning non-zero azimuths for smaller angles than before, which resulted in quite inaccurate non-zero reduced length and geodesic scale values too. To fix this, division by the cosine of the latitudes was changed to be done implicitly within atan2, which together with use of the haversine formula improves accuracy for all result values. --- .../geometry/formulas/andoyer_inverse.hpp | 42 ++++++++++--------- 1 file changed, 23 insertions(+), 19 deletions(-) diff --git a/include/boost/geometry/formulas/andoyer_inverse.hpp b/include/boost/geometry/formulas/andoyer_inverse.hpp index 708a574399..427db3c77e 100644 --- a/include/boost/geometry/formulas/andoyer_inverse.hpp +++ b/include/boost/geometry/formulas/andoyer_inverse.hpp @@ -73,28 +73,35 @@ class andoyer_inverse CT const c0 = CT(0); CT const c1 = CT(1); + CT const c2 = CT(2); CT const pi = math::pi(); CT const f = formula::flattening(spheroid); CT const dlon = lon2 - lon1; + CT const dlat = lat2 - lat1; CT const sin_dlon = sin(dlon); CT const cos_dlon = cos(dlon); + CT const hav_dlon = math::hav(dlon); CT const sin_lat1 = sin(lat1); CT const cos_lat1 = cos(lat1); CT const sin_lat2 = sin(lat2); CT const cos_lat2 = cos(lat2); + CT const hav_dlat = math::hav(dlat); - // H,G,T = infinity if cos_d = 1 or cos_d = -1 + // using the haversine formula instead of the cosine law for better accuracy around short distances + // H,G,T = infinity if hav_d = 0 or hav_d = 1 // lat1 == +-90 && lat2 == +-90 // lat1 == lat2 && lon1 == lon2 - CT cos_d = sin_lat1*sin_lat2 + cos_lat1*cos_lat2*cos_dlon; - // on some platforms cos_d may be outside valid range - if (cos_d < -c1) - cos_d = -c1; - else if (cos_d > c1) - cos_d = c1; - - CT const d = acos(cos_d); // [0, pi] + CT hav_d = hav_dlat + cos_lat1*cos_lat2*hav_dlon; + // on some platforms hav_d may be outside valid range + if (hav_d < c0) + hav_d = c0; + else if (hav_d > c1) + hav_d = c1; + + CT const sin_d_half = sqrt(hav_d); + CT const d_half = asin(sin_d_half); + CT const d = c2 * d_half; // [0, pi] CT const sin_d = sin(d); // [-1, 1] if BOOST_GEOMETRY_CONSTEXPR (EnableDistance) @@ -103,8 +110,8 @@ class andoyer_inverse CT const L = math::sqr(sin_lat1+sin_lat2); CT const three_sin_d = CT(3) * sin_d; - CT const one_minus_cos_d = c1 - cos_d; - CT const one_plus_cos_d = c1 + cos_d; + CT const one_minus_cos_d = c2 * hav_d; + CT const one_plus_cos_d = c2 - one_minus_cos_d; // cos_d = 1 means that the points are very close // cos_d = -1 means that the points are antipodal @@ -141,7 +148,7 @@ class andoyer_inverse // correctly and consistently across all formulas. // points very close - if (cos_d >= c0) + if (hav_d <= CT(0.5)) { result.azimuth = c0; result.reverse_azimuth = c0; @@ -164,7 +171,8 @@ class andoyer_inverse } else { - CT const c2 = CT(2); + CT const M = cos_lat1 * sin_lat2; + CT const N = sin_lat1 * cos_lat2; CT A = c0; CT U = c0; @@ -177,9 +185,7 @@ class andoyer_inverse } else { - CT const tan_lat2 = sin_lat2/cos_lat2; - CT const M = cos_lat1*tan_lat2-sin_lat1*cos_dlon; - A = atan2(sin_dlon, M); + A = atan2(sin_dlon*cos_lat2, M - N*cos_dlon); CT const sin_2A = sin(c2*A); U = (f/ c2)*math::sqr(cos_lat1)*sin_2A; } @@ -195,9 +201,7 @@ class andoyer_inverse } else { - CT const tan_lat1 = sin_lat1/cos_lat1; - CT const N = cos_lat2*tan_lat1-sin_lat2*cos_dlon; - B = atan2(sin_dlon, N); + B = atan2(sin_dlon*cos_lat1, N - M*cos_dlon); CT const sin_2B = sin(c2*B); V = (f/ c2)*math::sqr(cos_lat2)*sin_2B; }