1 // Boost.Geometry
2 
3 // Copyright (c) 2018 Oracle and/or its affiliates.
4 
5 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
6 // Contributed and/or modified by Adam Wulkiewicz, 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 #ifndef BOOST_GEOMETRY_FORMULAS_MERIDIAN_DIRECT_HPP
13 #define BOOST_GEOMETRY_FORMULAS_MERIDIAN_DIRECT_HPP
14 
15 #include <boost/math/constants/constants.hpp>
16 
17 #include <boost/geometry/core/radius.hpp>
18 
19 #include <boost/geometry/formulas/differential_quantities.hpp>
20 #include <boost/geometry/formulas/flattening.hpp>
21 #include <boost/geometry/formulas/meridian_inverse.hpp>
22 #include <boost/geometry/formulas/quarter_meridian.hpp>
23 #include <boost/geometry/formulas/result_direct.hpp>
24 
25 #include <boost/geometry/util/condition.hpp>
26 #include <boost/geometry/util/math.hpp>
27 
28 namespace boost { namespace geometry { namespace formula
29 {
30 
31 /*!
32 \brief Compute the direct geodesic problem on a meridian
33 */
34 
35 template <
36     typename CT,
37     bool EnableCoordinates = true,
38     bool EnableReverseAzimuth = false,
39     bool EnableReducedLength = false,
40     bool EnableGeodesicScale = false,
41     unsigned int Order = 4
42 >
43 class meridian_direct
44 {
45     static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale;
46     static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcQuantities;
47     static const bool CalcCoordinates = EnableCoordinates || CalcRevAzimuth;
48 
49 public:
50     typedef result_direct<CT> result_type;
51 
52     template <typename T, typename Dist, typename Spheroid>
apply(T const & lo1,T const & la1,Dist const & distance,bool north,Spheroid const & spheroid)53     static inline result_type apply(T const& lo1,
54                                     T const& la1,
55                                     Dist const& distance,
56                                     bool north,
57                                     Spheroid const& spheroid)
58     {
59         result_type result;
60 
61         CT const half_pi = math::half_pi<CT>();
62         CT const pi = math::pi<CT>();
63         CT const one_and_a_half_pi = pi + half_pi;
64         CT const c0 = 0;
65 
66         CT azimuth = north ? c0 : pi;
67 
68         if (BOOST_GEOMETRY_CONDITION(CalcCoordinates))
69         {
70             CT s0 = meridian_inverse<CT, Order>::apply(la1, spheroid);
71             int signed_distance = north ? distance : -distance;
72             result.lon2 = lo1;
73             result.lat2 = apply(s0 + signed_distance, spheroid);
74         }
75 
76         if (BOOST_GEOMETRY_CONDITION(CalcRevAzimuth))
77         {
78             result.reverse_azimuth = azimuth;
79 
80 
81             if (result.lat2 > half_pi &&
82                 result.lat2 < one_and_a_half_pi)
83             {
84                 result.reverse_azimuth =  pi;
85             }
86             else if (result.lat2 < -half_pi &&
87                      result.lat2 >  -one_and_a_half_pi)
88             {
89                 result.reverse_azimuth =  c0;
90             }
91 
92         }
93 
94         if (BOOST_GEOMETRY_CONDITION(CalcQuantities))
95         {
96             CT const b = CT(get_radius<2>(spheroid));
97             CT const f = formula::flattening<CT>(spheroid);
98 
99             boost::geometry::math::normalize_spheroidal_coordinates
100                 <
101                     boost::geometry::radian,
102                     double
103                 >(result.lon2, result.lat2);
104 
105             typedef differential_quantities
106             <
107                 CT,
108                 EnableReducedLength,
109                 EnableGeodesicScale,
110                 Order
111             > quantities;
112             quantities::apply(lo1, la1, result.lon2, result.lat2,
113                               azimuth, result.reverse_azimuth,
114                               b, f,
115                               result.reduced_length, result.geodesic_scale);
116         }
117         return result;
118     }
119 
120     // https://en.wikipedia.org/wiki/Meridian_arc#The_inverse_meridian_problem_for_the_ellipsoid
121     // latitudes are assumed to be in radians and in [-pi/2,pi/2]
122     template <typename T, typename Spheroid>
apply(T m,Spheroid const & spheroid)123     static CT apply(T m, Spheroid const& spheroid)
124     {
125         CT const f = formula::flattening<CT>(spheroid);
126         CT n = f / (CT(2) - f);
127         CT mp = formula::quarter_meridian<CT>(spheroid);
128         CT mu = geometry::math::pi<CT>()/CT(2) * m / mp;
129 
130         if (Order == 0)
131         {
132             return mu;
133         }
134 
135         CT H2 = 1.5 * n;
136 
137         if (Order == 1)
138         {
139             return mu + H2 * sin(2*mu);
140         }
141 
142         CT n2 = n * n;
143         CT H4 = 1.3125 * n2;
144 
145         if (Order == 2)
146         {
147             return mu + H2 * sin(2*mu) + H4 * sin(4*mu);
148         }
149 
150         CT n3 = n2 * n;
151         H2 -= 0.84375 * n3;
152         CT H6 = 1.572916667 * n3;
153 
154         if (Order == 3)
155         {
156             return mu + H2 * sin(2*mu) + H4 * sin(4*mu) + H6 * sin(6*mu);
157         }
158 
159         CT n4 = n2 * n2;
160         H4 -= 1.71875 * n4;
161         CT H8 = 2.142578125 * n4;
162 
163         // Order 4 or higher
164         return mu + H2 * sin(2*mu) + H4 * sin(4*mu) + H6 * sin(6*mu) + H8 * sin(8*mu);
165     }
166 };
167 
168 }}} // namespace boost::geometry::formula
169 
170 #endif // BOOST_GEOMETRY_FORMULAS_MERIDIAN_DIRECT_HPP
171