1 // Boost.Geometry 2 3 // Copyright (c) 2016-2019 Oracle and/or its affiliates. 4 5 // Contributed and/or modified by Adam Wulkiewicz, 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_INVERSE_DIFFERENTIAL_QUANTITIES_HPP 12 #define BOOST_GEOMETRY_FORMULAS_INVERSE_DIFFERENTIAL_QUANTITIES_HPP 13 14 #include <boost/geometry/core/assert.hpp> 15 16 #include <boost/geometry/util/condition.hpp> 17 #include <boost/geometry/util/math.hpp> 18 19 20 namespace boost { namespace geometry { namespace formula 21 { 22 23 /*! 24 \brief The solution of a part of the inverse problem - differential quantities. 25 \author See 26 - Charles F.F Karney, Algorithms for geodesics, 2011 27 https://arxiv.org/pdf/1109.4448.pdf 28 */ 29 template < 30 typename CT, 31 bool EnableReducedLength, 32 bool EnableGeodesicScale, 33 unsigned int Order = 2, 34 bool ApproxF = true 35 > 36 class differential_quantities 37 { 38 public: apply(CT const & lon1,CT const & lat1,CT const & lon2,CT const & lat2,CT const & azimuth,CT const & reverse_azimuth,CT const & b,CT const & f,CT & reduced_length,CT & geodesic_scale)39 static inline void apply(CT const& lon1, CT const& lat1, 40 CT const& lon2, CT const& lat2, 41 CT const& azimuth, CT const& reverse_azimuth, 42 CT const& b, CT const& f, 43 CT & reduced_length, CT & geodesic_scale) 44 { 45 CT const dlon = lon2 - lon1; 46 CT const sin_lat1 = sin(lat1); 47 CT const cos_lat1 = cos(lat1); 48 CT const sin_lat2 = sin(lat2); 49 CT const cos_lat2 = cos(lat2); 50 51 apply(dlon, sin_lat1, cos_lat1, sin_lat2, cos_lat2, 52 azimuth, reverse_azimuth, 53 b, f, 54 reduced_length, geodesic_scale); 55 } 56 apply(CT const & dlon,CT const & sin_lat1,CT const & cos_lat1,CT const & sin_lat2,CT const & cos_lat2,CT const & azimuth,CT const & reverse_azimuth,CT const & b,CT const & f,CT & reduced_length,CT & geodesic_scale)57 static inline void apply(CT const& dlon, 58 CT const& sin_lat1, CT const& cos_lat1, 59 CT const& sin_lat2, CT const& cos_lat2, 60 CT const& azimuth, CT const& reverse_azimuth, 61 CT const& b, CT const& f, 62 CT & reduced_length, CT & geodesic_scale) 63 { 64 CT const c0 = 0; 65 CT const c1 = 1; 66 CT const one_minus_f = c1 - f; 67 68 CT sin_bet1 = one_minus_f * sin_lat1; 69 CT sin_bet2 = one_minus_f * sin_lat2; 70 71 // equator 72 if (math::equals(sin_bet1, c0) && math::equals(sin_bet2, c0)) 73 { 74 CT const sig_12 = dlon / one_minus_f; 75 if (BOOST_GEOMETRY_CONDITION(EnableReducedLength)) 76 { 77 BOOST_GEOMETRY_ASSERT((-math::pi<CT>() <= azimuth && azimuth <= math::pi<CT>())); 78 79 int azi_sign = math::sign(azimuth) >= 0 ? 1 : -1; // for antipodal 80 CT m12 = azi_sign * sin(sig_12) * b; 81 reduced_length = m12; 82 } 83 84 if (BOOST_GEOMETRY_CONDITION(EnableGeodesicScale)) 85 { 86 CT M12 = cos(sig_12); 87 geodesic_scale = M12; 88 } 89 } 90 else 91 { 92 CT const c2 = 2; 93 CT const e2 = f * (c2 - f); 94 CT const ep2 = e2 / math::sqr(one_minus_f); 95 96 CT const sin_alp1 = sin(azimuth); 97 CT const cos_alp1 = cos(azimuth); 98 //CT const sin_alp2 = sin(reverse_azimuth); 99 CT const cos_alp2 = cos(reverse_azimuth); 100 101 CT cos_bet1 = cos_lat1; 102 CT cos_bet2 = cos_lat2; 103 104 normalize(sin_bet1, cos_bet1); 105 normalize(sin_bet2, cos_bet2); 106 107 CT sin_sig1 = sin_bet1; 108 CT cos_sig1 = cos_alp1 * cos_bet1; 109 CT sin_sig2 = sin_bet2; 110 CT cos_sig2 = cos_alp2 * cos_bet2; 111 112 normalize(sin_sig1, cos_sig1); 113 normalize(sin_sig2, cos_sig2); 114 115 CT const sin_alp0 = sin_alp1 * cos_bet1; 116 CT const cos_alp0_sqr = c1 - math::sqr(sin_alp0); 117 118 CT const J12 = BOOST_GEOMETRY_CONDITION(ApproxF) ? 119 J12_f(sin_sig1, cos_sig1, sin_sig2, cos_sig2, cos_alp0_sqr, f) : 120 J12_ep_sqr(sin_sig1, cos_sig1, sin_sig2, cos_sig2, cos_alp0_sqr, ep2) ; 121 122 CT const dn1 = math::sqrt(c1 + ep2 * math::sqr(sin_bet1)); 123 CT const dn2 = math::sqrt(c1 + ep2 * math::sqr(sin_bet2)); 124 125 if (BOOST_GEOMETRY_CONDITION(EnableReducedLength)) 126 { 127 CT const m12_b = dn2 * (cos_sig1 * sin_sig2) 128 - dn1 * (sin_sig1 * cos_sig2) 129 - cos_sig1 * cos_sig2 * J12; 130 CT const m12 = m12_b * b; 131 132 reduced_length = m12; 133 } 134 135 if (BOOST_GEOMETRY_CONDITION(EnableGeodesicScale)) 136 { 137 CT const cos_sig12 = cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2; 138 CT const t = ep2 * (cos_bet1 - cos_bet2) * (cos_bet1 + cos_bet2) / (dn1 + dn2); 139 CT const M12 = cos_sig12 + (t * sin_sig2 - cos_sig2 * J12) * sin_sig1 / dn1; 140 141 geodesic_scale = M12; 142 } 143 } 144 } 145 146 private: 147 /*! Approximation of J12, expanded into taylor series in f 148 Maxima script: 149 ep2: f * (2-f) / ((1-f)^2); 150 k2: ca02 * ep2; 151 assume(f < 1); 152 assume(sig > 0); 153 I1(sig):= integrate(sqrt(1 + k2 * sin(s)^2), s, 0, sig); 154 I2(sig):= integrate(1/sqrt(1 + k2 * sin(s)^2), s, 0, sig); 155 J(sig):= I1(sig) - I2(sig); 156 S: taylor(J(sig), f, 0, 3); 157 S1: factor( 2*integrate(sin(s)^2,s,0,sig)*ca02*f ); 158 S2: factor( ((integrate(-6*ca02^2*sin(s)^4+6*ca02*sin(s)^2,s,0,sig)+integrate(-2*ca02^2*sin(s)^4+6*ca02*sin(s)^2,s,0,sig))*f^2)/4 ); 159 S3: factor( ((integrate(30*ca02^3*sin(s)^6-54*ca02^2*sin(s)^4+24*ca02*sin(s)^2,s,0,sig)+integrate(6*ca02^3*sin(s)^6-18*ca02^2*sin(s)^4+24*ca02*sin(s)^2,s,0,sig))*f^3)/12 ); 160 */ J12_f(CT const & sin_sig1,CT const & cos_sig1,CT const & sin_sig2,CT const & cos_sig2,CT const & cos_alp0_sqr,CT const & f)161 static inline CT J12_f(CT const& sin_sig1, CT const& cos_sig1, 162 CT const& sin_sig2, CT const& cos_sig2, 163 CT const& cos_alp0_sqr, CT const& f) 164 { 165 if (Order == 0) 166 { 167 return 0; 168 } 169 170 CT const c2 = 2; 171 172 CT const sig_12 = atan2(cos_sig1 * sin_sig2 - sin_sig1 * cos_sig2, 173 cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2); 174 CT const sin_2sig1 = c2 * cos_sig1 * sin_sig1; // sin(2sig1) 175 CT const sin_2sig2 = c2 * cos_sig2 * sin_sig2; // sin(2sig2) 176 CT const sin_2sig_12 = sin_2sig2 - sin_2sig1; 177 CT const L1 = sig_12 - sin_2sig_12 / c2; 178 179 if (Order == 1) 180 { 181 return cos_alp0_sqr * f * L1; 182 } 183 184 CT const sin_4sig1 = c2 * sin_2sig1 * (math::sqr(cos_sig1) - math::sqr(sin_sig1)); // sin(4sig1) 185 CT const sin_4sig2 = c2 * sin_2sig2 * (math::sqr(cos_sig2) - math::sqr(sin_sig2)); // sin(4sig2) 186 CT const sin_4sig_12 = sin_4sig2 - sin_4sig1; 187 188 CT const c8 = 8; 189 CT const c12 = 12; 190 CT const c16 = 16; 191 CT const c24 = 24; 192 193 CT const L2 = -( cos_alp0_sqr * sin_4sig_12 194 + (-c8 * cos_alp0_sqr + c12) * sin_2sig_12 195 + (c12 * cos_alp0_sqr - c24) * sig_12) 196 / c16; 197 198 if (Order == 2) 199 { 200 return cos_alp0_sqr * f * (L1 + f * L2); 201 } 202 203 CT const c4 = 4; 204 CT const c9 = 9; 205 CT const c48 = 48; 206 CT const c60 = 60; 207 CT const c64 = 64; 208 CT const c96 = 96; 209 CT const c128 = 128; 210 CT const c144 = 144; 211 212 CT const cos_alp0_quad = math::sqr(cos_alp0_sqr); 213 CT const sin3_2sig1 = math::sqr(sin_2sig1) * sin_2sig1; 214 CT const sin3_2sig2 = math::sqr(sin_2sig2) * sin_2sig2; 215 CT const sin3_2sig_12 = sin3_2sig2 - sin3_2sig1; 216 217 CT const A = (c9 * cos_alp0_quad - c12 * cos_alp0_sqr) * sin_4sig_12; 218 CT const B = c4 * cos_alp0_quad * sin3_2sig_12; 219 CT const C = (-c48 * cos_alp0_quad + c96 * cos_alp0_sqr - c64) * sin_2sig_12; 220 CT const D = (c60 * cos_alp0_quad - c144 * cos_alp0_sqr + c128) * sig_12; 221 222 CT const L3 = (A + B + C + D) / c64; 223 224 // Order 3 and higher 225 return cos_alp0_sqr * f * (L1 + f * (L2 + f * L3)); 226 } 227 228 /*! Approximation of J12, expanded into taylor series in e'^2 229 Maxima script: 230 k2: ca02 * ep2; 231 assume(sig > 0); 232 I1(sig):= integrate(sqrt(1 + k2 * sin(s)^2), s, 0, sig); 233 I2(sig):= integrate(1/sqrt(1 + k2 * sin(s)^2), s, 0, sig); 234 J(sig):= I1(sig) - I2(sig); 235 S: taylor(J(sig), ep2, 0, 3); 236 S1: factor( integrate(sin(s)^2,s,0,sig)*ca02*ep2 ); 237 S2: factor( (integrate(sin(s)^4,s,0,sig)*ca02^2*ep2^2)/2 ); 238 S3: factor( (3*integrate(sin(s)^6,s,0,sig)*ca02^3*ep2^3)/8 ); 239 */ J12_ep_sqr(CT const & sin_sig1,CT const & cos_sig1,CT const & sin_sig2,CT const & cos_sig2,CT const & cos_alp0_sqr,CT const & ep_sqr)240 static inline CT J12_ep_sqr(CT const& sin_sig1, CT const& cos_sig1, 241 CT const& sin_sig2, CT const& cos_sig2, 242 CT const& cos_alp0_sqr, CT const& ep_sqr) 243 { 244 if (Order == 0) 245 { 246 return 0; 247 } 248 249 CT const c2 = 2; 250 CT const c4 = 4; 251 252 CT const c2a0ep2 = cos_alp0_sqr * ep_sqr; 253 254 CT const sig_12 = atan2(cos_sig1 * sin_sig2 - sin_sig1 * cos_sig2, 255 cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2); // sig2 - sig1 256 CT const sin_2sig1 = c2 * cos_sig1 * sin_sig1; // sin(2sig1) 257 CT const sin_2sig2 = c2 * cos_sig2 * sin_sig2; // sin(2sig2) 258 CT const sin_2sig_12 = sin_2sig2 - sin_2sig1; 259 260 CT const L1 = (c2 * sig_12 - sin_2sig_12) / c4; 261 262 if (Order == 1) 263 { 264 return c2a0ep2 * L1; 265 } 266 267 CT const c8 = 8; 268 CT const c64 = 64; 269 270 CT const sin_4sig1 = c2 * sin_2sig1 * (math::sqr(cos_sig1) - math::sqr(sin_sig1)); // sin(4sig1) 271 CT const sin_4sig2 = c2 * sin_2sig2 * (math::sqr(cos_sig2) - math::sqr(sin_sig2)); // sin(4sig2) 272 CT const sin_4sig_12 = sin_4sig2 - sin_4sig1; 273 274 CT const L2 = (sin_4sig_12 - c8 * sin_2sig_12 + 12 * sig_12) / c64; 275 276 if (Order == 2) 277 { 278 return c2a0ep2 * (L1 + c2a0ep2 * L2); 279 } 280 281 CT const sin3_2sig1 = math::sqr(sin_2sig1) * sin_2sig1; 282 CT const sin3_2sig2 = math::sqr(sin_2sig2) * sin_2sig2; 283 CT const sin3_2sig_12 = sin3_2sig2 - sin3_2sig1; 284 285 CT const c9 = 9; 286 CT const c48 = 48; 287 CT const c60 = 60; 288 CT const c512 = 512; 289 290 CT const L3 = (c9 * sin_4sig_12 + c4 * sin3_2sig_12 - c48 * sin_2sig_12 + c60 * sig_12) / c512; 291 292 // Order 3 and higher 293 return c2a0ep2 * (L1 + c2a0ep2 * (L2 + c2a0ep2 * L3)); 294 } 295 normalize(CT & x,CT & y)296 static inline void normalize(CT & x, CT & y) 297 { 298 CT const len = math::sqrt(math::sqr(x) + math::sqr(y)); 299 x /= len; 300 y /= len; 301 } 302 }; 303 304 }}} // namespace boost::geometry::formula 305 306 307 #endif // BOOST_GEOMETRY_FORMULAS_INVERSE_DIFFERENTIAL_QUANTITIES_HPP 308