Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:35:26 

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_SJOBERG_INTERSECTION_HPP
0013 #define BOOST_GEOMETRY_FORMULAS_SJOBERG_INTERSECTION_HPP
0014 
0015 
0016 #include <boost/math/constants/constants.hpp>
0017 
0018 #include <boost/geometry/core/radius.hpp>
0019 
0020 #include <boost/geometry/util/condition.hpp>
0021 #include <boost/geometry/util/math.hpp>
0022 #include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
0023 
0024 #include <boost/geometry/formulas/flattening.hpp>
0025 #include <boost/geometry/formulas/spherical.hpp>
0026 
0027 
0028 namespace boost { namespace geometry { namespace formula
0029 {
0030 
0031 /*!
0032 \brief The intersection of two great circles as proposed by Sjoberg.
0033 \see See
0034     - [Sjoberg02] Lars E. Sjoberg, Intersections on the sphere and ellipsoid, 2002
0035       http://link.springer.com/article/10.1007/s00190-001-0230-9
0036 */
0037 template <typename CT>
0038 struct sjoberg_intersection_spherical_02
0039 {
0040     // TODO: if it will be used as standalone formula
0041     //       support segments on equator and endpoints on poles
0042 
0043     static inline bool apply(CT const& lon1, CT const& lat1, CT const& lon_a2, CT const& lat_a2,
0044                              CT const& lon2, CT const& lat2, CT const& lon_b2, CT const& lat_b2,
0045                              CT & lon, CT & lat)
0046     {
0047         CT tan_lat = 0;
0048         bool res = apply_alt(lon1, lat1, lon_a2, lat_a2,
0049                              lon2, lat2, lon_b2, lat_b2,
0050                              lon, tan_lat);
0051 
0052         if (res)
0053         {
0054             lat = atan(tan_lat);
0055         }
0056 
0057         return res;
0058     }
0059 
0060     static inline bool apply_alt(CT const& lon1, CT const& lat1, CT const& lon_a2, CT const& lat_a2,
0061                                  CT const& lon2, CT const& lat2, CT const& lon_b2, CT const& lat_b2,
0062                                  CT & lon, CT & tan_lat)
0063     {
0064         CT const cos_lon1 = cos(lon1);
0065         CT const sin_lon1 = sin(lon1);
0066         CT const cos_lon2 = cos(lon2);
0067         CT const sin_lon2 = sin(lon2);
0068         CT const sin_lat1 = sin(lat1);
0069         CT const sin_lat2 = sin(lat2);
0070         CT const cos_lat1 = cos(lat1);
0071         CT const cos_lat2 = cos(lat2);
0072 
0073         CT const tan_lat_a2 = tan(lat_a2);
0074         CT const tan_lat_b2 = tan(lat_b2);
0075 
0076         return apply(lon1, lon_a2, lon2, lon_b2,
0077                      sin_lon1, cos_lon1, sin_lat1, cos_lat1,
0078                      sin_lon2, cos_lon2, sin_lat2, cos_lat2,
0079                      tan_lat_a2, tan_lat_b2,
0080                      lon, tan_lat);
0081     }
0082 
0083 private:
0084     static inline bool apply(CT const& lon1, CT const& lon_a2, CT const& lon2, CT const& lon_b2,
0085                              CT const& sin_lon1, CT const& cos_lon1, CT const& sin_lat1, CT const& cos_lat1,
0086                              CT const& sin_lon2, CT const& cos_lon2, CT const& sin_lat2, CT const& cos_lat2,
0087                              CT const& tan_lat_a2, CT const& tan_lat_b2,
0088                              CT & lon, CT & tan_lat)
0089     {
0090         // NOTE:
0091         // cos_lat_ = 0 <=> segment on equator
0092         // tan_alpha_ = 0 <=> segment vertical
0093 
0094         CT const tan_lat1 = sin_lat1 / cos_lat1; //tan(lat1);
0095         CT const tan_lat2 = sin_lat2 / cos_lat2; //tan(lat2);
0096 
0097         CT const dlon1 = lon_a2 - lon1;
0098         CT const sin_dlon1 = sin(dlon1);
0099         CT const dlon2 = lon_b2 - lon2;
0100         CT const sin_dlon2 = sin(dlon2);
0101 
0102         CT const cos_dlon1 = cos(dlon1);
0103         CT const cos_dlon2 = cos(dlon2);
0104 
0105         CT const tan_alpha1_x = cos_lat1 * tan_lat_a2 - sin_lat1 * cos_dlon1;
0106         CT const tan_alpha2_x = cos_lat2 * tan_lat_b2 - sin_lat2 * cos_dlon2;
0107 
0108         CT const c0 = 0;
0109         bool const is_vertical1 = math::equals(sin_dlon1, c0) || math::equals(tan_alpha1_x, c0);
0110         bool const is_vertical2 = math::equals(sin_dlon2, c0) || math::equals(tan_alpha2_x, c0);
0111 
0112         CT tan_alpha1 = 0;
0113         CT tan_alpha2 = 0;
0114 
0115         if (is_vertical1 && is_vertical2)
0116         {
0117             // circles intersect at one of the poles or are collinear
0118             return false;
0119         }
0120         else if (is_vertical1)
0121         {
0122             tan_alpha2 = sin_dlon2 / tan_alpha2_x;
0123 
0124             lon = lon1;
0125         }
0126         else if (is_vertical2)
0127         {
0128             tan_alpha1 = sin_dlon1 / tan_alpha1_x;
0129 
0130             lon = lon2;
0131         }
0132         else
0133         {
0134             tan_alpha1 = sin_dlon1 / tan_alpha1_x;
0135             tan_alpha2 = sin_dlon2 / tan_alpha2_x;
0136 
0137             CT const T1 = tan_alpha1 * cos_lat1;
0138             CT const T2 = tan_alpha2 * cos_lat2;
0139             CT const T1T2 = T1*T2;
0140             CT const tan_lon_y = T1 * sin_lon2 - T2 * sin_lon1 + T1T2 * (tan_lat1 * cos_lon1 - tan_lat2 * cos_lon2);
0141             CT const tan_lon_x = T1 * cos_lon2 - T2 * cos_lon1 - T1T2 * (tan_lat1 * sin_lon1 - tan_lat2 * sin_lon2);
0142 
0143             lon = atan2(tan_lon_y, tan_lon_x);
0144         }
0145 
0146         // choose closer result
0147         CT const pi = math::pi<CT>();
0148         CT const lon_2 = lon > c0 ? lon - pi : lon + pi;
0149         CT const lon_dist1 = (std::max)((std::min)(math::longitude_difference<radian>(lon1, lon),
0150                                                    math::longitude_difference<radian>(lon_a2, lon)),
0151                                         (std::min)(math::longitude_difference<radian>(lon2, lon),
0152                                                    math::longitude_difference<radian>(lon_b2, lon)));
0153         CT const lon_dist2 = (std::max)((std::min)(math::longitude_difference<radian>(lon1, lon_2),
0154                                                    math::longitude_difference<radian>(lon_a2, lon_2)),
0155                                         (std::min)(math::longitude_difference<radian>(lon2, lon_2),
0156                                                    math::longitude_difference<radian>(lon_b2, lon_2)));
0157         if (lon_dist2 < lon_dist1)
0158         {
0159             lon = lon_2;
0160         }
0161 
0162         CT const sin_lon = sin(lon);
0163         CT const cos_lon = cos(lon);
0164 
0165         if (math::abs(tan_alpha1) >= math::abs(tan_alpha2)) // pick less vertical segment
0166         {
0167             CT const sin_dlon_1 = sin_lon * cos_lon1 - cos_lon * sin_lon1;
0168             CT const cos_dlon_1 = cos_lon * cos_lon1 + sin_lon * sin_lon1;
0169             CT const lat_y_1 = sin_dlon_1 + tan_alpha1 * sin_lat1 * cos_dlon_1;
0170             CT const lat_x_1 = tan_alpha1 * cos_lat1;
0171             tan_lat = lat_y_1 / lat_x_1;
0172         }
0173         else
0174         {
0175             CT const sin_dlon_2 = sin_lon * cos_lon2 - cos_lon * sin_lon2;
0176             CT const cos_dlon_2 = cos_lon * cos_lon2 + sin_lon * sin_lon2;
0177             CT const lat_y_2 = sin_dlon_2 + tan_alpha2 * sin_lat2 * cos_dlon_2;
0178             CT const lat_x_2 = tan_alpha2 * cos_lat2;
0179             tan_lat = lat_y_2 / lat_x_2;
0180         }
0181 
0182         return true;
0183     }
0184 };
0185 
0186 
0187 /*! Approximation of dLambda_j [Sjoberg07], expanded into taylor series in e^2
0188     Maxima script:
0189     dLI_j(c_j, sinB_j, sinB) := integrate(1 / (sqrt(1 - c_j ^ 2 - x ^ 2)*(1 + sqrt(1 - e2*(1 - x ^ 2)))), x, sinB_j, sinB);
0190     dL_j(c_j, B_j, B) := -e2 * c_j * dLI_j(c_j, B_j, B);
0191     S: taylor(dLI_j(c_j, sinB_j, sinB), e2, 0, 3);
0192     assume(c_j < 1);
0193     assume(c_j > 0);
0194     L1: factor(integrate(sqrt(-x ^ 2 - c_j ^ 2 + 1) / (x ^ 2 + c_j ^ 2 - 1), x));
0195     L2: factor(integrate(((x ^ 2 - 1)*sqrt(-x ^ 2 - c_j ^ 2 + 1)) / (x ^ 2 + c_j ^ 2 - 1), x));
0196     L3: factor(integrate(((x ^ 4 - 2 * x ^ 2 + 1)*sqrt(-x ^ 2 - c_j ^ 2 + 1)) / (x ^ 2 + c_j ^ 2 - 1), x));
0197     L4: factor(integrate(((x ^ 6 - 3 * x ^ 4 + 3 * x ^ 2 - 1)*sqrt(-x ^ 2 - c_j ^ 2 + 1)) / (x ^ 2 + c_j ^ 2 - 1), x));
0198 
0199 \see See
0200     - [Sjoberg07] Lars E. Sjoberg, Geodetic intersection on the ellipsoid, 2007
0201       http://link.springer.com/article/10.1007/s00190-007-0204-7
0202 */
0203 template <unsigned int Order, typename CT>
0204 inline CT sjoberg_d_lambda_e_sqr(CT const& sin_betaj, CT const& sin_beta,
0205                                  CT const& Cj, CT const& sqrt_1_Cj_sqr,
0206                                  CT const& e_sqr)
0207 {
0208     using math::detail::bounded;
0209 
0210     if (BOOST_GEOMETRY_CONDITION(Order == 0))
0211     {
0212         return 0;
0213     }
0214 
0215     CT const c1 = 1;
0216     CT const c2 = 2;
0217 
0218     CT const asin_B = asin(bounded(sin_beta / sqrt_1_Cj_sqr, -c1, c1));
0219     CT const asin_Bj = asin(sin_betaj / sqrt_1_Cj_sqr);
0220     CT const L0 = (asin_B - asin_Bj) / c2;
0221 
0222     if (BOOST_GEOMETRY_CONDITION(Order == 1))
0223     {
0224         return -Cj * e_sqr * L0;
0225     }
0226 
0227     CT const c0 = 0;
0228     CT const c16 = 16;
0229 
0230     CT const X = sin_beta;
0231     CT const Xj = sin_betaj;
0232     CT const X_sqr = math::sqr(X);
0233     CT const Xj_sqr = math::sqr(Xj);
0234     CT const Cj_sqr = math::sqr(Cj);
0235     CT const Cj_sqr_plus_one = Cj_sqr + c1;
0236     CT const one_minus_Cj_sqr = c1 - Cj_sqr;
0237     CT const sqrt_Y = math::sqrt(bounded(-X_sqr + one_minus_Cj_sqr, c0));
0238     CT const sqrt_Yj = math::sqrt(-Xj_sqr + one_minus_Cj_sqr);
0239     CT const L1 = (Cj_sqr_plus_one * (asin_B - asin_Bj) + X * sqrt_Y - Xj * sqrt_Yj) / c16;
0240 
0241     if (BOOST_GEOMETRY_CONDITION(Order == 2))
0242     {
0243         return -Cj * e_sqr * (L0 + e_sqr * L1);
0244     }
0245 
0246     CT const c3 = 3;
0247     CT const c5 = 5;
0248     CT const c128 = 128;
0249 
0250     CT const E = Cj_sqr * (c3 * Cj_sqr + c2) + c3;
0251     CT const F = X * (-c2 * X_sqr + c3 * Cj_sqr + c5);
0252     CT const Fj = Xj * (-c2 * Xj_sqr + c3 * Cj_sqr + c5);
0253     CT const L2 = (E * (asin_B - asin_Bj) + F * sqrt_Y - Fj * sqrt_Yj) / c128;
0254 
0255     if (BOOST_GEOMETRY_CONDITION(Order == 3))
0256     {
0257         return -Cj * e_sqr * (L0 + e_sqr * (L1 + e_sqr * L2));
0258     }
0259 
0260     CT const c8 = 8;
0261     CT const c9 = 9;
0262     CT const c10 = 10;
0263     CT const c15 = 15;
0264     CT const c24 = 24;
0265     CT const c26 = 26;
0266     CT const c33 = 33;
0267     CT const c6144 = 6144;
0268 
0269     CT const G = Cj_sqr * (Cj_sqr * (Cj_sqr * c15 + c9) + c9) + c15;
0270     CT const H = -c10 * Cj_sqr - c26;
0271     CT const I = Cj_sqr * (Cj_sqr * c15 + c24) + c33;
0272     CT const J = X_sqr * (X * (c8 * X_sqr + H)) + X * I;
0273     CT const Jj = Xj_sqr * (Xj * (c8 * Xj_sqr + H)) + Xj * I;
0274     CT const L3 = (G * (asin_B - asin_Bj) + J * sqrt_Y - Jj * sqrt_Yj) / c6144;
0275 
0276     // Order 4 and higher
0277     return -Cj * e_sqr * (L0 + e_sqr * (L1 + e_sqr * (L2 + e_sqr * L3)));
0278 }
0279 
0280 /*!
0281 \brief The representation of geodesic as proposed by Sjoberg.
0282 \see See
0283     - [Sjoberg07] Lars E. Sjoberg, Geodetic intersection on the ellipsoid, 2007
0284       http://link.springer.com/article/10.1007/s00190-007-0204-7
0285     - [Sjoberg12] Lars E. Sjoberg, Solutions to the ellipsoidal Clairaut constant
0286       and the inverse geodetic problem by numerical integration, 2012
0287       https://www.degruyter.com/view/j/jogs.2012.2.issue-3/v10156-011-0037-4/v10156-011-0037-4.xml
0288 */
0289 template <typename CT, unsigned int Order>
0290 class sjoberg_geodesic
0291 {
0292     sjoberg_geodesic() {}
0293 
0294     static int sign_C(CT const& alphaj)
0295     {
0296         CT const c0 = 0;
0297         CT const c2 = 2;
0298         CT const pi = math::pi<CT>();
0299         CT const pi_half = pi / c2;
0300 
0301         return (pi_half < alphaj && alphaj < pi) || (-pi_half < alphaj && alphaj < c0) ? -1 : 1;
0302     }
0303 
0304 public:
0305     sjoberg_geodesic(CT const& lon, CT const& lat, CT const& alpha, CT const& f)
0306         : lonj(lon)
0307         , latj(lat)
0308         , alphaj(alpha)
0309     {
0310         CT const c0 = 0;
0311         CT const c1 = 1;
0312         CT const c2 = 2;
0313         //CT const pi = math::pi<CT>();
0314         //CT const pi_half = pi / c2;
0315 
0316         one_minus_f = c1 - f;
0317         e_sqr = f * (c2 - f);
0318 
0319         tan_latj = tan(lat);
0320         tan_betaj = one_minus_f * tan_latj;
0321         betaj = atan(tan_betaj);
0322         sin_betaj = sin(betaj);
0323 
0324         cos_betaj = cos(betaj);
0325         sin_alphaj = sin(alphaj);
0326         // Clairaut constant (lower-case in the paper)
0327         Cj = sign_C(alphaj) * cos_betaj * sin_alphaj;
0328         Cj_sqr = math::sqr(Cj);
0329         sqrt_1_Cj_sqr = math::sqrt(c1 - Cj_sqr);
0330 
0331         sign_lon_diff = alphaj >= 0 ? 1 : -1; // || alphaj == -pi ?
0332         //sign_lon_diff = 1;
0333 
0334         is_on_equator = math::equals(sqrt_1_Cj_sqr, c0);
0335         is_Cj_zero = math::equals(Cj, c0);
0336 
0337         t0j = c0;
0338         asin_tj_t0j = c0;
0339 
0340         if (! is_Cj_zero)
0341         {
0342             t0j = sqrt_1_Cj_sqr / Cj;
0343         }
0344 
0345         if (! is_on_equator)
0346         {
0347             //asin_tj_t0j = asin(tan_betaj / t0j);
0348             asin_tj_t0j = asin(tan_betaj * Cj / sqrt_1_Cj_sqr);
0349         }
0350     }
0351 
0352     struct vertex_data
0353     {
0354         //CT beta0j;
0355         CT sin_beta0j;
0356         CT dL0j;
0357         CT lon0j;
0358     };
0359 
0360     vertex_data get_vertex_data() const
0361     {
0362         CT const c2 = 2;
0363         CT const pi = math::pi<CT>();
0364         CT const pi_half = pi / c2;
0365 
0366         vertex_data res;
0367 
0368         if (! is_Cj_zero)
0369         {
0370             //res.beta0j = atan(t0j);
0371             //res.sin_beta0j = sin(res.beta0j);
0372             res.sin_beta0j = math::sign(t0j) * sqrt_1_Cj_sqr;
0373             res.dL0j = d_lambda(res.sin_beta0j);
0374             res.lon0j = lonj + sign_lon_diff * (pi_half - asin_tj_t0j + res.dL0j);
0375         }
0376         else
0377         {
0378             //res.beta0j = pi_half;
0379             //res.sin_beta0j = betaj >= 0 ? 1 : -1;
0380             res.sin_beta0j = 1;
0381             res.dL0j = 0;
0382             res.lon0j = lonj;
0383         }
0384 
0385         return res;
0386     }
0387 
0388     bool is_sin_beta_ok(CT const& sin_beta) const
0389     {
0390         CT const c1 = 1;
0391         return math::abs(sin_beta / sqrt_1_Cj_sqr) <= c1;
0392     }
0393 
0394     bool k_diff(CT const& sin_beta,
0395                 CT & delta_k) const
0396     {
0397         if (is_Cj_zero)
0398         {
0399             delta_k = 0;
0400             return true;
0401         }
0402 
0403         // beta out of bounds and not close
0404         if (! (is_sin_beta_ok(sin_beta)
0405                 || math::equals(math::abs(sin_beta), sqrt_1_Cj_sqr)) )
0406         {
0407             return false;
0408         }
0409 
0410         // NOTE: beta may be slightly out of bounds here but d_lambda handles that
0411         CT const dLj = d_lambda(sin_beta);
0412         delta_k = sign_lon_diff * (/*asin_t_t0j*/ - asin_tj_t0j + dLj);
0413 
0414         return true;
0415     }
0416 
0417     bool lon_diff(CT const& sin_beta, CT const& t,
0418                   CT & delta_lon) const
0419     {
0420         using math::detail::bounded;
0421         CT const c1 = 1;
0422 
0423         if (is_Cj_zero)
0424         {
0425             delta_lon = 0;
0426             return true;
0427         }
0428 
0429         CT delta_k = 0;
0430         if (! k_diff(sin_beta, delta_k))
0431         {
0432             return false;
0433         }
0434 
0435         CT const t_t0j = t / t0j;
0436         // NOTE: t may be slightly out of bounds here
0437         CT const asin_t_t0j = asin(bounded(t_t0j, -c1, c1));
0438         delta_lon = sign_lon_diff * asin_t_t0j + delta_k;
0439 
0440         return true;
0441     }
0442 
0443     bool k_diffs(CT const& sin_beta, vertex_data const& vd,
0444                  CT & delta_k_before, CT & delta_k_behind,
0445                  bool check_sin_beta = true) const
0446     {
0447         CT const pi = math::pi<CT>();
0448 
0449         if (is_Cj_zero)
0450         {
0451             delta_k_before = 0;
0452             delta_k_behind = sign_lon_diff * pi;
0453             return true;
0454         }
0455 
0456         // beta out of bounds and not close
0457         if (check_sin_beta
0458             && ! (is_sin_beta_ok(sin_beta)
0459                     || math::equals(math::abs(sin_beta), sqrt_1_Cj_sqr)) )
0460         {
0461             return false;
0462         }
0463 
0464         // NOTE: beta may be slightly out of bounds here but d_lambda handles that
0465         CT const dLj = d_lambda(sin_beta);
0466         delta_k_before = sign_lon_diff * (/*asin_t_t0j*/ - asin_tj_t0j + dLj);
0467 
0468         // This version require no additional dLj calculation
0469         delta_k_behind = sign_lon_diff * (pi /*- asin_t_t0j*/ - asin_tj_t0j + vd.dL0j + (vd.dL0j - dLj));
0470 
0471         // [Sjoberg12]
0472         //CT const dL101 = d_lambda(sin_betaj, vd.sin_beta0j);
0473         // WARNING: the following call might not work if beta was OoB because only the second argument is bounded
0474         //CT const dL_01 = d_lambda(sin_beta, vd.sin_beta0j);
0475         //delta_k_behind = sign_lon_diff * (pi /*- asin_t_t0j*/ - asin_tj_t0j + dL101 + dL_01);
0476 
0477         return true;
0478     }
0479 
0480     bool lon_diffs(CT const& sin_beta, CT const& t, vertex_data const& vd,
0481                    CT & delta_lon_before, CT & delta_lon_behind) const
0482     {
0483         using math::detail::bounded;
0484         CT const c1 = 1;
0485         CT const pi = math::pi<CT>();
0486 
0487         if (is_Cj_zero)
0488         {
0489             delta_lon_before = 0;
0490             delta_lon_behind = sign_lon_diff * pi;
0491             return true;
0492         }
0493 
0494         CT delta_k_before = 0, delta_k_behind = 0;
0495         if (! k_diffs(sin_beta, vd, delta_k_before, delta_k_behind))
0496         {
0497             return false;
0498         }
0499 
0500         CT const t_t0j = t / t0j;
0501         // NOTE: t may be slightly out of bounds here
0502         CT const asin_t_t0j = asin(bounded(t_t0j, -c1, c1));
0503         CT const sign_asin_t_t0j = sign_lon_diff * asin_t_t0j;
0504         delta_lon_before = sign_asin_t_t0j + delta_k_before;
0505         delta_lon_behind = -sign_asin_t_t0j + delta_k_behind;
0506 
0507         return true;
0508     }
0509 
0510     bool lon(CT const& sin_beta, CT const& t, vertex_data const& vd,
0511              CT & lon_before, CT & lon_behind) const
0512     {
0513         using math::detail::bounded;
0514         CT const c1 = 1;
0515         CT const pi = math::pi<CT>();
0516 
0517         if (is_Cj_zero)
0518         {
0519             lon_before = lonj;
0520             lon_behind = lonj + sign_lon_diff * pi;
0521             return true;
0522         }
0523 
0524         if (! (is_sin_beta_ok(sin_beta)
0525                 || math::equals(math::abs(sin_beta), sqrt_1_Cj_sqr)) )
0526         {
0527             return false;
0528         }
0529 
0530         CT const t_t0j = t / t0j;
0531         CT const asin_t_t0j = asin(bounded(t_t0j, -c1, c1));
0532         CT const dLj = d_lambda(sin_beta);
0533         lon_before = lonj + sign_lon_diff * (asin_t_t0j - asin_tj_t0j + dLj);
0534         lon_behind = vd.lon0j + (vd.lon0j - lon_before);
0535 
0536         return true;
0537     }
0538 
0539 
0540     CT lon(CT const& delta_lon) const
0541     {
0542         return lonj + delta_lon;
0543     }
0544 
0545     CT lat(CT const& t) const
0546     {
0547         // t = tan(beta) = (1-f)tan(lat)
0548         return atan(t / one_minus_f);
0549     }
0550 
0551     void vertex(CT & lon, CT & lat) const
0552     {
0553         lon = get_vertex_data().lon0j;
0554         if (! is_Cj_zero)
0555         {
0556             lat = sjoberg_geodesic::lat(t0j);
0557         }
0558         else
0559         {
0560             CT const c2 = 2;
0561             lat = math::pi<CT>() / c2;
0562         }
0563     }
0564 
0565     CT lon_of_equator_intersection() const
0566     {
0567         CT const c0 = 0;
0568         CT const dLj = d_lambda(c0);
0569         CT const asin_tj_t0j = asin(Cj * tan_betaj / sqrt_1_Cj_sqr);
0570         return lonj - asin_tj_t0j + dLj;
0571     }
0572 
0573     CT d_lambda(CT const& sin_beta) const
0574     {
0575         return sjoberg_d_lambda_e_sqr<Order>(sin_betaj, sin_beta, Cj, sqrt_1_Cj_sqr, e_sqr);
0576     }
0577 
0578     // [Sjoberg12]
0579     /*CT d_lambda(CT const& sin_beta1, CT const& sin_beta2) const
0580     {
0581         return sjoberg_d_lambda_e_sqr<Order>(sin_beta1, sin_beta2, Cj, sqrt_1_Cj_sqr, e_sqr);
0582     }*/
0583 
0584     CT lonj;
0585     CT latj;
0586     CT alphaj;
0587 
0588     CT one_minus_f;
0589     CT e_sqr;
0590 
0591     CT tan_latj;
0592     CT tan_betaj;
0593     CT betaj;
0594     CT sin_betaj;
0595     CT cos_betaj;
0596     CT sin_alphaj;
0597     CT Cj;
0598     CT Cj_sqr;
0599     CT sqrt_1_Cj_sqr;
0600 
0601     int sign_lon_diff;
0602 
0603     bool is_on_equator;
0604     bool is_Cj_zero;
0605 
0606     CT t0j;
0607     CT asin_tj_t0j;
0608 };
0609 
0610 
0611 /*!
0612 \brief The intersection of two geodesics as proposed by Sjoberg.
0613 \see See
0614     - [Sjoberg02] Lars E. Sjoberg, Intersections on the sphere and ellipsoid, 2002
0615       http://link.springer.com/article/10.1007/s00190-001-0230-9
0616     - [Sjoberg07] Lars E. Sjoberg, Geodetic intersection on the ellipsoid, 2007
0617       http://link.springer.com/article/10.1007/s00190-007-0204-7
0618     - [Sjoberg12] Lars E. Sjoberg, Solutions to the ellipsoidal Clairaut constant
0619       and the inverse geodetic problem by numerical integration, 2012
0620       https://www.degruyter.com/view/j/jogs.2012.2.issue-3/v10156-011-0037-4/v10156-011-0037-4.xml
0621 */
0622 template
0623 <
0624     typename CT,
0625     template <typename, bool, bool, bool, bool, bool> class Inverse,
0626     unsigned int Order = 4
0627 >
0628 class sjoberg_intersection
0629 {
0630     typedef sjoberg_geodesic<CT, Order> geodesic_type;
0631     typedef Inverse<CT, false, true, false, false, false> inverse_type;
0632     typedef typename inverse_type::result_type inverse_result;
0633 
0634     static bool const enable_02 = true;
0635     static int const max_iterations_02 = 10;
0636     static int const max_iterations_07 = 20;
0637 
0638 public:
0639     template <typename T1, typename T2, typename Spheroid>
0640     static inline bool apply(T1 const& lona1, T1 const& lata1,
0641                              T1 const& lona2, T1 const& lata2,
0642                              T2 const& lonb1, T2 const& latb1,
0643                              T2 const& lonb2, T2 const& latb2,
0644                              CT & lon, CT & lat,
0645                              Spheroid const& spheroid)
0646     {
0647         CT const lon_a1 = lona1;
0648         CT const lat_a1 = lata1;
0649         CT const lon_a2 = lona2;
0650         CT const lat_a2 = lata2;
0651         CT const lon_b1 = lonb1;
0652         CT const lat_b1 = latb1;
0653         CT const lon_b2 = lonb2;
0654         CT const lat_b2 = latb2;
0655 
0656         inverse_result const res1 = inverse_type::apply(lon_a1, lat_a1, lon_a2, lat_a2, spheroid);
0657         inverse_result const res2 = inverse_type::apply(lon_b1, lat_b1, lon_b2, lat_b2, spheroid);
0658 
0659         return apply(lon_a1, lat_a1, lon_a2, lat_a2, res1.azimuth,
0660                      lon_b1, lat_b1, lon_b2, lat_b2, res2.azimuth,
0661                      lon, lat, spheroid);
0662     }
0663 
0664     // TODO: Currently may not work correctly if one of the endpoints is the pole
0665     template <typename Spheroid>
0666     static inline bool apply(CT const& lon_a1, CT const& lat_a1, CT const& lon_a2, CT const& lat_a2, CT const& alpha_a1,
0667                              CT const& lon_b1, CT const& lat_b1, CT const& lon_b2, CT const& lat_b2, CT const& alpha_b1,
0668                              CT & lon, CT & lat,
0669                              Spheroid const& spheroid)
0670     {
0671         // coordinates in radians
0672 
0673         CT const c0 = 0;
0674         CT const c1 = 1;
0675 
0676         CT const f = formula::flattening<CT>(spheroid);
0677         CT const one_minus_f = c1 - f;
0678 
0679         geodesic_type geod1(lon_a1, lat_a1, alpha_a1, f);
0680         geodesic_type geod2(lon_b1, lat_b1, alpha_b1, f);
0681 
0682         // Cj = 1 if on equator <=> sqrt_1_Cj_sqr = 0
0683         // Cj = 0 if vertical <=> sqrt_1_Cj_sqr = 1
0684 
0685         if (geod1.is_on_equator && geod2.is_on_equator)
0686         {
0687             return false;
0688         }
0689         else if (geod1.is_on_equator)
0690         {
0691             lon = geod2.lon_of_equator_intersection();
0692             lat = c0;
0693             return true;
0694         }
0695         else if (geod2.is_on_equator)
0696         {
0697             lon = geod1.lon_of_equator_intersection();
0698             lat = c0;
0699             return true;
0700         }
0701 
0702         // (lon1 - lon2) normalized to (-180, 180]
0703         CT const lon1_minus_lon2 = math::longitude_distance_signed<radian>(geod2.lonj, geod1.lonj);
0704 
0705         // vertical segments
0706         if (geod1.is_Cj_zero && geod2.is_Cj_zero)
0707         {
0708             CT const pi = math::pi<CT>();
0709 
0710             // the geodesics are parallel, the intersection point cannot be calculated
0711             if ( math::equals(lon1_minus_lon2, c0)
0712               || math::equals(lon1_minus_lon2 + (lon1_minus_lon2 < c0 ? pi : -pi), c0) )
0713             {
0714                 return false;
0715             }
0716 
0717             lon = c0;
0718 
0719             // the geodesics intersect at one of the poles
0720             CT const pi_half = pi / CT(2);
0721             CT const abs_lat_a1 = math::abs(lat_a1);
0722             CT const abs_lat_a2 = math::abs(lat_a2);
0723             if (math::equals(abs_lat_a1, abs_lat_a2))
0724             {
0725                 lat = pi_half;
0726             }
0727             else
0728             {
0729                 // pick the pole closest to one of the points of the first segment
0730                 CT const& closer_lat = abs_lat_a1 > abs_lat_a2 ? lat_a1 : lat_a2;
0731                 lat = closer_lat >= 0 ? pi_half : -pi_half;
0732             }
0733 
0734             return true;
0735         }
0736 
0737         CT lon_sph = 0;
0738 
0739         // Starting tan(beta)
0740         CT t = 0;
0741 
0742         /*if (geod1.is_Cj_zero)
0743         {
0744             CT const k_base = lon1_minus_lon2 + geod2.sign_lon_diff * geod2.asin_tj_t0j;
0745             t = sin(k_base) * geod2.t0j;
0746             lon_sph = vertical_intersection_longitude(geod1.lonj, lon_b1, lon_b2);
0747         }
0748         else if (geod2.is_Cj_zero)
0749         {
0750             CT const k_base = lon1_minus_lon2 - geod1.sign_lon_diff * geod1.asin_tj_t0j;
0751             t = sin(-k_base) * geod1.t0j;
0752             lon_sph = vertical_intersection_longitude(geod2.lonj, lon_a1, lon_a2);
0753         }
0754         else*/
0755         {
0756             // TODO: Consider using betas instead of latitudes.
0757             //       Some function calls might be saved this way.
0758             CT tan_lat_sph = 0;
0759             sjoberg_intersection_spherical_02<CT>::apply_alt(lon_a1, lat_a1, lon_a2, lat_a2,
0760                                                              lon_b1, lat_b1, lon_b2, lat_b2,
0761                                                              lon_sph, tan_lat_sph);
0762 
0763             // Return for sphere
0764             if (math::equals(f, c0))
0765             {
0766                 lon = lon_sph;
0767                 lat = atan(tan_lat_sph);
0768                 return true;
0769             }
0770 
0771             t = one_minus_f * tan_lat_sph; // tan(beta)
0772         }
0773 
0774         // TODO: no need to calculate atan here if reduced latitudes were used
0775         //       instead of latitudes above, in sjoberg_intersection_spherical_02
0776         CT const beta = atan(t);
0777 
0778         if (enable_02 && newton_method(geod1, geod2, beta, t, lon1_minus_lon2, lon_sph, lon, lat))
0779         {
0780             // TODO: Newton's method may return wrong result in some specific cases
0781             // Detected for sphere and nearly sphere, e.g. A=6371228, B=6371227
0782             // and segments s1=(-121 -19,37 8) and s2=(-19 -15,-104 -58)
0783             // It's unclear if this is a bug or a characteristic of this method
0784             // so until this is investigated check if the resulting longitude is
0785             // between endpoints of the segments. It should be since before calling
0786             // this formula sides of endpoints WRT other segments are checked.
0787             if ( is_result_longitude_ok(geod1, lon_a1, lon_a2, lon)
0788               && is_result_longitude_ok(geod2, lon_b1, lon_b2, lon) )
0789             {
0790                 return true;
0791             }
0792         }
0793 
0794         return converge_07(geod1, geod2, beta, t, lon1_minus_lon2, lon_sph, lon, lat);
0795     }
0796 
0797 private:
0798     static inline bool newton_method(geodesic_type const& geod1, geodesic_type const& geod2, // in
0799                                      CT beta, CT t, CT const& lon1_minus_lon2, CT const& lon_sph, // in
0800                                      CT & lon, CT & lat) // out
0801     {
0802         CT const c0 = 0;
0803         CT const c1 = 1;
0804 
0805         CT const e_sqr = geod1.e_sqr;
0806 
0807         CT lon1_diff = 0;
0808         CT lon2_diff = 0;
0809 
0810         // The segment is vertical and intersection point is behind the vertex
0811         // this method is unable to calculate correct result
0812         if (geod1.is_Cj_zero && math::abs(geod1.lonj - lon_sph) > math::half_pi<CT>())
0813             return false;
0814         if (geod2.is_Cj_zero && math::abs(geod2.lonj - lon_sph) > math::half_pi<CT>())
0815             return false;
0816 
0817         CT abs_dbeta_last = 0;
0818 
0819         // [Sjoberg02] converges faster than solution in [Sjoberg07]
0820         // Newton-Raphson method
0821         for (int i = 0; i < max_iterations_02; ++i)
0822         {
0823             CT const cos_beta = cos(beta);
0824 
0825             if (math::equals(cos_beta, c0))
0826             {
0827                 return false;
0828             }
0829 
0830             CT const sin_beta = sin(beta);
0831             CT const cos_beta_sqr = math::sqr(cos_beta);
0832             CT const G = c1 - e_sqr * cos_beta_sqr;
0833 
0834             CT f1 = 0;
0835             CT f2 = 0;
0836 
0837             if (!geod1.is_Cj_zero)
0838             {
0839                 bool is_beta_ok = geod1.lon_diff(sin_beta, t, lon1_diff);
0840 
0841                 if (is_beta_ok)
0842                 {
0843                     CT const H = cos_beta_sqr - geod1.Cj_sqr;
0844                     if (math::equals(H, c0))
0845                     {
0846                         return false;
0847                     }
0848                     f1 = geod1.Cj / cos_beta * math::sqrt(G / H);
0849                 }
0850                 else
0851                 {
0852                     return false;
0853                 }
0854             }
0855 
0856             if (!geod2.is_Cj_zero)
0857             {
0858                 bool is_beta_ok = geod2.lon_diff(sin_beta, t, lon2_diff);
0859 
0860                 if (is_beta_ok)
0861                 {
0862                     CT const H = cos_beta_sqr - geod2.Cj_sqr;
0863                     if (math::equals(H, c0))
0864                     {
0865                         // NOTE: This may happen for segment nearly
0866                         // at the equator. Detected for (radian):
0867                         // (-0.0872665 -0.0872665, -0.0872665 0.0872665)
0868                         // x
0869                         // (0 1.57e-07, -0.392699 1.57e-07)
0870                         return false;
0871                     }
0872                     f2 = geod2.Cj / cos_beta * math::sqrt(G / H);
0873                 }
0874                 else
0875                 {
0876                     return false;
0877                 }
0878             }
0879 
0880             // NOTE: Things may go wrong if the IP is near the vertex
0881             //   1. May converge into the wrong direction (from the other way around).
0882             //      This happens when the starting point is on the other side than the vertex
0883             //   2. During converging may "jump" into the other side of the vertex.
0884             //      In this case sin_beta/sqrt_1_Cj_sqr and t/t0j is not in [-1, 1]
0885             //   3. f1-f2 may be 0 which means that the intermediate point is on the vertex
0886             //      In this case it's not possible to check if this is the correct result
0887             //   4. f1-f2 may also be 0 in other cases, e.g.
0888             //      geodesics are symmetrical wrt equator and longitude directions are different
0889 
0890             CT const dbeta_denom = f1 - f2;
0891             //CT const dbeta_denom = math::abs(f1) + math::abs(f2);
0892 
0893             if (math::equals(dbeta_denom, c0))
0894             {
0895                 return false;
0896             }
0897 
0898             // The sign of dbeta is changed WRT [Sjoberg02]
0899             CT const dbeta = (lon1_minus_lon2 + lon1_diff - lon2_diff) / dbeta_denom;
0900 
0901             CT const abs_dbeta = math::abs(dbeta);
0902             if (i > 0 && abs_dbeta > abs_dbeta_last)
0903             {
0904                 // The algorithm is not converging
0905                 // The intersection may be on the other side of the vertex
0906                 return false;
0907             }
0908             abs_dbeta_last = abs_dbeta;
0909 
0910             if (math::equals(dbeta, c0))
0911             {
0912                 // Result found
0913                 break;
0914             }
0915 
0916             // Because the sign of dbeta is changed WRT [Sjoberg02] dbeta is subtracted here
0917             beta = beta - dbeta;
0918 
0919             t = tan(beta);
0920         }
0921 
0922         lat = geod1.lat(t);
0923         // NOTE: if Cj is 0 then the result is lonj or lonj+180
0924         lon = ! geod1.is_Cj_zero
0925                 ? geod1.lon(lon1_diff)
0926                 : geod2.lon(lon2_diff);
0927 
0928         return true;
0929     }
0930 
0931     static inline bool is_result_longitude_ok(geodesic_type const& geod,
0932                                               CT const& lon1, CT const& lon2, CT const& lon)
0933     {
0934         CT const c0 = 0;
0935 
0936         if (geod.is_Cj_zero)
0937             return true; // don't check vertical segment
0938 
0939         CT dist1p = math::longitude_distance_signed<radian>(lon1, lon);
0940         CT dist12 = math::longitude_distance_signed<radian>(lon1, lon2);
0941 
0942         if (dist12 < c0)
0943         {
0944             dist1p = -dist1p;
0945             dist12 = -dist12;
0946         }
0947 
0948         return (c0 <= dist1p && dist1p <= dist12)
0949             || math::equals(dist1p, c0)
0950             || math::equals(dist1p, dist12);
0951     }
0952 
0953     struct geodesics_type
0954     {
0955         geodesics_type(geodesic_type const& g1, geodesic_type const& g2)
0956             : geod1(g1)
0957             , geod2(g2)
0958             , vertex1(geod1.get_vertex_data())
0959             , vertex2(geod2.get_vertex_data())
0960         {}
0961 
0962         geodesic_type const& geod1;
0963         geodesic_type const& geod2;
0964         typename geodesic_type::vertex_data vertex1;
0965         typename geodesic_type::vertex_data vertex2;
0966     };
0967 
0968     struct converge_07_result
0969     {
0970         converge_07_result()
0971             : lon1(0), lon2(0), k1_diff(0), k2_diff(0), t1(0), t2(0)
0972         {}
0973 
0974         CT lon1, lon2;
0975         CT k1_diff, k2_diff;
0976         CT t1, t2;
0977     };
0978 
0979     static inline bool converge_07(geodesic_type const& geod1, geodesic_type const& geod2,
0980                                    CT beta, CT t,
0981                                    CT const& lon1_minus_lon2, CT const& lon_sph,
0982                                    CT & lon, CT & lat)
0983     {
0984         //CT const c0 = 0;
0985         //CT const c1 = 1;
0986         //CT const c2 = 2;
0987         //CT const pi = math::pi<CT>();
0988 
0989         geodesics_type geodesics(geod1, geod2);
0990         converge_07_result result;
0991 
0992         // calculate first pair of longitudes
0993         if (!converge_07_step_one(CT(sin(beta)), t, lon1_minus_lon2, geodesics, lon_sph, result, false))
0994         {
0995             return false;
0996         }
0997 
0998         int t_direction = 0;
0999 
1000         CT lon_diff_prev = math::longitude_difference<radian>(result.lon1, result.lon2);
1001 
1002         // [Sjoberg07]
1003         for (int i = 2; i < max_iterations_07; ++i)
1004         {
1005             // pick t candidates from previous result based on dir
1006             CT t_cand1 = result.t1;
1007             CT t_cand2 = result.t2;
1008             // if direction is 0 the closer one is the first
1009             if (t_direction < 0)
1010             {
1011                 t_cand1 = (std::min)(result.t1, result.t2);
1012                 t_cand2 = (std::max)(result.t1, result.t2);
1013             }
1014             else if (t_direction > 0)
1015             {
1016                 t_cand1 = (std::max)(result.t1, result.t2);
1017                 t_cand2 = (std::min)(result.t1, result.t2);
1018             }
1019             else
1020             {
1021                 t_direction = t_cand1 < t_cand2 ? -1 : 1;
1022             }
1023 
1024             CT t1 = t;
1025             CT beta1 = beta;
1026             // check if the further calculation is needed
1027             if (converge_07_update(t1, beta1, t_cand1))
1028             {
1029                 break;
1030             }
1031 
1032             bool try_t2 = false;
1033             converge_07_result result_curr;
1034             if (converge_07_step_one(CT(sin(beta1)), t1, lon1_minus_lon2, geodesics, lon_sph, result_curr))
1035             {
1036                 CT const lon_diff1 = math::longitude_difference<radian>(result_curr.lon1, result_curr.lon2);
1037                 if (lon_diff_prev > lon_diff1)
1038                 {
1039                     t = t1;
1040                     beta = beta1;
1041                     lon_diff_prev = lon_diff1;
1042                     result = result_curr;
1043                 }
1044                 else if (t_cand1 != t_cand2)
1045                 {
1046                     try_t2 = true;
1047                 }
1048                 else
1049                 {
1050                     // the result is not fully correct but it won't be more accurate
1051                     break;
1052                 }
1053             }
1054             // ! converge_07_step_one
1055             else
1056             {
1057                 if (t_cand1 != t_cand2)
1058                 {
1059                     try_t2 = true;
1060                 }
1061                 else
1062                 {
1063                     return false;
1064                 }
1065             }
1066 
1067 
1068             if (try_t2)
1069             {
1070                 CT t2 = t;
1071                 CT beta2 = beta;
1072                 // check if the further calculation is needed
1073                 if (converge_07_update(t2, beta2, t_cand2))
1074                 {
1075                     break;
1076                 }
1077 
1078                 if (! converge_07_step_one(CT(sin(beta2)), t2, lon1_minus_lon2, geodesics, lon_sph, result_curr))
1079                 {
1080                     return false;
1081                 }
1082 
1083                 CT const lon_diff2 = math::longitude_difference<radian>(result_curr.lon1, result_curr.lon2);
1084                 if (lon_diff_prev > lon_diff2)
1085                 {
1086                     t_direction *= -1;
1087                     t = t2;
1088                     beta = beta2;
1089                     lon_diff_prev = lon_diff2;
1090                     result = result_curr;
1091                 }
1092                 else
1093                 {
1094                     // the result is not fully correct but it won't be more accurate
1095                     break;
1096                 }
1097             }
1098         }
1099 
1100         lat = geod1.lat(t);
1101         lon = ! geod1.is_Cj_zero ? result.lon1 : result.lon2;
1102         math::normalize_longitude<radian>(lon);
1103 
1104         return true;
1105     }
1106 
1107     static inline bool converge_07_update(CT & t, CT & beta, CT const& t_new)
1108     {
1109         CT const c0 = 0;
1110 
1111         CT const beta_new = atan(t_new);
1112         CT const dbeta = beta_new - beta;
1113         beta = beta_new;
1114         t = t_new;
1115 
1116         return math::equals(dbeta, c0);
1117     }
1118 
1119     static inline CT const& pick_t(CT const& t1, CT const& t2, int direction)
1120     {
1121         return direction < 0 ? (std::min)(t1, t2) : (std::max)(t1, t2);
1122     }
1123 
1124     static inline bool converge_07_step_one(CT const& sin_beta,
1125                                             CT const& t,
1126                                             CT const& lon1_minus_lon2,
1127                                             geodesics_type const& geodesics,
1128                                             CT const& lon_sph,
1129                                             converge_07_result & result,
1130                                             bool check_sin_beta = true)
1131     {
1132         bool ok = converge_07_one_geod(sin_beta, t, geodesics.geod1, geodesics.vertex1, lon_sph,
1133                                        result.lon1, result.k1_diff, check_sin_beta)
1134                && converge_07_one_geod(sin_beta, t, geodesics.geod2, geodesics.vertex2, lon_sph,
1135                                        result.lon2, result.k2_diff, check_sin_beta);
1136 
1137         if (!ok)
1138         {
1139             return false;
1140         }
1141 
1142         CT const k = lon1_minus_lon2 + result.k1_diff - result.k2_diff;
1143 
1144         // get 2 possible ts one lesser and one greater than t
1145         // t1 is the closer one
1146         calc_ts(t, k, geodesics.geod1, geodesics.geod2, result.t1, result.t2);
1147 
1148         return true;
1149     }
1150 
1151     static inline bool converge_07_one_geod(CT const& sin_beta, CT const& t,
1152                                             geodesic_type const& geod,
1153                                             typename geodesic_type::vertex_data const& vertex,
1154                                             CT const& lon_sph,
1155                                             CT & lon, CT & k_diff,
1156                                             bool check_sin_beta)
1157     {
1158         using math::detail::bounded;
1159         CT const c1 = 1;
1160 
1161         CT k_diff_before = 0;
1162         CT k_diff_behind = 0;
1163 
1164         bool is_beta_ok = geod.k_diffs(sin_beta, vertex, k_diff_before, k_diff_behind, check_sin_beta);
1165 
1166         if (! is_beta_ok)
1167         {
1168             return false;
1169         }
1170 
1171         CT const asin_t_t0j = ! geod.is_Cj_zero ? asin(bounded(t / geod.t0j, -c1, c1)) : 0;
1172         CT const sign_asin_t_t0j = geod.sign_lon_diff * asin_t_t0j;
1173 
1174         CT const lon_before = geod.lonj + sign_asin_t_t0j + k_diff_before;
1175         CT const lon_behind = geod.lonj - sign_asin_t_t0j + k_diff_behind;
1176 
1177         CT const lon_dist_before = math::longitude_distance_signed<radian>(lon_before, lon_sph);
1178         CT const lon_dist_behind = math::longitude_distance_signed<radian>(lon_behind, lon_sph);
1179         if (math::abs(lon_dist_before) <= math::abs(lon_dist_behind))
1180         {
1181             k_diff = k_diff_before;
1182             lon = lon_before;
1183         }
1184         else
1185         {
1186             k_diff = k_diff_behind;
1187             lon = lon_behind;
1188         }
1189 
1190         return true;
1191     }
1192 
1193     static inline void calc_ts(CT const& t, CT const& k,
1194                                geodesic_type const& geod1, geodesic_type const& geod2,
1195                                CT & t1, CT& t2)
1196     {
1197         CT const c0 = 0;
1198         CT const c1 = 1;
1199         CT const c2 = 2;
1200 
1201         CT const K = sin(k);
1202 
1203         BOOST_GEOMETRY_ASSERT(!geod1.is_Cj_zero || !geod2.is_Cj_zero);
1204         if (geod1.is_Cj_zero)
1205         {
1206             t1 = K * geod2.t0j;
1207             t2 = -t1;
1208         }
1209         else if (geod2.is_Cj_zero)
1210         {
1211             t1 = -K * geod1.t0j;
1212             t2 = -t1;
1213         }
1214         else
1215         {
1216             CT const A = math::sqr(geod1.t0j) + math::sqr(geod2.t0j);
1217             CT const B = c2 * geod1.t0j * geod2.t0j * math::sqrt(c1 - math::sqr(K));
1218 
1219             CT const K_t01_t02 = K * geod1.t0j * geod2.t0j;
1220             CT const D1 = math::sqrt(A + B);
1221             CT const D2 = math::sqrt(A - B);
1222             CT const t_new1 = math::equals(D1, c0) ? c0 : K_t01_t02 / D1;
1223             CT const t_new2 = math::equals(D2, c0) ? c0 : K_t01_t02 / D2;
1224             CT const t_new3 = -t_new1;
1225             CT const t_new4 = -t_new2;
1226 
1227             // Pick 2 nearest t_new, one greater and one lesser than current t
1228             CT const abs_t_new1 = math::abs(t_new1);
1229             CT const abs_t_new2 = math::abs(t_new2);
1230             CT const abs_t_max = (std::max)(abs_t_new1, abs_t_new2);
1231             t1 = -abs_t_max; // lesser
1232             t2 = abs_t_max; // greater
1233             if (t1 < t)
1234             {
1235                 if (t_new1 < t && t_new1 > t1)
1236                     t1 = t_new1;
1237                 if (t_new2 < t && t_new2 > t1)
1238                     t1 = t_new2;
1239                 if (t_new3 < t && t_new3 > t1)
1240                     t1 = t_new3;
1241                 if (t_new4 < t && t_new4 > t1)
1242                     t1 = t_new4;
1243             }
1244             if (t2 > t)
1245             {
1246                 if (t_new1 > t && t_new1 < t2)
1247                     t2 = t_new1;
1248                 if (t_new2 > t && t_new2 < t2)
1249                     t2 = t_new2;
1250                 if (t_new3 > t && t_new3 < t2)
1251                     t2 = t_new3;
1252                 if (t_new4 > t && t_new4 < t2)
1253                     t2 = t_new4;
1254             }
1255         }
1256 
1257         // the first one is the closer one
1258         if (math::abs(t - t2) < math::abs(t - t1))
1259         {
1260             std::swap(t2, t1);
1261         }
1262     }
1263 
1264     static inline CT fj(CT const& cos_beta, CT const& cos2_beta, CT const& Cj, CT const& e_sqr)
1265     {
1266         CT const c1 = 1;
1267         CT const Cj_sqr = math::sqr(Cj);
1268         return Cj / cos_beta * math::sqrt((c1 - e_sqr * cos2_beta) / (cos2_beta - Cj_sqr));
1269     }
1270 
1271     /*static inline CT vertical_intersection_longitude(CT const& ip_lon, CT const& seg_lon1, CT const& seg_lon2)
1272     {
1273         CT const c0 = 0;
1274         CT const lon_2 = ip_lon > c0 ? ip_lon - pi : ip_lon + pi;
1275 
1276         return (std::min)(math::longitude_difference<radian>(ip_lon, seg_lon1),
1277                           math::longitude_difference<radian>(ip_lon, seg_lon2))
1278             <=
1279                (std::min)(math::longitude_difference<radian>(lon_2, seg_lon1),
1280                           math::longitude_difference<radian>(lon_2, seg_lon2))
1281             ? ip_lon : lon_2;
1282     }*/
1283 };
1284 
1285 }}} // namespace boost::geometry::formula
1286 
1287 
1288 #endif // BOOST_GEOMETRY_FORMULAS_SJOBERG_INTERSECTION_HPP
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!