geometry/formulas/meridian_direct.hpp

0001 // Boost.Geometry
0002
0003 // Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
0004
0005 // Copyright (c) 2018 Oracle and/or its affiliates.
0006 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0007 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0008
0009 // Use, modification and distribution is subject to the Boost Software License,
0010 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0011 // http://www.boost.org/LICENSE_1_0.txt)
0012
0013 #ifndef BOOST_GEOMETRY_FORMULAS_MERIDIAN_DIRECT_HPP
0014 #define BOOST_GEOMETRY_FORMULAS_MERIDIAN_DIRECT_HPP
0015
0016 #include <boost/math/constants/constants.hpp>
0017
0018 #include <boost/geometry/core/radius.hpp>
0019
0020 #include <boost/geometry/formulas/differential_quantities.hpp>
0021 #include <boost/geometry/formulas/flattening.hpp>
0022 #include <boost/geometry/formulas/meridian_inverse.hpp>
0023 #include <boost/geometry/formulas/quarter_meridian.hpp>
0024 #include <boost/geometry/formulas/result_direct.hpp>
0025
0026 #include <boost/geometry/util/condition.hpp>
0027 #include <boost/geometry/util/math.hpp>
0028
0029 namespace boost { namespace geometry { namespace formula
0030 {
0031
0032 /*!
0033 \brief Compute the direct geodesic problem on a meridian
0034 */
0035
0036 template <
0037     typename CT,
0038     bool EnableCoordinates = true,
0039     bool EnableReverseAzimuth = false,
0040     bool EnableReducedLength = false,
0041     bool EnableGeodesicScale = false,
0042     unsigned int Order = 4
0043 >
0044 class meridian_direct
0045 {
0046     static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale;
0047     static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcQuantities;
0048     static const bool CalcCoordinates = EnableCoordinates || CalcRevAzimuth;
0049
0050 public:
0051     typedef result_direct<CT> result_type;
0052
0053     template <typename T, typename Dist, typename Spheroid>
0054     static inline result_type apply(T const& lo1,
0055                                     T const& la1,
0056                                     Dist const& distance,
0057                                     bool north,
0058                                     Spheroid const& spheroid)
0059     {
0060         result_type result;
0061
0062         CT const half_pi = math::half_pi<CT>();
0063         CT const pi = math::pi<CT>();
0064         CT const one_and_a_half_pi = pi + half_pi;
0065         CT const c0 = 0;
0066
0067         CT azimuth = north ? c0 : pi;
0068
0069         if (BOOST_GEOMETRY_CONDITION(CalcCoordinates))
0070         {
0071             CT s0 = meridian_inverse<CT, Order>::apply(la1, spheroid);
0072             int signed_distance = north ? distance : -distance;
0073             result.lon2 = lo1;
0074             result.lat2 = apply(s0 + signed_distance, spheroid);
0075         }
0076
0077         if (BOOST_GEOMETRY_CONDITION(CalcRevAzimuth))
0078         {
0079             result.reverse_azimuth = azimuth;
0080
0081
0082             if (result.lat2 > half_pi &&
0083                 result.lat2 < one_and_a_half_pi)
0084             {
0085                 result.reverse_azimuth =  pi;
0086             }
0087             else if (result.lat2 > -one_and_a_half_pi &&
0088                      result.lat2 < -half_pi)
0089             {
0090                 result.reverse_azimuth =  c0;
0091             }
0092
0093         }
0094
0095         if (BOOST_GEOMETRY_CONDITION(CalcQuantities))
0096         {
0097             CT const b = CT(get_radius<2>(spheroid));
0098             CT const f = formula::flattening<CT>(spheroid);
0099
0100             boost::geometry::math::normalize_spheroidal_coordinates
0101                 <
0102                     boost::geometry::radian,
0103                     double
0104                 >(result.lon2, result.lat2);
0105
0106             typedef differential_quantities
0107             <
0108                 CT,
0109                 EnableReducedLength,
0110                 EnableGeodesicScale,
0111                 Order
0112             > quantities;
0113             quantities::apply(lo1, la1, result.lon2, result.lat2,
0114                               azimuth, result.reverse_azimuth,
0115                               b, f,
0116                               result.reduced_length, result.geodesic_scale);
0117         }
0118         return result;
0119     }
0120
0121     // https://en.wikipedia.org/wiki/Meridian_arc#The_inverse_meridian_problem_for_the_ellipsoid
0122     // latitudes are assumed to be in radians and in [-pi/2,pi/2]
0123     template <typename T, typename Spheroid>
0124     static CT apply(T m, Spheroid const& spheroid)
0125     {
0126         CT const f = formula::flattening<CT>(spheroid);
0127         CT n = f / (CT(2) - f);
0128         CT mp = formula::quarter_meridian<CT>(spheroid);
0129         CT mu = geometry::math::pi<CT>()/CT(2) * m / mp;
0130
0131         if (BOOST_GEOMETRY_CONDITION(Order == 0))
0132         {
0133             return mu;
0134         }
0135
0136         CT H2 = 1.5 * n;
0137
0138         if (BOOST_GEOMETRY_CONDITION(Order == 1))
0139         {
0140             return mu + H2 * sin(2*mu);
0141         }
0142
0143         CT n2 = n * n;
0144         CT H4 = 1.3125 * n2;
0145
0146         if (BOOST_GEOMETRY_CONDITION(Order == 2))
0147         {
0148             return mu + H2 * sin(2*mu) + H4 * sin(4*mu);
0149         }
0150
0151         CT n3 = n2 * n;
0152         H2 -= 0.84375 * n3;
0153         CT H6 = 1.572916667 * n3;
0154
0155         if (BOOST_GEOMETRY_CONDITION(Order == 3))
0156         {
0157             return mu + H2 * sin(2*mu) + H4 * sin(4*mu) + H6 * sin(6*mu);
0158         }
0159
0160         CT n4 = n2 * n2;
0161         H4 -= 1.71875 * n4;
0162         CT H8 = 2.142578125 * n4;
0163
0164         // Order 4 or higher
0165         return mu + H2 * sin(2*mu) + H4 * sin(4*mu) + H6 * sin(6*mu) + H8 * sin(8*mu);
0166     }
0167 };
0168
0169 }}} // namespace boost::geometry::formula
0170
0171 #endif // BOOST_GEOMETRY_FORMULAS_MERIDIAN_DIRECT_HPP