geometry/formulas/spherical.hpp

0001 // Boost.Geometry
0002
0003 // Copyright (c) 2016-2020, Oracle and/or its affiliates.
0004 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0005 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0006
0007 // Use, modification and distribution is subject to the Boost Software License,
0008 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0009 // http://www.boost.org/LICENSE_1_0.txt)
0010
0011 #ifndef BOOST_GEOMETRY_FORMULAS_SPHERICAL_HPP
0012 #define BOOST_GEOMETRY_FORMULAS_SPHERICAL_HPP
0013
0014 #include <boost/geometry/core/coordinate_system.hpp>
0015 #include <boost/geometry/core/coordinate_type.hpp>
0016 #include <boost/geometry/core/cs.hpp>
0017 #include <boost/geometry/core/access.hpp>
0018 #include <boost/geometry/core/radian_access.hpp>
0019 #include <boost/geometry/core/radius.hpp>
0020
0021 //#include <boost/geometry/arithmetic/arithmetic.hpp>
0022 #include <boost/geometry/arithmetic/cross_product.hpp>
0023 #include <boost/geometry/arithmetic/dot_product.hpp>
0024
0025 #include <boost/geometry/util/math.hpp>
0026 #include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
0027 #include <boost/geometry/util/select_coordinate_type.hpp>
0028
0029 #include <boost/geometry/formulas/result_direct.hpp>
0030
0031 namespace boost { namespace geometry {
0032
0033 namespace formula {
0034
0035 template <typename T>
0036 struct result_spherical
0037 {
0038     result_spherical()
0039         : azimuth(0)
0040         , reverse_azimuth(0)
0041     {}
0042
0043     T azimuth;
0044     T reverse_azimuth;
0045 };
0046
0047 template <typename T>
0048 static inline void sph_to_cart3d(T const& lon, T const& lat, T & x, T & y, T & z)
0049 {
0050     T const cos_lat = cos(lat);
0051     x = cos_lat * cos(lon);
0052     y = cos_lat * sin(lon);
0053     z = sin(lat);
0054 }
0055
0056 template <typename Point3d, typename PointSph>
0057 static inline Point3d sph_to_cart3d(PointSph const& point_sph)
0058 {
0059     typedef typename coordinate_type<Point3d>::type calc_t;
0060
0061     calc_t const lon = get_as_radian<0>(point_sph);
0062     calc_t const lat = get_as_radian<1>(point_sph);
0063     calc_t x, y, z;
0064     sph_to_cart3d(lon, lat, x, y, z);
0065
0066     Point3d res;
0067     set<0>(res, x);
0068     set<1>(res, y);
0069     set<2>(res, z);
0070
0071     return res;
0072 }
0073
0074 template <typename T>
0075 static inline void cart3d_to_sph(T const& x, T const& y, T const& z, T & lon, T & lat)
0076 {
0077     lon = atan2(y, x);
0078     lat = asin(z);
0079 }
0080
0081 template <typename PointSph, typename Point3d>
0082 static inline PointSph cart3d_to_sph(Point3d const& point_3d)
0083 {
0084     typedef typename coordinate_type<PointSph>::type coord_t;
0085     typedef typename coordinate_type<Point3d>::type calc_t;
0086
0087     calc_t const x = get<0>(point_3d);
0088     calc_t const y = get<1>(point_3d);
0089     calc_t const z = get<2>(point_3d);
0090     calc_t lonr, latr;
0091     cart3d_to_sph(x, y, z, lonr, latr);
0092
0093     PointSph res;
0094     set_from_radian<0>(res, lonr);
0095     set_from_radian<1>(res, latr);
0096
0097     coord_t lon = get<0>(res);
0098     coord_t lat = get<1>(res);
0099
0100     math::normalize_spheroidal_coordinates
0101         <
0102             typename geometry::detail::cs_angular_units<PointSph>::type,
0103             coord_t
0104         >(lon, lat);
0105
0106     set<0>(res, lon);
0107     set<1>(res, lat);
0108
0109     return res;
0110 }
0111
0112 // -1 right
0113 // 1 left
0114 // 0 on
0115 template <typename Point3d1, typename Point3d2>
0116 static inline int sph_side_value(Point3d1 const& norm, Point3d2 const& pt)
0117 {
0118     typedef typename select_coordinate_type<Point3d1, Point3d2>::type calc_t;
0119     calc_t c0 = 0;
0120     calc_t d = dot_product(norm, pt);
0121     return math::equals(d, c0) ? 0
0122         : d > c0 ? 1
0123         : -1; // d < 0
0124 }
0125
0126 template <typename CT, bool ReverseAzimuth, typename T1, typename T2>
0127 static inline result_spherical<CT> spherical_azimuth(T1 const& lon1,
0128                                                      T1 const& lat1,
0129                                                      T2 const& lon2,
0130                                                      T2 const& lat2)
0131 {
0132     typedef result_spherical<CT> result_type;
0133     result_type result;
0134
0135     // http://williams.best.vwh.net/avform.htm#Crs
0136     // https://en.wikipedia.org/wiki/Great-circle_navigation
0137     CT dlon = lon2 - lon1;
0138
0139     // An optimization which should kick in often for Boxes
0140     //if ( math::equals(dlon, ReturnType(0)) )
0141     //if ( get<0>(p1) == get<0>(p2) )
0142     //{
0143     //    return - sin(get_as_radian<1>(p1)) * cos_p2lat);
0144     //}
0145
0146     CT const cos_dlon = cos(dlon);
0147     CT const sin_dlon = sin(dlon);
0148     CT const cos_lat1 = cos(lat1);
0149     CT const cos_lat2 = cos(lat2);
0150     CT const sin_lat1 = sin(lat1);
0151     CT const sin_lat2 = sin(lat2);
0152
0153     {
0154         // "An alternative formula, not requiring the pre-computation of d"
0155         // In the formula below dlon is used as "d"
0156         CT const y = sin_dlon * cos_lat2;
0157         CT const x = cos_lat1 * sin_lat2 - sin_lat1 * cos_lat2 * cos_dlon;
0158         result.azimuth = atan2(y, x);
0159     }
0160
0161     if (ReverseAzimuth)
0162     {
0163         CT const y = sin_dlon * cos_lat1;
0164         CT const x = sin_lat2 * cos_lat1 * cos_dlon - cos_lat2 * sin_lat1;
0165         result.reverse_azimuth = atan2(y, x);
0166     }
0167
0168     return result;
0169 }
0170
0171 template <typename ReturnType, typename T1, typename T2>
0172 inline ReturnType spherical_azimuth(T1 const& lon1, T1 const& lat1,
0173                                     T2 const& lon2, T2 const& lat2)
0174 {
0175     return spherical_azimuth<ReturnType, false>(lon1, lat1, lon2, lat2).azimuth;
0176 }
0177
0178 template <typename T>
0179 inline T spherical_azimuth(T const& lon1, T const& lat1, T const& lon2, T const& lat2)
0180 {
0181     return spherical_azimuth<T, false>(lon1, lat1, lon2, lat2).azimuth;
0182 }
0183
0184 template <typename T>
0185 inline int azimuth_side_value(T const& azi_a1_p, T const& azi_a1_a2)
0186 {
0187     T const c0 = 0;
0188     T const pi = math::pi<T>();
0189
0190     // instead of the formula from XTD
0191     //calc_t a_diff = asin(sin(azi_a1_p - azi_a1_a2));
0192
0193     T a_diff = azi_a1_p - azi_a1_a2;
0194     // normalize, angle in (-pi, pi]
0195     math::detail::normalize_angle_loop<radian>(a_diff);
0196
0197     // NOTE: in general it shouldn't be required to support the pi/-pi case
0198     // because in non-cartesian systems it makes sense to check the side
0199     // only "between" the endpoints.
0200     // However currently the winding strategy calls the side strategy
0201     // for vertical segments to check if the point is "between the endpoints.
0202     // This could be avoided since the side strategy is not required for that
0203     // because meridian is the shortest path. So a difference of
0204     // longitudes would be sufficient (of course normalized to (-pi, pi]).
0205
0206     // NOTE: with the above said, the pi/-pi check is temporary
0207     // however in case if this was required
0208     // the geodesics on ellipsoid aren't "symmetrical"
0209     // therefore instead of comparing a_diff to pi and -pi
0210     // one should probably use inverse azimuths and compare
0211     // the difference to 0 as well
0212
0213     // positive azimuth is on the right side
0214     return math::equals(a_diff, c0)
0215         || math::equals(a_diff, pi)
0216         || math::equals(a_diff, -pi) ? 0
0217         : a_diff > 0 ? -1 // right
0218         : 1; // left
0219 }
0220
0221 template
0222 <
0223     bool Coordinates,
0224     bool ReverseAzimuth,
0225     typename CT,
0226     typename Sphere
0227 >
0228 inline result_direct<CT> spherical_direct(CT const& lon1,
0229                                           CT const& lat1,
0230                                           CT const& sig12,
0231                                           CT const& alp1,
0232                                           Sphere const& sphere)
0233 {
0234     result_direct<CT> result;
0235
0236     CT const sin_alp1 = sin(alp1);
0237     CT const sin_lat1 = sin(lat1);
0238     CT const cos_alp1 = cos(alp1);
0239     CT const cos_lat1 = cos(lat1);
0240
0241     CT const norm = math::sqrt(cos_alp1 * cos_alp1 + sin_alp1 * sin_alp1
0242                                                    * sin_lat1 * sin_lat1);
0243     CT const alp0 = atan2(sin_alp1 * cos_lat1, norm);
0244     CT const sig1 = atan2(sin_lat1, cos_alp1 * cos_lat1);
0245     CT const sig2 = sig1 + sig12 / get_radius<0>(sphere);
0246
0247     CT const cos_sig2 = cos(sig2);
0248     CT const sin_alp0 = sin(alp0);
0249     CT const cos_alp0 = cos(alp0);
0250
0251     if (Coordinates)
0252     {
0253         CT const sin_sig2 = sin(sig2);
0254         CT const sin_sig1 = sin(sig1);
0255         CT const cos_sig1 = cos(sig1);
0256
0257         CT const norm2 = math::sqrt(cos_alp0 * cos_alp0 * cos_sig2 * cos_sig2
0258                                     + sin_alp0 * sin_alp0);
0259         CT const lat2 = atan2(cos_alp0 * sin_sig2, norm2);
0260
0261         CT const omg1 = atan2(sin_alp0 * sin_sig1, cos_sig1);
0262         CT const lon2 = atan2(sin_alp0 * sin_sig2, cos_sig2);
0263
0264         result.lon2 = lon1 + lon2 - omg1;
0265         result.lat2 = lat2;
0266
0267         // For longitudes close to the antimeridian the result can be out
0268         // of range. Therefore normalize.
0269         math::detail::normalize_angle_cond<radian>(result.lon2);
0270     }
0271
0272     if (ReverseAzimuth)
0273     {
0274         CT const alp2 = atan2(sin_alp0, cos_alp0 * cos_sig2);
0275         result.reverse_azimuth = alp2;
0276     }
0277
0278     return result;
0279 }
0280
0281 } // namespace formula
0282
0283 }} // namespace boost::geometry
0284
0285 #endif // BOOST_GEOMETRY_FORMULAS_SPHERICAL_HPP