geometry/formulas/meridian_inverse.hpp

0001 // Boost.Geometry
0002
0003 // Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
0004
0005 // Copyright (c) 2017-2018 Oracle and/or its affiliates.
0006 // Contributed and/or modified by Vissarion Fysikopoulos, 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_MERIDIAN_INVERSE_HPP
0013 #define BOOST_GEOMETRY_FORMULAS_MERIDIAN_INVERSE_HPP
0014
0015 #include <boost/math/constants/constants.hpp>
0016
0017 #include <boost/geometry/core/radius.hpp>
0018
0019 #include <boost/geometry/util/condition.hpp>
0020 #include <boost/geometry/util/math.hpp>
0021 #include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
0022
0023 #include <boost/geometry/formulas/flattening.hpp>
0024 #include <boost/geometry/formulas/meridian_segment.hpp>
0025
0026 namespace boost { namespace geometry { namespace formula
0027 {
0028
0029 /*!
0030 \brief Compute the arc length of an ellipse.
0031 */
0032
0033 template <typename CT, unsigned int Order = 1>
0034 class meridian_inverse
0035 {
0036
0037 public :
0038
0039     struct result
0040     {
0041         result()
0042             : distance(0)
0043             , meridian(false)
0044         {}
0045
0046         CT distance;
0047         bool meridian;
0048     };
0049
0050     template <typename T>
0051     static bool meridian_not_crossing_pole(T lat1, T lat2, CT diff)
0052     {
0053         CT half_pi = math::pi<CT>()/CT(2);
0054         return math::equals(diff, CT(0)) ||
0055                     (math::equals(lat2, half_pi) && math::equals(lat1, -half_pi));
0056     }
0057
0058     static bool meridian_crossing_pole(CT diff)
0059     {
0060         return math::equals(math::abs(diff), math::pi<CT>());
0061     }
0062
0063
0064     template <typename T, typename Spheroid>
0065     static CT meridian_not_crossing_pole_dist(T lat1, T lat2, Spheroid const& spheroid)
0066     {
0067         return math::abs(apply(lat2, spheroid) - apply(lat1, spheroid));
0068     }
0069
0070     template <typename T, typename Spheroid>
0071     static CT meridian_crossing_pole_dist(T lat1, T lat2, Spheroid const& spheroid)
0072     {
0073         CT c0 = 0;
0074         CT half_pi = math::pi<CT>()/CT(2);
0075         CT lat_sign = 1;
0076         if (lat1+lat2 < c0)
0077         {
0078             lat_sign = CT(-1);
0079         }
0080         return math::abs(lat_sign * CT(2) * apply(half_pi, spheroid)
0081                          - apply(lat1, spheroid) - apply(lat2, spheroid));
0082     }
0083
0084     template <typename T, typename Spheroid>
0085     static result apply(T lon1, T lat1, T lon2, T lat2, Spheroid const& spheroid)
0086     {
0087         result res;
0088
0089         CT diff = geometry::math::longitude_distance_signed<geometry::radian>(lon1, lon2);
0090
0091         if (lat1 > lat2)
0092         {
0093             std::swap(lat1, lat2);
0094         }
0095
0096         if ( meridian_not_crossing_pole(lat1, lat2, diff) )
0097         {
0098             res.distance = meridian_not_crossing_pole_dist(lat1, lat2, spheroid);
0099             res.meridian = true;
0100         }
0101         else if ( meridian_crossing_pole(diff) )
0102         {
0103             res.distance = meridian_crossing_pole_dist(lat1, lat2, spheroid);
0104             res.meridian = true;
0105         }
0106         return res;
0107     }
0108
0109     // Distance computation on meridians using series approximations
0110     // to elliptic integrals. Formula to compute distance from lattitude 0 to lat
0111     // https://en.wikipedia.org/wiki/Meridian_arc
0112     // latitudes are assumed to be in radians and in [-pi/2,pi/2]
0113     template <typename T, typename Spheroid>
0114     static CT apply(T lat, Spheroid const& spheroid)
0115     {
0116         CT const a = get_radius<0>(spheroid);
0117         CT const f = formula::flattening<CT>(spheroid);
0118         CT n = f / (CT(2) - f);
0119         CT M = a/(1+n);
0120         CT C0 = 1;
0121
0122         if (BOOST_GEOMETRY_CONDITION(Order == 0))
0123         {
0124            return M * C0 * lat;
0125         }
0126
0127         CT C2 = -1.5 * n;
0128
0129         if (BOOST_GEOMETRY_CONDITION(Order == 1))
0130         {
0131             return M * (C0 * lat + C2 * sin(2*lat));
0132         }
0133
0134         CT n2 = n * n;
0135         C0 += .25 * n2;
0136         CT C4 = 0.9375 * n2;
0137
0138         if (BOOST_GEOMETRY_CONDITION(Order == 2))
0139         {
0140             return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat));
0141         }
0142
0143         CT n3 = n2 * n;
0144         C2 += 0.1875 * n3;
0145         CT C6 = -0.729166667 * n3;
0146
0147         if (BOOST_GEOMETRY_CONDITION(Order == 3))
0148         {
0149             return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat)
0150                       + C6 * sin(6*lat));
0151         }
0152
0153         CT n4 = n2 * n2;
0154         C4 -= 0.234375 * n4;
0155         CT C8 = 0.615234375 * n4;
0156
0157         if (BOOST_GEOMETRY_CONDITION(Order == 4))
0158         {
0159             return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat)
0160                       + C6 * sin(6*lat) + C8 * sin(8*lat));
0161         }
0162
0163         CT n5 = n4 * n;
0164         C6 += 0.227864583 * n5;
0165         CT C10 = -0.54140625 * n5;
0166
0167         // Order 5 or higher
0168         return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat)
0169                   + C6 * sin(6*lat) + C8 * sin(8*lat) + C10 * sin(10*lat));
0170
0171     }
0172 };
0173
0174 }}} // namespace boost::geometry::formula
0175
0176
0177 #endif // BOOST_GEOMETRY_FORMULAS_MERIDIAN_INVERSE_HPP