1 // Boost.Geometry
2 // Unit Test
3
4 // Copyright (c) 2016-2018 Oracle and/or its affiliates.
5
6 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
7 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
8 // Contributed and/or modified by Adeel Ahmad, as part of Google Summer of Code 2018 program
9
10 // Use, modification and distribution is subject to the Boost Software License,
11 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
12 // http://www.boost.org/LICENSE_1_0.txt)
13
14 #include "test_formula.hpp"
15 #include "direct_cases.hpp"
16 #include "direct_cases_antipodal.hpp"
17
18 #include <boost/geometry/formulas/vincenty_direct.hpp>
19 #include <boost/geometry/formulas/thomas_direct.hpp>
20 #include <boost/geometry/formulas/karney_direct.hpp>
21 //#include <boost/geometry/formulas/series_expansion_direct.hpp>
22 #include <boost/geometry/formulas/spherical.hpp>
23
24 #include <boost/geometry/srs/srs.hpp>
25
26 #ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
27 #include <GeographicLib/Geodesic.hpp>
28 #include <GeographicLib/Constants.hpp>
29 #endif // BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
30
symmetrize_wrt_origin(expected_result & r)31 inline void symmetrize_wrt_origin(expected_result & r)
32 {
33 r.lon2 = -r.lon2;
34 r.lat2 = -r.lat2;
35 r.reduced_length = -r.reduced_length;
36 }
37
symmetric_wrt_origin(expected_results r)38 inline expected_results symmetric_wrt_origin(expected_results r)
39 {
40 r.distance = -r.distance;
41 symmetrize_wrt_origin(r.karney);
42 symmetrize_wrt_origin(r.series);
43 symmetrize_wrt_origin(r.spherical);
44 symmetrize_wrt_origin(r.thomas);
45 symmetrize_wrt_origin(r.thomas1st);
46 symmetrize_wrt_origin(r.vincenty);
47 return r;
48 }
49
50 template <typename Result>
check_direct(Result const & result,expected_result const & expected,expected_result const & reference,double reference_error,bool check_reference_only=false)51 void check_direct(Result const& result, expected_result const& expected, expected_result const& reference,
52 double reference_error, bool check_reference_only = false)
53 {
54 check_direct_sph(result, expected, reference, reference_error, check_reference_only);
55 check_one(result.reduced_length, expected.reduced_length, reference.reduced_length, reference_error);
56 check_one(result.geodesic_scale, expected.geodesic_scale, reference.geodesic_scale, reference_error);
57 }
58
59 template <typename Result>
check_direct_sph(Result const & result,expected_result const & expected,expected_result const & reference,double reference_error,bool check_reference_only=false)60 void check_direct_sph(Result const& result, expected_result const& expected, expected_result const& reference,
61 double reference_error, bool check_reference_only = false)
62 {
63 check_one(result.lon2, expected.lon2, reference.lon2, reference_error, true, check_reference_only);
64 check_one(result.lat2, expected.lat2, reference.lat2, reference_error, true, check_reference_only);
65 check_one(result.reverse_azimuth, expected.reverse_azimuth, reference.reverse_azimuth, reference_error, true, check_reference_only);
66 }
67
test_all(expected_results const & results)68 void test_all(expected_results const& results)
69 {
70 double const d2r = bg::math::d2r<double>();
71 double const r2d = bg::math::r2d<double>();
72
73 double lon1r = results.p1.lon * d2r;
74 double lat1r = results.p1.lat * d2r;
75 double distance = results.distance;
76 double azi12r = results.azimuth12 * d2r;
77
78 double lon1d = results.p1.lon;
79 double lat1d = results.p1.lat;
80 double azi12d = results.azimuth12;
81
82 // WGS84
83 bg::srs::spheroid<double> spheroid(6378137.0, 6356752.3142451793);
84 bg::srs::sphere<double> const sphere;
85
86 bg::formula::result_direct<double> result;
87
88 typedef bg::formula::vincenty_direct<double, true, true, true, true> vi_t;
89 result = vi_t::apply(lon1r, lat1r, distance, azi12r, spheroid);
90 result.lon2 *= r2d;
91 result.lat2 *= r2d;
92 result.reverse_azimuth *= r2d;
93 check_direct(result, results.vincenty, results.karney, 0.00000001);
94
95 typedef bg::formula::thomas_direct<double, true, true, true, true, true> th_t;
96 result = th_t::apply(lon1r, lat1r, distance, azi12r, spheroid);
97 result.lon2 *= r2d;
98 result.lat2 *= r2d;
99 result.reverse_azimuth *= r2d;
100 check_direct(result, results.thomas, results.karney, 0.0000001);
101
102 typedef bg::formula::thomas_direct<double, false, true, true, true, true> th_t1st;
103 result = th_t1st::apply(lon1r, lat1r, distance, azi12r, spheroid);
104 result.lon2 *= r2d;
105 result.lat2 *= r2d;
106 result.reverse_azimuth *= r2d;
107 check_direct(result, results.thomas1st, results.karney, 0.0000001);
108 /*
109 typedef bg::formula::series_expansion_direct<double, true, true, true, true, 4> series;
110 result = series::apply(lon1r, lat1r, distance, azi12r, spheroid);
111 result.lon2 *= r2d;
112 result.lat2 *= r2d;
113 result.reverse_azimuth *= r2d;
114 check_direct(result, results.series, results.karney, 0.0000001);
115 */
116 result = bg::formula::spherical_direct<true, true>(lon1r, lat1r, distance,
117 azi12r, sphere);
118 result.lon2 *= r2d;
119 result.lat2 *= r2d;
120 result.reverse_azimuth *= r2d;
121 check_direct_sph(result, results.spherical, results.karney, 0.1);
122
123 typedef bg::formula::karney_direct<double, true, true, true, true, 2> ka_t;
124 result = ka_t::apply(lon1d, lat1d, distance, azi12d, spheroid);
125 check_direct(result, results.karney, results.karney, 0.0000001);
126
127 #ifdef BOOST_GEOEMTRY_TEST_WITH_GEOGRAPHICLIB
128 {
129 using namespace GeographicLib;
130 Geodesic geod(Constants::WGS84_a(), Constants::WGS84_f());
131 double foo = 0;
132 geod.Direct(lat1d, lon1d, azi12d, distance,
133 result.lat2, result.lon2, result.reverse_azimuth,
134 result.reduced_length, result.geodesic_scale, foo);
135 boost::ignore_unused(foo);
136 check_direct(result, results.karney, results.karney, 0.0000001);
137 }
138 #endif
139 }
140
test_karney_antipodal(expected_results_antipodal const & results)141 void test_karney_antipodal(expected_results_antipodal const& results)
142 {
143 double lon1d = results.p1.lon;
144 double lat1d = results.p1.lat;
145 double distance = results.distance;
146 double azi12d = results.azimuth12;
147
148 // WGS84
149 bg::srs::spheroid<double> spheroid(6378137.0, 6356752.3142451793);
150
151 bg::formula::result_direct<double> result;
152
153 typedef bg::formula::karney_direct<double, true, true, true, true, 8> ka_t;
154 result = ka_t::apply(lon1d, lat1d, distance, azi12d, spheroid);
155 check_direct(result, results.karney, results.karney, 0.0000001, true);
156 }
157
test_main(int,char * [])158 int test_main(int, char*[])
159 {
160 for (size_t i = 0; i < expected_size; ++i)
161 {
162 test_all(expected[i]);
163
164 if (expected[i].p1.lon == 0 && expected[i].p1.lat == 0)
165 {
166 test_all(symmetric_wrt_origin(expected[i]));
167 }
168 }
169
170 for (size_t i = 0; i < expected_size_antipodal; ++i)
171 {
172 test_karney_antipodal(expected_antipodal[i]);
173 }
174
175 return 0;
176 }
177