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