geometry/formulas/differential_quantities.hpp

0001 // Boost.Geometry
0002
0003 // Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
0004
0005 // Copyright (c) 2016-2019 Oracle and/or its affiliates.
0006 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0007
0008 // Use, modification and distribution is subject to the Boost Software License,
0009 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0010 // http://www.boost.org/LICENSE_1_0.txt)
0011
0012 #ifndef BOOST_GEOMETRY_FORMULAS_INVERSE_DIFFERENTIAL_QUANTITIES_HPP
0013 #define BOOST_GEOMETRY_FORMULAS_INVERSE_DIFFERENTIAL_QUANTITIES_HPP
0014
0015 #include <boost/geometry/core/assert.hpp>
0016
0017 #include <boost/geometry/util/condition.hpp>
0018 #include <boost/geometry/util/math.hpp>
0019
0020
0021 namespace boost { namespace geometry { namespace formula
0022 {
0023
0024 /*!
0025 \brief The solution of a part of the inverse problem - differential quantities.
0026 \author See
0027 - Charles F.F Karney, Algorithms for geodesics, 2011
0028 https://arxiv.org/pdf/1109.4448.pdf
0029 */
0030 template <
0031     typename CT,
0032     bool EnableReducedLength,
0033     bool EnableGeodesicScale,
0034     unsigned int Order = 2,
0035     bool ApproxF = true
0036 >
0037 class differential_quantities
0038 {
0039 public:
0040     static inline void apply(CT const& lon1, CT const& lat1,
0041                              CT const& lon2, CT const& lat2,
0042                              CT const& azimuth, CT const& reverse_azimuth,
0043                              CT const& b, CT const& f,
0044                              CT & reduced_length, CT & geodesic_scale)
0045     {
0046         CT const dlon = lon2 - lon1;
0047         CT const sin_lat1 = sin(lat1);
0048         CT const cos_lat1 = cos(lat1);
0049         CT const sin_lat2 = sin(lat2);
0050         CT const cos_lat2 = cos(lat2);
0051
0052         apply(dlon, sin_lat1, cos_lat1, sin_lat2, cos_lat2,
0053               azimuth, reverse_azimuth,
0054               b, f,
0055               reduced_length, geodesic_scale);
0056     }
0057
0058     static inline void apply(CT const& dlon,
0059                              CT const& sin_lat1, CT const& cos_lat1,
0060                              CT const& sin_lat2, CT const& cos_lat2,
0061                              CT const& azimuth, CT const& reverse_azimuth,
0062                              CT const& b, CT const& f,
0063                              CT & reduced_length, CT & geodesic_scale)
0064     {
0065         CT const c0 = 0;
0066         CT const c1 = 1;
0067         CT const one_minus_f = c1 - f;
0068
0069         CT sin_bet1 = one_minus_f * sin_lat1;
0070         CT sin_bet2 = one_minus_f * sin_lat2;
0071
0072         // equator
0073         if (math::equals(sin_bet1, c0) && math::equals(sin_bet2, c0))
0074         {
0075             CT const sig_12 = dlon / one_minus_f;
0076             if (BOOST_GEOMETRY_CONDITION(EnableReducedLength))
0077             {
0078                 BOOST_GEOMETRY_ASSERT((-math::pi<CT>() <= azimuth && azimuth <= math::pi<CT>()));
0079
0080                 int azi_sign = math::sign(azimuth) >= 0 ? 1 : -1; // for antipodal
0081                 CT m12 = azi_sign * sin(sig_12) * b;
0082                 reduced_length = m12;
0083             }
0084
0085             if (BOOST_GEOMETRY_CONDITION(EnableGeodesicScale))
0086             {
0087                 CT M12 = cos(sig_12);
0088                 geodesic_scale = M12;
0089             }
0090         }
0091         else
0092         {
0093             CT const c2 = 2;
0094             CT const e2 = f * (c2 - f);
0095             CT const ep2 = e2 / math::sqr(one_minus_f);
0096
0097             CT const sin_alp1 = sin(azimuth);
0098             CT const cos_alp1 = cos(azimuth);
0099             //CT const sin_alp2 = sin(reverse_azimuth);
0100             CT const cos_alp2 = cos(reverse_azimuth);
0101
0102             CT cos_bet1 = cos_lat1;
0103             CT cos_bet2 = cos_lat2;
0104
0105             normalize(sin_bet1, cos_bet1);
0106             normalize(sin_bet2, cos_bet2);
0107
0108             CT sin_sig1 = sin_bet1;
0109             CT cos_sig1 = cos_alp1 * cos_bet1;
0110             CT sin_sig2 = sin_bet2;
0111             CT cos_sig2 = cos_alp2 * cos_bet2;
0112
0113             normalize(sin_sig1, cos_sig1);
0114             normalize(sin_sig2, cos_sig2);
0115
0116             CT const sin_alp0 = sin_alp1 * cos_bet1;
0117             CT const cos_alp0_sqr = c1 - math::sqr(sin_alp0);
0118
0119             CT const J12 = BOOST_GEOMETRY_CONDITION(ApproxF) ?
0120                            J12_f(sin_sig1, cos_sig1, sin_sig2, cos_sig2, cos_alp0_sqr, f) :
0121                            J12_ep_sqr(sin_sig1, cos_sig1, sin_sig2, cos_sig2, cos_alp0_sqr, ep2) ;
0122
0123             CT const dn1 = math::sqrt(c1 + ep2 * math::sqr(sin_bet1));
0124             CT const dn2 = math::sqrt(c1 + ep2 * math::sqr(sin_bet2));
0125
0126             if (BOOST_GEOMETRY_CONDITION(EnableReducedLength))
0127             {
0128                 CT const m12_b = dn2 * (cos_sig1 * sin_sig2)
0129                                - dn1 * (sin_sig1 * cos_sig2)
0130                                - cos_sig1 * cos_sig2 * J12;
0131                 CT const m12 = m12_b * b;
0132
0133                 reduced_length = m12;
0134             }
0135
0136             if (BOOST_GEOMETRY_CONDITION(EnableGeodesicScale))
0137             {
0138                 CT const cos_sig12 = cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2;
0139                 CT const t = ep2 * (cos_bet1 - cos_bet2) * (cos_bet1 + cos_bet2) / (dn1 + dn2);
0140                 CT const M12 = cos_sig12 + (t * sin_sig2 - cos_sig2 * J12) * sin_sig1 / dn1;
0141
0142                 geodesic_scale = M12;
0143             }
0144         }
0145     }
0146
0147 private:
0148     /*! Approximation of J12, expanded into taylor series in f
0149         Maxima script:
0150         ep2: f * (2-f) / ((1-f)^2);
0151         k2: ca02 * ep2;
0152         assume(f < 1);
0153         assume(sig > 0);
0154         I1(sig):= integrate(sqrt(1 + k2 * sin(s)^2), s, 0, sig);
0155         I2(sig):= integrate(1/sqrt(1 + k2 * sin(s)^2), s, 0, sig);
0156         J(sig):= I1(sig) - I2(sig);
0157         S: taylor(J(sig), f, 0, 3);
0158         S1: factor( 2*integrate(sin(s)^2,s,0,sig)*ca02*f );
0159         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 );
0160         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 );
0161     */
0162     static inline CT J12_f(CT const& sin_sig1, CT const& cos_sig1,
0163                            CT const& sin_sig2, CT const& cos_sig2,
0164                            CT const& cos_alp0_sqr, CT const& f)
0165     {
0166         if (BOOST_GEOMETRY_CONDITION(Order == 0))
0167         {
0168             return 0;
0169         }
0170
0171         CT const c2 = 2;
0172
0173         CT const sig_12 = atan2(cos_sig1 * sin_sig2 - sin_sig1 * cos_sig2,
0174                                 cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2);
0175         CT const sin_2sig1 = c2 * cos_sig1 * sin_sig1; // sin(2sig1)
0176         CT const sin_2sig2 = c2 * cos_sig2 * sin_sig2; // sin(2sig2)
0177         CT const sin_2sig_12 = sin_2sig2 - sin_2sig1;
0178         CT const L1 = sig_12 - sin_2sig_12 / c2;
0179
0180         if (BOOST_GEOMETRY_CONDITION(Order == 1))
0181         {
0182             return cos_alp0_sqr * f * L1;
0183         }
0184
0185         CT const sin_4sig1 = c2 * sin_2sig1 * (math::sqr(cos_sig1) - math::sqr(sin_sig1)); // sin(4sig1)
0186         CT const sin_4sig2 = c2 * sin_2sig2 * (math::sqr(cos_sig2) - math::sqr(sin_sig2)); // sin(4sig2)
0187         CT const sin_4sig_12 = sin_4sig2 - sin_4sig1;
0188
0189         CT const c8 = 8;
0190         CT const c12 = 12;
0191         CT const c16 = 16;
0192         CT const c24 = 24;
0193
0194         CT const L2 = -( cos_alp0_sqr * sin_4sig_12
0195                          + (-c8 * cos_alp0_sqr + c12) * sin_2sig_12
0196                          + (c12 * cos_alp0_sqr - c24) * sig_12)
0197                        / c16;
0198
0199         if (BOOST_GEOMETRY_CONDITION(Order == 2))
0200         {
0201             return cos_alp0_sqr * f * (L1 + f * L2);
0202         }
0203
0204         CT const c4 = 4;
0205         CT const c9 = 9;
0206         CT const c48 = 48;
0207         CT const c60 = 60;
0208         CT const c64 = 64;
0209         CT const c96 = 96;
0210         CT const c128 = 128;
0211         CT const c144 = 144;
0212
0213         CT const cos_alp0_quad = math::sqr(cos_alp0_sqr);
0214         CT const sin3_2sig1 = math::sqr(sin_2sig1) * sin_2sig1;
0215         CT const sin3_2sig2 = math::sqr(sin_2sig2) * sin_2sig2;
0216         CT const sin3_2sig_12 = sin3_2sig2 - sin3_2sig1;
0217
0218         CT const A = (c9 * cos_alp0_quad - c12 * cos_alp0_sqr) * sin_4sig_12;
0219         CT const B = c4 * cos_alp0_quad * sin3_2sig_12;
0220         CT const C = (-c48 * cos_alp0_quad + c96 * cos_alp0_sqr - c64) * sin_2sig_12;
0221         CT const D = (c60 * cos_alp0_quad - c144 * cos_alp0_sqr + c128) * sig_12;
0222
0223         CT const L3 = (A + B + C + D) / c64;
0224
0225         // Order 3 and higher
0226         return cos_alp0_sqr * f * (L1 + f * (L2 + f * L3));
0227     }
0228
0229     /*! Approximation of J12, expanded into taylor series in e'^2
0230         Maxima script:
0231         k2: ca02 * ep2;
0232         assume(sig > 0);
0233         I1(sig):= integrate(sqrt(1 + k2 * sin(s)^2), s, 0, sig);
0234         I2(sig):= integrate(1/sqrt(1 + k2 * sin(s)^2), s, 0, sig);
0235         J(sig):= I1(sig) - I2(sig);
0236         S: taylor(J(sig), ep2, 0, 3);
0237         S1: factor( integrate(sin(s)^2,s,0,sig)*ca02*ep2 );
0238         S2: factor( (integrate(sin(s)^4,s,0,sig)*ca02^2*ep2^2)/2 );
0239         S3: factor( (3*integrate(sin(s)^6,s,0,sig)*ca02^3*ep2^3)/8 );
0240     */
0241     static inline CT J12_ep_sqr(CT const& sin_sig1, CT const& cos_sig1,
0242                                 CT const& sin_sig2, CT const& cos_sig2,
0243                                 CT const& cos_alp0_sqr, CT const& ep_sqr)
0244     {
0245         if (BOOST_GEOMETRY_CONDITION(Order == 0))
0246         {
0247             return 0;
0248         }
0249
0250         CT const c2 = 2;
0251         CT const c4 = 4;
0252
0253         CT const c2a0ep2 = cos_alp0_sqr * ep_sqr;
0254
0255         CT const sig_12 = atan2(cos_sig1 * sin_sig2 - sin_sig1 * cos_sig2,
0256                                 cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2); // sig2 - sig1
0257         CT const sin_2sig1 = c2 * cos_sig1 * sin_sig1; // sin(2sig1)
0258         CT const sin_2sig2 = c2 * cos_sig2 * sin_sig2; // sin(2sig2)
0259         CT const sin_2sig_12 = sin_2sig2 - sin_2sig1;
0260
0261         CT const L1 = (c2 * sig_12 - sin_2sig_12) / c4;
0262
0263         if (BOOST_GEOMETRY_CONDITION(Order == 1))
0264         {
0265             return c2a0ep2 * L1;
0266         }
0267
0268         CT const c8 = 8;
0269         CT const c64 = 64;
0270
0271         CT const sin_4sig1 = c2 * sin_2sig1 * (math::sqr(cos_sig1) - math::sqr(sin_sig1)); // sin(4sig1)
0272         CT const sin_4sig2 = c2 * sin_2sig2 * (math::sqr(cos_sig2) - math::sqr(sin_sig2)); // sin(4sig2)
0273         CT const sin_4sig_12 = sin_4sig2 - sin_4sig1;
0274
0275         CT const L2 = (sin_4sig_12 - c8 * sin_2sig_12 + 12 * sig_12) / c64;
0276
0277         if (BOOST_GEOMETRY_CONDITION(Order == 2))
0278         {
0279             return c2a0ep2 * (L1 + c2a0ep2 * L2);
0280         }
0281
0282         CT const sin3_2sig1 = math::sqr(sin_2sig1) * sin_2sig1;
0283         CT const sin3_2sig2 = math::sqr(sin_2sig2) * sin_2sig2;
0284         CT const sin3_2sig_12 = sin3_2sig2 - sin3_2sig1;
0285
0286         CT const c9 = 9;
0287         CT const c48 = 48;
0288         CT const c60 = 60;
0289         CT const c512 = 512;
0290
0291         CT const L3 = (c9 * sin_4sig_12 + c4 * sin3_2sig_12 - c48 * sin_2sig_12 + c60 * sig_12) / c512;
0292
0293         // Order 3 and higher
0294         return c2a0ep2 * (L1 + c2a0ep2 * (L2 + c2a0ep2 * L3));
0295     }
0296
0297     static inline void normalize(CT & x, CT & y)
0298     {
0299         CT const len = math::sqrt(math::sqr(x) + math::sqr(y));
0300         x /= len;
0301         y /= len;
0302     }
0303 };
0304
0305 }}} // namespace boost::geometry::formula
0306
0307
0308 #endif // BOOST_GEOMETRY_FORMULAS_INVERSE_DIFFERENTIAL_QUANTITIES_HPP