1 // Boost.Geometry
2 // Unit Test
3
4 // Copyright (c) 2017 Oracle and/or its affiliates.
5
6 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
7
8 // Use, modification and distribution is subject to the Boost Software License,
9 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
10 // http://www.boost.org/LICENSE_1_0.txt)
11
12 #include <geometry_test_common.hpp>
13
14 #include "test_formula.hpp"
15 #include "vertex_longitude_cases.hpp"
16
17 #include <boost/geometry/formulas/vertex_latitude.hpp>
18 #include <boost/geometry/formulas/vertex_longitude.hpp>
19 #include <boost/geometry/formulas/vincenty_inverse.hpp>
20 #include <boost/geometry/formulas/thomas_inverse.hpp>
21 #include <boost/geometry/formulas/andoyer_inverse.hpp>
22
23 #include <boost/geometry/strategies/strategies.hpp>
24
25 #include <boost/geometry/util/math.hpp>
26
27 #define BOOST_GEOMETRY_TEST_DEBUG
28
29 namespace bg = boost::geometry;
30
31 template<typename CT>
test_vrt_lon_sph(CT lon1r,CT lat1r,CT lon2r,CT lat2r)32 CT test_vrt_lon_sph(CT lon1r,
33 CT lat1r,
34 CT lon2r,
35 CT lat2r)
36 {
37 CT a1 = bg::formula::spherical_azimuth(lon1r, lat1r, lon2r, lat2r);
38
39 typedef bg::model::point<CT, 2,
40 bg::cs::spherical_equatorial<bg::radian> > point;
41
42 bg::model::segment<point> segment(point(lon1r, lat1r),
43 point(lon2r, lat2r));
44 bg::model::box<point> box;
45 bg::envelope(segment, box);
46
47 CT vertex_lat;
48 CT lat_sum = lat1r + lat2r;
49 if (lat_sum > CT(0))
50 {
51 vertex_lat = bg::get_as_radian<bg::max_corner, 1>(box);
52 } else {
53 vertex_lat = bg::get_as_radian<bg::min_corner, 1>(box);
54 }
55
56 bg::strategy::azimuth::spherical<> azimuth;
57
58 return bg::formula::vertex_longitude
59 <CT, bg::spherical_equatorial_tag>::
60 apply(lon1r, lat1r,
61 lon2r, lat2r,
62 vertex_lat,
63 a1,
64 azimuth);
65 }
66
67 template
68 <
69 template <typename, bool, bool, bool, bool, bool> class FormulaPolicy,
70 typename CT
71 >
test_vrt_lon_geo(CT lon1r,CT lat1r,CT lon2r,CT lat2r)72 CT test_vrt_lon_geo(CT lon1r,
73 CT lat1r,
74 CT lon2r,
75 CT lat2r)
76 {
77 // WGS84
78 bg::srs::spheroid<CT> spheroid(6378137.0, 6356752.3142451793);
79
80 typedef FormulaPolicy<CT, false, true, false, false, false> formula;
81 CT a1 = formula::apply(lon1r, lat1r, lon2r, lat2r, spheroid).azimuth;
82
83 typedef bg::model::point<CT, 2, bg::cs::geographic<bg::radian> > geo_point;
84
85 bg::model::segment<geo_point> segment(geo_point(lon1r, lat1r),
86 geo_point(lon2r, lat2r));
87 bg::model::box<geo_point> box;
88 bg::envelope(segment, box);
89
90 CT vertex_lat;
91 CT lat_sum = lat1r + lat2r;
92 if (lat_sum > CT(0))
93 {
94 vertex_lat = bg::get_as_radian<bg::max_corner, 1>(box);
95 } else {
96 vertex_lat = bg::get_as_radian<bg::min_corner, 1>(box);
97 }
98
99 bg::strategy::azimuth::geographic<> azimuth_geographic;
100
101 return bg::formula::vertex_longitude
102 <CT, bg::geographic_tag>::apply(lon1r, lat1r,
103 lon2r, lat2r,
104 vertex_lat,
105 a1,
106 azimuth_geographic);
107
108 }
109
test_all(expected_results const & results)110 void test_all(expected_results const& results)
111 {
112 double const d2r = bg::math::d2r<double>();
113
114 double lon1r = results.p1.lon * d2r;
115 double lat1r = results.p1.lat * d2r;
116 double lon2r = results.p2.lon * d2r;
117 double lat2r = results.p2.lat * d2r;
118
119 if(lon1r > lon2r)
120 {
121 std::swap(lon1r, lon2r);
122 std::swap(lat1r, lat2r);
123 }
124
125 double res_an = test_vrt_lon_geo<bg::formula::andoyer_inverse>
126 (lon1r, lat1r, lon2r, lat2r);
127 double res_th = test_vrt_lon_geo<bg::formula::thomas_inverse>
128 (lon1r, lat1r, lon2r, lat2r);
129 double res_vi = test_vrt_lon_geo<bg::formula::vincenty_inverse>
130 (lon1r, lat1r, lon2r, lat2r);
131 double res_sh = test_vrt_lon_sph(lon1r, lat1r, lon2r, lat2r);
132
133 bg::math::normalize_longitude<bg::radian, double>(res_an);
134 bg::math::normalize_longitude<bg::radian, double>(res_th);
135 bg::math::normalize_longitude<bg::radian, double>(res_vi);
136 bg::math::normalize_longitude<bg::radian, double>(res_sh);
137
138 check_one(res_an, results.andoyer * d2r, res_vi, 0.001);
139 check_one(res_th, results.thomas * d2r, res_vi, 0.01);//in some tests thomas gives low accuracy
140 check_one(res_vi, results.vincenty * d2r, res_vi, 0.0000001);
141 check_one(res_sh, results.spherical * d2r, res_vi, 1);
142 }
143
test_main(int,char * [])144 int test_main(int, char*[])
145 {
146
147 for (size_t i = 0; i < expected_size; ++i)
148 {
149 test_all(expected[i]);
150 }
151
152 return 0;
153 }
154
155