boost/geometry
2
0
Fork 0
mirror of https://github.com/boostorg/geometry.git synced 2026-01-27 19:02:12 +00:00
Code Issues Projects Releases Wiki Activity

Merge pull request #1087 from vissarion/fix/kanrey_direct_units2

Browse Source 
Make all direct formulas consistent
This commit is contained in:
Vissarion Fisikopoulos
2022-11-16 15:04:55 +02:00
committed by GitHub
parent 75742b0d24 9a4335a57d
commit 1663d2bb0b
2 changed files with 29 additions and 52 deletions
Download Patch File Download Diff File Expand all files Collapse all files

46
include/boost/geometry/formulas/karney_direct.hpp
View File

@@ -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>();
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 +152,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 +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<CT>();
}
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<CT>();
// 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<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 +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<degree, CT>(result.lon2);
result.lon2 *= math::d2r<CT>();
}
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;

35
test/formulas/direct.cpp
View File

@@ -48,7 +48,7 @@ inline expected_results symmetric_wrt_origin(expected_results r)
}
template <typename Result>
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 <typename Result>
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<double>();
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>();
double const r2d = bg::math::r2d<double>();
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<double> spheroid(6378137.0, 6356752.3142451793);
bg::srs::sphere<double> const sphere;
@@ -87,23 +86,14 @@ void test_all(expected_results const& results)
typedef bg::formula::vincenty_direct<double, true, true, true, true> 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<double, true, true, true, true, true> 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<double, false, true, true, true, true> 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<double, true, true, true, true, 4> series;
@@ -115,13 +105,10 @@ void test_all(expected_results const& results)
*/
result = bg::formula::spherical_direct<true, true>(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<double, true, true, true, true, 2> 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>();
double lat1r = results.p1.lat * bg::math::d2r<double>();
double distance = results.distance;
double azi12d = results.azimuth12;
double azi12r = results.azimuth12 * bg::math::d2r<double>();
// WGS84
bg::srs::spheroid<double> spheroid(6378137.0, 6356752.3142451793);
@@ -151,7 +138,7 @@ void test_karney_antipodal(expected_results_antipodal const& results)
bg::formula::result_direct<double> result;
typedef bg::formula::karney_direct<double, true, true, true, true, 8> 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);
}