Skip to content

Improve short-distance accuracy of Andoyer inverse #1461

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
mjacobse wants to merge 6 commits into
base: develop
Choose a base branch
from
+30 −20
Open

Improve short-distance accuracy of Andoyer inverse #1461

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
07c112c
Improve short-distance accuracy of Andoyer inverse
mjacobse Apr 27, 2026
7c18419
Qualify sqrt call with math namespace
mjacobse May 5, 2026
d3a4f79
Add curly braces for if statement
mjacobse May 5, 2026
38ec1bd
Improve readability of formulas
mjacobse May 5, 2026
14b6c8e
Update test reference value
mjacobse May 5, 2026
905ea42
Add comment for haversine - cosine relation
mjacobse May 7, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
48 changes: 29 additions & 19 deletions include/boost/geometry/formulas/andoyer_inverse.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,28 +73,39 @@ class andoyer_inverse

CT const c0 = CT(0);
CT const c1 = CT(1);
CT const c2 = CT(2);
CT const pi = math::pi<CT>();
CT const f = formula::flattening<CT>(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;
Comment thread
mjacobse marked this conversation as resolved.
}

CT const sin_d_half = math::sqrt(hav_d);
CT const d_half = asin(sin_d_half);
Comment thread
mjacobse marked this conversation as resolved.
CT const d = c2 * d_half; // [0, pi]
CT const sin_d = sin(d); // [-1, 1]

if BOOST_GEOMETRY_CONSTEXPR (EnableDistance)
Expand All @@ -103,8 +114,10 @@ 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;
// follows from definition of haversine and trigonometric power reduction formula:
// hav(d) = sin^2(d/2) = (1 - cos(d))/2
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
Comment thread
mjacobse marked this conversation as resolved.

Expand Down Expand Up @@ -141,7 +154,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;
Expand All @@ -164,7 +177,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;
Expand All @@ -177,9 +191,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;
}
Expand All @@ -195,9 +207,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;
}
Expand Down
2 changes: 1 addition & 1 deletion test/algorithms/area/area_sph_geo.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -501,7 +501,7 @@ void test_spherical_geo()

bg::read_wkt(wkt, geometry_geo_d);
area = bg::area(geometry_geo_d, area_a);
BOOST_CHECK_CLOSE(area, -25.47837, 0.001);
BOOST_CHECK_CLOSE(area, -25.55885, 0.001);
area = bg::area(geometry_geo_d, area_t);
BOOST_CHECK_CLOSE(area, -25.57355, 0.001);
area = bg::area(geometry_geo_d, area_v);
Expand Down