summaryrefslogtreecommitdiff
path: root/boost/geometry/formulas
diff options
context:
space:
mode:
Diffstat (limited to 'boost/geometry/formulas')
-rw-r--r--boost/geometry/formulas/andoyer_inverse.hpp27
-rw-r--r--boost/geometry/formulas/area_formulas.hpp578
-rw-r--r--boost/geometry/formulas/differential_quantities.hpp19
-rw-r--r--boost/geometry/formulas/eccentricity_sqr.hpp70
-rw-r--r--boost/geometry/formulas/flattening.hpp70
-rw-r--r--boost/geometry/formulas/geographic.hpp457
-rw-r--r--boost/geometry/formulas/gnomonic_spheroid.hpp3
-rw-r--r--boost/geometry/formulas/sjoberg_intersection.hpp1310
-rw-r--r--boost/geometry/formulas/spherical.hpp152
-rw-r--r--boost/geometry/formulas/thomas_direct.hpp5
-rw-r--r--boost/geometry/formulas/thomas_inverse.hpp5
-rw-r--r--boost/geometry/formulas/vertex_latitude.hpp148
-rw-r--r--boost/geometry/formulas/vincenty_direct.hpp5
-rw-r--r--boost/geometry/formulas/vincenty_inverse.hpp38
14 files changed, 2594 insertions, 293 deletions
diff --git a/boost/geometry/formulas/andoyer_inverse.hpp b/boost/geometry/formulas/andoyer_inverse.hpp
index 57b5ab5384..902fd7d8f6 100644
--- a/boost/geometry/formulas/andoyer_inverse.hpp
+++ b/boost/geometry/formulas/andoyer_inverse.hpp
@@ -1,6 +1,6 @@
// Boost.Geometry
-// Copyright (c) 2015-2016 Oracle and/or its affiliates.
+// Copyright (c) 2015-2017 Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
@@ -20,9 +20,8 @@
#include <boost/geometry/util/condition.hpp>
#include <boost/geometry/util/math.hpp>
-#include <boost/geometry/algorithms/detail/flattening.hpp>
-
#include <boost/geometry/formulas/differential_quantities.hpp>
+#include <boost/geometry/formulas/flattening.hpp>
#include <boost/geometry/formulas/result_inverse.hpp>
@@ -75,7 +74,7 @@ public:
CT const c0 = CT(0);
CT const c1 = CT(1);
CT const pi = math::pi<CT>();
- CT const f = detail::flattening<CT>(spheroid);
+ CT const f = formula::flattening<CT>(spheroid);
CT const dlon = lon2 - lon1;
CT const sin_dlon = sin(dlon);
@@ -97,7 +96,7 @@ public:
CT const d = acos(cos_d); // [0, pi]
CT const sin_d = sin(d); // [-1, 1]
-
+
if ( BOOST_GEOMETRY_CONDITION(EnableDistance) )
{
CT const K = math::sqr(sin_lat1-sin_lat2);
@@ -138,7 +137,14 @@ public:
CT A = c0;
CT U = c0;
- if ( ! math::equals(cos_lat2, c0) )
+ if (math::equals(cos_lat2, c0))
+ {
+ if (sin_lat2 < c0)
+ {
+ A = pi;
+ }
+ }
+ else
{
CT const tan_lat2 = sin_lat2/cos_lat2;
CT const M = cos_lat1*tan_lat2-sin_lat1*cos_dlon;
@@ -149,7 +155,14 @@ public:
CT B = c0;
CT V = c0;
- if ( ! math::equals(cos_lat1, c0) )
+ if (math::equals(cos_lat1, c0))
+ {
+ if (sin_lat1 < c0)
+ {
+ B = pi;
+ }
+ }
+ else
{
CT const tan_lat1 = sin_lat1/cos_lat1;
CT const N = cos_lat2*tan_lat1-sin_lat2*cos_dlon;
diff --git a/boost/geometry/formulas/area_formulas.hpp b/boost/geometry/formulas/area_formulas.hpp
new file mode 100644
index 0000000000..6a0b525e25
--- /dev/null
+++ b/boost/geometry/formulas/area_formulas.hpp
@@ -0,0 +1,578 @@
+// Boost.Geometry
+
+// Copyright (c) 2015-2016 Oracle and/or its affiliates.
+
+// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_GEOMETRY_FORMULAS_AREA_FORMULAS_HPP
+#define BOOST_GEOMETRY_FORMULAS_AREA_FORMULAS_HPP
+
+#include <boost/geometry/formulas/flattening.hpp>
+#include <boost/math/special_functions/hypot.hpp>
+
+namespace boost { namespace geometry { namespace formula
+{
+
+/*!
+\brief Formulas for computing spherical and ellipsoidal polygon area.
+ The current class computes the area of the trapezoid defined by a segment
+ the two meridians passing by the endpoints and the equator.
+\author See
+- Danielsen JS, The area under the geodesic. Surv Rev 30(232):
+61–66, 1989
+- Charles F.F Karney, Algorithms for geodesics, 2011
+https://arxiv.org/pdf/1109.4448.pdf
+*/
+
+template <
+ typename CT,
+ std::size_t SeriesOrder = 2,
+ bool ExpandEpsN = true
+>
+class area_formulas
+{
+
+public:
+
+ //TODO: move the following to a more general space to be used by other
+ // classes as well
+ /*
+ Evaluate the polynomial in x using Horner's method.
+ */
+ template <typename NT, typename IteratorType>
+ static inline NT horner_evaluate(NT x,
+ IteratorType begin,
+ IteratorType end)
+ {
+ NT result(0);
+ IteratorType it = end;
+ do
+ {
+ result = result * x + *--it;
+ }
+ while (it != begin);
+ return result;
+ }
+
+ /*
+ Clenshaw algorithm for summing trigonometric series
+ https://en.wikipedia.org/wiki/Clenshaw_algorithm
+ */
+ template <typename NT, typename IteratorType>
+ static inline NT clenshaw_sum(NT cosx,
+ IteratorType begin,
+ IteratorType end)
+ {
+ IteratorType it = end;
+ bool odd = true;
+ CT b_k, b_k1(0), b_k2(0);
+ do
+ {
+ CT c_k = odd ? *--it : NT(0);
+ b_k = c_k + NT(2) * cosx * b_k1 - b_k2;
+ b_k2 = b_k1;
+ b_k1 = b_k;
+ odd = !odd;
+ }
+ while (it != begin);
+
+ return *begin + b_k1 * cosx - b_k2;
+ }
+
+ template<typename T>
+ static inline void normalize(T& x, T& y)
+ {
+ T h = boost::math::hypot(x, y);
+ x /= h;
+ y /= h;
+ }
+
+ /*
+ Generate and evaluate the series expansion of the following integral
+
+ I4 = -integrate( (t(ep2) - t(k2*sin(sigma1)^2)) / (ep2 - k2*sin(sigma1)^2)
+ * sin(sigma1)/2, sigma1, pi/2, sigma )
+ where
+
+ t(x) = sqrt(1+1/x)*asinh(sqrt(x)) + x
+
+ valid for ep2 and k2 small. We substitute k2 = 4 * eps / (1 - eps)^2
+ and ep2 = 4 * n / (1 - n)^2 and expand in eps and n.
+
+ The resulting sum of the series is of the form
+
+ sum(C4[l] * cos((2*l+1)*sigma), l, 0, maxpow-1) )
+
+ The above expansion is performed in Computer Algebra System Maxima.
+ The C++ code (that yields the function evaluate_coeffs_n below) is generated
+ by the following Maxima script and is based on script:
+ http://geographiclib.sourceforge.net/html/geod.mac
+
+ // Maxima script begin
+ taylordepth:5$
+ ataylor(expr,var,ord):=expand(ratdisrep(taylor(expr,var,0,ord)))$
+ jtaylor(expr,var1,var2,ord):=block([zz],expand(subst([zz=1],
+ ratdisrep(taylor(subst([var1=zz*var1,var2=zz*var2],expr),zz,0,ord)))))$
+
+ compute(maxpow):=block([int,t,intexp,area, x,ep2,k2],
+ maxpow:maxpow-1,
+ t : sqrt(1+1/x) * asinh(sqrt(x)) + x,
+ int:-(tf(ep2) - tf(k2*sin(sigma)^2)) / (ep2 - k2*sin(sigma)^2)
+ * sin(sigma)/2,
+ int:subst([tf(ep2)=subst([x=ep2],t),
+ tf(k2*sin(sigma)^2)=subst([x=k2*sin(sigma)^2],t)],
+ int),
+ int:subst([abs(sin(sigma))=sin(sigma)],int),
+ int:subst([k2=4*eps/(1-eps)^2,ep2=4*n/(1-n)^2],int),
+ intexp:jtaylor(int,n,eps,maxpow),
+ area:trigreduce(integrate(intexp,sigma)),
+ area:expand(area-subst(sigma=%pi/2,area)),
+ for i:0 thru maxpow do C4[i]:coeff(area,cos((2*i+1)*sigma)),
+ if expand(area-sum(C4[i]*cos((2*i+1)*sigma),i,0,maxpow)) # 0
+ then error("left over terms in I4"),
+ 'done)$
+
+ printcode(maxpow):=
+ block([tab2:" ",tab3:" "],
+ print(" switch (SeriesOrder) {"),
+ for nn:1 thru maxpow do block([c],
+ print(concat(tab2,"case ",string(nn-1),":")),
+ c:0,
+ for m:0 thru nn-1 do block(
+ [q:jtaylor(subst([n=n],C4[m]),n,eps,nn-1),
+ linel:1200],
+ for j:m thru nn-1 do (
+ print(concat(tab3,"coeffs_n[",c,"] = ",
+ string(horner(coeff(q,eps,j))),";")),
+ c:c+1)
+ ),
+ print(concat(tab3,"break;"))),
+ print(" }"),
+ 'done)$
+
+ maxpow:6$
+ compute(maxpow)$
+ printcode(maxpow)$
+ // Maxima script end
+
+ In the resulting code we should replace each number x by CT(x)
+ e.g. using the following scirpt:
+ sed -e 's/[0-9]\+/CT(&)/g; s/\[CT(/\[/g; s/)\]/\]/g;
+ s/case\sCT(/case /g; s/):/:/g'
+ */
+
+ static inline void evaluate_coeffs_n(CT n, CT coeffs_n[])
+ {
+
+ switch (SeriesOrder) {
+ case 0:
+ coeffs_n[0] = CT(2)/CT(3);
+ break;
+ case 1:
+ coeffs_n[0] = (CT(10)-CT(4)*n)/CT(15);
+ coeffs_n[1] = -CT(1)/CT(5);
+ coeffs_n[2] = CT(1)/CT(45);
+ break;
+ case 2:
+ coeffs_n[0] = (n*(CT(8)*n-CT(28))+CT(70))/CT(105);
+ coeffs_n[1] = (CT(16)*n-CT(7))/CT(35);
+ coeffs_n[2] = -CT(2)/CT(105);
+ coeffs_n[3] = (CT(7)-CT(16)*n)/CT(315);
+ coeffs_n[4] = -CT(2)/CT(105);
+ coeffs_n[5] = CT(4)/CT(525);
+ break;
+ case 3:
+ coeffs_n[0] = (n*(n*(CT(4)*n+CT(24))-CT(84))+CT(210))/CT(315);
+ coeffs_n[1] = ((CT(48)-CT(32)*n)*n-CT(21))/CT(105);
+ coeffs_n[2] = (-CT(32)*n-CT(6))/CT(315);
+ coeffs_n[3] = CT(11)/CT(315);
+ coeffs_n[4] = (n*(CT(32)*n-CT(48))+CT(21))/CT(945);
+ coeffs_n[5] = (CT(64)*n-CT(18))/CT(945);
+ coeffs_n[6] = -CT(1)/CT(105);
+ coeffs_n[7] = (CT(12)-CT(32)*n)/CT(1575);
+ coeffs_n[8] = -CT(8)/CT(1575);
+ coeffs_n[9] = CT(8)/CT(2205);
+ break;
+ case 4:
+ coeffs_n[0] = (n*(n*(n*(CT(16)*n+CT(44))+CT(264))-CT(924))+CT(2310))/CT(3465);
+ coeffs_n[1] = (n*(n*(CT(48)*n-CT(352))+CT(528))-CT(231))/CT(1155);
+ coeffs_n[2] = (n*(CT(1088)*n-CT(352))-CT(66))/CT(3465);
+ coeffs_n[3] = (CT(121)-CT(368)*n)/CT(3465);
+ coeffs_n[4] = CT(4)/CT(1155);
+ coeffs_n[5] = (n*((CT(352)-CT(48)*n)*n-CT(528))+CT(231))/CT(10395);
+ coeffs_n[6] = ((CT(704)-CT(896)*n)*n-CT(198))/CT(10395);
+ coeffs_n[7] = (CT(80)*n-CT(99))/CT(10395);
+ coeffs_n[8] = CT(4)/CT(1155);
+ coeffs_n[9] = (n*(CT(320)*n-CT(352))+CT(132))/CT(17325);
+ coeffs_n[10] = (CT(384)*n-CT(88))/CT(17325);
+ coeffs_n[11] = -CT(8)/CT(1925);
+ coeffs_n[12] = (CT(88)-CT(256)*n)/CT(24255);
+ coeffs_n[13] = -CT(16)/CT(8085);
+ coeffs_n[14] = CT(64)/CT(31185);
+ break;
+ case 5:
+ coeffs_n[0] = (n*(n*(n*(n*(CT(100)*n+CT(208))+CT(572))+CT(3432))-CT(12012))+CT(30030))
+ /CT(45045);
+ coeffs_n[1] = (n*(n*(n*(CT(64)*n+CT(624))-CT(4576))+CT(6864))-CT(3003))/CT(15015);
+ coeffs_n[2] = (n*((CT(14144)-CT(10656)*n)*n-CT(4576))-CT(858))/CT(45045);
+ coeffs_n[3] = ((-CT(224)*n-CT(4784))*n+CT(1573))/CT(45045);
+ coeffs_n[4] = (CT(1088)*n+CT(156))/CT(45045);
+ coeffs_n[5] = CT(97)/CT(15015);
+ coeffs_n[6] = (n*(n*((-CT(64)*n-CT(624))*n+CT(4576))-CT(6864))+CT(3003))/CT(135135);
+ coeffs_n[7] = (n*(n*(CT(5952)*n-CT(11648))+CT(9152))-CT(2574))/CT(135135);
+ coeffs_n[8] = (n*(CT(5792)*n+CT(1040))-CT(1287))/CT(135135);
+ coeffs_n[9] = (CT(468)-CT(2944)*n)/CT(135135);
+ coeffs_n[10] = CT(1)/CT(9009);
+ coeffs_n[11] = (n*((CT(4160)-CT(1440)*n)*n-CT(4576))+CT(1716))/CT(225225);
+ coeffs_n[12] = ((CT(4992)-CT(8448)*n)*n-CT(1144))/CT(225225);
+ coeffs_n[13] = (CT(1856)*n-CT(936))/CT(225225);
+ coeffs_n[14] = CT(8)/CT(10725);
+ coeffs_n[15] = (n*(CT(3584)*n-CT(3328))+CT(1144))/CT(315315);
+ coeffs_n[16] = (CT(1024)*n-CT(208))/CT(105105);
+ coeffs_n[17] = -CT(136)/CT(63063);
+ coeffs_n[18] = (CT(832)-CT(2560)*n)/CT(405405);
+ coeffs_n[19] = -CT(128)/CT(135135);
+ coeffs_n[20] = CT(128)/CT(99099);
+ break;
+ }
+ }
+
+ /*
+ Expand in k2 and ep2.
+ */
+ static inline void evaluate_coeffs_ep(CT ep, CT coeffs_n[])
+ {
+ switch (SeriesOrder) {
+ case 0:
+ coeffs_n[0] = CT(2)/CT(3);
+ break;
+ case 1:
+ coeffs_n[0] = (CT(10)-ep)/CT(15);
+ coeffs_n[1] = -CT(1)/CT(20);
+ coeffs_n[2] = CT(1)/CT(180);
+ break;
+ case 2:
+ coeffs_n[0] = (ep*(CT(4)*ep-CT(7))+CT(70))/CT(105);
+ coeffs_n[1] = (CT(4)*ep-CT(7))/CT(140);
+ coeffs_n[2] = CT(1)/CT(42);
+ coeffs_n[3] = (CT(7)-CT(4)*ep)/CT(1260);
+ coeffs_n[4] = -CT(1)/CT(252);
+ coeffs_n[5] = CT(1)/CT(2100);
+ break;
+ case 3:
+ coeffs_n[0] = (ep*((CT(12)-CT(8)*ep)*ep-CT(21))+CT(210))/CT(315);
+ coeffs_n[1] = ((CT(12)-CT(8)*ep)*ep-CT(21))/CT(420);
+ coeffs_n[2] = (CT(3)-CT(2)*ep)/CT(126);
+ coeffs_n[3] = -CT(1)/CT(72);
+ coeffs_n[4] = (ep*(CT(8)*ep-CT(12))+CT(21))/CT(3780);
+ coeffs_n[5] = (CT(2)*ep-CT(3))/CT(756);
+ coeffs_n[6] = CT(1)/CT(360);
+ coeffs_n[7] = (CT(3)-CT(2)*ep)/CT(6300);
+ coeffs_n[8] = -CT(1)/CT(1800);
+ coeffs_n[9] = CT(1)/CT(17640);
+ break;
+ case 4:
+ coeffs_n[0] = (ep*(ep*(ep*(CT(64)*ep-CT(88))+CT(132))-CT(231))+CT(2310))/CT(3465);
+ coeffs_n[1] = (ep*(ep*(CT(64)*ep-CT(88))+CT(132))-CT(231))/CT(4620);
+ coeffs_n[2] = (ep*(CT(16)*ep-CT(22))+CT(33))/CT(1386);
+ coeffs_n[3] = (CT(8)*ep-CT(11))/CT(792);
+ coeffs_n[4] = CT(1)/CT(110);
+ coeffs_n[5] = (ep*((CT(88)-CT(64)*ep)*ep-CT(132))+CT(231))/CT(41580);
+ coeffs_n[6] = ((CT(22)-CT(16)*ep)*ep-CT(33))/CT(8316);
+ coeffs_n[7] = (CT(11)-CT(8)*ep)/CT(3960);
+ coeffs_n[8] = -CT(1)/CT(495);
+ coeffs_n[9] = (ep*(CT(16)*ep-CT(22))+CT(33))/CT(69300);
+ coeffs_n[10] = (CT(8)*ep-CT(11))/CT(19800);
+ coeffs_n[11] = CT(1)/CT(1925);
+ coeffs_n[12] = (CT(11)-CT(8)*ep)/CT(194040);
+ coeffs_n[13] = -CT(1)/CT(10780);
+ coeffs_n[14] = CT(1)/CT(124740);
+ break;
+ case 5:
+ coeffs_n[0] = (ep*(ep*(ep*((CT(832)-CT(640)*ep)*ep-CT(1144))+CT(1716))-CT(3003))+CT(30030))/CT(45045);
+ coeffs_n[1] = (ep*(ep*((CT(832)-CT(640)*ep)*ep-CT(1144))+CT(1716))-CT(3003))/CT(60060);
+ coeffs_n[2] = (ep*((CT(208)-CT(160)*ep)*ep-CT(286))+CT(429))/CT(18018);
+ coeffs_n[3] = ((CT(104)-CT(80)*ep)*ep-CT(143))/CT(10296);
+ coeffs_n[4] = (CT(13)-CT(10)*ep)/CT(1430);
+ coeffs_n[5] = -CT(1)/CT(156);
+ coeffs_n[6] = (ep*(ep*(ep*(CT(640)*ep-CT(832))+CT(1144))-CT(1716))+CT(3003))/CT(540540);
+ coeffs_n[7] = (ep*(ep*(CT(160)*ep-CT(208))+CT(286))-CT(429))/CT(108108);
+ coeffs_n[8] = (ep*(CT(80)*ep-CT(104))+CT(143))/CT(51480);
+ coeffs_n[9] = (CT(10)*ep-CT(13))/CT(6435);
+ coeffs_n[10] = CT(5)/CT(3276);
+ coeffs_n[11] = (ep*((CT(208)-CT(160)*ep)*ep-CT(286))+CT(429))/CT(900900);
+ coeffs_n[12] = ((CT(104)-CT(80)*ep)*ep-CT(143))/CT(257400);
+ coeffs_n[13] = (CT(13)-CT(10)*ep)/CT(25025);
+ coeffs_n[14] = -CT(1)/CT(2184);
+ coeffs_n[15] = (ep*(CT(80)*ep-CT(104))+CT(143))/CT(2522520);
+ coeffs_n[16] = (CT(10)*ep-CT(13))/CT(140140);
+ coeffs_n[17] = CT(5)/CT(45864);
+ coeffs_n[18] = (CT(13)-CT(10)*ep)/CT(1621620);
+ coeffs_n[19] = -CT(1)/CT(58968);
+ coeffs_n[20] = CT(1)/CT(792792);
+ break;
+ }
+ }
+
+ /*
+ Given the set of coefficients coeffs1[] evaluate on var2 and return
+ the set of coefficients coeffs2[]
+ */
+ static inline void evaluate_coeffs_var2(CT var2,
+ CT coeffs1[],
+ CT coeffs2[]){
+ std::size_t begin(0), end(0);
+ for(std::size_t i = 0; i <= SeriesOrder; i++){
+ end = begin + SeriesOrder + 1 - i;
+ coeffs2[i] = ((i==0) ? CT(1) : pow(var2,int(i)))
+ * horner_evaluate(var2, coeffs1 + begin, coeffs1 + end);
+ begin = end;
+ }
+ }
+
+ /*
+ Compute the spherical excess of a geodesic (or shperical) segment
+ */
+ template <
+ bool LongSegment,
+ typename PointOfSegment
+ >
+ static inline CT spherical(PointOfSegment const& p1,
+ PointOfSegment const& p2)
+ {
+ CT excess;
+
+ if(LongSegment) // not for segments parallel to equator
+ {
+ CT cbet1 = cos(geometry::get_as_radian<1>(p1));
+ CT sbet1 = sin(geometry::get_as_radian<1>(p1));
+ CT cbet2 = cos(geometry::get_as_radian<1>(p2));
+ CT sbet2 = sin(geometry::get_as_radian<1>(p2));
+
+ CT omg12 = geometry::get_as_radian<0>(p1)
+ - geometry::get_as_radian<0>(p2);
+ CT comg12 = cos(omg12);
+ CT somg12 = sin(omg12);
+
+ CT alp1 = atan2(cbet1 * sbet2
+ - sbet1 * cbet2 * comg12,
+ cbet2 * somg12);
+
+ CT alp2 = atan2(cbet1 * sbet2 * comg12
+ - sbet1 * cbet2,
+ cbet1 * somg12);
+
+ excess = alp2 - alp1;
+
+ } else {
+
+ // Trapezoidal formula
+
+ CT tan_lat1 =
+ tan(geometry::get_as_radian<1>(p1) / 2.0);
+ CT tan_lat2 =
+ tan(geometry::get_as_radian<1>(p2) / 2.0);
+
+ excess = CT(2.0)
+ * atan(((tan_lat1 + tan_lat2) / (CT(1) + tan_lat1 * tan_lat2))
+ * tan((geometry::get_as_radian<0>(p2)
+ - geometry::get_as_radian<0>(p1)) / 2));
+ }
+
+ return excess;
+ }
+
+ struct return_type_ellipsoidal
+ {
+ return_type_ellipsoidal()
+ : spherical_term(0),
+ ellipsoidal_term(0)
+ {}
+
+ CT spherical_term;
+ CT ellipsoidal_term;
+ };
+
+ /*
+ Compute the ellipsoidal correction of a geodesic (or shperical) segment
+ */
+ template <
+ template <typename, bool, bool, bool, bool, bool> class Inverse,
+ //typename AzimuthStrategy,
+ typename PointOfSegment,
+ typename SpheroidConst
+ >
+ static inline return_type_ellipsoidal ellipsoidal(PointOfSegment const& p1,
+ PointOfSegment const& p2,
+ SpheroidConst spheroid_const)
+ {
+ return_type_ellipsoidal result;
+
+ // Azimuth Approximation
+
+ typedef Inverse<CT, false, true, true, false, false> inverse_type;
+ typedef typename inverse_type::result_type inverse_result;
+
+ inverse_result i_res = inverse_type::apply(get_as_radian<0>(p1),
+ get_as_radian<1>(p1),
+ get_as_radian<0>(p2),
+ get_as_radian<1>(p2),
+ spheroid_const.m_spheroid);
+
+ CT alp1 = i_res.azimuth;
+ CT alp2 = i_res.reverse_azimuth;
+
+ // Constants
+
+ CT const ep = spheroid_const.m_ep;
+ CT const f = formula::flattening<CT>(spheroid_const.m_spheroid);
+ CT const one_minus_f = CT(1) - f;
+ std::size_t const series_order_plus_one = SeriesOrder + 1;
+ std::size_t const series_order_plus_two = SeriesOrder + 2;
+
+ // Basic trigonometric computations
+
+ CT tan_bet1 = tan(get_as_radian<1>(p1)) * one_minus_f;
+ CT tan_bet2 = tan(get_as_radian<1>(p2)) * one_minus_f;
+ CT cos_bet1 = cos(atan(tan_bet1));
+ CT cos_bet2 = cos(atan(tan_bet2));
+ CT sin_bet1 = tan_bet1 * cos_bet1;
+ CT sin_bet2 = tan_bet2 * cos_bet2;
+ CT sin_alp1 = sin(alp1);
+ CT cos_alp1 = cos(alp1);
+ CT cos_alp2 = cos(alp2);
+ CT sin_alp0 = sin_alp1 * cos_bet1;
+
+ // Spherical term computation
+
+ CT sin_omg1 = sin_alp0 * sin_bet1;
+ CT cos_omg1 = cos_alp1 * cos_bet1;
+ CT sin_omg2 = sin_alp0 * sin_bet2;
+ CT cos_omg2 = cos_alp2 * cos_bet2;
+ CT cos_omg12 = cos_omg1 * cos_omg2 + sin_omg1 * sin_omg2;
+ CT excess;
+
+ bool meridian = get<0>(p2) - get<0>(p1) == CT(0)
+ || get<1>(p1) == CT(90) || get<1>(p1) == -CT(90)
+ || get<1>(p2) == CT(90) || get<1>(p2) == -CT(90);
+
+ if (!meridian && cos_omg12 > -CT(0.7)
+ && sin_bet2 - sin_bet1 < CT(1.75)) // short segment
+ {
+ CT sin_omg12 = cos_omg1 * sin_omg2 - sin_omg1 * cos_omg2;
+ normalize(sin_omg12, cos_omg12);
+
+ CT cos_omg12p1 = CT(1) + cos_omg12;
+ CT cos_bet1p1 = CT(1) + cos_bet1;
+ CT cos_bet2p1 = CT(1) + cos_bet2;
+ excess = CT(2) * atan2(sin_omg12 * (sin_bet1 * cos_bet2p1 + sin_bet2 * cos_bet1p1),
+ cos_omg12p1 * (sin_bet1 * sin_bet2 + cos_bet1p1 * cos_bet2p1));
+ }
+ else
+ {
+ /*
+ CT sin_alp2 = sin(alp2);
+ CT sin_alp12 = sin_alp2 * cos_alp1 - cos_alp2 * sin_alp1;
+ CT cos_alp12 = cos_alp2 * cos_alp1 + sin_alp2 * sin_alp1;
+ excess = atan2(sin_alp12, cos_alp12);
+ */
+ excess = alp2 - alp1;
+ }
+
+ result.spherical_term = excess;
+
+ // Ellipsoidal term computation (uses integral approximation)
+
+ CT cos_alp0 = math::sqrt(CT(1) - math::sqr(sin_alp0));
+ CT cos_sig1 = cos_alp1 * cos_bet1;
+ CT cos_sig2 = cos_alp2 * cos_bet2;
+ CT sin_sig1 = sin_bet1;
+ CT sin_sig2 = sin_bet2;
+
+ normalize(sin_sig1, cos_sig1);
+ normalize(sin_sig2, cos_sig2);
+
+ CT coeffs[SeriesOrder + 1];
+ const std::size_t coeffs_var_size = (series_order_plus_two
+ * series_order_plus_one) / 2;
+ CT coeffs_var[coeffs_var_size];
+
+ if(ExpandEpsN){ // expand by eps and n
+
+ CT k2 = math::sqr(ep * cos_alp0);
+ CT sqrt_k2_plus_one = math::sqrt(CT(1) + k2);
+ CT eps = (sqrt_k2_plus_one - CT(1)) / (sqrt_k2_plus_one + CT(1));
+ CT n = f / (CT(2) - f);
+
+ // Generate and evaluate the polynomials on n
+ // to get the series coefficients (that depend on eps)
+ evaluate_coeffs_n(n, coeffs_var);
+
+ // Generate and evaluate the polynomials on eps (i.e. var2 = eps)
+ // to get the final series coefficients
+ evaluate_coeffs_var2(eps, coeffs_var, coeffs);
+
+ }else{ // expand by k2 and ep
+
+ CT k2 = math::sqr(ep * cos_alp0);
+ CT ep2 = math::sqr(ep);
+
+ // Generate and evaluate the polynomials on ep2
+ evaluate_coeffs_ep(ep2, coeffs_var);
+
+ // Generate and evaluate the polynomials on k2 (i.e. var2 = k2)
+ evaluate_coeffs_var2(k2, coeffs_var, coeffs);
+
+ }
+
+ // Evaluate the trigonometric sum
+ CT I12 = clenshaw_sum(cos_sig2, coeffs, coeffs + series_order_plus_one)
+ - clenshaw_sum(cos_sig1, coeffs, coeffs + series_order_plus_one);
+
+ // The part of the ellipsodal correction that depends on
+ // point coordinates
+ result.ellipsoidal_term = cos_alp0 * sin_alp0 * I12;
+
+ return result;
+ }
+
+ // Keep track whenever a segment crosses the prime meridian
+ // First normalize to [0,360)
+ template <typename PointOfSegment, typename StateType>
+ static inline int crosses_prime_meridian(PointOfSegment const& p1,
+ PointOfSegment const& p2,
+ StateType& state)
+ {
+ CT const pi
+ = geometry::math::pi<CT>();
+ CT const two_pi
+ = geometry::math::two_pi<CT>();
+
+ CT p1_lon = get_as_radian<0>(p1)
+ - ( floor( get_as_radian<0>(p1) / two_pi )
+ * two_pi );
+ CT p2_lon = get_as_radian<0>(p2)
+ - ( floor( get_as_radian<0>(p2) / two_pi )
+ * two_pi );
+
+ CT max_lon = (std::max)(p1_lon, p2_lon);
+ CT min_lon = (std::min)(p1_lon, p2_lon);
+
+ if(max_lon > pi && min_lon < pi && max_lon - min_lon > pi)
+ {
+ state.m_crosses_prime_meridian++;
+ }
+
+ return state.m_crosses_prime_meridian;
+ }
+
+};
+
+}}} // namespace boost::geometry::formula
+
+
+#endif // BOOST_GEOMETRY_FORMULAS_AREA_FORMULAS_HPP
diff --git a/boost/geometry/formulas/differential_quantities.hpp b/boost/geometry/formulas/differential_quantities.hpp
index 9a92f14e18..ff2ec539db 100644
--- a/boost/geometry/formulas/differential_quantities.hpp
+++ b/boost/geometry/formulas/differential_quantities.hpp
@@ -1,6 +1,6 @@
// Boost.Geometry
-// Copyright (c) 2016 Oracle and/or its affiliates.
+// Copyright (c) 2016-2017 Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
@@ -64,8 +64,8 @@ public:
CT const c1 = 1;
CT const one_minus_f = c1 - f;
- CT const sin_bet1 = one_minus_f * sin_lat1;
- CT const sin_bet2 = one_minus_f * sin_lat2;
+ CT sin_bet1 = one_minus_f * sin_lat1;
+ CT sin_bet2 = one_minus_f * sin_lat2;
// equator
if (math::equals(sin_bet1, c0) && math::equals(sin_bet2, c0))
@@ -89,14 +89,17 @@ public:
CT const e2 = f * (c2 - f);
CT const ep2 = e2 / math::sqr(one_minus_f);
- CT const cos_bet1 = cos_lat1;
- CT const cos_bet2 = cos_lat2;
-
CT const sin_alp1 = sin(azimuth);
CT const cos_alp1 = cos(azimuth);
//CT const sin_alp2 = sin(reverse_azimuth);
CT const cos_alp2 = cos(reverse_azimuth);
+ CT cos_bet1 = cos_lat1;
+ CT cos_bet2 = cos_lat2;
+
+ normalize(sin_bet1, cos_bet1);
+ normalize(sin_bet2, cos_bet2);
+
CT sin_sig1 = sin_bet1;
CT cos_sig1 = cos_alp1 * cos_bet1;
CT sin_sig2 = sin_bet2;
@@ -112,8 +115,8 @@ public:
J12_f(sin_sig1, cos_sig1, sin_sig2, cos_sig2, cos_alp0_sqr, f) :
J12_ep_sqr(sin_sig1, cos_sig1, sin_sig2, cos_sig2, cos_alp0_sqr, ep2) ;
- CT const dn1 = math::sqrt(c1 + e2 * math::sqr(sin_lat1));
- CT const dn2 = math::sqrt(c1 + e2 * math::sqr(sin_lat2));
+ CT const dn1 = math::sqrt(c1 + ep2 * math::sqr(sin_bet1));
+ CT const dn2 = math::sqrt(c1 + ep2 * math::sqr(sin_bet2));
if (BOOST_GEOMETRY_CONDITION(EnableReducedLength))
{
diff --git a/boost/geometry/formulas/eccentricity_sqr.hpp b/boost/geometry/formulas/eccentricity_sqr.hpp
new file mode 100644
index 0000000000..01a9beacb9
--- /dev/null
+++ b/boost/geometry/formulas/eccentricity_sqr.hpp
@@ -0,0 +1,70 @@
+// Boost.Geometry
+
+// Copyright (c) 2016 Oracle and/or its affiliates.
+
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_GEOMETRY_FORMULAS_ECCENCRICITY_SQR_HPP
+#define BOOST_GEOMETRY_FORMULAS_ECCENCRICITY_SQR_HPP
+
+#include <boost/geometry/core/radius.hpp>
+#include <boost/geometry/core/tag.hpp>
+#include <boost/geometry/core/tags.hpp>
+
+#include <boost/geometry/algorithms/not_implemented.hpp>
+
+namespace boost { namespace geometry
+{
+
+#ifndef DOXYGEN_NO_DISPATCH
+namespace formula_dispatch
+{
+
+template <typename ResultType, typename Geometry, typename Tag = typename tag<Geometry>::type>
+struct eccentricity_sqr
+ : not_implemented<Tag>
+{};
+
+template <typename ResultType, typename Geometry>
+struct eccentricity_sqr<ResultType, Geometry, srs_sphere_tag>
+{
+ static inline ResultType apply(Geometry const& /*geometry*/)
+ {
+ return ResultType(0);
+ }
+};
+
+template <typename ResultType, typename Geometry>
+struct eccentricity_sqr<ResultType, Geometry, srs_spheroid_tag>
+{
+ static inline ResultType apply(Geometry const& geometry)
+ {
+ // 1 - (b / a)^2
+ return ResultType(1) - math::sqr(ResultType(get_radius<2>(geometry))
+ / ResultType(get_radius<0>(geometry)));
+ }
+};
+
+} // namespace formula_dispatch
+#endif // DOXYGEN_NO_DISPATCH
+
+#ifndef DOXYGEN_NO_DETAIL
+namespace formula
+{
+
+template <typename ResultType, typename Geometry>
+ResultType eccentricity_sqr(Geometry const& geometry)
+{
+ return formula_dispatch::eccentricity_sqr<ResultType, Geometry>::apply(geometry);
+}
+
+} // namespace formula
+#endif // DOXYGEN_NO_DETAIL
+
+}} // namespace boost::geometry
+
+#endif // BOOST_GEOMETRY_FORMULAS_ECCENCRICITY_SQR_HPP
diff --git a/boost/geometry/formulas/flattening.hpp b/boost/geometry/formulas/flattening.hpp
new file mode 100644
index 0000000000..f94ead65b0
--- /dev/null
+++ b/boost/geometry/formulas/flattening.hpp
@@ -0,0 +1,70 @@
+// Boost.Geometry
+
+// Copyright (c) 2014-2016 Oracle and/or its affiliates.
+
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_GEOMETRY_FORMULAS_FLATTENING_HPP
+#define BOOST_GEOMETRY_FORMULAS_FLATTENING_HPP
+
+#include <boost/geometry/core/radius.hpp>
+#include <boost/geometry/core/tag.hpp>
+#include <boost/geometry/core/tags.hpp>
+
+#include <boost/geometry/algorithms/not_implemented.hpp>
+
+namespace boost { namespace geometry
+{
+
+#ifndef DOXYGEN_NO_DISPATCH
+namespace formula_dispatch
+{
+
+template <typename ResultType, typename Geometry, typename Tag = typename tag<Geometry>::type>
+struct flattening
+ : not_implemented<Tag>
+{};
+
+template <typename ResultType, typename Geometry>
+struct flattening<ResultType, Geometry, srs_sphere_tag>
+{
+ static inline ResultType apply(Geometry const& /*geometry*/)
+ {
+ return ResultType(0);
+ }
+};
+
+template <typename ResultType, typename Geometry>
+struct flattening<ResultType, Geometry, srs_spheroid_tag>
+{
+ static inline ResultType apply(Geometry const& geometry)
+ {
+ // (a - b) / a
+ return ResultType(get_radius<0>(geometry) - get_radius<2>(geometry))
+ / ResultType(get_radius<0>(geometry));
+ }
+};
+
+} // namespace formula_dispatch
+#endif // DOXYGEN_NO_DISPATCH
+
+#ifndef DOXYGEN_NO_DETAIL
+namespace formula
+{
+
+template <typename ResultType, typename Geometry>
+ResultType flattening(Geometry const& geometry)
+{
+ return formula_dispatch::flattening<ResultType, Geometry>::apply(geometry);
+}
+
+} // namespace formula
+#endif // DOXYGEN_NO_DETAIL
+
+}} // namespace boost::geometry
+
+#endif // BOOST_GEOMETRY_FORMULAS_FLATTENING_HPP
diff --git a/boost/geometry/formulas/geographic.hpp b/boost/geometry/formulas/geographic.hpp
new file mode 100644
index 0000000000..f6feb66633
--- /dev/null
+++ b/boost/geometry/formulas/geographic.hpp
@@ -0,0 +1,457 @@
+// Boost.Geometry
+
+// Copyright (c) 2016, Oracle and/or its affiliates.
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_GEOMETRY_FORMULAS_GEOGRAPHIC_HPP
+#define BOOST_GEOMETRY_FORMULAS_GEOGRAPHIC_HPP
+
+#include <boost/geometry/core/coordinate_system.hpp>
+#include <boost/geometry/core/coordinate_type.hpp>
+#include <boost/geometry/core/access.hpp>
+#include <boost/geometry/core/radian_access.hpp>
+
+#include <boost/geometry/arithmetic/arithmetic.hpp>
+#include <boost/geometry/arithmetic/cross_product.hpp>
+#include <boost/geometry/arithmetic/dot_product.hpp>
+#include <boost/geometry/arithmetic/normalize.hpp>
+
+#include <boost/geometry/formulas/eccentricity_sqr.hpp>
+#include <boost/geometry/formulas/flattening.hpp>
+
+#include <boost/geometry/util/math.hpp>
+#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
+#include <boost/geometry/util/select_coordinate_type.hpp>
+
+namespace boost { namespace geometry {
+
+namespace formula {
+
+template <typename Point3d, typename PointGeo, typename Spheroid>
+inline Point3d geo_to_cart3d(PointGeo const& point_geo, Spheroid const& spheroid)
+{
+ typedef typename coordinate_type<Point3d>::type calc_t;
+
+ calc_t const c1 = 1;
+ calc_t const e_sqr = eccentricity_sqr<calc_t>(spheroid);
+
+ calc_t const lon = get_as_radian<0>(point_geo);
+ calc_t const lat = get_as_radian<1>(point_geo);
+
+ Point3d res;
+
+ calc_t const sin_lat = sin(lat);
+
+ // "unit" spheroid, a = 1
+ calc_t const N = c1 / math::sqrt(c1 - e_sqr * math::sqr(sin_lat));
+ calc_t const N_cos_lat = N * cos(lat);
+
+ set<0>(res, N_cos_lat * cos(lon));
+ set<1>(res, N_cos_lat * sin(lon));
+ set<2>(res, N * (c1 - e_sqr) * sin_lat);
+
+ return res;
+}
+
+template <typename PointGeo, typename Spheroid, typename Point3d>
+inline void geo_to_cart3d(PointGeo const& point_geo, Point3d & result, Point3d & north, Point3d & east, Spheroid const& spheroid)
+{
+ typedef typename coordinate_type<Point3d>::type calc_t;
+
+ calc_t const c1 = 1;
+ calc_t const e_sqr = eccentricity_sqr<calc_t>(spheroid);
+
+ calc_t const lon = get_as_radian<0>(point_geo);
+ calc_t const lat = get_as_radian<1>(point_geo);
+
+ calc_t const sin_lon = sin(lon);
+ calc_t const cos_lon = cos(lon);
+ calc_t const sin_lat = sin(lat);
+ calc_t const cos_lat = cos(lat);
+
+ // "unit" spheroid, a = 1
+ calc_t const N = c1 / math::sqrt(c1 - e_sqr * math::sqr(sin_lat));
+ calc_t const N_cos_lat = N * cos_lat;
+
+ set<0>(result, N_cos_lat * cos_lon);
+ set<1>(result, N_cos_lat * sin_lon);
+ set<2>(result, N * (c1 - e_sqr) * sin_lat);
+
+ set<0>(east, -sin_lon);
+ set<1>(east, cos_lon);
+ set<2>(east, 0);
+
+ set<0>(north, -sin_lat * cos_lon);
+ set<1>(north, -sin_lat * sin_lon);
+ set<2>(north, cos_lat);
+}
+
+template <typename PointGeo, typename Point3d, typename Spheroid>
+inline PointGeo cart3d_to_geo(Point3d const& point_3d, Spheroid const& spheroid)
+{
+ typedef typename coordinate_type<PointGeo>::type coord_t;
+ typedef typename coordinate_type<Point3d>::type calc_t;
+
+ calc_t const c1 = 1;
+ //calc_t const c2 = 2;
+ calc_t const e_sqr = eccentricity_sqr<calc_t>(spheroid);
+
+ calc_t const x = get<0>(point_3d);
+ calc_t const y = get<1>(point_3d);
+ calc_t const z = get<2>(point_3d);
+ calc_t const xy_l = math::sqrt(math::sqr(x) + math::sqr(y));
+
+ calc_t const lonr = atan2(y, x);
+
+ // NOTE: Alternative version
+ // http://www.iag-aig.org/attach/989c8e501d9c5b5e2736955baf2632f5/V60N2_5FT.pdf
+ // calc_t const lonr = c2 * atan2(y, x + xy_l);
+
+ calc_t const latr = atan2(z, (c1 - e_sqr) * xy_l);
+
+ // NOTE: If h is equal to 0 then there is no need to improve value of latitude
+ // because then N_i / (N_i + h_i) = 1
+ // http://www.navipedia.net/index.php/Ellipsoidal_and_Cartesian_Coordinates_Conversion
+
+ PointGeo res;
+
+ set_from_radian<0>(res, lonr);
+ set_from_radian<1>(res, latr);
+
+ coord_t lon = get<0>(res);
+ coord_t lat = get<1>(res);
+
+ math::normalize_spheroidal_coordinates
+ <
+ typename coordinate_system<PointGeo>::type::units,
+ coord_t
+ >(lon, lat);
+
+ set<0>(res, lon);
+ set<1>(res, lat);
+
+ return res;
+}
+
+template <typename Point3d, typename Spheroid>
+inline Point3d projected_to_xy(Point3d const& point_3d, Spheroid const& spheroid)
+{
+ typedef typename coordinate_type<Point3d>::type coord_t;
+
+ // len_xy = sqrt(x^2 + y^2)
+ // r = len_xy - |z / tan(lat)|
+ // assuming h = 0
+ // lat = atan2(z, (1 - e^2) * len_xy);
+ // |z / tan(lat)| = (1 - e^2) * len_xy
+ // r = e^2 * len_xy
+ // x_res = r * cos(lon) = e^2 * len_xy * x / len_xy = e^2 * x
+ // y_res = r * sin(lon) = e^2 * len_xy * y / len_xy = e^2 * y
+
+ coord_t const c0 = 0;
+ coord_t const e_sqr = formula::eccentricity_sqr<coord_t>(spheroid);
+
+ Point3d res;
+
+ set<0>(res, e_sqr * get<0>(point_3d));
+ set<1>(res, e_sqr * get<1>(point_3d));
+ set<2>(res, c0);
+
+ return res;
+}
+
+template <typename Point3d, typename Spheroid>
+inline Point3d projected_to_surface(Point3d const& direction, Spheroid const& spheroid)
+{
+ typedef typename coordinate_type<Point3d>::type coord_t;
+
+ //coord_t const c0 = 0;
+ coord_t const c2 = 2;
+ coord_t const c4 = 4;
+
+ // calculate the point of intersection of a ray and spheroid's surface
+ // the origin is the origin of the coordinate system
+ //(x*x+y*y)/(a*a) + z*z/(b*b) = 1
+ // x = d.x * t
+ // y = d.y * t
+ // z = d.z * t
+ coord_t const dx = get<0>(direction);
+ coord_t const dy = get<1>(direction);
+ coord_t const dz = get<2>(direction);
+
+ //coord_t const a_sqr = math::sqr(get_radius<0>(spheroid));
+ //coord_t const b_sqr = math::sqr(get_radius<2>(spheroid));
+ // "unit" spheroid, a = 1
+ coord_t const a_sqr = 1;
+ coord_t const b_sqr = math::sqr(get_radius<2>(spheroid) / get_radius<0>(spheroid));
+
+ coord_t const param_a = (dx*dx + dy*dy) / a_sqr + dz*dz / b_sqr;
+ coord_t const delta = c4 * param_a;
+ // delta >= 0
+ coord_t const t = math::sqrt(delta) / (c2 * param_a);
+
+ // result = direction * t
+ Point3d result = direction;
+ multiply_value(result, t);
+
+ return result;
+}
+
+template <typename Point3d, typename Spheroid>
+inline bool projected_to_surface(Point3d const& origin, Point3d const& direction, Point3d & result1, Point3d & result2, Spheroid const& spheroid)
+{
+ typedef typename coordinate_type<Point3d>::type coord_t;
+
+ coord_t const c0 = 0;
+ coord_t const c1 = 1;
+ coord_t const c2 = 2;
+ coord_t const c4 = 4;
+
+ // calculate the point of intersection of a ray and spheroid's surface
+ //(x*x+y*y)/(a*a) + z*z/(b*b) = 1
+ // x = o.x + d.x * t
+ // y = o.y + d.y * t
+ // z = o.z + d.z * t
+ coord_t const ox = get<0>(origin);
+ coord_t const oy = get<1>(origin);
+ coord_t const oz = get<2>(origin);
+ coord_t const dx = get<0>(direction);
+ coord_t const dy = get<1>(direction);
+ coord_t const dz = get<2>(direction);
+
+ //coord_t const a_sqr = math::sqr(get_radius<0>(spheroid));
+ //coord_t const b_sqr = math::sqr(get_radius<2>(spheroid));
+ // "unit" spheroid, a = 1
+ coord_t const a_sqr = 1;
+ coord_t const b_sqr = math::sqr(get_radius<2>(spheroid) / get_radius<0>(spheroid));
+
+ coord_t const param_a = (dx*dx + dy*dy) / a_sqr + dz*dz / b_sqr;
+ coord_t const param_b = c2 * ((ox*dx + oy*dy) / a_sqr + oz*dz / b_sqr);
+ coord_t const param_c = (ox*ox + oy*oy) / a_sqr + oz*oz / b_sqr - c1;
+
+ coord_t const delta = math::sqr(param_b) - c4 * param_a*param_c;
+
+ // equals() ?
+ if (delta < c0 || param_a == 0)
+ {
+ return false;
+ }
+
+ // result = origin + direction * t
+
+ coord_t const sqrt_delta = math::sqrt(delta);
+ coord_t const two_a = c2 * param_a;
+
+ coord_t const t1 = (-param_b + sqrt_delta) / two_a;
+ result1 = direction;
+ multiply_value(result1, t1);
+ add_point(result1, origin);
+
+ coord_t const t2 = (-param_b - sqrt_delta) / two_a;
+ result2 = direction;
+ multiply_value(result2, t2);
+ add_point(result2, origin);
+
+ return true;
+}
+
+template <typename Point3d, typename Spheroid>
+inline bool great_elliptic_intersection(Point3d const& a1, Point3d const& a2,
+ Point3d const& b1, Point3d const& b2,
+ Point3d & result,
+ Spheroid const& spheroid)
+{
+ typedef typename coordinate_type<Point3d>::type coord_t;
+
+ coord_t c0 = 0;
+ coord_t c1 = 1;
+
+ Point3d n1 = cross_product(a1, a2);
+ Point3d n2 = cross_product(b1, b2);
+
+ // intersection direction
+ Point3d id = cross_product(n1, n2);
+ coord_t id_len_sqr = dot_product(id, id);
+
+ if (math::equals(id_len_sqr, c0))
+ {
+ return false;
+ }
+
+ // no need to normalize a1 and a2 because the intersection point on
+ // the opposite side of the globe is at the same distance from the origin
+ coord_t cos_a1i = dot_product(a1, id);
+ coord_t cos_a2i = dot_product(a2, id);
+ coord_t gri = math::detail::greatest(cos_a1i, cos_a2i);
+ Point3d neg_id = id;
+ multiply_value(neg_id, -c1);
+ coord_t cos_a1ni = dot_product(a1, neg_id);
+ coord_t cos_a2ni = dot_product(a2, neg_id);
+ coord_t grni = math::detail::greatest(cos_a1ni, cos_a2ni);
+
+ if (gri >= grni)
+ {
+ result = projected_to_surface(id, spheroid);
+ }
+ else
+ {
+ result = projected_to_surface(neg_id, spheroid);
+ }
+
+ return true;
+}
+
+template <typename Point3d1, typename Point3d2>
+static inline int elliptic_side_value(Point3d1 const& origin, Point3d1 const& norm, Point3d2 const& pt)
+{
+ typedef typename coordinate_type<Point3d1>::type calc_t;
+ calc_t c0 = 0;
+
+ // vector oposite to pt - origin
+ // only for the purpose of assigning origin
+ Point3d1 vec = origin;
+ subtract_point(vec, pt);
+
+ calc_t d = dot_product(norm, vec);
+
+ // since the vector is opposite the signs are opposite
+ return math::equals(d, c0) ? 0
+ : d < c0 ? 1
+ : -1; // d > 0
+}
+
+template <typename Point3d, typename Spheroid>
+inline bool planes_spheroid_intersection(Point3d const& o1, Point3d const& n1,
+ Point3d const& o2, Point3d const& n2,
+ Point3d & ip1, Point3d & ip2,
+ Spheroid const& spheroid)
+{
+ typedef typename coordinate_type<Point3d>::type coord_t;
+
+ coord_t c0 = 0;
+ coord_t c1 = 1;
+
+ // Below
+ // n . (p - o) = 0
+ // n . p - n . o = 0
+ // n . p + d = 0
+ // n . p = h
+
+ // intersection direction
+ Point3d id = cross_product(n1, n2);
+
+ if (math::equals(dot_product(id, id), c0))
+ {
+ return false;
+ }
+
+ coord_t dot_n1_n2 = dot_product(n1, n2);
+ coord_t dot_n1_n2_sqr = math::sqr(dot_n1_n2);
+
+ coord_t h1 = dot_product(n1, o1);
+ coord_t h2 = dot_product(n2, o2);
+
+ coord_t denom = c1 - dot_n1_n2_sqr;
+ coord_t C1 = (h1 - h2 * dot_n1_n2) / denom;
+ coord_t C2 = (h2 - h1 * dot_n1_n2) / denom;
+
+ // C1 * n1 + C2 * n2
+ Point3d C1_n1 = n1;
+ multiply_value(C1_n1, C1);
+ Point3d C2_n2 = n2;
+ multiply_value(C2_n2, C2);
+ Point3d io = C1_n1;
+ add_point(io, C2_n2);
+
+ if (! projected_to_surface(io, id, ip1, ip2, spheroid))
+ {
+ return false;
+ }
+
+ return true;
+}
+
+
+template <typename Point3d, typename Spheroid>
+inline void experimental_elliptic_plane(Point3d const& p1, Point3d const& p2,
+ Point3d & v1, Point3d & v2,
+ Point3d & origin, Point3d & normal,
+ Spheroid const& spheroid)
+{
+ typedef typename coordinate_type<Point3d>::type coord_t;
+
+ Point3d xy1 = projected_to_xy(p1, spheroid);
+ Point3d xy2 = projected_to_xy(p2, spheroid);
+
+ // origin = (xy1 + xy2) / 2
+ origin = xy1;
+ add_point(origin, xy2);
+ multiply_value(origin, coord_t(0.5));
+
+ // v1 = p1 - origin
+ v1 = p1;
+ subtract_point(v1, origin);
+ // v2 = p2 - origin
+ v2 = p2;
+ subtract_point(v2, origin);
+
+ normal = cross_product(v1, v2);
+}
+
+template <typename Point3d, typename Spheroid>
+inline void experimental_elliptic_plane(Point3d const& p1, Point3d const& p2,
+ Point3d & origin, Point3d & normal,
+ Spheroid const& spheroid)
+{
+ Point3d v1, v2;
+ experimental_elliptic_plane(p1, p2, v1, v2, origin, normal, spheroid);
+}
+
+template <typename Point3d, typename Spheroid>
+inline bool experimental_elliptic_intersection(Point3d const& a1, Point3d const& a2,
+ Point3d const& b1, Point3d const& b2,
+ Point3d & result,
+ Spheroid const& spheroid)
+{
+ typedef typename coordinate_type<Point3d>::type coord_t;
+
+ coord_t c0 = 0;
+ coord_t c1 = 1;
+
+ Point3d a1v, a2v, o1, n1;
+ experimental_elliptic_plane(a1, a2, a1v, a2v, o1, n1, spheroid);
+ Point3d b1v, b2v, o2, n2;
+ experimental_elliptic_plane(b1, b2, b1v, b2v, o2, n2, spheroid);
+
+ if (! detail::vec_normalize(n1) || ! detail::vec_normalize(n2))
+ {
+ return false;
+ }
+
+ Point3d ip1_s, ip2_s;
+ if (! planes_spheroid_intersection(o1, n1, o2, n2, ip1_s, ip2_s, spheroid))
+ {
+ return false;
+ }
+
+ // NOTE: simplified test, may not work in all cases
+ coord_t dot_a1i1 = dot_product(a1, ip1_s);
+ coord_t dot_a2i1 = dot_product(a2, ip1_s);
+ coord_t gri1 = math::detail::greatest(dot_a1i1, dot_a2i1);
+ coord_t dot_a1i2 = dot_product(a1, ip2_s);
+ coord_t dot_a2i2 = dot_product(a2, ip2_s);
+ coord_t gri2 = math::detail::greatest(dot_a1i2, dot_a2i2);
+
+ result = gri1 >= gri2 ? ip1_s : ip2_s;
+
+ return true;
+}
+
+} // namespace formula
+
+}} // namespace boost::geometry
+
+#endif // BOOST_GEOMETRY_FORMULAS_GEOGRAPHIC_HPP
diff --git a/boost/geometry/formulas/gnomonic_spheroid.hpp b/boost/geometry/formulas/gnomonic_spheroid.hpp
index 3457397b0f..4710c6c063 100644
--- a/boost/geometry/formulas/gnomonic_spheroid.hpp
+++ b/boost/geometry/formulas/gnomonic_spheroid.hpp
@@ -14,12 +14,11 @@
#include <boost/geometry/core/radius.hpp>
-#include <boost/geometry/algorithms/detail/flattening.hpp>
-
#include <boost/geometry/util/condition.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/geometry/formulas/andoyer_inverse.hpp>
+#include <boost/geometry/formulas/flattening.hpp>
#include <boost/geometry/formulas/thomas_inverse.hpp>
#include <boost/geometry/formulas/vincenty_direct.hpp>
#include <boost/geometry/formulas/vincenty_inverse.hpp>
diff --git a/boost/geometry/formulas/sjoberg_intersection.hpp b/boost/geometry/formulas/sjoberg_intersection.hpp
index 03bd4bc97e..92f9e8e78e 100644
--- a/boost/geometry/formulas/sjoberg_intersection.hpp
+++ b/boost/geometry/formulas/sjoberg_intersection.hpp
@@ -20,19 +20,606 @@
#include <boost/geometry/util/condition.hpp>
#include <boost/geometry/util/math.hpp>
-#include <boost/geometry/algorithms/detail/flattening.hpp>
+#include <boost/geometry/formulas/flattening.hpp>
+#include <boost/geometry/formulas/spherical.hpp>
namespace boost { namespace geometry { namespace formula
{
/*!
+\brief The intersection of two great circles as proposed by Sjoberg.
+\see See
+ - [Sjoberg02] Lars E. Sjoberg, Intersections on the sphere and ellipsoid, 2002
+ http://link.springer.com/article/10.1007/s00190-001-0230-9
+*/
+template <typename CT>
+struct sjoberg_intersection_spherical_02
+{
+ // TODO: if it will be used as standalone formula
+ // support segments on equator and endpoints on poles
+
+ static inline bool apply(CT const& lon1, CT const& lat1, CT const& lon_a2, CT const& lat_a2,
+ CT const& lon2, CT const& lat2, CT const& lon_b2, CT const& lat_b2,
+ CT & lon, CT & lat)
+ {
+ CT tan_lat = 0;
+ bool res = apply_alt(lon1, lat1, lon_a2, lat_a2,
+ lon2, lat2, lon_b2, lat_b2,
+ lon, tan_lat);
+
+ if (res)
+ {
+ lat = atan(tan_lat);
+ }
+
+ return res;
+ }
+
+ static inline bool apply_alt(CT const& lon1, CT const& lat1, CT const& lon_a2, CT const& lat_a2,
+ CT const& lon2, CT const& lat2, CT const& lon_b2, CT const& lat_b2,
+ CT & lon, CT & tan_lat)
+ {
+ CT const cos_lon1 = cos(lon1);
+ CT const sin_lon1 = sin(lon1);
+ CT const cos_lon2 = cos(lon2);
+ CT const sin_lon2 = sin(lon2);
+ CT const sin_lat1 = sin(lat1);
+ CT const sin_lat2 = sin(lat2);
+ CT const cos_lat1 = cos(lat1);
+ CT const cos_lat2 = cos(lat2);
+
+ CT const tan_lat_a2 = tan(lat_a2);
+ CT const tan_lat_b2 = tan(lat_b2);
+
+ return apply(lon1, lon_a2, lon2, lon_b2,
+ sin_lon1, cos_lon1, sin_lat1, cos_lat1,
+ sin_lon2, cos_lon2, sin_lat2, cos_lat2,
+ tan_lat_a2, tan_lat_b2,
+ lon, tan_lat);
+ }
+
+private:
+ static inline bool apply(CT const& lon1, CT const& lon_a2, CT const& lon2, CT const& lon_b2,
+ CT const& sin_lon1, CT const& cos_lon1, CT const& sin_lat1, CT const& cos_lat1,
+ CT const& sin_lon2, CT const& cos_lon2, CT const& sin_lat2, CT const& cos_lat2,
+ CT const& tan_lat_a2, CT const& tan_lat_b2,
+ CT & lon, CT & tan_lat)
+ {
+ // NOTE:
+ // cos_lat_ = 0 <=> segment on equator
+ // tan_alpha_ = 0 <=> segment vertical
+
+ CT const tan_lat1 = sin_lat1 / cos_lat1; //tan(lat1);
+ CT const tan_lat2 = sin_lat2 / cos_lat2; //tan(lat2);
+
+ CT const dlon1 = lon_a2 - lon1;
+ CT const sin_dlon1 = sin(dlon1);
+ CT const dlon2 = lon_b2 - lon2;
+ CT const sin_dlon2 = sin(dlon2);
+
+ CT const c0 = 0;
+ bool const is_vertical1 = math::equals(sin_dlon1, c0);
+ bool const is_vertical2 = math::equals(sin_dlon2, c0);
+
+ CT tan_alpha1 = 0;
+ CT tan_alpha2 = 0;
+
+ if (is_vertical1 && is_vertical2)
+ {
+ // circles intersect at one of the poles or are collinear
+ return false;
+ }
+ else if (is_vertical1)
+ {
+ CT const cos_dlon2 = cos(dlon2);
+ CT const tan_alpha2_x = cos_lat2 * tan_lat_b2 - sin_lat2 * cos_dlon2;
+ tan_alpha2 = sin_dlon2 / tan_alpha2_x;
+
+ lon = lon1;
+ }
+ else if (is_vertical2)
+ {
+ CT const cos_dlon1 = cos(dlon1);
+ CT const tan_alpha1_x = cos_lat1 * tan_lat_a2 - sin_lat1 * cos_dlon1;
+ tan_alpha1 = sin_dlon1 / tan_alpha1_x;
+
+ lon = lon2;
+ }
+ else
+ {
+ CT const cos_dlon1 = cos(dlon1);
+ CT const cos_dlon2 = cos(dlon2);
+
+ CT const tan_alpha1_x = cos_lat1 * tan_lat_a2 - sin_lat1 * cos_dlon1;
+ CT const tan_alpha2_x = cos_lat2 * tan_lat_b2 - sin_lat2 * cos_dlon2;
+ tan_alpha1 = sin_dlon1 / tan_alpha1_x;
+ tan_alpha2 = sin_dlon2 / tan_alpha2_x;
+
+ CT const T1 = tan_alpha1 * cos_lat1;
+ CT const T2 = tan_alpha2 * cos_lat2;
+ CT const T1T2 = T1*T2;
+ CT const tan_lon_y = T1 * sin_lon2 - T2 * sin_lon1 + T1T2 * (tan_lat1 * cos_lon1 - tan_lat2 * cos_lon2);
+ CT const tan_lon_x = T1 * cos_lon2 - T2 * cos_lon1 - T1T2 * (tan_lat1 * sin_lon1 - tan_lat2 * sin_lon2);
+
+ lon = atan2(tan_lon_y, tan_lon_x);
+ }
+
+ // choose closer result
+ CT const pi = math::pi<CT>();
+ CT const lon_2 = lon > c0 ? lon - pi : lon + pi;
+ CT const lon_dist1 = (std::max)((std::min)(math::longitude_difference<radian>(lon1, lon),
+ math::longitude_difference<radian>(lon_a2, lon)),
+ (std::min)(math::longitude_difference<radian>(lon2, lon),
+ math::longitude_difference<radian>(lon_b2, lon)));
+ CT const lon_dist2 = (std::max)((std::min)(math::longitude_difference<radian>(lon1, lon_2),
+ math::longitude_difference<radian>(lon_a2, lon_2)),
+ (std::min)(math::longitude_difference<radian>(lon2, lon_2),
+ math::longitude_difference<radian>(lon_b2, lon_2)));
+ if (lon_dist2 < lon_dist1)
+ {
+ lon = lon_2;
+ }
+
+ CT const sin_lon = sin(lon);
+ CT const cos_lon = cos(lon);
+
+ if (math::abs(tan_alpha1) >= math::abs(tan_alpha2)) // pick less vertical segment
+ {
+ CT const sin_dlon_1 = sin_lon * cos_lon1 - cos_lon * sin_lon1;
+ CT const cos_dlon_1 = cos_lon * cos_lon1 + sin_lon * sin_lon1;
+ CT const lat_y_1 = sin_dlon_1 + tan_alpha1 * sin_lat1 * cos_dlon_1;
+ CT const lat_x_1 = tan_alpha1 * cos_lat1;
+ tan_lat = lat_y_1 / lat_x_1;
+ }
+ else
+ {
+ CT const sin_dlon_2 = sin_lon * cos_lon2 - cos_lon * sin_lon2;
+ CT const cos_dlon_2 = cos_lon * cos_lon2 + sin_lon * sin_lon2;
+ CT const lat_y_2 = sin_dlon_2 + tan_alpha2 * sin_lat2 * cos_dlon_2;
+ CT const lat_x_2 = tan_alpha2 * cos_lat2;
+ tan_lat = lat_y_2 / lat_x_2;
+ }
+
+ return true;
+ }
+};
+
+
+/*! Approximation of dLambda_j [Sjoberg07], expanded into taylor series in e^2
+ Maxima script:
+ 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);
+ dL_j(c_j, B_j, B) := -e2 * c_j * dLI_j(c_j, B_j, B);
+ S: taylor(dLI_j(c_j, sinB_j, sinB), e2, 0, 3);
+ assume(c_j < 1);
+ assume(c_j > 0);
+ L1: factor(integrate(sqrt(-x ^ 2 - c_j ^ 2 + 1) / (x ^ 2 + c_j ^ 2 - 1), x));
+ L2: factor(integrate(((x ^ 2 - 1)*sqrt(-x ^ 2 - c_j ^ 2 + 1)) / (x ^ 2 + c_j ^ 2 - 1), x));
+ L3: factor(integrate(((x ^ 4 - 2 * x ^ 2 + 1)*sqrt(-x ^ 2 - c_j ^ 2 + 1)) / (x ^ 2 + c_j ^ 2 - 1), x));
+ 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));
+
+\see See
+ - [Sjoberg07] Lars E. Sjoberg, Geodetic intersection on the ellipsoid, 2007
+ http://link.springer.com/article/10.1007/s00190-007-0204-7
+*/
+template <unsigned int Order, typename CT>
+inline CT sjoberg_d_lambda_e_sqr(CT const& sin_betaj, CT const& sin_beta,
+ CT const& Cj, CT const& sqrt_1_Cj_sqr,
+ CT const& e_sqr)
+{
+ using math::detail::bounded;
+
+ if (Order == 0)
+ {
+ return 0;
+ }
+
+ CT const c1 = 1;
+ CT const c2 = 2;
+
+ CT const asin_B = asin(bounded(sin_beta / sqrt_1_Cj_sqr, -c1, c1));
+ CT const asin_Bj = asin(sin_betaj / sqrt_1_Cj_sqr);
+ CT const L0 = (asin_B - asin_Bj) / c2;
+
+ if (Order == 1)
+ {
+ return -Cj * e_sqr * L0;
+ }
+
+ CT const c0 = 0;
+ CT const c16 = 16;
+
+ CT const X = sin_beta;
+ CT const Xj = sin_betaj;
+ CT const X_sqr = math::sqr(X);
+ CT const Xj_sqr = math::sqr(Xj);
+ CT const Cj_sqr = math::sqr(Cj);
+ CT const Cj_sqr_plus_one = Cj_sqr + c1;
+ CT const one_minus_Cj_sqr = c1 - Cj_sqr;
+ CT const sqrt_Y = math::sqrt(bounded(-X_sqr + one_minus_Cj_sqr, c0));
+ CT const sqrt_Yj = math::sqrt(-Xj_sqr + one_minus_Cj_sqr);
+ CT const L1 = (Cj_sqr_plus_one * (asin_B - asin_Bj) + X * sqrt_Y - Xj * sqrt_Yj) / c16;
+
+ if (Order == 2)
+ {
+ return -Cj * e_sqr * (L0 + e_sqr * L1);
+ }
+
+ CT const c3 = 3;
+ CT const c5 = 5;
+ CT const c128 = 128;
+
+ CT const E = Cj_sqr * (c3 * Cj_sqr + c2) + c3;
+ CT const F = X * (-c2 * X_sqr + c3 * Cj_sqr + c5);
+ CT const Fj = Xj * (-c2 * Xj_sqr + c3 * Cj_sqr + c5);
+ CT const L2 = (E * (asin_B - asin_Bj) + F * sqrt_Y - Fj * sqrt_Yj) / c128;
+
+ if (Order == 3)
+ {
+ return -Cj * e_sqr * (L0 + e_sqr * (L1 + e_sqr * L2));
+ }
+
+ CT const c8 = 8;
+ CT const c9 = 9;
+ CT const c10 = 10;
+ CT const c15 = 15;
+ CT const c24 = 24;
+ CT const c26 = 26;
+ CT const c33 = 33;
+ CT const c6144 = 6144;
+
+ CT const G = Cj_sqr * (Cj_sqr * (Cj_sqr * c15 + c9) + c9) + c15;
+ CT const H = -c10 * Cj_sqr - c26;
+ CT const I = Cj_sqr * (Cj_sqr * c15 + c24) + c33;
+ CT const J = X_sqr * (X * (c8 * X_sqr + H)) + X * I;
+ CT const Jj = Xj_sqr * (Xj * (c8 * Xj_sqr + H)) + Xj * I;
+ CT const L3 = (G * (asin_B - asin_Bj) + J * sqrt_Y - Jj * sqrt_Yj) / c6144;
+
+ // Order 4 and higher
+ return -Cj * e_sqr * (L0 + e_sqr * (L1 + e_sqr * (L2 + e_sqr * L3)));
+}
+
+/*!
+\brief The representation of geodesic as proposed by Sjoberg.
+\see See
+ - [Sjoberg07] Lars E. Sjoberg, Geodetic intersection on the ellipsoid, 2007
+ http://link.springer.com/article/10.1007/s00190-007-0204-7
+ - [Sjoberg12] Lars E. Sjoberg, Solutions to the ellipsoidal Clairaut constant
+ and the inverse geodetic problem by numerical integration, 2012
+ https://www.degruyter.com/view/j/jogs.2012.2.issue-3/v10156-011-0037-4/v10156-011-0037-4.xml
+*/
+template <typename CT, unsigned int Order>
+class sjoberg_geodesic
+{
+ sjoberg_geodesic() {}
+
+ static int sign_C(CT const& alphaj)
+ {
+ CT const c0 = 0;
+ CT const c2 = 2;
+ CT const pi = math::pi<CT>();
+ CT const pi_half = pi / c2;
+
+ return (pi_half < alphaj && alphaj < pi) || (-pi_half < alphaj && alphaj < c0) ? -1 : 1;
+ }
+
+public:
+ sjoberg_geodesic(CT const& lon, CT const& lat, CT const& alpha, CT const& f)
+ : lonj(lon)
+ , latj(lat)
+ , alphaj(alpha)
+ {
+ CT const c0 = 0;
+ CT const c1 = 1;
+ CT const c2 = 2;
+ //CT const pi = math::pi<CT>();
+ //CT const pi_half = pi / c2;
+
+ one_minus_f = c1 - f;
+ e_sqr = f * (c2 - f);
+
+ tan_latj = tan(lat);
+ tan_betaj = one_minus_f * tan_latj;
+ betaj = atan(tan_betaj);
+ sin_betaj = sin(betaj);
+
+ cos_betaj = cos(betaj);
+ sin_alphaj = sin(alphaj);
+ // Clairaut constant (lower-case in the paper)
+ Cj = sign_C(alphaj) * cos_betaj * sin_alphaj;
+ Cj_sqr = math::sqr(Cj);
+ sqrt_1_Cj_sqr = math::sqrt(c1 - Cj_sqr);
+
+ sign_lon_diff = alphaj >= 0 ? 1 : -1; // || alphaj == -pi ?
+ //sign_lon_diff = 1;
+
+ is_on_equator = math::equals(sqrt_1_Cj_sqr, c0);
+ is_Cj_zero = math::equals(Cj, c0);
+
+ t0j = c0;
+ asin_tj_t0j = c0;
+
+ if (! is_Cj_zero)
+ {
+ t0j = sqrt_1_Cj_sqr / Cj;
+ }
+
+ if (! is_on_equator)
+ {
+ //asin_tj_t0j = asin(tan_betaj / t0j);
+ asin_tj_t0j = asin(tan_betaj * Cj / sqrt_1_Cj_sqr);
+ }
+ }
+
+ struct vertex_data
+ {
+ //CT beta0j;
+ CT sin_beta0j;
+ CT dL0j;
+ CT lon0j;
+ };
+
+ vertex_data get_vertex_data() const
+ {
+ CT const c2 = 2;
+ CT const pi = math::pi<CT>();
+ CT const pi_half = pi / c2;
+
+ vertex_data res;
+
+ if (! is_Cj_zero)
+ {
+ //res.beta0j = atan(t0j);
+ //res.sin_beta0j = sin(res.beta0j);
+ res.sin_beta0j = math::sign(t0j) * sqrt_1_Cj_sqr;
+ res.dL0j = d_lambda(res.sin_beta0j);
+ res.lon0j = lonj + sign_lon_diff * (pi_half - asin_tj_t0j + res.dL0j);
+ }
+ else
+ {
+ //res.beta0j = pi_half;
+ //res.sin_beta0j = betaj >= 0 ? 1 : -1;
+ res.sin_beta0j = 1;
+ res.dL0j = 0;
+ res.lon0j = lonj;
+ }
+
+ return res;
+ }
+
+ bool is_sin_beta_ok(CT const& sin_beta) const
+ {
+ CT const c1 = 1;
+ return math::abs(sin_beta / sqrt_1_Cj_sqr) <= c1;
+ }
+
+ bool k_diff(CT const& sin_beta,
+ CT & delta_k) const
+ {
+ if (is_Cj_zero)
+ {
+ delta_k = 0;
+ return true;
+ }
+
+ // beta out of bounds and not close
+ if (! (is_sin_beta_ok(sin_beta)
+ || math::equals(math::abs(sin_beta), sqrt_1_Cj_sqr)) )
+ {
+ return false;
+ }
+
+ // NOTE: beta may be slightly out of bounds here but d_lambda handles that
+ CT const dLj = d_lambda(sin_beta);
+ delta_k = sign_lon_diff * (/*asin_t_t0j*/ - asin_tj_t0j + dLj);
+
+ return true;
+ }
+
+ bool lon_diff(CT const& sin_beta, CT const& t,
+ CT & delta_lon) const
+ {
+ using math::detail::bounded;
+ CT const c1 = 1;
+
+ if (is_Cj_zero)
+ {
+ delta_lon = 0;
+ return true;
+ }
+
+ CT delta_k = 0;
+ if (! k_diff(sin_beta, delta_k))
+ {
+ return false;
+ }
+
+ CT const t_t0j = t / t0j;
+ // NOTE: t may be slightly out of bounds here
+ CT const asin_t_t0j = asin(bounded(t_t0j, -c1, c1));
+ delta_lon = sign_lon_diff * asin_t_t0j + delta_k;
+
+ return true;
+ }
+
+ bool k_diffs(CT const& sin_beta, vertex_data const& vd,
+ CT & delta_k_before, CT & delta_k_behind,
+ bool check_sin_beta = true) const
+ {
+ CT const pi = math::pi<CT>();
+
+ if (is_Cj_zero)
+ {
+ delta_k_before = 0;
+ delta_k_behind = sign_lon_diff * pi;
+ return true;
+ }
+
+ // beta out of bounds and not close
+ if (check_sin_beta
+ && ! (is_sin_beta_ok(sin_beta)
+ || math::equals(math::abs(sin_beta), sqrt_1_Cj_sqr)) )
+ {
+ return false;
+ }
+
+ // NOTE: beta may be slightly out of bounds here but d_lambda handles that
+ CT const dLj = d_lambda(sin_beta);
+ delta_k_before = sign_lon_diff * (/*asin_t_t0j*/ - asin_tj_t0j + dLj);
+
+ // This version require no additional dLj calculation
+ delta_k_behind = sign_lon_diff * (pi /*- asin_t_t0j*/ - asin_tj_t0j + vd.dL0j + (vd.dL0j - dLj));
+
+ // [Sjoberg12]
+ //CT const dL101 = d_lambda(sin_betaj, vd.sin_beta0j);
+ // WARNING: the following call might not work if beta was OoB because only the second argument is bounded
+ //CT const dL_01 = d_lambda(sin_beta, vd.sin_beta0j);
+ //delta_k_behind = sign_lon_diff * (pi /*- asin_t_t0j*/ - asin_tj_t0j + dL101 + dL_01);
+
+ return true;
+ }
+
+ bool lon_diffs(CT const& sin_beta, CT const& t, vertex_data const& vd,
+ CT & delta_lon_before, CT & delta_lon_behind) const
+ {
+ using math::detail::bounded;
+ CT const c1 = 1;
+ CT const pi = math::pi<CT>();
+
+ if (is_Cj_zero)
+ {
+ delta_lon_before = 0;
+ delta_lon_behind = sign_lon_diff * pi;
+ return true;
+ }
+
+ CT delta_k_before = 0, delta_k_behind = 0;
+ if (! k_diffs(sin_beta, vd, delta_k_before, delta_k_behind))
+ {
+ return false;
+ }
+
+ CT const t_t0j = t / t0j;
+ // NOTE: t may be slightly out of bounds here
+ CT const asin_t_t0j = asin(bounded(t_t0j, -c1, c1));
+ CT const sign_asin_t_t0j = sign_lon_diff * asin_t_t0j;
+ delta_lon_before = sign_asin_t_t0j + delta_k_before;
+ delta_lon_behind = -sign_asin_t_t0j + delta_k_behind;
+
+ return true;
+ }
+
+ bool lon(CT const& sin_beta, CT const& t, vertex_data const& vd,
+ CT & lon_before, CT & lon_behind) const
+ {
+ using math::detail::bounded;
+ CT const c1 = 1;
+ CT const pi = math::pi<CT>();
+
+ if (is_Cj_zero)
+ {
+ lon_before = lonj;
+ lon_behind = lonj + sign_lon_diff * pi;
+ return true;
+ }
+
+ if (! (is_sin_beta_ok(sin_beta)
+ || math::equals(math::abs(sin_beta), sqrt_1_Cj_sqr)) )
+ {
+ return false;
+ }
+
+ CT const t_t0j = t / t0j;
+ CT const asin_t_t0j = asin(bounded(t_t0j, -c1, c1));
+ CT const dLj = d_lambda(sin_beta);
+ lon_before = lonj + sign_lon_diff * (asin_t_t0j - asin_tj_t0j + dLj);
+ lon_behind = vd.lon0j + (vd.lon0j - lon_before);
+
+ return true;
+ }
+
+
+ CT lon(CT const& delta_lon) const
+ {
+ return lonj + delta_lon;
+ }
+
+ CT lat(CT const& t) const
+ {
+ // t = tan(beta) = (1-f)tan(lat)
+ return atan(t / one_minus_f);
+ }
+
+ void vertex(CT & lon, CT & lat) const
+ {
+ lon = get_vertex_data().lon0j;
+ if (! is_Cj_zero)
+ {
+ lat = sjoberg_geodesic::lat(t0j);
+ }
+ else
+ {
+ CT const c2 = 2;
+ lat = math::pi<CT>() / c2;
+ }
+ }
+
+ CT lon_of_equator_intersection() const
+ {
+ CT const c0 = 0;
+ CT const dLj = d_lambda(c0);
+ CT const asin_tj_t0j = asin(Cj * tan_betaj / sqrt_1_Cj_sqr);
+ return lonj - asin_tj_t0j + dLj;
+ }
+
+ CT d_lambda(CT const& sin_beta) const
+ {
+ return sjoberg_d_lambda_e_sqr<Order>(sin_betaj, sin_beta, Cj, sqrt_1_Cj_sqr, e_sqr);
+ }
+
+ // [Sjoberg12]
+ /*CT d_lambda(CT const& sin_beta1, CT const& sin_beta2) const
+ {
+ return sjoberg_d_lambda_e_sqr<Order>(sin_beta1, sin_beta2, Cj, sqrt_1_Cj_sqr, e_sqr);
+ }*/
+
+ CT lonj;
+ CT latj;
+ CT alphaj;
+
+ CT one_minus_f;
+ CT e_sqr;
+
+ CT tan_latj;
+ CT tan_betaj;
+ CT betaj;
+ CT sin_betaj;
+ CT cos_betaj;
+ CT sin_alphaj;
+ CT Cj;
+ CT Cj_sqr;
+ CT sqrt_1_Cj_sqr;
+
+ int sign_lon_diff;
+
+ bool is_on_equator;
+ bool is_Cj_zero;
+
+ CT t0j;
+ CT asin_tj_t0j;
+};
+
+
+/*!
\brief The intersection of two geodesics as proposed by Sjoberg.
-\author See
+\see See
- [Sjoberg02] Lars E. Sjoberg, Intersections on the sphere and ellipsoid, 2002
http://link.springer.com/article/10.1007/s00190-001-0230-9
- [Sjoberg07] Lars E. Sjoberg, Geodetic intersection on the ellipsoid, 2007
http://link.springer.com/article/10.1007/s00190-007-0204-7
+ - [Sjoberg12] Lars E. Sjoberg, Solutions to the ellipsoidal Clairaut constant
+ and the inverse geodetic problem by numerical integration, 2012
+ https://www.degruyter.com/view/j/jogs.2012.2.issue-3/v10156-011-0037-4/v10156-011-0037-4.xml
*/
template
<
@@ -42,9 +629,14 @@ template
>
class sjoberg_intersection
{
+ typedef sjoberg_geodesic<CT, Order> geodesic_type;
typedef Inverse<CT, false, true, false, false, false> inverse_type;
typedef typename inverse_type::result_type inverse_result;
+ static bool const enable_02 = true;
+ static int const max_iterations_02 = 10;
+ static int const max_iterations_07 = 20;
+
public:
template <typename T1, typename T2, typename Spheroid>
static inline bool apply(T1 const& lona1, T1 const& lata1,
@@ -63,320 +655,562 @@ public:
CT const lon_b2 = lonb2;
CT const lat_b2 = latb2;
- CT const alpha1 = inverse_type::apply(lon_a1, lat_a1, lon_a2, lat_a2, spheroid).azimuth;
- CT const alpha2 = inverse_type::apply(lon_b1, lat_b1, lon_b2, lat_b2, spheroid).azimuth;
+ inverse_result const res1 = inverse_type::apply(lon_a1, lat_a1, lon_a2, lat_a2, spheroid);
+ inverse_result const res2 = inverse_type::apply(lon_b1, lat_b1, lon_b2, lat_b2, spheroid);
- return apply(lon_a1, lat_a1, alpha1, lon_b1, lat_b1, alpha2, lon, lat, spheroid);
+ return apply(lon_a1, lat_a1, lon_a2, lat_a2, res1.azimuth,
+ lon_b1, lat_b1, lon_b2, lat_b2, res2.azimuth,
+ lon, lat, spheroid);
}
-
+
+ // TODO: Currently may not work correctly if one of the endpoints is the pole
template <typename Spheroid>
- static inline bool apply(CT const& lon1, CT const& lat1, CT const& alpha1,
- CT const& lon2, CT const& lat2, CT const& alpha2,
+ static inline bool apply(CT const& lon_a1, CT const& lat_a1, CT const& lon_a2, CT const& lat_a2, CT const& alpha_a1,
+ CT const& lon_b1, CT const& lat_b1, CT const& lon_b2, CT const& lat_b2, CT const& alpha_b1,
CT & lon, CT & lat,
Spheroid const& spheroid)
{
// coordinates in radians
- // TODO - handle special cases like degenerated segments, equator, poles, etc.
-
CT const c0 = 0;
CT const c1 = 1;
- CT const c2 = 2;
- CT const pi = math::pi<CT>();
- CT const pi_half = pi / c2;
- CT const f = detail::flattening<CT>(spheroid);
+ CT const f = formula::flattening<CT>(spheroid);
CT const one_minus_f = c1 - f;
- CT const e_sqr = f * (c2 - f);
-
- CT const sin_alpha1 = sin(alpha1);
- CT const sin_alpha2 = sin(alpha2);
-
- CT const tan_beta1 = one_minus_f * tan(lat1);
- CT const tan_beta2 = one_minus_f * tan(lat2);
- CT const beta1 = atan(tan_beta1);
- CT const beta2 = atan(tan_beta2);
- CT const cos_beta1 = cos(beta1);
- CT const cos_beta2 = cos(beta2);
- CT const sin_beta1 = sin(beta1);
- CT const sin_beta2 = sin(beta2);
-
- // Clairaut constants (lower-case in the paper)
- int const sign_C1 = math::abs(alpha1) <= pi_half ? 1 : -1;
- int const sign_C2 = math::abs(alpha2) <= pi_half ? 1 : -1;
- // Cj = 1 if on equator
- CT const C1 = sign_C1 * cos_beta1 * sin_alpha1;
- CT const C2 = sign_C2 * cos_beta2 * sin_alpha2;
-
- CT const sqrt_1_C1_sqr = math::sqrt(c1 - math::sqr(C1));
- CT const sqrt_1_C2_sqr = math::sqrt(c1 - math::sqr(C2));
-
- // handle special case: segments on the equator
- bool const on_equator1 = math::equals(sqrt_1_C1_sqr, c0);
- bool const on_equator2 = math::equals(sqrt_1_C2_sqr, c0);
- if (on_equator1 && on_equator2)
+
+ geodesic_type geod1(lon_a1, lat_a1, alpha_a1, f);
+ geodesic_type geod2(lon_b1, lat_b1, alpha_b1, f);
+
+ // Cj = 1 if on equator <=> sqrt_1_Cj_sqr = 0
+ // Cj = 0 if vertical <=> sqrt_1_Cj_sqr = 1
+
+ if (geod1.is_on_equator && geod2.is_on_equator)
{
return false;
}
- else if (on_equator1)
+ else if (geod1.is_on_equator)
{
- CT const dL2 = d_lambda_e_sqr(sin_beta2, c0, C2, sqrt_1_C2_sqr, e_sqr);
- CT const asin_t2_t02 = asin(C2 * tan_beta2 / sqrt_1_C2_sqr);
+ lon = geod2.lon_of_equator_intersection();
lat = c0;
- lon = lon2 - asin_t2_t02 + dL2;
return true;
}
- else if (on_equator2)
+ else if (geod2.is_on_equator)
{
- CT const dL1 = d_lambda_e_sqr(sin_beta1, c0, C1, sqrt_1_C1_sqr, e_sqr);
- CT const asin_t1_t01 = asin(C1 * tan_beta1 / sqrt_1_C1_sqr);
+ lon = geod1.lon_of_equator_intersection();
lat = c0;
- lon = lon1 - asin_t1_t01 + dL1;
return true;
}
- CT const t01 = sqrt_1_C1_sqr / C1;
- CT const t02 = sqrt_1_C2_sqr / C2;
+ // (lon1 - lon2) normalized to (-180, 180]
+ CT const lon1_minus_lon2 = math::longitude_distance_signed<radian>(geod2.lonj, geod1.lonj);
- CT const asin_t1_t01 = asin(tan_beta1 / t01);
- CT const asin_t2_t02 = asin(tan_beta2 / t02);
- CT const t01_t02 = t01 * t02;
- CT const t01_t02_2 = c2 * t01_t02;
- CT const sqr_t01_sqr_t02 = math::sqr(t01) + math::sqr(t02);
+ // vertical segments
+ if (geod1.is_Cj_zero && geod2.is_Cj_zero)
+ {
+ CT const pi = math::pi<CT>();
- CT t = tan_beta1;
- int t_id = 0;
+ // the geodesics are parallel, the intersection point cannot be calculated
+ if ( math::equals(lon1_minus_lon2, c0)
+ || math::equals(lon1_minus_lon2 + (lon1_minus_lon2 < c0 ? pi : -pi), c0) )
+ {
+ return false;
+ }
- // find the initial t using simplified spherical solution
- // though not entirely since the reduced latitudes and azimuths are spheroidal
- // [Sjoberg07]
- CT const k_base = lon1 - lon2 + asin_t2_t02 - asin_t1_t01;
-
+ lon = c0;
+
+ // the geodesics intersect at one of the poles
+ CT const pi_half = pi / CT(2);
+ CT const abs_lat_a1 = math::abs(lat_a1);
+ CT const abs_lat_a2 = math::abs(lat_a2);
+ if (math::equals(abs_lat_a1, abs_lat_a2))
+ {
+ lat = pi_half;
+ }
+ else
+ {
+ // pick the pole closest to one of the points of the first segment
+ CT const& closer_lat = abs_lat_a1 > abs_lat_a2 ? lat_a1 : lat_a2;
+ lat = closer_lat >= 0 ? pi_half : -pi_half;
+ }
+
+ return true;
+ }
+
+ CT lon_sph = 0;
+
+ // Starting tan(beta)
+ CT t = 0;
+
+ /*if (geod1.is_Cj_zero)
{
- CT const K = sin(k_base);
- CT const d1 = sqr_t01_sqr_t02;
- //CT const d2 = t01_t02_2 * math::sqrt(c1 - math::sqr(K));
- CT const d2 = t01_t02_2 * cos(k_base);
- CT const D1 = math::sqrt(d1 - d2);
- CT const D2 = math::sqrt(d1 + d2);
- CT const K_t01_t02 = K * t01_t02;
-
- CT const T1 = K_t01_t02 / D1;
- CT const T2 = K_t01_t02 / D2;
- CT asin_T1_t01 = 0;
- CT asin_T1_t02 = 0;
- CT asin_T2_t01 = 0;
- CT asin_T2_t02 = 0;
-
- // test 4 possible results
- CT l1 = 0, l2 = 0, dl = 0;
- bool found = check_t<0>( T1,
- lon1, asin_T1_t01 = asin(T1 / t01), asin_t1_t01,
- lon2, asin_T1_t02 = asin(T1 / t02), asin_t2_t02,
- t, l1, l2, dl, t_id)
- || check_t<1>(-T1,
- lon1, -asin_T1_t01 , asin_t1_t01,
- lon2, -asin_T1_t02 , asin_t2_t02,
- t, l1, l2, dl, t_id)
- || check_t<2>( T2,
- lon1, asin_T2_t01 = asin(T2 / t01), asin_t1_t01,
- lon2, asin_T2_t02 = asin(T2 / t02), asin_t2_t02,
- t, l1, l2, dl, t_id)
- || check_t<3>(-T2,
- lon1, -asin_T2_t01 , asin_t1_t01,
- lon2, -asin_T2_t02 , asin_t2_t02,
- t, l1, l2, dl, t_id);
-
- boost::ignore_unused(found);
+ CT const k_base = lon1_minus_lon2 + geod2.sign_lon_diff * geod2.asin_tj_t0j;
+ t = sin(k_base) * geod2.t0j;
+ lon_sph = vertical_intersection_longitude(geod1.lonj, lon_b1, lon_b2);
}
+ else if (geod2.is_Cj_zero)
+ {
+ CT const k_base = lon1_minus_lon2 - geod1.sign_lon_diff * geod1.asin_tj_t0j;
+ t = sin(-k_base) * geod1.t0j;
+ lon_sph = vertical_intersection_longitude(geod2.lonj, lon_a1, lon_a2);
+ }
+ else*/
+ {
+ // TODO: Consider using betas instead of latitudes.
+ // Some function calls might be saved this way.
+ CT tan_lat_sph = 0;
+ sjoberg_intersection_spherical_02<CT>::apply_alt(lon_a1, lat_a1, lon_a2, lat_a2,
+ lon_b1, lat_b1, lon_b2, lat_b2,
+ lon_sph, tan_lat_sph);
+ t = one_minus_f * tan_lat_sph; // tan(beta)
+ }
+
+ // TODO: no need to calculate atan here if reduced latitudes were used
+ // instead of latitudes above, in sjoberg_intersection_spherical_02
+ CT const beta = atan(t);
+
+ if (enable_02 && newton_method(geod1, geod2, beta, t, lon1_minus_lon2, lon, lat))
+ {
+ return true;
+ }
+
+ return converge_07(geod1, geod2, beta, t, lon1_minus_lon2, lon_sph, lon, lat);
+ }
+
+private:
+ static inline bool newton_method(geodesic_type const& geod1, geodesic_type const& geod2, // in
+ CT beta, CT t, CT const& lon1_minus_lon2, // in
+ CT & lon, CT & lat) // out
+ {
+ CT const c0 = 0;
+ CT const c1 = 1;
+
+ CT const e_sqr = geod1.e_sqr;
- // [Sjoberg07]
- //int const d2_sign = t_id < 2 ? -1 : 1;
- int const t_sign = (t_id % 2) ? -1 : 1;
- // [Sjoberg02]
- CT const C1_sqr = math::sqr(C1);
- CT const C2_sqr = math::sqr(C2);
-
- CT beta = atan(t);
- CT dL1 = 0, dL2 = 0;
- CT asin_t_t01 = 0;
- CT asin_t_t02 = 0;
+ CT lon1_diff = 0;
+ CT lon2_diff = 0;
- for (int i = 0; i < 10; ++i)
+ CT abs_dbeta_last = 0;
+
+ // [Sjoberg02] converges faster than solution in [Sjoberg07]
+ // Newton-Raphson method
+ for (int i = 0; i < max_iterations_02; ++i)
{
CT const sin_beta = sin(beta);
-
- // integrals approximation
- dL1 = d_lambda_e_sqr(sin_beta1, sin_beta, C1, sqrt_1_C1_sqr, e_sqr);
- dL2 = d_lambda_e_sqr(sin_beta2, sin_beta, C2, sqrt_1_C2_sqr, e_sqr);
-
- // [Sjoberg07]
- /*CT const k = k_base + dL1 - dL2;
- CT const K = sin(k);
- CT const d1 = sqr_t01_sqr_t02;
- //CT const d2 = t01_t02_2 * math::sqrt(c1 - math::sqr(K));
- CT const d2 = t01_t02_2 * cos(k);
- CT const D = math::sqrt(d1 + d2_sign * d2);
- CT const t_new = t_sign * K * t01_t02 / D;
- CT const dt = math::abs(t_new - t);
- t = t_new;
- CT const new_beta = atan(t);
- CT const dbeta = math::abs(new_beta - beta);
- beta = new_beta;*/
-
- // [Sjoberg02] - it converges faster
- // Newton–Raphson method
- asin_t_t01 = asin(t / t01);
- asin_t_t02 = asin(t / t02);
- CT const R1 = asin_t_t01 + dL1;
- CT const R2 = asin_t_t02 + dL2;
CT const cos_beta = cos(beta);
CT const cos_beta_sqr = math::sqr(cos_beta);
CT const G = c1 - e_sqr * cos_beta_sqr;
- CT const f1 = C1 / cos_beta * math::sqrt(G / (cos_beta_sqr - C1_sqr));
- CT const f2 = C2 / cos_beta * math::sqrt(G / (cos_beta_sqr - C2_sqr));
- CT const abs_f1 = math::abs(f1);
- CT const abs_f2 = math::abs(f2);
- CT const dbeta = t_sign * (k_base - R2 + R1) / (abs_f1 + abs_f2);
-
- if (math::equals(dbeta, CT(0)))
+
+ CT f1 = 0;
+ CT f2 = 0;
+
+ if (!geod1.is_Cj_zero)
+ {
+ bool is_beta_ok = geod1.lon_diff(sin_beta, t, lon1_diff);
+
+ if (is_beta_ok)
+ {
+ CT const H = cos_beta_sqr - geod1.Cj_sqr;
+ f1 = geod1.Cj / cos_beta * math::sqrt(G / H);
+ }
+ else
+ {
+ return false;
+ }
+ }
+
+ if (!geod2.is_Cj_zero)
+ {
+ bool is_beta_ok = geod2.lon_diff(sin_beta, t, lon2_diff);
+
+ if (is_beta_ok)
+ {
+ CT const H = cos_beta_sqr - geod2.Cj_sqr;
+ f2 = geod2.Cj / cos_beta * math::sqrt(G / H);
+ }
+ else
+ {
+ return false;
+ }
+ }
+
+ // NOTE: Things may go wrong if the IP is near the vertex
+ // 1. May converge into the wrong direction (from the other way around).
+ // This happens when the starting point is on the other side than the vertex
+ // 2. During converging may "jump" into the other side of the vertex.
+ // In this case sin_beta/sqrt_1_Cj_sqr and t/t0j is not in [-1, 1]
+ // 3. f1-f2 may be 0 which means that the intermediate point is on the vertex
+ // In this case it's not possible to check if this is the correct result
+
+ CT const dbeta_denom = f1 - f2;
+ //CT const dbeta_denom = math::abs(f1) + math::abs(f2);
+
+ if (math::equals(dbeta_denom, c0))
+ {
+ return false;
+ }
+
+ // The sign of dbeta is changed WRT [Sjoberg02]
+ CT const dbeta = (lon1_minus_lon2 + lon1_diff - lon2_diff) / dbeta_denom;
+
+ CT const abs_dbeta = math::abs(dbeta);
+ if (i > 0 && abs_dbeta > abs_dbeta_last)
{
+ // The algorithm is not converging
+ // The intersection may be on the other side of the vertex
+ return false;
+ }
+ abs_dbeta_last = abs_dbeta;
+
+ if (math::equals(dbeta, c0))
+ {
+ // Result found
break;
}
+ // Because the sign of dbeta is changed WRT [Sjoberg02] dbeta is subtracted here
beta = beta - dbeta;
+
t = tan(beta);
}
-
- // t = tan(beta) = (1-f)tan(lat)
- lat = atan(t / one_minus_f);
- CT const l1 = lon1 + asin_t_t01 - asin_t1_t01 + dL1;
- //CT const l2 = lon2 + asin_t_t02 - asin_t2_t02 + dL2;
- lon = l1;
+ lat = geod1.lat(t);
+ // NOTE: if Cj is 0 then the result is lonj or lonj+180
+ lon = ! geod1.is_Cj_zero
+ ? geod1.lon(lon1_diff)
+ : geod2.lon(lon2_diff);
return true;
}
-private:
- /*! Approximation of dLambda_j [Sjoberg07], expanded into taylor series in e^2
- Maxima script:
- 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);
- dL_j(c_j, B_j, B) := -e2 * c_j * dLI_j(c_j, B_j, B);
- S: taylor(dLI_j(c_j, sinB_j, sinB), e2, 0, 3);
- assume(c_j < 1);
- assume(c_j > 0);
- L1: factor(integrate(sqrt(-x ^ 2 - c_j ^ 2 + 1) / (x ^ 2 + c_j ^ 2 - 1), x));
- L2: factor(integrate(((x ^ 2 - 1)*sqrt(-x ^ 2 - c_j ^ 2 + 1)) / (x ^ 2 + c_j ^ 2 - 1), x));
- L3: factor(integrate(((x ^ 4 - 2 * x ^ 2 + 1)*sqrt(-x ^ 2 - c_j ^ 2 + 1)) / (x ^ 2 + c_j ^ 2 - 1), x));
- 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));
- */
- static inline CT d_lambda_e_sqr(CT const& sin_betaj, CT const& sin_beta,
- CT const& Cj, CT const& sqrt_1_Cj_sqr,
- CT const& e_sqr)
- {
- if (Order == 0)
- {
- return 0;
- }
+ struct geodesics_type
+ {
+ geodesics_type(geodesic_type const& g1, geodesic_type const& g2)
+ : geod1(g1)
+ , geod2(g2)
+ , vertex1(geod1.get_vertex_data())
+ , vertex2(geod2.get_vertex_data())
+ {}
- CT const c2 = 2;
-
- CT const asin_B = asin(sin_beta / sqrt_1_Cj_sqr);
- CT const asin_Bj = asin(sin_betaj / sqrt_1_Cj_sqr);
- CT const L0 = (asin_B - asin_Bj) / c2;
+ geodesic_type const& geod1;
+ geodesic_type const& geod2;
+ typename geodesic_type::vertex_data vertex1;
+ typename geodesic_type::vertex_data vertex2;
+ };
- if (Order == 1)
+ struct converge_07_result
+ {
+ converge_07_result()
+ : lon1(0), lon2(0), k1_diff(0), k2_diff(0), t1(0), t2(0)
+ {}
+
+ CT lon1, lon2;
+ CT k1_diff, k2_diff;
+ CT t1, t2;
+ };
+
+ static inline bool converge_07(geodesic_type const& geod1, geodesic_type const& geod2,
+ CT beta, CT t,
+ CT const& lon1_minus_lon2, CT const& lon_sph,
+ CT & lon, CT & lat)
+ {
+ //CT const c0 = 0;
+ //CT const c1 = 1;
+ //CT const c2 = 2;
+ //CT const pi = math::pi<CT>();
+
+ geodesics_type geodesics(geod1, geod2);
+ converge_07_result result;
+
+ // calculate first pair of longitudes
+ if (!converge_07_step_one(CT(sin(beta)), t, lon1_minus_lon2, geodesics, lon_sph, result, false))
{
- return -Cj * e_sqr * L0;
+ return false;
}
- CT const c1 = 1;
- CT const c16 = 16;
+ int t_direction = 0;
- CT const X = sin_beta;
- CT const Xj = sin_betaj;
- CT const Cj_sqr = math::sqr(Cj);
- CT const Cj_sqr_plus_one = Cj_sqr + c1;
- CT const one_minus_Cj_sqr = c1 - Cj_sqr;
- CT const sqrt_Y = math::sqrt(-math::sqr(X) + one_minus_Cj_sqr);
- CT const sqrt_Yj = math::sqrt(-math::sqr(Xj) + one_minus_Cj_sqr);
- CT const L1 = (Cj_sqr_plus_one * (asin_B - asin_Bj) + X * sqrt_Y - Xj * sqrt_Yj) / c16;
+ CT lon_diff_prev = math::longitude_difference<radian>(result.lon1, result.lon2);
- if (Order == 2)
+ // [Sjoberg07]
+ for (int i = 2; i < max_iterations_07; ++i)
{
- return -Cj * e_sqr * (L0 + e_sqr * L1);
+ // pick t candidates from previous result based on dir
+ CT t_cand1 = result.t1;
+ CT t_cand2 = result.t2;
+ // if direction is 0 the closer one is the first
+ if (t_direction < 0)
+ {
+ t_cand1 = (std::min)(result.t1, result.t2);
+ t_cand2 = (std::max)(result.t1, result.t2);
+ }
+ else if (t_direction > 0)
+ {
+ t_cand1 = (std::max)(result.t1, result.t2);
+ t_cand2 = (std::min)(result.t1, result.t2);
+ }
+ else
+ {
+ t_direction = t_cand1 < t_cand2 ? -1 : 1;
+ }
+
+ CT t1 = t;
+ CT beta1 = beta;
+ // check if the further calculation is needed
+ if (converge_07_update(t1, beta1, t_cand1))
+ {
+ break;
+ }
+
+ bool try_t2 = false;
+ converge_07_result result_curr;
+ if (converge_07_step_one(CT(sin(beta1)), t1, lon1_minus_lon2, geodesics, lon_sph, result_curr))
+ {
+ CT const lon_diff1 = math::longitude_difference<radian>(result_curr.lon1, result_curr.lon2);
+ if (lon_diff_prev > lon_diff1)
+ {
+ t = t1;
+ beta = beta1;
+ lon_diff_prev = lon_diff1;
+ result = result_curr;
+ }
+ else if (t_cand1 != t_cand2)
+ {
+ try_t2 = true;
+ }
+ else
+ {
+ // the result is not fully correct but it won't be more accurate
+ break;
+ }
+ }
+ // ! converge_07_step_one
+ else
+ {
+ if (t_cand1 != t_cand2)
+ {
+ try_t2 = true;
+ }
+ else
+ {
+ return false;
+ }
+ }
+
+
+ if (try_t2)
+ {
+ CT t2 = t;
+ CT beta2 = beta;
+ // check if the further calculation is needed
+ if (converge_07_update(t2, beta2, t_cand2))
+ {
+ break;
+ }
+
+ if (! converge_07_step_one(CT(sin(beta2)), t2, lon1_minus_lon2, geodesics, lon_sph, result_curr))
+ {
+ return false;
+ }
+
+ CT const lon_diff2 = math::longitude_difference<radian>(result_curr.lon1, result_curr.lon2);
+ if (lon_diff_prev > lon_diff2)
+ {
+ t_direction *= -1;
+ t = t2;
+ beta = beta2;
+ lon_diff_prev = lon_diff2;
+ result = result_curr;
+ }
+ else
+ {
+ // the result is not fully correct but it won't be more accurate
+ break;
+ }
+ }
}
- CT const c3 = 3;
- CT const c5 = 5;
- CT const c128 = 128;
+ lat = geod1.lat(t);
+ lon = ! geod1.is_Cj_zero ? result.lon1 : result.lon2;
+ math::normalize_longitude<radian>(lon);
+
+ return true;
+ }
+
+ static inline bool converge_07_update(CT & t, CT & beta, CT const& t_new)
+ {
+ CT const c0 = 0;
+
+ CT const beta_new = atan(t_new);
+ CT const dbeta = beta_new - beta;
+ beta = beta_new;
+ t = t_new;
- CT const E = Cj_sqr * (c3 * Cj_sqr + c2) + c3;
- CT const X_sqr = math::sqr(X);
- CT const Xj_sqr = math::sqr(Xj);
- CT const F = X * (-c2 * X_sqr + c3 * Cj_sqr + c5);
- CT const Fj = Xj * (-c2 * Xj_sqr + c3 * Cj_sqr + c5);
- CT const L2 = (E * (asin_B - asin_Bj) + F * sqrt_Y - Fj * sqrt_Yj) / c128;
+ return math::equals(dbeta, c0);
+ }
- if (Order == 3)
+ static inline CT const& pick_t(CT const& t1, CT const& t2, int direction)
+ {
+ return direction < 0 ? (std::min)(t1, t2) : (std::max)(t1, t2);
+ }
+
+ static inline bool converge_07_step_one(CT const& sin_beta,
+ CT const& t,
+ CT const& lon1_minus_lon2,
+ geodesics_type const& geodesics,
+ CT const& lon_sph,
+ converge_07_result & result,
+ bool check_sin_beta = true)
+ {
+ bool ok = converge_07_one_geod(sin_beta, t, geodesics.geod1, geodesics.vertex1, lon_sph,
+ result.lon1, result.k1_diff, check_sin_beta)
+ && converge_07_one_geod(sin_beta, t, geodesics.geod2, geodesics.vertex2, lon_sph,
+ result.lon2, result.k2_diff, check_sin_beta);
+
+ if (!ok)
{
- return -Cj * e_sqr * (L0 + e_sqr * (L1 + e_sqr * L2));
+ return false;
}
- CT const c8 = 8;
- CT const c9 = 9;
- CT const c10 = 10;
- CT const c15 = 15;
- CT const c24 = 24;
- CT const c26 = 26;
- CT const c33 = 33;
- CT const c6144 = 6144;
+ CT const k = lon1_minus_lon2 + result.k1_diff - result.k2_diff;
- CT const G = Cj_sqr * (Cj_sqr * (Cj_sqr * c15 + c9) + c9) + c15;
- CT const H = -c10 * Cj_sqr - c26;
- CT const I = Cj_sqr * (Cj_sqr * c15 + c24) + c33;
- CT const J = X_sqr * (X * (c8 * X_sqr + H)) + X * I;
- CT const Jj = Xj_sqr * (Xj * (c8 * Xj_sqr + H)) + Xj * I;
- CT const L3 = (G * (asin_B - asin_Bj) + J * sqrt_Y - Jj * sqrt_Yj) / c6144;
+ // get 2 possible ts one lesser and one greater than t
+ // t1 is the closer one
+ calc_ts(t, k, geodesics.geod1, geodesics.geod2, result.t1, result.t2);
- // Order 4 and higher
- return -Cj * e_sqr * (L0 + e_sqr * (L1 + e_sqr * (L2 + e_sqr * L3)));
+ return true;
}
- static inline CT fj(CT const& cos_beta, CT const& cos2_beta, CT const& Cj, CT const& e_sqr)
+ static inline bool converge_07_one_geod(CT const& sin_beta, CT const& t,
+ geodesic_type const& geod,
+ typename geodesic_type::vertex_data const& vertex,
+ CT const& lon_sph,
+ CT & lon, CT & k_diff,
+ bool check_sin_beta)
{
+ using math::detail::bounded;
CT const c1 = 1;
- CT const Cj_sqr = math::sqr(Cj);
- return Cj / cos_beta * math::sqrt((c1 - e_sqr * cos2_beta) / (cos2_beta - Cj_sqr));
+
+ CT k_diff_before = 0;
+ CT k_diff_behind = 0;
+
+ bool is_beta_ok = geod.k_diffs(sin_beta, vertex, k_diff_before, k_diff_behind, check_sin_beta);
+
+ if (! is_beta_ok)
+ {
+ return false;
+ }
+
+ CT const asin_t_t0j = ! geod.is_Cj_zero ? asin(bounded(t / geod.t0j, -c1, c1)) : 0;
+ CT const sign_asin_t_t0j = geod.sign_lon_diff * asin_t_t0j;
+
+ CT const lon_before = geod.lonj + sign_asin_t_t0j + k_diff_before;
+ CT const lon_behind = geod.lonj - sign_asin_t_t0j + k_diff_behind;
+
+ CT const lon_dist_before = math::longitude_distance_signed<radian>(lon_before, lon_sph);
+ CT const lon_dist_behind = math::longitude_distance_signed<radian>(lon_behind, lon_sph);
+ if (math::abs(lon_dist_before) <= math::abs(lon_dist_behind))
+ {
+ k_diff = k_diff_before;
+ lon = lon_before;
+ }
+ else
+ {
+ k_diff = k_diff_behind;
+ lon = lon_behind;
+ }
+
+ return true;
}
- template <int TId>
- static inline bool check_t(CT const& t,
- CT const& lon_a1, CT const& asin_t_t01, CT const& asin_t1_t01,
- CT const& lon_b1, CT const& asin_t_t02, CT const& asin_t2_t02,
- CT & current_t, CT & current_lon1, CT & current_lon2, CT & current_dlon,
- int & t_id)
+ static inline void calc_ts(CT const& t, CT const& k,
+ geodesic_type const& geod1, geodesic_type const& geod2,
+ CT & t1, CT& t2)
{
- CT const lon1 = lon_a1 + asin_t_t01 - asin_t1_t01;
- CT const lon2 = lon_b1 + asin_t_t02 - asin_t2_t02;
+ CT const c1 = 1;
+ CT const c2 = 2;
- // TODO - true angle difference
- CT const dlon = math::abs(lon2 - lon1);
+ CT const K = sin(k);
- bool are_equal = math::equals(dlon, CT(0));
-
- if ((TId == 0) || are_equal || dlon < current_dlon)
+ BOOST_GEOMETRY_ASSERT(!geod1.is_Cj_zero || !geod2.is_Cj_zero);
+ if (geod1.is_Cj_zero)
+ {
+ t1 = K * geod2.t0j;
+ t2 = -t1;
+ }
+ else if (geod2.is_Cj_zero)
{
- current_t = t;
- current_lon1 = lon1;
- current_lon2 = lon2;
- current_dlon = dlon;
- t_id = TId;
+ t1 = -K * geod1.t0j;
+ t2 = -t1;
}
+ else
+ {
+ CT const A = math::sqr(geod1.t0j) + math::sqr(geod2.t0j);
+ CT const B = c2 * geod1.t0j * geod2.t0j * math::sqrt(c1 - math::sqr(K));
+
+ CT const K_t01_t02 = K * geod1.t0j * geod2.t0j;
+ CT const D1 = math::sqrt(A + B);
+ CT const D2 = math::sqrt(A - B);
+ CT const t_new1 = K_t01_t02 / D1;
+ CT const t_new2 = K_t01_t02 / D2;
+ CT const t_new3 = -t_new1;
+ CT const t_new4 = -t_new2;
- return are_equal;
+ // Pick 2 nearest t_new, one greater and one lesser than current t
+ CT const abs_t_new1 = math::abs(t_new1);
+ CT const abs_t_new2 = math::abs(t_new2);
+ CT const abs_t_max = (std::max)(abs_t_new1, abs_t_new2);
+ t1 = -abs_t_max; // lesser
+ t2 = abs_t_max; // greater
+ if (t1 < t)
+ {
+ if (t_new1 < t && t_new1 > t1)
+ t1 = t_new1;
+ if (t_new2 < t && t_new2 > t1)
+ t1 = t_new2;
+ if (t_new3 < t && t_new3 > t1)
+ t1 = t_new3;
+ if (t_new4 < t && t_new4 > t1)
+ t1 = t_new4;
+ }
+ if (t2 > t)
+ {
+ if (t_new1 > t && t_new1 < t2)
+ t2 = t_new1;
+ if (t_new2 > t && t_new2 < t2)
+ t2 = t_new2;
+ if (t_new3 > t && t_new3 < t2)
+ t2 = t_new3;
+ if (t_new4 > t && t_new4 < t2)
+ t2 = t_new4;
+ }
+ }
+
+ // the first one is the closer one
+ if (math::abs(t - t2) < math::abs(t - t1))
+ {
+ std::swap(t2, t1);
+ }
+ }
+
+ static inline CT fj(CT const& cos_beta, CT const& cos2_beta, CT const& Cj, CT const& e_sqr)
+ {
+ CT const c1 = 1;
+ CT const Cj_sqr = math::sqr(Cj);
+ return Cj / cos_beta * math::sqrt((c1 - e_sqr * cos2_beta) / (cos2_beta - Cj_sqr));
}
+
+ /*static inline CT vertical_intersection_longitude(CT const& ip_lon, CT const& seg_lon1, CT const& seg_lon2)
+ {
+ CT const c0 = 0;
+ CT const lon_2 = ip_lon > c0 ? ip_lon - pi : ip_lon + pi;
+
+ return (std::min)(math::longitude_difference<radian>(ip_lon, seg_lon1),
+ math::longitude_difference<radian>(ip_lon, seg_lon2))
+ <=
+ (std::min)(math::longitude_difference<radian>(lon_2, seg_lon1),
+ math::longitude_difference<radian>(lon_2, seg_lon2))
+ ? ip_lon : lon_2;
+ }*/
};
}}} // namespace boost::geometry::formula
diff --git a/boost/geometry/formulas/spherical.hpp b/boost/geometry/formulas/spherical.hpp
index 2195bbbe10..ff24c51a88 100644
--- a/boost/geometry/formulas/spherical.hpp
+++ b/boost/geometry/formulas/spherical.hpp
@@ -1,6 +1,7 @@
// Boost.Geometry
// Copyright (c) 2016, Oracle and/or its affiliates.
+// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
// Use, modification and distribution is subject to the Boost Software License,
@@ -27,38 +28,67 @@ namespace boost { namespace geometry {
namespace formula {
+template <typename T>
+struct result_spherical
+{
+ result_spherical()
+ : azimuth(0)
+ , reverse_azimuth(0)
+ {}
+
+ T azimuth;
+ T reverse_azimuth;
+};
+
+template <typename T>
+static inline void sph_to_cart3d(T const& lon, T const& lat, T & x, T & y, T & z)
+{
+ T const cos_lat = cos(lat);
+ x = cos_lat * cos(lon);
+ y = cos_lat * sin(lon);
+ z = sin(lat);
+}
+
template <typename Point3d, typename PointSph>
static inline Point3d sph_to_cart3d(PointSph const& point_sph)
{
typedef typename coordinate_type<Point3d>::type calc_t;
- Point3d res;
-
- calc_t lon = get_as_radian<0>(point_sph);
- calc_t lat = get_as_radian<1>(point_sph);
+ calc_t const lon = get_as_radian<0>(point_sph);
+ calc_t const lat = get_as_radian<1>(point_sph);
+ calc_t x, y, z;
+ sph_to_cart3d(lon, lat, x, y, z);
- calc_t const cos_lat = cos(lat);
- set<0>(res, cos_lat * cos(lon));
- set<1>(res, cos_lat * sin(lon));
- set<2>(res, sin(lat));
+ Point3d res;
+ set<0>(res, x);
+ set<1>(res, y);
+ set<2>(res, z);
return res;
}
+template <typename T>
+static inline void cart3d_to_sph(T const& x, T const& y, T const& z, T & lon, T & lat)
+{
+ lon = atan2(y, x);
+ lat = asin(z);
+}
+
template <typename PointSph, typename Point3d>
static inline PointSph cart3d_to_sph(Point3d const& point_3d)
{
typedef typename coordinate_type<PointSph>::type coord_t;
typedef typename coordinate_type<Point3d>::type calc_t;
- PointSph res;
-
calc_t const x = get<0>(point_3d);
calc_t const y = get<1>(point_3d);
calc_t const z = get<2>(point_3d);
+ calc_t lonr, latr;
+ cart3d_to_sph(x, y, z, lonr, latr);
- set_from_radian<0>(res, atan2(y, x));
- set_from_radian<1>(res, asin(z));
+ PointSph res;
+ set_from_radian<0>(res, lonr);
+ set_from_radian<1>(res, latr);
coord_t lon = get<0>(res);
coord_t lat = get<1>(res);
@@ -89,6 +119,104 @@ static inline int sph_side_value(Point3d1 const& norm, Point3d2 const& pt)
: -1; // d < 0
}
+template <typename CT, bool ReverseAzimuth, typename T1, typename T2>
+static inline result_spherical<CT> spherical_azimuth(T1 const& lon1,
+ T1 const& lat1,
+ T2 const& lon2,
+ T2 const& lat2)
+{
+ typedef result_spherical<CT> result_type;
+ result_type result;
+
+ // http://williams.best.vwh.net/avform.htm#Crs
+ // https://en.wikipedia.org/wiki/Great-circle_navigation
+ CT dlon = lon2 - lon1;
+
+ // An optimization which should kick in often for Boxes
+ //if ( math::equals(dlon, ReturnType(0)) )
+ //if ( get<0>(p1) == get<0>(p2) )
+ //{
+ // return - sin(get_as_radian<1>(p1)) * cos_p2lat);
+ //}
+
+ CT const cos_dlon = cos(dlon);
+ CT const sin_dlon = sin(dlon);
+ CT const cos_lat1 = cos(lat1);
+ CT const cos_lat2 = cos(lat2);
+ CT const sin_lat1 = sin(lat1);
+ CT const sin_lat2 = sin(lat2);
+
+ {
+ // "An alternative formula, not requiring the pre-computation of d"
+ // In the formula below dlon is used as "d"
+ CT const y = sin_dlon * cos_lat2;
+ CT const x = cos_lat1 * sin_lat2 - sin_lat1 * cos_lat2 * cos_dlon;
+ result.azimuth = atan2(y, x);
+ }
+
+ if (ReverseAzimuth)
+ {
+ CT const y = sin_dlon * cos_lat1;
+ CT const x = sin_lat2 * cos_lat1 * cos_dlon - cos_lat2 * sin_lat1;
+ result.reverse_azimuth = atan2(y, x);
+ }
+
+ return result;
+}
+
+template <typename ReturnType, typename T1, typename T2>
+inline ReturnType spherical_azimuth(T1 const& lon1, T1 const& lat1,
+ T2 const& lon2, T2 const& lat2)
+{
+ return spherical_azimuth<ReturnType, false>(lon1, lat1, lon2, lat2).azimuth;
+}
+
+template <typename T>
+inline T spherical_azimuth(T const& lon1, T const& lat1, T const& lon2, T const& lat2)
+{
+ return spherical_azimuth<T, false>(lon1, lat1, lon2, lat2).azimuth;
+}
+
+template <typename T>
+inline int azimuth_side_value(T const& azi_a1_p, T const& azi_a1_a2)
+{
+ T const pi = math::pi<T>();
+ T const two_pi = math::two_pi<T>();
+
+ // instead of the formula from XTD
+ //calc_t a_diff = asin(sin(azi_a1_p - azi_a1_a2));
+
+ T a_diff = azi_a1_p - azi_a1_a2;
+ // normalize, angle in [-pi, pi]
+ while (a_diff > pi)
+ a_diff -= two_pi;
+ while (a_diff < -pi)
+ a_diff += two_pi;
+
+ // NOTE: in general it shouldn't be required to support the pi/-pi case
+ // because in non-cartesian systems it makes sense to check the side
+ // only "between" the endpoints.
+ // However currently the winding strategy calls the side strategy
+ // for vertical segments to check if the point is "between the endpoints.
+ // This could be avoided since the side strategy is not required for that
+ // because meridian is the shortest path. So a difference of
+ // longitudes would be sufficient (of course normalized to [-pi, pi]).
+
+ // NOTE: with the above said, the pi/-pi check is temporary
+ // however in case if this was required
+ // the geodesics on ellipsoid aren't "symmetrical"
+ // therefore instead of comparing a_diff to pi and -pi
+ // one should probably use inverse azimuths and compare
+ // the difference to 0 as well
+
+ // positive azimuth is on the right side
+ return math::equals(a_diff, 0)
+ || math::equals(a_diff, pi)
+ || math::equals(a_diff, -pi) ? 0
+ : a_diff > 0 ? -1 // right
+ : 1; // left
+}
+
} // namespace formula
}} // namespace boost::geometry
diff --git a/boost/geometry/formulas/thomas_direct.hpp b/boost/geometry/formulas/thomas_direct.hpp
index f8a7f83943..f208167cf5 100644
--- a/boost/geometry/formulas/thomas_direct.hpp
+++ b/boost/geometry/formulas/thomas_direct.hpp
@@ -20,9 +20,8 @@
#include <boost/geometry/util/condition.hpp>
#include <boost/geometry/util/math.hpp>
-#include <boost/geometry/algorithms/detail/flattening.hpp>
-
#include <boost/geometry/formulas/differential_quantities.hpp>
+#include <boost/geometry/formulas/flattening.hpp>
#include <boost/geometry/formulas/result_direct.hpp>
@@ -82,7 +81,7 @@ public:
CT const a = CT(get_radius<0>(spheroid));
CT const b = CT(get_radius<2>(spheroid));
- CT const f = detail::flattening<CT>(spheroid);
+ CT const f = formula::flattening<CT>(spheroid);
CT const one_minus_f = c1 - f;
CT const pi = math::pi<CT>();
diff --git a/boost/geometry/formulas/thomas_inverse.hpp b/boost/geometry/formulas/thomas_inverse.hpp
index d68c9de054..0853a36980 100644
--- a/boost/geometry/formulas/thomas_inverse.hpp
+++ b/boost/geometry/formulas/thomas_inverse.hpp
@@ -20,9 +20,8 @@
#include <boost/geometry/util/condition.hpp>
#include <boost/geometry/util/math.hpp>
-#include <boost/geometry/algorithms/detail/flattening.hpp>
-
#include <boost/geometry/formulas/differential_quantities.hpp>
+#include <boost/geometry/formulas/flattening.hpp>
#include <boost/geometry/formulas/result_inverse.hpp>
@@ -78,7 +77,7 @@ public:
CT const c4 = 4;
CT const pi_half = math::pi<CT>() / c2;
- CT const f = detail::flattening<CT>(spheroid);
+ CT const f = formula::flattening<CT>(spheroid);
CT const one_minus_f = c1 - f;
// CT const tan_theta1 = one_minus_f * tan(lat1);
diff --git a/boost/geometry/formulas/vertex_latitude.hpp b/boost/geometry/formulas/vertex_latitude.hpp
new file mode 100644
index 0000000000..755067b08d
--- /dev/null
+++ b/boost/geometry/formulas/vertex_latitude.hpp
@@ -0,0 +1,148 @@
+// Boost.Geometry
+
+// Copyright (c) 2016-2017 Oracle and/or its affiliates.
+
+// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_GEOMETRY_FORMULAS_MAXIMUM_LATITUDE_HPP
+#define BOOST_GEOMETRY_FORMULAS_MAXIMUM_LATITUDE_HPP
+
+#include <boost/geometry/core/srs.hpp>
+#include <boost/geometry/formulas/flattening.hpp>
+#include <boost/geometry/formulas/spherical.hpp>
+#include <boost/mpl/assert.hpp>
+
+namespace boost { namespace geometry { namespace formula
+{
+
+/*!
+\brief Algorithm to compute the vertex latitude of a geodesic segment. Vertex is
+a point on the geodesic that maximizes (or minimizes) the latitude.
+\author See
+ [Wood96] Wood - Vertex Latitudes on Ellipsoid Geodesics, SIAM Rev., 38(4),
+ 637–644, 1996
+*/
+
+template <typename CT>
+class vertex_latitude_on_sphere
+{
+
+public:
+ template<typename T1, typename T2>
+ static inline CT apply(T1 const& lat1,
+ T2 const& alp1)
+ {
+ return std::acos( math::abs(cos(lat1) * sin(alp1)) );
+ }
+};
+
+template <typename CT>
+class vertex_latitude_on_spheroid
+{
+
+public:
+/*
+ * formula based on paper
+ * [Wood96] Wood - Vertex Latitudes on Ellipsoid Geodesics, SIAM Rev., 38(4),
+ * 637–644, 1996
+ template <typename T1, typename T2, typename Spheroid>
+ static inline CT apply(T1 const& lat1,
+ T2 const& alp1,
+ Spheroid const& spheroid)
+ {
+ CT const f = formula::flattening<CT>(spheroid);
+
+ CT const e2 = f * (CT(2) - f);
+ CT const sin_alp1 = sin(alp1);
+ CT const sin2_lat1 = math::sqr(sin(lat1));
+ CT const cos2_lat1 = CT(1) - sin2_lat1;
+
+ CT const e2_sin2 = CT(1) - e2 * sin2_lat1;
+ CT const cos2_sin2 = cos2_lat1 * math::sqr(sin_alp1);
+ CT const vertex_lat = std::asin( math::sqrt((e2_sin2 - cos2_sin2)
+ / (e2_sin2 - e2 * cos2_sin2)));
+ return vertex_lat;
+ }
+*/
+
+ // simpler formula based on Clairaut relation for spheroids
+ template <typename T1, typename T2, typename Spheroid>
+ static inline CT apply(T1 const& lat1,
+ T2 const& alp1,
+ Spheroid const& spheroid)
+ {
+ CT const f = formula::flattening<CT>(spheroid);
+
+ CT const one_minus_f = (CT(1) - f);
+
+ //get the reduced latitude
+ CT const bet1 = atan( one_minus_f * tan(lat1) );
+
+ //apply Clairaut relation
+ CT const betv = vertex_latitude_on_sphere<CT>::apply(bet1, alp1);
+
+ //return the spheroid latitude
+ return atan( tan(betv) / one_minus_f );
+ }
+
+ /*
+ template <typename T>
+ inline static void sign_adjustment(CT lat1, CT lat2, CT vertex_lat, T& vrt_result)
+ {
+ // signbit returns a non-zero value (true) if the sign is negative;
+ // and zero (false) otherwise.
+ bool sign = std::signbit(std::abs(lat1) > std::abs(lat2) ? lat1 : lat2);
+
+ vrt_result.north = sign ? std::max(lat1, lat2) : vertex_lat;
+ vrt_result.south = sign ? vertex_lat * CT(-1) : std::min(lat1, lat2);
+ }
+
+ template <typename T>
+ inline static bool vertex_on_segment(CT alp1, CT alp2, CT lat1, CT lat2, T& vrt_result)
+ {
+ CT const half_pi = math::pi<CT>() / CT(2);
+
+ // if the segment does not contain the vertex of the geodesic
+ // then return the endpoint of max (min) latitude
+ if ((alp1 < half_pi && alp2 < half_pi)
+ || (alp1 > half_pi && alp2 > half_pi))
+ {
+ vrt_result.north = std::max(lat1, lat2);
+ vrt_result.south = std::min(lat1, lat2);
+ return false;
+ }
+ return true;
+ }
+ */
+};
+
+
+template <typename CT, typename CS_Tag>
+struct vertex_latitude
+{
+ BOOST_MPL_ASSERT_MSG
+ (
+ false, NOT_IMPLEMENTED_FOR_THIS_COORDINATE_SYSTEM, (types<CS_Tag>)
+ );
+
+};
+
+template <typename CT>
+struct vertex_latitude<CT, spherical_equatorial_tag>
+ : vertex_latitude_on_sphere<CT>
+{};
+
+template <typename CT>
+struct vertex_latitude<CT, geographic_tag>
+ : vertex_latitude_on_spheroid<CT>
+{};
+
+
+}}} // namespace boost::geometry::formula
+
+#endif // BOOST_GEOMETRY_FORMULAS_MAXIMUM_LATITUDE_HPP
diff --git a/boost/geometry/formulas/vincenty_direct.hpp b/boost/geometry/formulas/vincenty_direct.hpp
index f3647ff4e6..1697e5fb63 100644
--- a/boost/geometry/formulas/vincenty_direct.hpp
+++ b/boost/geometry/formulas/vincenty_direct.hpp
@@ -23,9 +23,8 @@
#include <boost/geometry/util/condition.hpp>
#include <boost/geometry/util/math.hpp>
-#include <boost/geometry/algorithms/detail/flattening.hpp>
-
#include <boost/geometry/formulas/differential_quantities.hpp>
+#include <boost/geometry/formulas/flattening.hpp>
#include <boost/geometry/formulas/result_direct.hpp>
@@ -85,7 +84,7 @@ public:
CT const radius_a = CT(get_radius<0>(spheroid));
CT const radius_b = CT(get_radius<2>(spheroid));
- CT const flattening = geometry::detail::flattening<CT>(spheroid);
+ CT const flattening = formula::flattening<CT>(spheroid);
CT const sin_azimuth12 = sin(azimuth12);
CT const cos_azimuth12 = cos(azimuth12);
diff --git a/boost/geometry/formulas/vincenty_inverse.hpp b/boost/geometry/formulas/vincenty_inverse.hpp
index bbda00036b..032e16e291 100644
--- a/boost/geometry/formulas/vincenty_inverse.hpp
+++ b/boost/geometry/formulas/vincenty_inverse.hpp
@@ -2,8 +2,8 @@
// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
-// This file was modified by Oracle on 2014, 2016.
-// Modifications copyright (c) 2014-2016 Oracle and/or its affiliates.
+// This file was modified by Oracle on 2014, 2016, 2017.
+// Modifications copyright (c) 2014-2017 Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
@@ -23,9 +23,8 @@
#include <boost/geometry/util/condition.hpp>
#include <boost/geometry/util/math.hpp>
-#include <boost/geometry/algorithms/detail/flattening.hpp>
-
#include <boost/geometry/formulas/differential_quantities.hpp>
+#include <boost/geometry/formulas/flattening.hpp>
#include <boost/geometry/formulas/result_inverse.hpp>
@@ -41,7 +40,7 @@ namespace boost { namespace geometry { namespace formula
\brief The solution of the inverse problem of geodesics on latlong coordinates, after Vincenty, 1975
\author See
- http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
- - http://www.icsm.gov.au/gda/gdav2.3.pdf
+ - http://www.icsm.gov.au/gda/gda-v_2.4.pdf
\author Adapted from various implementations to get it close to the original document
- http://www.movable-type.co.uk/scripts/LatLongVincenty.html
- http://exogen.case.edu/projects/geopy/source/geopy.distance.html
@@ -99,10 +98,10 @@ public:
CT const radius_a = CT(get_radius<0>(spheroid));
CT const radius_b = CT(get_radius<2>(spheroid));
- CT const flattening = geometry::detail::flattening<CT>(spheroid);
+ CT const f = formula::flattening<CT>(spheroid);
// U: reduced latitude, defined by tan U = (1-f) tan phi
- CT const one_min_f = c1 - flattening;
+ CT const one_min_f = c1 - f;
CT const tan_U1 = one_min_f * tan(lat1); // above (1)
CT const tan_U2 = one_min_f * tan(lat2); // above (1)
@@ -113,8 +112,9 @@ public:
CT const cos_U1 = c1 / temp_den_U1;
CT const cos_U2 = c1 / temp_den_U2;
// sin = tan / sqrt(1 + tan^2)
- CT const sin_U1 = tan_U1 / temp_den_U1;
- CT const sin_U2 = tan_U2 / temp_den_U2;
+ // sin = tan * cos
+ CT const sin_U1 = tan_U1 * cos_U1;
+ CT const sin_U2 = tan_U2 * cos_U2;
// calculate sin U and cos U directly
//CT const U1 = atan(tan_U1);
@@ -130,7 +130,8 @@ public:
CT sin_sigma;
CT sin_alpha;
CT cos2_alpha;
- CT cos2_sigma_m;
+ CT cos_2sigma_m;
+ CT cos2_2sigma_m;
CT sigma;
int counter = 0; // robustness
@@ -144,12 +145,13 @@ public:
CT cos_sigma = sin_U1 * sin_U2 + cos_U1 * cos_U2 * cos_lambda; // (15)
sin_alpha = cos_U1 * cos_U2 * sin_lambda / sin_sigma; // (17)
cos2_alpha = c1 - math::sqr(sin_alpha);
- cos2_sigma_m = math::equals(cos2_alpha, 0) ? 0 : cos_sigma - c2 * sin_U1 * sin_U2 / cos2_alpha; // (18)
+ cos_2sigma_m = math::equals(cos2_alpha, 0) ? 0 : cos_sigma - c2 * sin_U1 * sin_U2 / cos2_alpha; // (18)
+ cos2_2sigma_m = math::sqr(cos_2sigma_m);
- CT C = flattening/c16 * cos2_alpha * (c4 + flattening * (c4 - c3 * cos2_alpha)); // (10)
+ CT C = f/c16 * cos2_alpha * (c4 + f * (c4 - c3 * cos2_alpha)); // (10)
sigma = atan2(sin_sigma, cos_sigma); // (16)
- lambda = L + (c1 - C) * flattening * sin_alpha *
- (sigma + C * sin_sigma * ( cos2_sigma_m + C * cos_sigma * (-c1 + c2 * math::sqr(cos2_sigma_m)))); // (11)
+ lambda = L + (c1 - C) * f * sin_alpha *
+ (sigma + C * sin_sigma * (cos_2sigma_m + C * cos_sigma * (-c1 + c2 * cos2_2sigma_m))); // (11)
++counter; // robustness
@@ -182,8 +184,10 @@ public:
CT A = c1 + sqr_u/c16384 * (c4096 + sqr_u * (-c768 + sqr_u * (c320 - c175 * sqr_u))); // (3)
CT B = sqr_u/c1024 * (c256 + sqr_u * ( -c128 + sqr_u * (c74 - c47 * sqr_u))); // (4)
- CT delta_sigma = B * sin_sigma * ( cos2_sigma_m + (B/c4) * (cos(sigma)* (-c1 + c2 * cos2_sigma_m)
- - (B/c6) * cos2_sigma_m * (-c3 + c4 * math::sqr(sin_sigma)) * (-c3 + c4 * cos2_sigma_m))); // (6)
+ CT const cos_sigma = cos(sigma);
+ CT const sin2_sigma = math::sqr(sin_sigma);
+ CT delta_sigma = B * sin_sigma * (cos_2sigma_m + (B/c4) * (cos_sigma* (-c1 + c2 * cos2_2sigma_m)
+ - (B/c6) * cos_2sigma_m * (-c3 + c4 * sin2_sigma) * (-c3 + c4 * cos2_2sigma_m))); // (6)
result.distance = radius_b * A * (sigma - delta_sigma); // (19)
}
@@ -206,7 +210,7 @@ public:
typedef differential_quantities<CT, EnableReducedLength, EnableGeodesicScale, 2> quantities;
quantities::apply(lon1, lat1, lon2, lat2,
result.azimuth, result.reverse_azimuth,
- radius_b, flattening,
+ radius_b, f,
result.reduced_length, result.geodesic_scale);
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-01 10:11:16 +0000