1 // Boost.Geometry
2 
3 // Copyright (c) 2017-2018 Oracle and/or its affiliates.
4 
5 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
6 
7 // Use, modification and distribution is subject to the Boost Software License,
8 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
9 // http://www.boost.org/LICENSE_1_0.txt)
10 
11 #ifndef BOOST_GEOMETRY_FORMULAS_MERIDIAN_INVERSE_HPP
12 #define BOOST_GEOMETRY_FORMULAS_MERIDIAN_INVERSE_HPP
13 
14 #include <boost/math/constants/constants.hpp>
15 
16 #include <boost/geometry/core/radius.hpp>
17 
18 #include <boost/geometry/util/condition.hpp>
19 #include <boost/geometry/util/math.hpp>
20 #include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
21 
22 #include <boost/geometry/formulas/flattening.hpp>
23 #include <boost/geometry/formulas/meridian_segment.hpp>
24 
25 namespace boost { namespace geometry { namespace formula
26 {
27 
28 /*!
29 \brief Compute the arc length of an ellipse.
30 */
31 
32 template <typename CT, unsigned int Order = 1>
33 class meridian_inverse
34 {
35 
36 public :
37 
38     struct result
39     {
resultboost::geometry::formula::meridian_inverse::result40         result()
41             : distance(0)
42             , meridian(false)
43         {}
44 
45         CT distance;
46         bool meridian;
47     };
48 
49     template <typename T>
meridian_not_crossing_pole(T lat1,T lat2,CT diff)50     static bool meridian_not_crossing_pole(T lat1, T lat2, CT diff)
51     {
52         CT half_pi = math::pi<CT>()/CT(2);
53         return math::equals(diff, CT(0)) ||
54                     (math::equals(lat2, half_pi) && math::equals(lat1, -half_pi));
55     }
56 
meridian_crossing_pole(CT diff)57     static bool meridian_crossing_pole(CT diff)
58     {
59         return math::equals(math::abs(diff), math::pi<CT>());
60     }
61 
62 
63     template <typename T, typename Spheroid>
meridian_not_crossing_pole_dist(T lat1,T lat2,Spheroid const & spheroid)64     static CT meridian_not_crossing_pole_dist(T lat1, T lat2, Spheroid const& spheroid)
65     {
66         return math::abs(apply(lat2, spheroid) - apply(lat1, spheroid));
67     }
68 
69     template <typename T, typename Spheroid>
meridian_crossing_pole_dist(T lat1,T lat2,Spheroid const & spheroid)70     static CT meridian_crossing_pole_dist(T lat1, T lat2, Spheroid const& spheroid)
71     {
72         CT c0 = 0;
73         CT half_pi = math::pi<CT>()/CT(2);
74         CT lat_sign = 1;
75         if (lat1+lat2 < c0)
76         {
77             lat_sign = CT(-1);
78         }
79         return math::abs(lat_sign * CT(2) * apply(half_pi, spheroid)
80                          - apply(lat1, spheroid) - apply(lat2, spheroid));
81     }
82 
83     template <typename T, typename Spheroid>
apply(T lon1,T lat1,T lon2,T lat2,Spheroid const & spheroid)84     static result apply(T lon1, T lat1, T lon2, T lat2, Spheroid const& spheroid)
85     {
86         result res;
87 
88         CT diff = geometry::math::longitude_distance_signed<geometry::radian>(lon1, lon2);
89 
90         if (lat1 > lat2)
91         {
92             std::swap(lat1, lat2);
93         }
94 
95         if ( meridian_not_crossing_pole(lat1, lat2, diff) )
96         {
97             res.distance = meridian_not_crossing_pole_dist(lat1, lat2, spheroid);
98             res.meridian = true;
99         }
100         else if ( meridian_crossing_pole(diff) )
101         {
102             res.distance = meridian_crossing_pole_dist(lat1, lat2, spheroid);
103             res.meridian = true;
104         }
105         return res;
106     }
107 
108     // Distance computation on meridians using series approximations
109     // to elliptic integrals. Formula to compute distance from lattitude 0 to lat
110     // https://en.wikipedia.org/wiki/Meridian_arc
111     // latitudes are assumed to be in radians and in [-pi/2,pi/2]
112     template <typename T, typename Spheroid>
apply(T lat,Spheroid const & spheroid)113     static CT apply(T lat, Spheroid const& spheroid)
114     {
115         CT const a = get_radius<0>(spheroid);
116         CT const f = formula::flattening<CT>(spheroid);
117         CT n = f / (CT(2) - f);
118         CT M = a/(1+n);
119         CT C0 = 1;
120 
121         if (Order == 0)
122         {
123            return M * C0 * lat;
124         }
125 
126         CT C2 = -1.5 * n;
127 
128         if (Order == 1)
129         {
130             return M * (C0 * lat + C2 * sin(2*lat));
131         }
132 
133         CT n2 = n * n;
134         C0 += .25 * n2;
135         CT C4 = 0.9375 * n2;
136 
137         if (Order == 2)
138         {
139             return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat));
140         }
141 
142         CT n3 = n2 * n;
143         C2 += 0.1875 * n3;
144         CT C6 = -0.729166667 * n3;
145 
146         if (Order == 3)
147         {
148             return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat)
149                       + C6 * sin(6*lat));
150         }
151 
152         CT n4 = n2 * n2;
153         C4 -= 0.234375 * n4;
154         CT C8 = 0.615234375 * n4;
155 
156         if (Order == 4)
157         {
158             return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat)
159                       + C6 * sin(6*lat) + C8 * sin(8*lat));
160         }
161 
162         CT n5 = n4 * n;
163         C6 += 0.227864583 * n5;
164         CT C10 = -0.54140625 * n5;
165 
166         // Order 5 or higher
167         return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat)
168                   + C6 * sin(6*lat) + C8 * sin(8*lat) + C10 * sin(10*lat));
169 
170     }
171 };
172 
173 }}} // namespace boost::geometry::formula
174 
175 
176 #endif // BOOST_GEOMETRY_FORMULAS_MERIDIAN_INVERSE_HPP
177