Initial commit

winx64/include/boost/geometry/formulas/andoyer_inverse.hpp

@@ -0,0 +1,260 @@
+// Boost.Geometry
+
+// Copyright (c) 2015-2017 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_ANDOYER_INVERSE_HPP
+#define BOOST_GEOMETRY_FORMULAS_ANDOYER_INVERSE_HPP
+
+
+#include <boost/math/constants/constants.hpp>
+
+#include <boost/geometry/core/radius.hpp>
+#include <boost/geometry/core/srs.hpp>
+
+#include <boost/geometry/util/condition.hpp>
+#include <boost/geometry/util/math.hpp>
+
+#include <boost/geometry/formulas/differential_quantities.hpp>
+#include <boost/geometry/formulas/flattening.hpp>
+#include <boost/geometry/formulas/result_inverse.hpp>
+
+
+namespace boost { namespace geometry { namespace formula
+{
+
+/*!
+\brief The solution of the inverse problem of geodesics on latlong coordinates,
+       Forsyth-Andoyer-Lambert type approximation with first order terms.
+\author See
+    - Technical Report: PAUL D. THOMAS, MATHEMATICAL MODELS FOR NAVIGATION SYSTEMS, 1965
+      http://www.dtic.mil/docs/citations/AD0627893
+    - Technical Report: PAUL D. THOMAS, SPHEROIDAL GEODESICS, REFERENCE SYSTEMS, AND LOCAL GEOMETRY, 1970
+      http://www.dtic.mil/docs/citations/AD703541
+*/
+template <
+    typename CT,
+    bool EnableDistance,
+    bool EnableAzimuth,
+    bool EnableReverseAzimuth = false,
+    bool EnableReducedLength = false,
+    bool EnableGeodesicScale = false
+>
+class andoyer_inverse
+{
+    static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale;
+    static const bool CalcAzimuths = EnableAzimuth || EnableReverseAzimuth || CalcQuantities;
+    static const bool CalcFwdAzimuth = EnableAzimuth || CalcQuantities;
+    static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcQuantities;
+
+public:
+    typedef result_inverse<CT> result_type;
+
+    template <typename T1, typename T2, typename Spheroid>
+    static inline result_type apply(T1 const& lon1,
+                                    T1 const& lat1,
+                                    T2 const& lon2,
+                                    T2 const& lat2,
+                                    Spheroid const& spheroid)
+    {
+        result_type result;
+
+        // coordinates in radians
+
+        if ( math::equals(lon1, lon2) && math::equals(lat1, lat2) )
+        {
+            return result;
+        }
+
+        CT const c0 = CT(0);
+        CT const c1 = CT(1);
+        CT const pi = math::pi<CT>();
+        CT const f = formula::flattening<CT>(spheroid);
+
+        CT const dlon = lon2 - lon1;
+        CT const sin_dlon = sin(dlon);
+        CT const cos_dlon = cos(dlon);
+        CT const sin_lat1 = sin(lat1);
+        CT const cos_lat1 = cos(lat1);
+        CT const sin_lat2 = sin(lat2);
+        CT const cos_lat2 = cos(lat2);
+
+        // H,G,T = infinity if cos_d = 1 or cos_d = -1
+        // lat1 == +-90 && lat2 == +-90
+        // lat1 == lat2 && lon1 == lon2
+        CT cos_d = sin_lat1*sin_lat2 + cos_lat1*cos_lat2*cos_dlon;
+        // on some platforms cos_d may be outside valid range
+        if (cos_d < -c1)
+            cos_d = -c1;
+        else if (cos_d > c1)
+            cos_d = c1;
+
+        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);
+            CT const L = math::sqr(sin_lat1+sin_lat2);
+            CT const three_sin_d = CT(3) * sin_d;
+
+            CT const one_minus_cos_d = c1 - cos_d;
+            CT const one_plus_cos_d = c1 + cos_d;
+            // cos_d = 1 or cos_d = -1 means that the points are antipodal
+
+            CT const H = math::equals(one_minus_cos_d, c0) ?
+                            c0 :
+                            (d + three_sin_d) / one_minus_cos_d;
+            CT const G = math::equals(one_plus_cos_d, c0) ?
+                            c0 :
+                            (d - three_sin_d) / one_plus_cos_d;
+
+            CT const dd = -(f/CT(4))*(H*K+G*L);
+
+            CT const a = get_radius<0>(spheroid);
+
+            result.distance = a * (d + dd);
+        }
+
+        if ( BOOST_GEOMETRY_CONDITION(CalcAzimuths) )
+        {
+            // sin_d = 0 <=> antipodal points
+            if (math::equals(sin_d, c0))
+            {
+                // T = inf
+                // dA = inf
+                // azimuth = -inf
+                result.azimuth = lat1 <= lat2 ? c0 : pi;
+            }
+            else
+            {
+                CT const c2 = CT(2);
+
+                CT A = c0;
+                CT U = 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;
+                    A = atan2(sin_dlon, M);
+                    CT const sin_2A = sin(c2*A);
+                    U = (f/ c2)*math::sqr(cos_lat1)*sin_2A;
+                }
+
+                CT B = c0;
+                CT V = 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;
+                    B = atan2(sin_dlon, N);
+                    CT const sin_2B = sin(c2*B);
+                    V = (f/ c2)*math::sqr(cos_lat2)*sin_2B;
+                }
+
+                CT const T = d / sin_d;
+
+                // even with sin_d == 0 checked above if the second point
+                // is somewhere in the antipodal area T may still be great
+                // therefore dA and dB may be great and the resulting azimuths
+                // may be some more or less arbitrary angles
+
+                if (BOOST_GEOMETRY_CONDITION(CalcFwdAzimuth))
+                {
+                    CT const dA = V*T - U;
+                    result.azimuth = A - dA;
+                    normalize_azimuth(result.azimuth, A, dA);
+                }
+
+                if (BOOST_GEOMETRY_CONDITION(CalcRevAzimuth))
+                {
+                    CT const dB = -U*T + V;
+                    result.reverse_azimuth = pi - B - dB;
+                    if (result.reverse_azimuth > pi)
+                    {
+                        result.reverse_azimuth -= 2 * pi;
+                    }
+                    normalize_azimuth(result.reverse_azimuth, B, dB);
+                }
+            }
+        }
+
+        if (BOOST_GEOMETRY_CONDITION(CalcQuantities))
+        {
+            typedef differential_quantities<CT, EnableReducedLength, EnableGeodesicScale, 1> quantities;
+            quantities::apply(dlon, sin_lat1, cos_lat1, sin_lat2, cos_lat2,
+                              result.azimuth, result.reverse_azimuth,
+                              get_radius<2>(spheroid), f,
+                              result.reduced_length, result.geodesic_scale);
+        }
+
+        return result;
+    }
+
+private:
+    static inline void normalize_azimuth(CT & azimuth, CT const& A, CT const& dA)
+    {
+        CT const c0 = 0;
+        
+        if (A >= c0) // A indicates Eastern hemisphere
+        {
+            if (dA >= c0) // A altered towards 0
+            {
+                if (azimuth < c0)
+                {
+                    azimuth = c0;
+                }
+            }
+            else // dA < 0, A altered towards pi
+            {
+                CT const pi = math::pi<CT>();
+                if (azimuth > pi)
+                {
+                    azimuth = pi;
+                }
+            }
+        }
+        else // A indicates Western hemisphere
+        {
+            if (dA <= c0) // A altered towards 0
+            {
+                if (azimuth > c0)
+                {
+                    azimuth = c0;
+                }
+            }
+            else // dA > 0, A altered towards -pi
+            {
+                CT const minus_pi = -math::pi<CT>();
+                if (azimuth < minus_pi)
+                {
+                    azimuth = minus_pi;
+                }
+            }
+        }
+    }
+};
+
+}}} // namespace boost::geometry::formula
+
+
+#endif // BOOST_GEOMETRY_FORMULAS_ANDOYER_INVERSE_HPP

winx64/include/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

winx64/include/boost/geometry/formulas/differential_quantities.hpp

@@ -0,0 +1,303 @@
+// Boost.Geometry
+
+// Copyright (c) 2016-2017 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_INVERSE_DIFFERENTIAL_QUANTITIES_HPP
+#define BOOST_GEOMETRY_FORMULAS_INVERSE_DIFFERENTIAL_QUANTITIES_HPP
+
+
+#include <boost/geometry/util/condition.hpp>
+#include <boost/geometry/util/math.hpp>
+
+
+namespace boost { namespace geometry { namespace formula
+{
+
+/*!
+\brief The solution of a part of the inverse problem - differential quantities.
+\author See
+- Charles F.F Karney, Algorithms for geodesics, 2011
+https://arxiv.org/pdf/1109.4448.pdf
+*/
+template <
+    typename CT,
+    bool EnableReducedLength,
+    bool EnableGeodesicScale,
+    unsigned int Order = 2,
+    bool ApproxF = true
+>
+class differential_quantities
+{
+public:
+    static inline void apply(CT const& lon1, CT const& lat1,
+                             CT const& lon2, CT const& lat2,
+                             CT const& azimuth, CT const& reverse_azimuth,
+                             CT const& b, CT const& f,
+                             CT & reduced_length, CT & geodesic_scale)
+    {
+        CT const dlon = lon2 - lon1;
+        CT const sin_lat1 = sin(lat1);
+        CT const cos_lat1 = cos(lat1);
+        CT const sin_lat2 = sin(lat2);
+        CT const cos_lat2 = cos(lat2);
+
+        apply(dlon, sin_lat1, cos_lat1, sin_lat2, cos_lat2,
+              azimuth, reverse_azimuth,
+              b, f,
+              reduced_length, geodesic_scale);
+    }
+
+    static inline void apply(CT const& dlon,
+                             CT const& sin_lat1, CT const& cos_lat1,
+                             CT const& sin_lat2, CT const& cos_lat2,
+                             CT const& azimuth, CT const& reverse_azimuth,
+                             CT const& b, CT const& f,
+                             CT & reduced_length, CT & geodesic_scale)
+    {
+        CT const c0 = 0;
+        CT const c1 = 1;
+        CT const one_minus_f = c1 - f;
+
+        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))
+        {
+            CT const sig_12 = math::abs(dlon) / one_minus_f;
+            if (BOOST_GEOMETRY_CONDITION(EnableReducedLength))
+            {
+                CT m12 = sin(sig_12) * b;
+                reduced_length = m12;
+            }
+                
+            if (BOOST_GEOMETRY_CONDITION(EnableGeodesicScale))
+            {
+                CT M12 = cos(sig_12);
+                geodesic_scale = M12;
+            }
+        }
+        else
+        {
+            CT const c2 = 2;
+            CT const e2 = f * (c2 - f);
+            CT const ep2 = e2 / math::sqr(one_minus_f);
+
+            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;
+            CT cos_sig2 = cos_alp2 * cos_bet2;
+
+            normalize(sin_sig1, cos_sig1);
+            normalize(sin_sig2, cos_sig2);
+
+            CT const sin_alp0 = sin_alp1 * cos_bet1;
+            CT const cos_alp0_sqr = c1 - math::sqr(sin_alp0);
+
+            CT const J12 = BOOST_GEOMETRY_CONDITION(ApproxF) ?
+                           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 + ep2 * math::sqr(sin_bet1));
+            CT const dn2 = math::sqrt(c1 + ep2 * math::sqr(sin_bet2));
+
+            if (BOOST_GEOMETRY_CONDITION(EnableReducedLength))
+            {
+                CT const m12_b = dn2 * (cos_sig1 * sin_sig2)
+                               - dn1 * (sin_sig1 * cos_sig2)
+                               - cos_sig1 * cos_sig2 * J12;
+                CT const m12 = m12_b * b;
+
+                reduced_length = m12;
+            }
+
+            if (BOOST_GEOMETRY_CONDITION(EnableGeodesicScale))
+            {
+                CT const cos_sig12 = cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2;
+                CT const t = ep2 * (cos_bet1 - cos_bet2) * (cos_bet1 + cos_bet2) / (dn1 + dn2);
+                CT const M12 = cos_sig12 + (t * sin_sig2 - cos_sig2 * J12) * sin_sig1 / dn1;
+
+                geodesic_scale = M12;
+            }
+        }
+    }
+
+private:
+    /*! Approximation of J12, expanded into taylor series in f
+        Maxima script:
+        ep2: f * (2-f) / ((1-f)^2);
+        k2: ca02 * ep2;
+        assume(f < 1);
+        assume(sig > 0);
+        I1(sig):= integrate(sqrt(1 + k2 * sin(s)^2), s, 0, sig);
+        I2(sig):= integrate(1/sqrt(1 + k2 * sin(s)^2), s, 0, sig);
+        J(sig):= I1(sig) - I2(sig);
+        S: taylor(J(sig), f, 0, 3);
+        S1: factor( 2*integrate(sin(s)^2,s,0,sig)*ca02*f );
+        S2: factor( ((integrate(-6*ca02^2*sin(s)^4+6*ca02*sin(s)^2,s,0,sig)+integrate(-2*ca02^2*sin(s)^4+6*ca02*sin(s)^2,s,0,sig))*f^2)/4 );
+        S3: factor( ((integrate(30*ca02^3*sin(s)^6-54*ca02^2*sin(s)^4+24*ca02*sin(s)^2,s,0,sig)+integrate(6*ca02^3*sin(s)^6-18*ca02^2*sin(s)^4+24*ca02*sin(s)^2,s,0,sig))*f^3)/12 );
+    */
+    static inline CT J12_f(CT const& sin_sig1, CT const& cos_sig1,
+                           CT const& sin_sig2, CT const& cos_sig2,
+                           CT const& cos_alp0_sqr, CT const& f)
+    {
+        if (Order == 0)
+        {
+            return 0;
+        }
+
+        CT const c2 = 2;
+
+        CT const sig_12 = atan2(cos_sig1 * sin_sig2 - sin_sig1 * cos_sig2,
+                                cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2);
+        CT const sin_2sig1 = c2 * cos_sig1 * sin_sig1; // sin(2sig1)
+        CT const sin_2sig2 = c2 * cos_sig2 * sin_sig2; // sin(2sig2)
+        CT const sin_2sig_12 = sin_2sig2 - sin_2sig1;
+        CT const L1 = sig_12 - sin_2sig_12 / c2;
+
+        if (Order == 1)
+        {
+            return cos_alp0_sqr * f * L1;
+        }
+        
+        CT const sin_4sig1 = c2 * sin_2sig1 * (math::sqr(cos_sig1) - math::sqr(sin_sig1)); // sin(4sig1)
+        CT const sin_4sig2 = c2 * sin_2sig2 * (math::sqr(cos_sig2) - math::sqr(sin_sig2)); // sin(4sig2)
+        CT const sin_4sig_12 = sin_4sig2 - sin_4sig1;
+        
+        CT const c8 = 8;
+        CT const c12 = 12;
+        CT const c16 = 16;
+        CT const c24 = 24;
+
+        CT const L2 = -( cos_alp0_sqr * sin_4sig_12
+                         + (-c8 * cos_alp0_sqr + c12) * sin_2sig_12
+                         + (c12 * cos_alp0_sqr - c24) * sig_12)
+                       / c16;
+
+        if (Order == 2)
+        {
+            return cos_alp0_sqr * f * (L1 + f * L2);
+        }
+
+        CT const c4 = 4;
+        CT const c9 = 9;
+        CT const c48 = 48;
+        CT const c60 = 60;
+        CT const c64 = 64;
+        CT const c96 = 96;
+        CT const c128 = 128;
+        CT const c144 = 144;
+
+        CT const cos_alp0_quad = math::sqr(cos_alp0_sqr);
+        CT const sin3_2sig1 = math::sqr(sin_2sig1) * sin_2sig1;
+        CT const sin3_2sig2 = math::sqr(sin_2sig2) * sin_2sig2;
+        CT const sin3_2sig_12 = sin3_2sig2 - sin3_2sig1;
+
+        CT const A = (c9 * cos_alp0_quad - c12 * cos_alp0_sqr) * sin_4sig_12;
+        CT const B = c4 * cos_alp0_quad * sin3_2sig_12;
+        CT const C = (-c48 * cos_alp0_quad + c96 * cos_alp0_sqr - c64) * sin_2sig_12;
+        CT const D = (c60 * cos_alp0_quad - c144 * cos_alp0_sqr + c128) * sig_12;
+
+        CT const L3 = (A + B + C + D) / c64;
+
+        // Order 3 and higher
+        return cos_alp0_sqr * f * (L1 + f * (L2 + f * L3));
+    }
+
+    /*! Approximation of J12, expanded into taylor series in e'^2
+        Maxima script:
+        k2: ca02 * ep2;
+        assume(sig > 0);
+        I1(sig):= integrate(sqrt(1 + k2 * sin(s)^2), s, 0, sig);
+        I2(sig):= integrate(1/sqrt(1 + k2 * sin(s)^2), s, 0, sig);
+        J(sig):= I1(sig) - I2(sig);
+        S: taylor(J(sig), ep2, 0, 3);
+        S1: factor( integrate(sin(s)^2,s,0,sig)*ca02*ep2 );
+        S2: factor( (integrate(sin(s)^4,s,0,sig)*ca02^2*ep2^2)/2 );
+        S3: factor( (3*integrate(sin(s)^6,s,0,sig)*ca02^3*ep2^3)/8 );
+    */
+    static inline CT J12_ep_sqr(CT const& sin_sig1, CT const& cos_sig1,
+                                CT const& sin_sig2, CT const& cos_sig2,
+                                CT const& cos_alp0_sqr, CT const& ep_sqr)
+    {
+        if (Order == 0)
+        {
+            return 0;
+        }
+
+        CT const c2 = 2;
+        CT const c4 = 4;
+
+        CT const c2a0ep2 = cos_alp0_sqr * ep_sqr;
+
+        CT const sig_12 = atan2(cos_sig1 * sin_sig2 - sin_sig1 * cos_sig2,
+                                cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2); // sig2 - sig1
+        CT const sin_2sig1 = c2 * cos_sig1 * sin_sig1; // sin(2sig1)
+        CT const sin_2sig2 = c2 * cos_sig2 * sin_sig2; // sin(2sig2)
+        CT const sin_2sig_12 = sin_2sig2 - sin_2sig1;
+
+        CT const L1 = (c2 * sig_12 - sin_2sig_12) / c4;
+
+        if (Order == 1)
+        {
+            return c2a0ep2 * L1;
+        }
+
+        CT const c8 = 8;
+        CT const c64 = 64;
+        
+        CT const sin_4sig1 = c2 * sin_2sig1 * (math::sqr(cos_sig1) - math::sqr(sin_sig1)); // sin(4sig1)
+        CT const sin_4sig2 = c2 * sin_2sig2 * (math::sqr(cos_sig2) - math::sqr(sin_sig2)); // sin(4sig2)
+        CT const sin_4sig_12 = sin_4sig2 - sin_4sig1;
+        
+        CT const L2 = (sin_4sig_12 - c8 * sin_2sig_12 + 12 * sig_12) / c64;
+
+        if (Order == 2)
+        {
+            return c2a0ep2 * (L1 + c2a0ep2 * L2);
+        }
+
+        CT const sin3_2sig1 = math::sqr(sin_2sig1) * sin_2sig1;
+        CT const sin3_2sig2 = math::sqr(sin_2sig2) * sin_2sig2;
+        CT const sin3_2sig_12 = sin3_2sig2 - sin3_2sig1;
+
+        CT const c9 = 9;
+        CT const c48 = 48;
+        CT const c60 = 60;
+        CT const c512 = 512;
+
+        CT const L3 = (c9 * sin_4sig_12 + c4 * sin3_2sig_12 - c48 * sin_2sig_12 + c60 * sig_12) / c512;
+
+        // Order 3 and higher
+        return c2a0ep2 * (L1 + c2a0ep2 * (L2 + c2a0ep2 * L3));
+    }
+
+    static inline void normalize(CT & x, CT & y)
+    {
+        CT const len = math::sqrt(math::sqr(x) + math::sqr(y));
+        x /= len;
+        y /= len;
+    }
+};
+
+}}} // namespace boost::geometry::formula
+
+
+#endif // BOOST_GEOMETRY_FORMULAS_INVERSE_DIFFERENTIAL_QUANTITIES_HPP

winx64/include/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

winx64/include/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

winx64/include/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

winx64/include/boost/geometry/formulas/gnomonic_intersection.hpp

@@ -0,0 +1,148 @@
+// 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_GNOMONIC_INTERSECTION_HPP
+#define BOOST_GEOMETRY_FORMULAS_GNOMONIC_INTERSECTION_HPP
+
+#include <boost/geometry/core/access.hpp>
+#include <boost/geometry/core/cs.hpp>
+
+#include <boost/geometry/arithmetic/cross_product.hpp>
+#include <boost/geometry/formulas/gnomonic_spheroid.hpp>
+#include <boost/geometry/geometries/point.hpp>
+#include <boost/geometry/util/math.hpp>
+
+
+namespace boost { namespace geometry { namespace formula
+{
+
+/*!
+\brief The intersection of two geodesics using spheroidal gnomonic projection
+       as proposed by Karney.
+\author See
+    - Charles F.F Karney, Algorithms for geodesics, 2011
+      https://arxiv.org/pdf/1109.4448.pdf
+    - GeographicLib forum thread: Intersection between two geodesic lines
+      https://sourceforge.net/p/geographiclib/discussion/1026621/thread/21aaff9f/
+*/
+template
+<
+    typename CT,
+    template <typename, bool, bool, bool, bool, bool> class Inverse,
+    template <typename, bool, bool, bool, bool> class Direct
+>
+class gnomonic_intersection
+{
+public:
+    template <typename T1, typename T2, typename Spheroid>
+    static inline bool apply(T1 const& lona1, T1 const& lata1,
+                             T1 const& lona2, T1 const& lata2,
+                             T2 const& lonb1, T2 const& latb1,
+                             T2 const& lonb2, T2 const& latb2,
+                             CT & lon, CT & lat,
+                             Spheroid const& spheroid)
+    {
+        CT const lon_a1 = lona1;
+        CT const lat_a1 = lata1;
+        CT const lon_a2 = lona2;
+        CT const lat_a2 = lata2;
+        CT const lon_b1 = lonb1;
+        CT const lat_b1 = latb1;
+        CT const lon_b2 = lonb2;
+        CT const lat_b2 = latb2;
+
+        return apply(lon_a1, lat_a1, lon_a2, lat_a2, lon_b1, lat_b1, lon_b2, lat_b2, lon, lat, spheroid);
+    }
+    
+    template <typename Spheroid>
+    static inline bool apply(CT const& lona1, CT const& lata1,
+                             CT const& lona2, CT const& lata2,
+                             CT const& lonb1, CT const& latb1,
+                             CT const& lonb2, CT const& latb2,
+                             CT & lon, CT & lat,
+                             Spheroid const& spheroid)
+    {
+        typedef gnomonic_spheroid<CT, Inverse, Direct> gnom_t;
+
+        lon = (lona1 + lona2 + lonb1 + lonb2) / 4;
+        lat = (lata1 + lata2 + latb1 + latb2) / 4;
+        // TODO: consider normalizing lon
+
+        for (int i = 0; i < 10; ++i)
+        {
+            CT xa1, ya1, xa2, ya2;
+            CT xb1, yb1, xb2, yb2;
+            CT x, y;
+            double lat1, lon1;
+
+            bool ok = gnom_t::forward(lon, lat, lona1, lata1, xa1, ya1, spheroid)
+                   && gnom_t::forward(lon, lat, lona2, lata2, xa2, ya2, spheroid)
+                   && gnom_t::forward(lon, lat, lonb1, latb1, xb1, yb1, spheroid)
+                   && gnom_t::forward(lon, lat, lonb2, latb2, xb2, yb2, spheroid)
+                   && intersect(xa1, ya1, xa2, ya2, xb1, yb1, xb2, yb2, x, y)
+                   && gnom_t::inverse(lon, lat, x, y, lon1, lat1, spheroid);
+
+            if (! ok)
+            {
+                return false;
+            }
+            
+            if (math::equals(lat1, lat) && math::equals(lon1, lon))
+            {
+                break;
+            }
+
+            lat = lat1;
+            lon = lon1;
+        }
+
+        // NOTE: true is also returned if the number of iterations is too great
+        //       which means that the accuracy of the result is low
+        return true;
+    }
+
+private:
+    static inline bool intersect(CT const& xa1, CT const& ya1, CT const& xa2, CT const& ya2,
+                                 CT const& xb1, CT const& yb1, CT const& xb2, CT const& yb2,
+                                 CT & x, CT & y)
+    {
+        typedef model::point<CT, 3, cs::cartesian> v3d_t;
+
+        CT const c0 = 0;
+        CT const c1 = 1;
+
+        v3d_t const va1(xa1, ya1, c1);
+        v3d_t const va2(xa2, ya2, c1);
+        v3d_t const vb1(xb1, yb1, c1);
+        v3d_t const vb2(xb2, yb2, c1);
+
+        v3d_t const la = cross_product(va1, va2);
+        v3d_t const lb = cross_product(vb1, vb2);
+        v3d_t const p = cross_product(la, lb);
+
+        CT const z = get<2>(p);
+
+        if (math::equals(z, c0))
+        {
+            // degenerated or collinear segments
+            return false;
+        }
+
+        x = get<0>(p) / z;
+        y = get<1>(p) / z;
+        
+        return true;
+    }
+};
+
+}}} // namespace boost::geometry::formula
+
+
+#endif // BOOST_GEOMETRY_FORMULAS_GNOMONIC_INTERSECTION_HPP

winx64/include/boost/geometry/formulas/gnomonic_spheroid.hpp

@@ -0,0 +1,125 @@
+// 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_GNOMONIC_SPHEROID_HPP
+#define BOOST_GEOMETRY_FORMULAS_GNOMONIC_SPHEROID_HPP
+
+
+#include <boost/geometry/core/radius.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>
+
+
+namespace boost { namespace geometry { namespace formula
+{
+
+/*!
+\brief Gnomonic projection on spheroid (ellipsoid of revolution).
+\author See
+- Charles F.F Karney, Algorithms for geodesics, 2011
+https://arxiv.org/pdf/1109.4448.pdf
+*/
+template <
+    typename CT,
+    template <typename, bool, bool, bool, bool ,bool> class Inverse,
+    template <typename, bool, bool, bool, bool> class Direct
+>
+class gnomonic_spheroid
+{
+    typedef Inverse<CT, false, true, true, true, true> inverse_type;
+    typedef typename inverse_type::result_type inverse_result;
+
+    typedef Direct<CT, false, false, true, true> direct_quantities_type;
+    typedef Direct<CT, true, false, false, false> direct_coordinates_type;
+    typedef typename direct_coordinates_type::result_type direct_result;
+
+public:
+    template <typename Spheroid>
+    static inline bool forward(CT const& lon0, CT const& lat0,
+                               CT const& lon, CT const& lat,
+                               CT & x, CT & y,
+                               Spheroid const& spheroid)
+    {
+        inverse_result i_res = inverse_type::apply(lon0, lat0, lon, lat, spheroid);
+        CT const& m = i_res.reduced_length;
+        CT const& M = i_res.geodesic_scale;
+
+        if (math::smaller_or_equals(M, CT(0)))
+        {
+            return false;
+        }
+
+        CT rho = m / M;
+        x = sin(i_res.azimuth) * rho;
+        y = cos(i_res.azimuth) * rho;
+        
+        return true;
+    }
+
+    template <typename Spheroid>
+    static inline bool inverse(CT const& lon0, CT const& lat0,
+                               CT const& x, CT const& y,
+                               CT & lon, CT & lat,
+                               Spheroid const& spheroid)
+    {
+        CT const a = get_radius<0>(spheroid);
+        CT const ds_threshold = a * std::numeric_limits<CT>::epsilon(); // TODO: 0 for non-fundamental type
+        
+        CT const azimuth = atan2(x, y);
+        CT const rho = math::sqrt(math::sqr(x) + math::sqr(y)); // use hypot?
+        CT distance = a * atan(rho / a);
+
+        bool found = false;
+        for (int i = 0 ; i < 10 ; ++i)
+        {
+            direct_result d_res = direct_quantities_type::apply(lon0, lat0, distance, azimuth, spheroid);
+            CT const& m = d_res.reduced_length;
+            CT const& M = d_res.geodesic_scale;
+
+            if (math::smaller_or_equals(M, CT(0)))
+            {
+                // found = false;
+                return found;
+            }
+            
+            CT const drho = m / M - rho; // rho = m / M
+            CT const ds = drho * math::sqr(M); // drho/ds = 1/M^2
+            distance -= ds;
+
+            // ds_threshold may be 0
+            if (math::abs(ds) <= ds_threshold)
+            {
+                found = true;
+                break;
+            }
+        }
+
+        if (found)
+        {
+            direct_result d_res = direct_coordinates_type::apply(lon0, lat0, distance, azimuth, spheroid);
+            lon = d_res.lon2;
+            lat = d_res.lat2;
+        }
+
+        return found;
+    }
+};
+
+}}} // namespace boost::geometry::formula
+
+
+#endif // BOOST_GEOMETRY_FORMULAS_GNOMONIC_SPHEROID_HPP

winx64/include/boost/geometry/formulas/result_direct.hpp

@@ -0,0 +1,39 @@
+// 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_RESULT_DIRECT_HPP
+#define BOOST_GEOMETRY_FORMULAS_RESULT_DIRECT_HPP
+
+
+namespace boost { namespace geometry { namespace formula
+{
+
+template <typename T>
+struct result_direct
+{
+    result_direct()
+        : lon2(0)
+        , lat2(0)
+        , reverse_azimuth(0)
+        , reduced_length(0)
+        , geodesic_scale(1)
+    {}
+
+    T lon2;
+    T lat2;
+    T reverse_azimuth;
+    T reduced_length;
+    T geodesic_scale;
+};
+
+}}} // namespace boost::geometry::formula
+
+
+#endif // BOOST_GEOMETRY_FORMULAS_RESULT_DIRECT_HPP

winx64/include/boost/geometry/formulas/result_inverse.hpp

@@ -0,0 +1,39 @@
+// Boost.Geometry
+
+// Copyright (c) 2015-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_RESULT_INVERSE_HPP
+#define BOOST_GEOMETRY_FORMULAS_RESULT_INVERSE_HPP
+
+
+namespace boost { namespace geometry { namespace formula
+{
+
+template <typename T>
+struct result_inverse
+{
+    result_inverse()
+        : distance(0)
+        , azimuth(0)
+        , reverse_azimuth(0)
+        , reduced_length(0)
+        , geodesic_scale(1)
+    {}
+
+    T distance;
+    T azimuth;
+    T reverse_azimuth;
+    T reduced_length;
+    T geodesic_scale;
+};
+
+}}} // namespace boost::geometry::formula
+
+
+#endif // BOOST_GEOMETRY_FORMULAS_RESULT_INVERSE_HPP

winx64/include/boost/geometry/formulas/spherical.hpp

@@ -0,0 +1,224 @@
+// 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,
+// 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_SPHERICAL_HPP
+#define BOOST_GEOMETRY_FORMULAS_SPHERICAL_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/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 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;
+
+    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);
+
+    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;
+
+    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);
+
+    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);
+
+    math::normalize_spheroidal_coordinates
+        <
+            typename coordinate_system<PointSph>::type::units,
+            coord_t
+        >(lon, lat);
+
+    set<0>(res, lon);
+    set<1>(res, lat);
+
+    return res;
+}
+
+// -1 right
+// 1 left
+// 0 on
+template <typename Point3d1, typename Point3d2>
+static inline int sph_side_value(Point3d1 const& norm, Point3d2 const& pt)
+{
+    typedef typename select_coordinate_type<Point3d1, Point3d2>::type calc_t;
+    calc_t c0 = 0;
+    calc_t d = dot_product(norm, pt);
+    return math::equals(d, c0) ? 0
+        : d > c0 ? 1
+        : -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
+
+#endif // BOOST_GEOMETRY_FORMULAS_SPHERICAL_HPP

winx64/include/boost/geometry/formulas/thomas_direct.hpp

@@ -0,0 +1,248 @@
+// 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_THOMAS_DIRECT_HPP
+#define BOOST_GEOMETRY_FORMULAS_THOMAS_DIRECT_HPP
+
+
+#include <boost/math/constants/constants.hpp>
+
+#include <boost/geometry/core/radius.hpp>
+#include <boost/geometry/core/srs.hpp>
+
+#include <boost/geometry/util/condition.hpp>
+#include <boost/geometry/util/math.hpp>
+
+#include <boost/geometry/formulas/differential_quantities.hpp>
+#include <boost/geometry/formulas/flattening.hpp>
+#include <boost/geometry/formulas/result_direct.hpp>
+
+
+namespace boost { namespace geometry { namespace formula
+{
+
+
+/*!
+\brief The solution of the direct problem of geodesics on latlong coordinates,
+       Forsyth-Andoyer-Lambert type approximation with second order terms.
+\author See
+    - Technical Report: PAUL D. THOMAS, MATHEMATICAL MODELS FOR NAVIGATION SYSTEMS, 1965
+      http://www.dtic.mil/docs/citations/AD0627893
+    - Technical Report: PAUL D. THOMAS, SPHEROIDAL GEODESICS, REFERENCE SYSTEMS, AND LOCAL GEOMETRY, 1970
+      http://www.dtic.mil/docs/citations/AD0703541
+
+*/
+template <
+    typename CT,
+    bool EnableCoordinates = true,
+    bool EnableReverseAzimuth = false,
+    bool EnableReducedLength = false,
+    bool EnableGeodesicScale = false
+>
+class thomas_direct
+{
+    static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale;
+    static const bool CalcCoordinates = EnableCoordinates || CalcQuantities;
+    static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcCoordinates || CalcQuantities;
+
+public:
+    typedef result_direct<CT> result_type;
+
+    template <typename T, typename Dist, typename Azi, typename Spheroid>
+    static inline result_type apply(T const& lo1,
+                                    T const& la1,
+                                    Dist const& distance,
+                                    Azi const& azimuth12,
+                                    Spheroid const& spheroid)
+    {
+        result_type result;
+
+        CT const lon1 = lo1;
+        CT const lat1 = la1;
+
+        if ( math::equals(distance, Dist(0)) || distance < Dist(0) )
+        {
+            result.lon2 = lon1;
+            result.lat2 = lat1;
+            return result;
+        }
+
+        CT const c0 = 0;
+        CT const c1 = 1;
+        CT const c2 = 2;
+        CT const c4 = 4;
+
+        CT const a = CT(get_radius<0>(spheroid));
+        CT const b = CT(get_radius<2>(spheroid));
+        CT const f = formula::flattening<CT>(spheroid);
+        CT const one_minus_f = c1 - f;
+
+        CT const pi = math::pi<CT>();
+        CT const pi_half = pi / c2;
+
+        // keep azimuth small - experiments show low accuracy
+        // if the azimuth is closer to (+-)180 deg.
+        CT azi12_alt = azimuth12;
+        CT lat1_alt = lat1;
+        bool alter_result = vflip_if_south(lat1, azimuth12, lat1_alt, azi12_alt);
+        
+        CT const theta1 = math::equals(lat1_alt, pi_half) ? lat1_alt :
+                          math::equals(lat1_alt, -pi_half) ? lat1_alt :
+                          atan(one_minus_f * tan(lat1_alt));
+        CT const sin_theta1 = sin(theta1);
+        CT const cos_theta1 = cos(theta1);
+
+        CT const sin_a12 = sin(azi12_alt);
+        CT const cos_a12 = cos(azi12_alt);
+
+        CT const M = cos_theta1 * sin_a12; // cos_theta0
+        CT const theta0 = acos(M);
+        CT const sin_theta0 = sin(theta0);
+
+        CT const N = cos_theta1 * cos_a12;
+        CT const C1 = f * M; // lower-case c1 in the technical report
+        CT const C2 = f * (c1 - math::sqr(M)) / c4; // lower-case c2 in the technical report
+        CT const D = (c1 - C2) * (c1 - C2 - C1 * M);
+        CT const P = C2 * (c1 + C1 * M / c2) / D;
+
+        // special case for equator:
+        // sin_theta0 = 0 <=> lat1 = 0 ^ |azimuth12| = pi/2
+        // NOTE: in this case it doesn't matter what's the value of cos_sigma1 because
+        //       theta1=0, theta0=0, M=1|-1, C2=0 so X=0 and Y=0 so d_sigma=d
+        //       cos_a12=0 so N=0, therefore
+        //       lat2=0, azi21=pi/2|-pi/2
+        //       d_eta = atan2(sin_d_sigma, cos_d_sigma)
+        //       H = C1 * d_sigma
+        CT const cos_sigma1 = math::equals(sin_theta0, c0)
+                                ? c1
+                                : normalized1_1(sin_theta1 / sin_theta0);
+        CT const sigma1 = acos(cos_sigma1);
+        CT const d = distance / (a * D);
+        CT const u = 2 * (sigma1 - d);
+        CT const cos_d = cos(d);
+        CT const sin_d = sin(d);
+        CT const cos_u = cos(u);
+        CT const sin_u = sin(u);
+
+        CT const W = c1 - c2 * P * cos_u;
+        CT const V = cos_u * cos_d - sin_u * sin_d;
+        CT const X = math::sqr(C2) * sin_d * cos_d * (2 * math::sqr(V) - c1);
+        CT const Y = c2 * P * V * W * sin_d;
+        CT const d_sigma = d + X - Y;
+        CT const sin_d_sigma = sin(d_sigma);
+        CT const cos_d_sigma = cos(d_sigma);
+
+        if (BOOST_GEOMETRY_CONDITION(CalcRevAzimuth))
+        {
+            result.reverse_azimuth = atan2(M, N * cos_d_sigma - sin_theta1 * sin_d_sigma);
+
+            if (alter_result)
+            {
+                vflip_rev_azi(result.reverse_azimuth, azimuth12);
+            }
+        }
+
+        if (BOOST_GEOMETRY_CONDITION(CalcCoordinates))
+        {
+            CT const S_sigma = c2 * sigma1 - d_sigma;
+            CT const cos_S_sigma = cos(S_sigma);
+            CT const d_eta = atan2(sin_d_sigma * sin_a12, cos_theta1 * cos_d_sigma - sin_theta1 * sin_d_sigma * cos_a12);
+            CT const H = C1 * (c1 - C2) * d_sigma - C1 * C2 * sin_d_sigma * cos_S_sigma;
+            CT const d_lambda = d_eta - H;
+            
+            result.lon2 = lon1 + d_lambda;
+
+            if (! math::equals(M, c0))
+            {
+                CT const sin_a21 = sin(result.reverse_azimuth);
+                CT const tan_theta2 = (sin_theta1 * cos_d_sigma + N * sin_d_sigma) * sin_a21 / M;
+                result.lat2 = atan(tan_theta2 / one_minus_f);
+            }
+            else
+            {
+                CT const sigma2 = S_sigma - sigma1;
+                //theta2 = asin(cos(sigma2)) <=> sin_theta0 = 1
+                CT const tan_theta2 = cos(sigma2) / sin(sigma2);
+                result.lat2 = atan(tan_theta2 / one_minus_f);
+            }
+
+            if (alter_result)
+            {
+                result.lat2 = -result.lat2;
+            }
+        }        
+
+        if (BOOST_GEOMETRY_CONDITION(CalcQuantities))
+        {
+            typedef differential_quantities<CT, EnableReducedLength, EnableGeodesicScale, 2> quantities;
+            quantities::apply(lon1, lat1, result.lon2, result.lat2,
+                              azimuth12, result.reverse_azimuth,
+                              b, f,
+                              result.reduced_length, result.geodesic_scale);
+        }
+
+        return result;
+    }
+
+private:
+    static inline bool vflip_if_south(CT const& lat1, CT const& azi12, CT & lat1_alt, CT & azi12_alt)
+    {
+        CT const c2 = 2;
+        CT const pi = math::pi<CT>();
+        CT const pi_half = pi / c2;
+
+        if (azi12 > pi_half)
+        {
+            azi12_alt = pi - azi12;
+            lat1_alt = -lat1;
+            return true;
+        }
+        else if (azi12 < -pi_half)
+        {
+            azi12_alt = -pi - azi12;
+            lat1_alt = -lat1;
+            return true;
+        }
+
+        return false;
+    }
+
+    static inline void vflip_rev_azi(CT & rev_azi, CT const& azimuth12)
+    {
+        CT const c0 = 0;
+        CT const pi = math::pi<CT>();
+
+        if (rev_azi == c0)
+        {
+            rev_azi = azimuth12 >= 0 ? pi : -pi;
+        }
+        else if (rev_azi > c0)
+        {
+            rev_azi = pi - rev_azi;
+        }
+        else
+        {
+            rev_azi = -pi - rev_azi;
+        }
+    }
+
+    static inline CT normalized1_1(CT const& value)
+    {
+        CT const c1 = 1;
+        return value > c1 ? c1 :
+               value < -c1 ? -c1 :
+               value;
+    }
+};
+
+}}} // namespace boost::geometry::formula
+
+
+#endif // BOOST_GEOMETRY_FORMULAS_THOMAS_DIRECT_HPP

winx64/include/boost/geometry/formulas/thomas_inverse.hpp

@@ -0,0 +1,220 @@
+// Boost.Geometry
+
+// Copyright (c) 2015-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_THOMAS_INVERSE_HPP
+#define BOOST_GEOMETRY_FORMULAS_THOMAS_INVERSE_HPP
+
+
+#include <boost/math/constants/constants.hpp>
+
+#include <boost/geometry/core/radius.hpp>
+#include <boost/geometry/core/srs.hpp>
+
+#include <boost/geometry/util/condition.hpp>
+#include <boost/geometry/util/math.hpp>
+
+#include <boost/geometry/formulas/differential_quantities.hpp>
+#include <boost/geometry/formulas/flattening.hpp>
+#include <boost/geometry/formulas/result_inverse.hpp>
+
+
+namespace boost { namespace geometry { namespace formula
+{
+
+/*!
+\brief The solution of the inverse problem of geodesics on latlong coordinates,
+       Forsyth-Andoyer-Lambert type approximation with second order terms.
+\author See
+    - Technical Report: PAUL D. THOMAS, MATHEMATICAL MODELS FOR NAVIGATION SYSTEMS, 1965
+      http://www.dtic.mil/docs/citations/AD0627893
+    - Technical Report: PAUL D. THOMAS, SPHEROIDAL GEODESICS, REFERENCE SYSTEMS, AND LOCAL GEOMETRY, 1970
+      http://www.dtic.mil/docs/citations/AD0703541
+*/
+template <
+    typename CT,
+    bool EnableDistance,
+    bool EnableAzimuth,
+    bool EnableReverseAzimuth = false,
+    bool EnableReducedLength = false,
+    bool EnableGeodesicScale = false
+>
+class thomas_inverse
+{
+    static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale;
+    static const bool CalcAzimuths = EnableAzimuth || EnableReverseAzimuth || CalcQuantities;
+    static const bool CalcFwdAzimuth = EnableAzimuth || CalcQuantities;
+    static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcQuantities;
+
+public:
+    typedef result_inverse<CT> result_type;
+
+    template <typename T1, typename T2, typename Spheroid>
+    static inline result_type apply(T1 const& lon1,
+                                    T1 const& lat1,
+                                    T2 const& lon2,
+                                    T2 const& lat2,
+                                    Spheroid const& spheroid)
+    {
+        result_type result;
+
+        // coordinates in radians
+
+        if ( math::equals(lon1, lon2) && math::equals(lat1, lat2) )
+        {
+            return result;
+        }
+
+        CT const c0 = 0;
+        CT const c1 = 1;
+        CT const c2 = 2;
+        CT const c4 = 4;
+
+        CT const pi_half = math::pi<CT>() / c2;
+        CT const f = formula::flattening<CT>(spheroid);
+        CT const one_minus_f = c1 - f;
+
+//        CT const tan_theta1 = one_minus_f * tan(lat1);
+//        CT const tan_theta2 = one_minus_f * tan(lat2);
+//        CT const theta1 = atan(tan_theta1);
+//        CT const theta2 = atan(tan_theta2);
+
+        CT const theta1 = math::equals(lat1, pi_half) ? lat1 :
+                          math::equals(lat1, -pi_half) ? lat1 :
+                          atan(one_minus_f * tan(lat1));
+        CT const theta2 = math::equals(lat2, pi_half) ? lat2 :
+                          math::equals(lat2, -pi_half) ? lat2 :
+                          atan(one_minus_f * tan(lat2));
+
+        CT const theta_m = (theta1 + theta2) / c2;
+        CT const d_theta_m = (theta2 - theta1) / c2;
+        CT const d_lambda = lon2 - lon1;
+        CT const d_lambda_m = d_lambda / c2;
+
+        CT const sin_theta_m = sin(theta_m);
+        CT const cos_theta_m = cos(theta_m);
+        CT const sin_d_theta_m = sin(d_theta_m);
+        CT const cos_d_theta_m = cos(d_theta_m);
+        CT const sin2_theta_m = math::sqr(sin_theta_m);
+        CT const cos2_theta_m = math::sqr(cos_theta_m);
+        CT const sin2_d_theta_m = math::sqr(sin_d_theta_m);
+        CT const cos2_d_theta_m = math::sqr(cos_d_theta_m);
+        CT const sin_d_lambda_m = sin(d_lambda_m);
+        CT const sin2_d_lambda_m = math::sqr(sin_d_lambda_m);
+
+        CT const H = cos2_theta_m - sin2_d_theta_m;
+        CT const L = sin2_d_theta_m + H * sin2_d_lambda_m;
+        CT const cos_d = c1 - c2 * L;
+        CT const d = acos(cos_d);
+        CT const sin_d = sin(d);
+
+        CT const one_minus_L = c1 - L;
+
+        if ( math::equals(sin_d, c0)
+          || math::equals(L, c0)
+          || math::equals(one_minus_L, c0) )
+        {
+            return result;
+        }
+
+        CT const U = c2 * sin2_theta_m * cos2_d_theta_m / one_minus_L;
+        CT const V = c2 * sin2_d_theta_m * cos2_theta_m / L;
+        CT const X = U + V;
+        CT const Y = U - V;
+        CT const T = d / sin_d;
+        CT const D = c4 * math::sqr(T);
+        CT const E = c2 * cos_d;
+        CT const A = D * E;
+        CT const B = c2 * D;
+        CT const C = T - (A - E) / c2;
+
+        CT const f_sqr = math::sqr(f);
+        CT const f_sqr_per_64 = f_sqr / CT(64);
+    
+        if ( BOOST_GEOMETRY_CONDITION(EnableDistance) )
+        {
+            CT const n1 = X * (A + C*X);
+            CT const n2 = Y * (B + E*Y);
+            CT const n3 = D*X*Y;
+
+            CT const delta1d = f * (T*X-Y) / c4;
+            CT const delta2d = f_sqr_per_64 * (n1 - n2 + n3);
+
+            CT const a = get_radius<0>(spheroid);
+
+            //result.distance = a * sin_d * (T - delta1d);
+            result.distance = a * sin_d * (T - delta1d + delta2d);
+        }
+    
+        if ( BOOST_GEOMETRY_CONDITION(CalcAzimuths) )
+        {
+            // NOTE: if both cos_latX == 0 then below we'd have 0 * INF
+            // it's a situation when the endpoints are on the poles +-90 deg
+            // in this case the azimuth could either be 0 or +-pi
+            // but above always 0 is returned
+
+            CT const F = c2*Y-E*(c4-X);
+            CT const M = CT(32)*T-(CT(20)*T-A)*X-(B+c4)*Y;
+            CT const G = f*T/c2 + f_sqr_per_64 * M;
+            
+            // TODO:
+            // If d_lambda is close to 90 or -90 deg then tan(d_lambda) is big
+            // and F is small. The result is not accurate.
+            // In the edge case the result may be 2 orders of magnitude less
+            // accurate than Andoyer's.
+            CT const tan_d_lambda = tan(d_lambda);
+            CT const Q = -(F*G*tan_d_lambda) / c4;
+            CT const d_lambda_m_p = (d_lambda + Q) / c2;
+            CT const tan_d_lambda_m_p = tan(d_lambda_m_p);
+
+            CT const v = atan2(cos_d_theta_m, sin_theta_m * tan_d_lambda_m_p);
+            CT const u = atan2(-sin_d_theta_m, cos_theta_m * tan_d_lambda_m_p);
+
+            CT const pi = math::pi<CT>();
+
+            if (BOOST_GEOMETRY_CONDITION(EnableAzimuth))
+            {
+                CT alpha1 = v + u;
+                if (alpha1 > pi)
+                {
+                    alpha1 -= c2 * pi;
+                }
+
+                result.azimuth = alpha1;
+            }
+
+            if (BOOST_GEOMETRY_CONDITION(EnableReverseAzimuth))
+            {
+                CT alpha2 = pi - (v - u);
+                if (alpha2 > pi)
+                {
+                    alpha2 -= c2 * pi;
+                }
+
+                result.reverse_azimuth = alpha2;
+            }
+        }
+
+        if (BOOST_GEOMETRY_CONDITION(CalcQuantities))
+        {
+            typedef differential_quantities<CT, EnableReducedLength, EnableGeodesicScale, 2> quantities;
+            quantities::apply(lon1, lat1, lon2, lat2,
+                              result.azimuth, result.reverse_azimuth,
+                              get_radius<2>(spheroid), f,
+                              result.reduced_length, result.geodesic_scale);
+        }
+
+        return result;
+    }
+};
+
+}}} // namespace boost::geometry::formula
+
+
+#endif // BOOST_GEOMETRY_FORMULAS_THOMAS_INVERSE_HPP

winx64/include/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

winx64/include/boost/geometry/formulas/vincenty_direct.hpp

@@ -0,0 +1,184 @@
+// Boost.Geometry
+
+// 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.
+
+// 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_VINCENTY_DIRECT_HPP
+#define BOOST_GEOMETRY_FORMULAS_VINCENTY_DIRECT_HPP
+
+
+#include <boost/math/constants/constants.hpp>
+
+#include <boost/geometry/core/radius.hpp>
+#include <boost/geometry/core/srs.hpp>
+
+#include <boost/geometry/util/condition.hpp>
+#include <boost/geometry/util/math.hpp>
+
+#include <boost/geometry/formulas/differential_quantities.hpp>
+#include <boost/geometry/formulas/flattening.hpp>
+#include <boost/geometry/formulas/result_direct.hpp>
+
+
+#ifndef BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS
+#define BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS 1000
+#endif
+
+
+namespace boost { namespace geometry { namespace formula
+{
+
+/*!
+\brief The solution of the direct 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
+\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
+    - http://futureboy.homeip.net/fsp/colorize.fsp?fileName=navigation.frink
+
+*/
+template <
+    typename CT,
+    bool EnableCoordinates = true,
+    bool EnableReverseAzimuth = false,
+    bool EnableReducedLength = false,
+    bool EnableGeodesicScale = false
+>
+class vincenty_direct
+{
+    static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale;
+    static const bool CalcCoordinates = EnableCoordinates || CalcQuantities;
+    static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcQuantities;
+
+public:
+    typedef result_direct<CT> result_type;
+
+    template <typename T, typename Dist, typename Azi, typename Spheroid>
+    static inline result_type apply(T const& lo1,
+                                    T const& la1,
+                                    Dist const& distance,
+                                    Azi const& azimuth12,
+                                    Spheroid const& spheroid)
+    {
+        result_type result;
+
+        CT const lon1 = lo1;
+        CT const lat1 = la1;
+
+        if ( math::equals(distance, Dist(0)) || distance < Dist(0) )
+        {
+            result.lon2 = lon1;
+            result.lat2 = lat1;
+            return result;
+        }
+
+        CT const radius_a = CT(get_radius<0>(spheroid));
+        CT const radius_b = CT(get_radius<2>(spheroid));
+        CT const flattening = formula::flattening<CT>(spheroid);
+
+        CT const sin_azimuth12 = sin(azimuth12);
+        CT const cos_azimuth12 = cos(azimuth12);
+
+        // U: reduced latitude, defined by tan U = (1-f) tan phi
+        CT const one_min_f = CT(1) - flattening;
+        CT const tan_U1 = one_min_f * tan(lat1);
+        CT const sigma1 = atan2(tan_U1, cos_azimuth12); // (1)
+
+        // may be calculated from tan using 1 sqrt()
+        CT const U1 = atan(tan_U1);
+        CT const sin_U1 = sin(U1);
+        CT const cos_U1 = cos(U1);
+
+        CT const sin_alpha = cos_U1 * sin_azimuth12; // (2)
+        CT const sin_alpha_sqr = math::sqr(sin_alpha);
+        CT const cos_alpha_sqr = CT(1) - sin_alpha_sqr;
+
+        CT const b_sqr = radius_b * radius_b;
+        CT const u_sqr = cos_alpha_sqr * (radius_a * radius_a - b_sqr) / b_sqr;
+        CT const A = CT(1) + (u_sqr/CT(16384)) * (CT(4096) + u_sqr*(CT(-768) + u_sqr*(CT(320) - u_sqr*CT(175)))); // (3)
+        CT const B = (u_sqr/CT(1024))*(CT(256) + u_sqr*(CT(-128) + u_sqr*(CT(74) - u_sqr*CT(47)))); // (4)
+
+        CT s_div_bA = distance / (radius_b * A);
+        CT sigma = s_div_bA; // (7)
+
+        CT previous_sigma;
+        CT sin_sigma;
+        CT cos_sigma;
+        CT cos_2sigma_m;
+        CT cos_2sigma_m_sqr;
+
+        int counter = 0; // robustness
+
+        do
+        {
+            previous_sigma = sigma;
+
+            CT const two_sigma_m = CT(2) * sigma1 + sigma; // (5)
+
+            sin_sigma = sin(sigma);
+            cos_sigma = cos(sigma);
+            CT const sin_sigma_sqr = math::sqr(sin_sigma);
+            cos_2sigma_m = cos(two_sigma_m);
+            cos_2sigma_m_sqr = math::sqr(cos_2sigma_m);
+
+            CT const delta_sigma = B * sin_sigma * (cos_2sigma_m
+                                        + (B/CT(4)) * ( cos_sigma * (CT(-1) + CT(2)*cos_2sigma_m_sqr)
+                                            - (B/CT(6) * cos_2sigma_m * (CT(-3)+CT(4)*sin_sigma_sqr) * (CT(-3)+CT(4)*cos_2sigma_m_sqr)) )); // (6)
+
+            sigma = s_div_bA + delta_sigma; // (7)
+
+            ++counter; // robustness
+
+        } while ( geometry::math::abs(previous_sigma - sigma) > CT(1e-12)
+               //&& geometry::math::abs(sigma) < pi
+               && counter < BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS ); // robustness
+
+        if (BOOST_GEOMETRY_CONDITION(CalcCoordinates))
+        {
+            result.lat2
+                = atan2( sin_U1 * cos_sigma + cos_U1 * sin_sigma * cos_azimuth12,
+                         one_min_f * math::sqrt(sin_alpha_sqr + math::sqr(sin_U1 * sin_sigma - cos_U1 * cos_sigma * cos_azimuth12))); // (8)
+            
+            CT const lambda = atan2( sin_sigma * sin_azimuth12,
+                                     cos_U1 * cos_sigma - sin_U1 * sin_sigma * cos_azimuth12); // (9)
+            CT const C = (flattening/CT(16)) * cos_alpha_sqr * ( CT(4) + flattening * ( CT(4) - CT(3) * cos_alpha_sqr ) ); // (10)
+            CT const L = lambda - (CT(1) - C) * flattening * sin_alpha
+                            * ( sigma + C * sin_sigma * ( cos_2sigma_m + C * cos_sigma * ( CT(-1) + CT(2) * cos_2sigma_m_sqr ) ) ); // (11)
+
+            result.lon2 = lon1 + L;
+        }
+
+        if (BOOST_GEOMETRY_CONDITION(CalcRevAzimuth))
+        {
+            result.reverse_azimuth
+                = atan2(sin_alpha, -sin_U1 * sin_sigma + cos_U1 * cos_sigma * cos_azimuth12); // (12)
+        }
+
+        if (BOOST_GEOMETRY_CONDITION(CalcQuantities))
+        {
+            typedef differential_quantities<CT, EnableReducedLength, EnableGeodesicScale, 2> quantities;
+            quantities::apply(lon1, lat1, result.lon2, result.lat2,
+                              azimuth12, result.reverse_azimuth,
+                              radius_b, flattening,
+                              result.reduced_length, result.geodesic_scale);
+        }
+
+        return result;
+    }
+
+};
+
+}}} // namespace boost::geometry::formula
+
+
+#endif // BOOST_GEOMETRY_FORMULAS_VINCENTY_DIRECT_HPP

winx64/include/boost/geometry/formulas/vincenty_inverse.hpp

@@ -0,0 +1,224 @@
+// Boost.Geometry
+
+// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
+
+// 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
+
+// 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_VINCENTY_INVERSE_HPP
+#define BOOST_GEOMETRY_FORMULAS_VINCENTY_INVERSE_HPP
+
+
+#include <boost/math/constants/constants.hpp>
+
+#include <boost/geometry/core/radius.hpp>
+#include <boost/geometry/core/srs.hpp>
+
+#include <boost/geometry/util/condition.hpp>
+#include <boost/geometry/util/math.hpp>
+
+#include <boost/geometry/formulas/differential_quantities.hpp>
+#include <boost/geometry/formulas/flattening.hpp>
+#include <boost/geometry/formulas/result_inverse.hpp>
+
+
+#ifndef BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS
+#define BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS 1000
+#endif
+
+
+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/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
+    - http://futureboy.homeip.net/fsp/colorize.fsp?fileName=navigation.frink
+
+*/
+template <
+    typename CT,
+    bool EnableDistance,
+    bool EnableAzimuth,
+    bool EnableReverseAzimuth = false,
+    bool EnableReducedLength = false,
+    bool EnableGeodesicScale = false
+>
+struct vincenty_inverse
+{
+    static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale;
+    static const bool CalcAzimuths = EnableAzimuth || EnableReverseAzimuth || CalcQuantities;
+    static const bool CalcFwdAzimuth = EnableAzimuth || CalcQuantities;
+    static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcQuantities;
+
+public:
+    typedef result_inverse<CT> result_type;
+
+    template <typename T1, typename T2, typename Spheroid>
+    static inline result_type apply(T1 const& lon1,
+                                    T1 const& lat1,
+                                    T2 const& lon2,
+                                    T2 const& lat2,
+                                    Spheroid const& spheroid)
+    {
+        result_type result;
+
+        if (math::equals(lat1, lat2) && math::equals(lon1, lon2))
+        {
+            return result;
+        }
+
+        CT const c1 = 1;
+        CT const c2 = 2;
+        CT const c3 = 3;
+        CT const c4 = 4;
+        CT const c16 = 16;
+        CT const c_e_12 = CT(1e-12);
+
+        CT const pi = geometry::math::pi<CT>();
+        CT const two_pi = c2 * pi;
+
+        // lambda: difference in longitude on an auxiliary sphere
+        CT L = lon2 - lon1;
+        CT lambda = L;
+
+        if (L < -pi) L += two_pi;
+        if (L > pi) L -= two_pi;
+
+        CT const radius_a = CT(get_radius<0>(spheroid));
+        CT const radius_b = CT(get_radius<2>(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 - f;
+        CT const tan_U1 = one_min_f * tan(lat1); // above (1)
+        CT const tan_U2 = one_min_f * tan(lat2); // above (1)
+
+        // calculate sin U and cos U using trigonometric identities
+        CT const temp_den_U1 = math::sqrt(c1 + math::sqr(tan_U1));
+        CT const temp_den_U2 = math::sqrt(c1 + math::sqr(tan_U2));
+        // cos = 1 / sqrt(1 + tan^2)
+        CT const cos_U1 = c1 / temp_den_U1;
+        CT const cos_U2 = c1 / temp_den_U2;
+        // sin = tan / sqrt(1 + tan^2)
+        // 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);
+        //CT const U2 = atan(tan_U2);
+        //cos_U1 = cos(U1);
+        //cos_U2 = cos(U2);
+        //sin_U1 = tan_U1 * cos_U1; // sin(U1);
+        //sin_U2 = tan_U2 * cos_U2; // sin(U2);
+
+        CT previous_lambda;
+        CT sin_lambda;
+        CT cos_lambda;
+        CT sin_sigma;
+        CT sin_alpha;
+        CT cos2_alpha;
+        CT cos_2sigma_m;
+        CT cos2_2sigma_m;
+        CT sigma;
+
+        int counter = 0; // robustness
+
+        do
+        {
+            previous_lambda = lambda; // (13)
+            sin_lambda = sin(lambda);
+            cos_lambda = cos(lambda);
+            sin_sigma = math::sqrt(math::sqr(cos_U2 * sin_lambda) + math::sqr(cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_lambda)); // (14)
+            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);
+            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 = f/c16 * cos2_alpha * (c4 + f * (c4 - c3 * cos2_alpha)); // (10)
+            sigma = atan2(sin_sigma, cos_sigma); // (16)
+            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
+
+        } while ( geometry::math::abs(previous_lambda - lambda) > c_e_12
+               && geometry::math::abs(lambda) < pi
+               && counter < BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS ); // robustness
+    
+        if ( BOOST_GEOMETRY_CONDITION(EnableDistance) )
+        {
+            // Oops getting hard here
+            // (again, problem is that ttmath cannot divide by doubles, which is OK)
+            CT const c1 = 1;
+            CT const c2 = 2;
+            CT const c3 = 3;
+            CT const c4 = 4;
+            CT const c6 = 6;
+            CT const c47 = 47;
+            CT const c74 = 74;
+            CT const c128 = 128;
+            CT const c256 = 256;
+            CT const c175 = 175;
+            CT const c320 = 320;
+            CT const c768 = 768;
+            CT const c1024 = 1024;
+            CT const c4096 = 4096;
+            CT const c16384 = 16384;
+
+            //CT sqr_u = cos2_alpha * (math::sqr(radius_a) - math::sqr(radius_b)) / math::sqr(radius_b); // above (1)
+            CT sqr_u = cos2_alpha * ( math::sqr(radius_a / radius_b) - c1 ); // above (1)
+
+            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 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)
+        }
+    
+        if ( BOOST_GEOMETRY_CONDITION(CalcAzimuths) )
+        {
+            if (BOOST_GEOMETRY_CONDITION(CalcFwdAzimuth))
+            {
+                result.azimuth = atan2(cos_U2 * sin_lambda, cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_lambda); // (20)
+            }
+
+            if (BOOST_GEOMETRY_CONDITION(CalcRevAzimuth))
+            {
+                result.reverse_azimuth = atan2(cos_U1 * sin_lambda, -sin_U1 * cos_U2 + cos_U1 * sin_U2 * cos_lambda); // (21)
+            }
+        }
+
+        if (BOOST_GEOMETRY_CONDITION(CalcQuantities))
+        {
+            typedef differential_quantities<CT, EnableReducedLength, EnableGeodesicScale, 2> quantities;
+            quantities::apply(lon1, lat1, lon2, lat2,
+                              result.azimuth, result.reverse_azimuth,
+                              radius_b, f,
+                              result.reduced_length, result.geodesic_scale);
+        }
+
+        return result;
+    }
+};
+
+}}} // namespace boost::geometry::formula
+
+
+#endif // BOOST_GEOMETRY_FORMULAS_VINCENTY_INVERSE_HPP