bgirgis/cpp-thirdparty
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

update boost on linux

Browse Source
This commit is contained in:
Bassem Girgis
2019-08-10 16:06:25 -05:00
parent 76ad52be58
commit 861b918727
5363 changed files with 483306 additions and 116507 deletions
Download Patch File Download Diff File Expand all files Collapse all files

44
linx64/include/boost/geometry/formulas/andoyer_inverse.hpp
View File

@@ -1,5 +1,7 @@
// Boost.Geometry
// Copyright (c) 2018 Adam Wulkiewicz, Lodz, Poland.
// Copyright (c) 2015-2017 Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
@@ -15,7 +17,6 @@
#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>
@@ -116,20 +117,40 @@ public:
CT const dd = -(f/CT(4))*(H*K+G*L);
CT const a = get_radius<0>(spheroid);
CT const a = CT(get_radius<0>(spheroid));
result.distance = a * (d + dd);
}
if ( BOOST_GEOMETRY_CONDITION(CalcAzimuths) )
{
// sin_d = 0 <=> antipodal points
// sin_d = 0 <=> antipodal points (incl. poles)
if (math::equals(sin_d, c0))
{
// T = inf
// dA = inf
// azimuth = -inf
result.azimuth = lat1 <= lat2 ? c0 : pi;
// TODO: The following azimuths are inconsistent with distance
// i.e. according to azimuths below a segment with antipodal endpoints
// travels through the north pole, however the distance returned above
// is the length of a segment traveling along the equator.
// Furthermore, this special case handling is only done in andoyer
// formula.
// The most correct way of fixing it is to handle antipodal regions
// correctly and consistently across all formulas.
// Set azimuth to 0 unless the first endpoint is the north pole
if (! math::equals(sin_lat1, c1))
{
result.azimuth = c0;
result.reverse_azimuth = pi;
}
else
{
result.azimuth = pi;
result.reverse_azimuth = 0;
}
}
else
{
@@ -188,11 +209,10 @@ public:
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;
}
if (B >= 0)
result.reverse_azimuth = pi - B - dB;
else
result.reverse_azimuth = -pi - B - dB;
normalize_azimuth(result.reverse_azimuth, B, dB);
}
}
@@ -200,10 +220,12 @@ public:
if (BOOST_GEOMETRY_CONDITION(CalcQuantities))
{
CT const b = CT(get_radius<2>(spheroid));
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,
b, f,
result.reduced_length, result.geodesic_scale);
}
@@ -214,7 +236,7 @@ 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

46
linx64/include/boost/geometry/formulas/area_formulas.hpp
View File

@@ -1,8 +1,9 @@
// Boost.Geometry
// Copyright (c) 2015-2016 Oracle and/or its affiliates.
// Copyright (c) 2015-2018 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
@@ -11,7 +12,9 @@
#ifndef BOOST_GEOMETRY_FORMULAS_AREA_FORMULAS_HPP
#define BOOST_GEOMETRY_FORMULAS_AREA_FORMULAS_HPP
#include <boost/geometry/core/radian_access.hpp>
#include <boost/geometry/formulas/flattening.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/math/special_functions/hypot.hpp>
namespace boost { namespace geometry { namespace formula
@@ -44,7 +47,7 @@ public:
Evaluate the polynomial in x using Horner's method.
*/
template <typename NT, typename IteratorType>
static inline NT horner_evaluate(NT x,
static inline NT horner_evaluate(NT const& x,
IteratorType begin,
IteratorType end)
{
@@ -63,7 +66,7 @@ public:
https://en.wikipedia.org/wiki/Clenshaw_algorithm
*/
template <typename NT, typename IteratorType>
static inline NT clenshaw_sum(NT cosx,
static inline NT clenshaw_sum(NT const& cosx,
IteratorType begin,
IteratorType end)
{
@@ -165,7 +168,7 @@ public:
s/case\sCT(/case /g; s/):/:/g'
*/
static inline void evaluate_coeffs_n(CT n, CT coeffs_n[])
static inline void evaluate_coeffs_n(CT const& n, CT coeffs_n[])
{
switch (SeriesOrder) {
@@ -244,7 +247,7 @@ public:
/*
Expand in k2 and ep2.
*/
static inline void evaluate_coeffs_ep(CT ep, CT coeffs_n[])
static inline void evaluate_coeffs_ep(CT const& ep, CT coeffs_n[])
{
switch (SeriesOrder) {
case 0:
@@ -322,18 +325,22 @@ public:
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[]){
static inline void evaluate_coeffs_var2(CT const& var2,
CT const coeffs1[],
CT coeffs2[])
{
std::size_t begin(0), end(0);
for(std::size_t i = 0; i <= SeriesOrder; i++){
for(std::size_t i = 0; i <= SeriesOrder; i++)
{
end = begin + SeriesOrder + 1 - i;
coeffs2[i] = ((i==0) ? CT(1) : pow(var2,int(i)))
coeffs2[i] = ((i==0) ? CT(1) : math::pow(var2, int(i)))
* horner_evaluate(var2, coeffs1 + begin, coeffs1 + end);
begin = end;
}
}
/*
Compute the spherical excess of a geodesic (or shperical) segment
*/
@@ -402,13 +409,12 @@ public:
*/
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)
SpheroidConst const& spheroid_const)
{
return_type_ellipsoidal result;
@@ -540,12 +546,11 @@ public:
return result;
}
// Keep track whenever a segment crosses the prime meridian
// Check 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)
template <typename PointOfSegment>
static inline bool crosses_prime_meridian(PointOfSegment const& p1,
PointOfSegment const& p2)
{
CT const pi
= geometry::math::pi<CT>();
@@ -562,12 +567,7 @@ public:
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;
return max_lon > pi && min_lon < pi && max_lon - min_lon > pi;
}
};

96
linx64/include/boost/geometry/formulas/authalic_radius_sqr.hpp Normal file
View File

@@ -0,0 +1,96 @@
// Boost.Geometry
// Copyright (c) 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_AUTHALIC_RADIUS_SQR_HPP
#define BOOST_GEOMETRY_FORMULAS_AUTHALIC_RADIUS_SQR_HPP
#include <boost/geometry/core/radius.hpp>
#include <boost/geometry/core/tag.hpp>
#include <boost/geometry/core/tags.hpp>
#include <boost/geometry/formulas/eccentricity_sqr.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/geometry/algorithms/not_implemented.hpp>
#include <boost/math/special_functions/atanh.hpp>
namespace boost { namespace geometry
{
#ifndef DOXYGEN_NO_DISPATCH
namespace formula_dispatch
{
template <typename ResultType, typename Geometry, typename Tag = typename tag<Geometry>::type>
struct authalic_radius_sqr
: not_implemented<Tag>
{};
template <typename ResultType, typename Geometry>
struct authalic_radius_sqr<ResultType, Geometry, srs_sphere_tag>
{
static inline ResultType apply(Geometry const& geometry)
{
return math::sqr<ResultType>(get_radius<0>(geometry));
}
};
template <typename ResultType, typename Geometry>
struct authalic_radius_sqr<ResultType, Geometry, srs_spheroid_tag>
{
static inline ResultType apply(Geometry const& geometry)
{
ResultType const a2 = math::sqr<ResultType>(get_radius<0>(geometry));
ResultType const e2 = formula::eccentricity_sqr<ResultType>(geometry);
return apply(a2, e2);
}
static inline ResultType apply(ResultType const& a2, ResultType const& e2)
{
ResultType const c0 = 0;
if (math::equals(e2, c0))
{
return a2;
}
ResultType const e = math::sqrt(e2);
ResultType const c2 = 2;
//ResultType const b2 = math::sqr(get_radius<2>(geometry));
//return a2 / c2 + b2 * boost::math::atanh(e) / (c2 * e);
ResultType const c1 = 1;
return (a2 / c2) * ( c1 + (c1 - e2) * boost::math::atanh(e) / e );
}
};
} // namespace formula_dispatch
#endif // DOXYGEN_NO_DISPATCH
#ifndef DOXYGEN_NO_DETAIL
namespace formula
{
template <typename ResultType, typename Geometry>
inline ResultType authalic_radius_sqr(Geometry const& geometry)
{
return formula_dispatch::authalic_radius_sqr<ResultType, Geometry>::apply(geometry);
}
} // namespace formula
#endif // DOXYGEN_NO_DETAIL
}} // namespace boost::geometry
#endif // BOOST_GEOMETRY_FORMULAS_AUTHALIC_RADIUS_SQR_HPP

7
linx64/include/boost/geometry/formulas/differential_quantities.hpp
View File

@@ -1,6 +1,6 @@
// Boost.Geometry
// Copyright (c) 2016-2017 Oracle and/or its affiliates.
// Copyright (c) 2016-2019 Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
@@ -70,10 +70,11 @@ public:
// equator
if (math::equals(sin_bet1, c0) && math::equals(sin_bet2, c0))
{
CT const sig_12 = math::abs(dlon) / one_minus_f;
CT const sig_12 = dlon / one_minus_f;
if (BOOST_GEOMETRY_CONDITION(EnableReducedLength))
{
CT m12 = sin(sig_12) * b;
int azi_sign = math::sign(azimuth) >= 0 ? 1 : -1; // for antipodal
CT m12 = azi_sign * sin(sig_12) * b;
reduced_length = m12;
}

6
linx64/include/boost/geometry/formulas/eccentricity_sqr.hpp
View File

@@ -1,6 +1,6 @@
// Boost.Geometry
// Copyright (c) 2016 Oracle and/or its affiliates.
// Copyright (c) 2016, 2018 Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
@@ -11,11 +11,13 @@
#ifndef BOOST_GEOMETRY_FORMULAS_ECCENCRICITY_SQR_HPP
#define BOOST_GEOMETRY_FORMULAS_ECCENCRICITY_SQR_HPP
#include <boost/geometry/algorithms/not_implemented.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>
#include <boost/geometry/util/math.hpp>
namespace boost { namespace geometry
{

7
linx64/include/boost/geometry/formulas/geographic.hpp
View File

@@ -1,6 +1,6 @@
// Boost.Geometry
// Copyright (c) 2016, Oracle and/or its affiliates.
// Copyright (c) 2016-2017, Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
// Use, modification and distribution is subject to the Boost Software License,
@@ -22,6 +22,7 @@
#include <boost/geometry/formulas/eccentricity_sqr.hpp>
#include <boost/geometry/formulas/flattening.hpp>
#include <boost/geometry/formulas/unit_spheroid.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
@@ -186,7 +187,7 @@ inline Point3d projected_to_surface(Point3d const& direction, Spheroid const& sp
//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 b_sqr = math::sqr(formula::unit_spheroid_b<coord_t>(spheroid));
coord_t const param_a = (dx*dx + dy*dy) / a_sqr + dz*dz / b_sqr;
coord_t const delta = c4 * param_a;
@@ -226,7 +227,7 @@ inline bool projected_to_surface(Point3d const& origin, Point3d const& direction
//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 b_sqr = math::sqr(formula::unit_spheroid_b<coord_t>(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);

107
linx64/include/boost/geometry/formulas/interpolate_point_spherical.hpp Normal file
View File

@@ -0,0 +1,107 @@
// Boost.Geometry
// Copyright (c) 2019 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_INTERPOLATE_POINT_SPHERICAL_HPP
#define BOOST_GEOMETRY_FORMULAS_INTERPOLATE_POINT_SPHERICAL_HPP
namespace boost { namespace geometry { namespace formula
{
template <typename CalculationType>
class interpolate_point_spherical
{
typedef model::point<CalculationType, 3, cs::cartesian> point3d_t;
public :
template <typename Point>
void compute_angle(Point const& p0,
Point const& p1,
CalculationType& angle01)
{
m_xyz0 = formula::sph_to_cart3d<point3d_t>(p0);
m_xyz1 = formula::sph_to_cart3d<point3d_t>(p1);
CalculationType const dot01 = geometry::dot_product(m_xyz0, m_xyz1);
angle01 = acos(dot01);
}
template <typename Point>
void compute_axis(Point const& p0,
CalculationType const& angle01)
{
CalculationType const c0 = 0, c1 = 1;
CalculationType const pi = math::pi<CalculationType>();
if (! math::equals(angle01, pi))
{
m_axis = geometry::cross_product(m_xyz0, m_xyz1);
geometry::detail::vec_normalize(m_axis);
}
else // antipodal
{
CalculationType const half_pi = math::half_pi<CalculationType>();
CalculationType const lat = geometry::get_as_radian<1>(p0);
if (math::equals(lat, half_pi))
{
// pointing east, segment lies on prime meridian, going south
m_axis = point3d_t(c0, c1, c0);
}
else if (math::equals(lat, -half_pi))
{
// pointing west, segment lies on prime meridian, going north
m_axis = point3d_t(c0, -c1, c0);
}
else
{
// lon rotated west by pi/2 at equator
CalculationType const lon = geometry::get_as_radian<0>(p0);
m_axis = point3d_t(sin(lon), -cos(lon), c0);
}
}
}
template <typename Point>
void compute_point(CalculationType const& a, Point& p)
{
CalculationType const c1 = 1;
// Axis-Angle rotation
// see: https://en.wikipedia.org/wiki/Axis-angle_representation
CalculationType const cos_a = cos(a);
CalculationType const sin_a = sin(a);
// cos_a * v
point3d_t s1 = m_xyz0;
geometry::multiply_value(s1, cos_a);
// sin_a * (n x v)
point3d_t s2 = geometry::cross_product(m_axis, m_xyz0);
geometry::multiply_value(s2, sin_a);
// (1 - cos_a)(n.v) * n
point3d_t s3 = m_axis;
geometry::multiply_value(s3, (c1 - cos_a) *
geometry::dot_product(m_axis, m_xyz0));
// v_rot = cos_a * v + sin_a * (n x v) + (1 - cos_a)(n.v) * e
point3d_t v_rot = s1;
geometry::add_point(v_rot, s2);
geometry::add_point(v_rot, s3);
p = formula::cart3d_to_sph<Point>(v_rot);
}
private :
point3d_t m_xyz0;
point3d_t m_xyz1;
point3d_t m_axis;
};
}}} // namespace boost::geometry::formula
#endif // BOOST_GEOMETRY_FORMULAS_INTERPOLATE_POINT_SPHERICAL_HPP

263
linx64/include/boost/geometry/formulas/karney_direct.hpp Normal file
View File

@@ -0,0 +1,263 @@
// Boost.Geometry
// Copyright (c) 2018 Adeel Ahmad, Islamabad, Pakistan.
// Contributed and/or modified by Adeel Ahmad,
// as part of Google Summer of Code 2018 program.
// This file was modified by Oracle on 2018.
// Modifications copyright (c) 2018 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)
// This file is converted from GeographicLib, https://geographiclib.sourceforge.io
// GeographicLib is originally written by Charles Karney.
// Author: Charles Karney (2008-2017)
// Last updated version of GeographicLib: 1.49
// Original copyright notice:
// Copyright (c) Charles Karney (2008-2017) <charles@karney.com> and licensed
// under the MIT/X11 License. For more information, see
// https://geographiclib.sourceforge.io
#ifndef BOOST_GEOMETRY_FORMULAS_KARNEY_DIRECT_HPP
#define BOOST_GEOMETRY_FORMULAS_KARNEY_DIRECT_HPP
#include <boost/array.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/hypot.hpp>
#include <boost/geometry/formulas/flattening.hpp>
#include <boost/geometry/formulas/result_direct.hpp>
#include <boost/geometry/util/condition.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
#include <boost/geometry/util/series_expansion.hpp>
namespace boost { namespace geometry { namespace formula
{
namespace se = series_expansion;
/*!
\brief The solution of the direct problem of geodesics on latlong coordinates,
after Karney (2011).
\author See
- Charles F.F Karney, Algorithms for geodesics, 2011
https://arxiv.org/pdf/1109.4448.pdf
*/
template <
typename CT,
bool EnableCoordinates = true,
bool EnableReverseAzimuth = false,
bool EnableReducedLength = false,
bool EnableGeodesicScale = false,
size_t SeriesOrder = 8
>
class karney_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 lon1 = lo1;
CT const lat1 = la1;
Azi azi12 = azimuth12;
math::normalize_azimuth<degree, Azi>(azi12);
CT const c0 = 0;
CT const c1 = 1;
CT const c2 = 2;
CT const b = CT(get_radius<2>(spheroid));
CT const f = formula::flattening<CT>(spheroid);
CT const one_minus_f = c1 - f;
CT const two_minus_f = c2 - f;
CT const n = f / two_minus_f;
CT const e2 = f * two_minus_f;
CT const ep2 = e2 / math::sqr(one_minus_f);
CT sin_alpha1, cos_alpha1;
math::sin_cos_degrees<CT>(math::round_angle<CT>(azi12), sin_alpha1, cos_alpha1);
// Find the reduced latitude.
CT sin_beta1, cos_beta1;
math::sin_cos_degrees<CT>(math::round_angle<CT>(lat1), sin_beta1, cos_beta1);
sin_beta1 *= one_minus_f;
math::normalize_unit_vector<CT>(sin_beta1, cos_beta1);
cos_beta1 = (std::max)(c0, cos_beta1);
// Obtain alpha 0 by solving the spherical triangle.
CT const sin_alpha0 = sin_alpha1 * cos_beta1;
CT const cos_alpha0 = boost::math::hypot(cos_alpha1, sin_alpha1 * sin_beta1);
CT const k2 = math::sqr(cos_alpha0) * ep2;
CT const epsilon = k2 / (c2 * (c1 + math::sqrt(c1 + k2)) + k2);
// Find the coefficients for A1 by computing the
// series expansion using Horner scehme.
CT const expansion_A1 = se::evaluate_A1<SeriesOrder>(epsilon);
// Index zero element of coeffs_C1 is unused.
se::coeffs_C1<SeriesOrder, CT> const coeffs_C1(epsilon);
// Tau is an integration variable.
CT const tau12 = distance / (b * (c1 + expansion_A1));
CT const sin_tau12 = sin(tau12);
CT const cos_tau12 = cos(tau12);
CT sin_sigma1 = sin_beta1;
CT sin_omega1 = sin_alpha0 * sin_beta1;
CT cos_sigma1, cos_omega1;
cos_sigma1 = cos_omega1 = sin_beta1 != c0 || cos_alpha1 != c0 ? cos_beta1 * cos_alpha1 : c1;
math::normalize_unit_vector<CT>(sin_sigma1, cos_sigma1);
CT const B11 = se::sin_cos_series(sin_sigma1, cos_sigma1, coeffs_C1);
CT const sin_B11 = sin(B11);
CT const cos_B11 = cos(B11);
CT const sin_tau1 = sin_sigma1 * cos_B11 + cos_sigma1 * sin_B11;
CT const cos_tau1 = cos_sigma1 * cos_B11 - sin_sigma1 * sin_B11;
// Index zero element of coeffs_C1p is unused.
se::coeffs_C1p<SeriesOrder, CT> const coeffs_C1p(epsilon);
CT const B12 = - se::sin_cos_series
(sin_tau1 * cos_tau12 + cos_tau1 * sin_tau12,
cos_tau1 * cos_tau12 - sin_tau1 * sin_tau12,
coeffs_C1p);
CT const sigma12 = tau12 - (B12 - B11);
CT const sin_sigma12 = sin(sigma12);
CT const cos_sigma12 = cos(sigma12);
CT const sin_sigma2 = sin_sigma1 * cos_sigma12 + cos_sigma1 * sin_sigma12;
CT const cos_sigma2 = cos_sigma1 * cos_sigma12 - sin_sigma1 * sin_sigma12;
if (BOOST_GEOMETRY_CONDITION(CalcRevAzimuth))
{
CT const sin_alpha2 = sin_alpha0;
CT const cos_alpha2 = cos_alpha0 * cos_sigma2;
result.reverse_azimuth = atan2(sin_alpha2, cos_alpha2);
// Convert the angle to radians.
result.reverse_azimuth /= math::d2r<CT>();
}
if (BOOST_GEOMETRY_CONDITION(CalcCoordinates))
{
// Find the latitude at the second point.
CT const sin_beta2 = cos_alpha0 * sin_sigma2;
CT const cos_beta2 = boost::math::hypot(sin_alpha0, cos_alpha0 * cos_sigma2);
result.lat2 = atan2(sin_beta2, one_minus_f * cos_beta2);
// Convert the coordinate to radians.
result.lat2 /= math::d2r<CT>();
// Find the longitude at the second point.
CT const sin_omega2 = sin_alpha0 * sin_sigma2;
CT const cos_omega2 = cos_sigma2;
CT const omega12 = atan2(sin_omega2 * cos_omega1 - cos_omega2 * sin_omega1,
cos_omega2 * cos_omega1 + sin_omega2 * sin_omega1);
se::coeffs_A3<SeriesOrder, CT> const coeffs_A3(n);
CT const A3 = math::horner_evaluate(epsilon, coeffs_A3.begin(), coeffs_A3.end());
CT const A3c = -f * sin_alpha0 * A3;
se::coeffs_C3<SeriesOrder, CT> const coeffs_C3(n, epsilon);
CT const B31 = se::sin_cos_series(sin_sigma1, cos_sigma1, coeffs_C3);
CT const lam12 = omega12 + A3c *
(sigma12 + (se::sin_cos_series
(sin_sigma2,
cos_sigma2,
coeffs_C3) - B31));
// Convert to radians to get the
// longitudinal difference.
CT lon12 = lam12 / math::d2r<CT>();
// Add the longitude at first point to the longitudinal
// difference and normalize the result.
math::normalize_longitude<degree, CT>(lon1);
math::normalize_longitude<degree, CT>(lon12);
result.lon2 = lon1 + lon12;
}
if (BOOST_GEOMETRY_CONDITION(CalcQuantities))
{
// Evaluate the coefficients for C2.
// Index zero element of coeffs_C2 is unused.
se::coeffs_C2<SeriesOrder, CT> const coeffs_C2(epsilon);
CT const B21 = se::sin_cos_series(sin_sigma1, cos_sigma1, coeffs_C2);
CT const B22 = se::sin_cos_series(sin_sigma2, cos_sigma2, coeffs_C2);
// Find the coefficients for A2 by computing the
// series expansion using Horner scehme.
CT const expansion_A2 = se::evaluate_A2<SeriesOrder>(epsilon);
CT const AB1 = (c1 + expansion_A1) * (B12 - B11);
CT const AB2 = (c1 + expansion_A2) * (B22 - B21);
CT const J12 = (expansion_A1 - expansion_A2) * sigma12 + (AB1 - AB2);
CT const dn1 = math::sqrt(c1 + ep2 * math::sqr(sin_beta1));
CT const dn2 = math::sqrt(c1 + k2 * math::sqr(sin_sigma2));
// Find the reduced length.
result.reduced_length = b * ((dn2 * (cos_sigma1 * sin_sigma2) -
dn1 * (sin_sigma1 * cos_sigma2)) -
cos_sigma1 * cos_sigma2 * J12);
// Find the geodesic scale.
CT const t = k2 * (sin_sigma2 - sin_sigma1) *
(sin_sigma2 + sin_sigma1) / (dn1 + dn2);
result.geodesic_scale = cos_sigma12 +
(t * sin_sigma2 - cos_sigma2 * J12) *
sin_sigma1 / dn1;
}
return result;
}
};
}}} // namespace boost::geometry::formula
#endif // BOOST_GEOMETRY_FORMULAS_KARNEY_DIRECT_HPP

71
linx64/include/boost/geometry/formulas/mean_radius.hpp Normal file
View File

@@ -0,0 +1,71 @@
// Boost.Geometry
// Copyright (c) 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_MEAN_RADIUS_HPP
#define BOOST_GEOMETRY_FORMULAS_MEAN_RADIUS_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 mean_radius
: not_implemented<Tag>
{};
template <typename ResultType, typename Geometry>
struct mean_radius<ResultType, Geometry, srs_sphere_tag>
{
static inline ResultType apply(Geometry const& geometry)
{
return ResultType(get_radius<0>(geometry));
}
};
template <typename ResultType, typename Geometry>
struct mean_radius<ResultType, Geometry, srs_spheroid_tag>
{
static inline ResultType apply(Geometry const& geometry)
{
// (2*a + b) / 3
return (ResultType(2) * ResultType(get_radius<0>(geometry))
+ ResultType(get_radius<2>(geometry)))
/ ResultType(3);
}
};
} // namespace formula_dispatch
#endif // DOXYGEN_NO_DISPATCH
#ifndef DOXYGEN_NO_DETAIL
namespace formula
{
template <typename ResultType, typename Geometry>
inline ResultType mean_radius(Geometry const& geometry)
{
return formula_dispatch::mean_radius<ResultType, Geometry>::apply(geometry);
}
} // namespace formula
#endif // DOXYGEN_NO_DETAIL
}} // namespace boost::geometry
#endif // BOOST_GEOMETRY_FORMULAS_MEAN_RADIUS_HPP

170
linx64/include/boost/geometry/formulas/meridian_direct.hpp Normal file
View File

@@ -0,0 +1,170 @@
// Boost.Geometry
// Copyright (c) 2018 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_MERIDIAN_DIRECT_HPP
#define BOOST_GEOMETRY_FORMULAS_MERIDIAN_DIRECT_HPP
#include <boost/math/constants/constants.hpp>
#include <boost/geometry/core/radius.hpp>
#include <boost/geometry/formulas/differential_quantities.hpp>
#include <boost/geometry/formulas/flattening.hpp>
#include <boost/geometry/formulas/meridian_inverse.hpp>
#include <boost/geometry/formulas/quarter_meridian.hpp>
#include <boost/geometry/formulas/result_direct.hpp>
#include <boost/geometry/util/condition.hpp>
#include <boost/geometry/util/math.hpp>
namespace boost { namespace geometry { namespace formula
{
/*!
\brief Compute the direct geodesic problem on a meridian
*/
template <
typename CT,
bool EnableCoordinates = true,
bool EnableReverseAzimuth = false,
bool EnableReducedLength = false,
bool EnableGeodesicScale = false,
unsigned int Order = 4
>
class meridian_direct
{
static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale;
static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcQuantities;
static const bool CalcCoordinates = EnableCoordinates || CalcRevAzimuth;
public:
typedef result_direct<CT> result_type;
template <typename T, typename Dist, typename Spheroid>
static inline result_type apply(T const& lo1,
T const& la1,
Dist const& distance,
bool north,
Spheroid const& spheroid)
{
result_type result;
CT const half_pi = math::half_pi<CT>();
CT const pi = math::pi<CT>();
CT const one_and_a_half_pi = pi + half_pi;
CT const c0 = 0;
CT azimuth = north ? c0 : pi;
if (BOOST_GEOMETRY_CONDITION(CalcCoordinates))
{
CT s0 = meridian_inverse<CT, Order>::apply(la1, spheroid);
int signed_distance = north ? distance : -distance;
result.lon2 = lo1;
result.lat2 = apply(s0 + signed_distance, spheroid);
}
if (BOOST_GEOMETRY_CONDITION(CalcRevAzimuth))
{
result.reverse_azimuth = azimuth;
if (result.lat2 > half_pi &&
result.lat2 < one_and_a_half_pi)
{
result.reverse_azimuth = pi;
}
else if (result.lat2 < -half_pi &&
result.lat2 > -one_and_a_half_pi)
{
result.reverse_azimuth = c0;
}
}
if (BOOST_GEOMETRY_CONDITION(CalcQuantities))
{
CT const b = CT(get_radius<2>(spheroid));
CT const f = formula::flattening<CT>(spheroid);
boost::geometry::math::normalize_spheroidal_coordinates
<
boost::geometry::radian,
double
>(result.lon2, result.lat2);
typedef differential_quantities
<
CT,
EnableReducedLength,
EnableGeodesicScale,
Order
> quantities;
quantities::apply(lo1, la1, result.lon2, result.lat2,
azimuth, result.reverse_azimuth,
b, f,
result.reduced_length, result.geodesic_scale);
}
return result;
}
// https://en.wikipedia.org/wiki/Meridian_arc#The_inverse_meridian_problem_for_the_ellipsoid
// latitudes are assumed to be in radians and in [-pi/2,pi/2]
template <typename T, typename Spheroid>
static CT apply(T m, Spheroid const& spheroid)
{
CT const f = formula::flattening<CT>(spheroid);
CT n = f / (CT(2) - f);
CT mp = formula::quarter_meridian<CT>(spheroid);
CT mu = geometry::math::pi<CT>()/CT(2) * m / mp;
if (Order == 0)
{
return mu;
}
CT H2 = 1.5 * n;
if (Order == 1)
{
return mu + H2 * sin(2*mu);
}
CT n2 = n * n;
CT H4 = 1.3125 * n2;
if (Order == 2)
{
return mu + H2 * sin(2*mu) + H4 * sin(4*mu);
}
CT n3 = n2 * n;
H2 -= 0.84375 * n3;
CT H6 = 1.572916667 * n3;
if (Order == 3)
{
return mu + H2 * sin(2*mu) + H4 * sin(4*mu) + H6 * sin(6*mu);
}
CT n4 = n2 * n2;
H4 -= 1.71875 * n4;
CT H8 = 2.142578125 * n4;
// Order 4 or higher
return mu + H2 * sin(2*mu) + H4 * sin(4*mu) + H6 * sin(6*mu) + H8 * sin(8*mu);
}
};
}}} // namespace boost::geometry::formula
#endif // BOOST_GEOMETRY_FORMULAS_MERIDIAN_DIRECT_HPP

176
linx64/include/boost/geometry/formulas/meridian_inverse.hpp Normal file
View File

@@ -0,0 +1,176 @@
// Boost.Geometry
// Copyright (c) 2017-2018 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_MERIDIAN_INVERSE_HPP
#define BOOST_GEOMETRY_FORMULAS_MERIDIAN_INVERSE_HPP
#include <boost/math/constants/constants.hpp>
#include <boost/geometry/core/radius.hpp>
#include <boost/geometry/util/condition.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
#include <boost/geometry/formulas/flattening.hpp>
#include <boost/geometry/formulas/meridian_segment.hpp>
namespace boost { namespace geometry { namespace formula
{
/*!
\brief Compute the arc length of an ellipse.
*/
template <typename CT, unsigned int Order = 1>
class meridian_inverse
{
public :
struct result
{
result()
: distance(0)
, meridian(false)
{}
CT distance;
bool meridian;
};
template <typename T>
static bool meridian_not_crossing_pole(T lat1, T lat2, CT diff)
{
CT half_pi = math::pi<CT>()/CT(2);
return math::equals(diff, CT(0)) ||
(math::equals(lat2, half_pi) && math::equals(lat1, -half_pi));
}
static bool meridian_crossing_pole(CT diff)
{
return math::equals(math::abs(diff), math::pi<CT>());
}
template <typename T, typename Spheroid>
static CT meridian_not_crossing_pole_dist(T lat1, T lat2, Spheroid const& spheroid)
{
return math::abs(apply(lat2, spheroid) - apply(lat1, spheroid));
}
template <typename T, typename Spheroid>
static CT meridian_crossing_pole_dist(T lat1, T lat2, Spheroid const& spheroid)
{
CT c0 = 0;
CT half_pi = math::pi<CT>()/CT(2);
CT lat_sign = 1;
if (lat1+lat2 < c0)
{
lat_sign = CT(-1);
}
return math::abs(lat_sign * CT(2) * apply(half_pi, spheroid)
- apply(lat1, spheroid) - apply(lat2, spheroid));
}
template <typename T, typename Spheroid>
static result apply(T lon1, T lat1, T lon2, T lat2, Spheroid const& spheroid)
{
result res;
CT diff = geometry::math::longitude_distance_signed<geometry::radian>(lon1, lon2);
if (lat1 > lat2)
{
std::swap(lat1, lat2);
}
if ( meridian_not_crossing_pole(lat1, lat2, diff) )
{
res.distance = meridian_not_crossing_pole_dist(lat1, lat2, spheroid);
res.meridian = true;
}
else if ( meridian_crossing_pole(diff) )
{
res.distance = meridian_crossing_pole_dist(lat1, lat2, spheroid);
res.meridian = true;
}
return res;
}
// Distance computation on meridians using series approximations
// to elliptic integrals. Formula to compute distance from lattitude 0 to lat
// https://en.wikipedia.org/wiki/Meridian_arc
// latitudes are assumed to be in radians and in [-pi/2,pi/2]
template <typename T, typename Spheroid>
static CT apply(T lat, Spheroid const& spheroid)
{
CT const a = get_radius<0>(spheroid);
CT const f = formula::flattening<CT>(spheroid);
CT n = f / (CT(2) - f);
CT M = a/(1+n);
CT C0 = 1;
if (Order == 0)
{
return M * C0 * lat;
}
CT C2 = -1.5 * n;
if (Order == 1)
{
return M * (C0 * lat + C2 * sin(2*lat));
}
CT n2 = n * n;
C0 += .25 * n2;
CT C4 = 0.9375 * n2;
if (Order == 2)
{
return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat));
}
CT n3 = n2 * n;
C2 += 0.1875 * n3;
CT C6 = -0.729166667 * n3;
if (Order == 3)
{
return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat)
+ C6 * sin(6*lat));
}
CT n4 = n2 * n2;
C4 -= 0.234375 * n4;
CT C8 = 0.615234375 * n4;
if (Order == 4)
{
return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat)
+ C6 * sin(6*lat) + C8 * sin(8*lat));
}
CT n5 = n4 * n;
C6 += 0.227864583 * n5;
CT C10 = -0.54140625 * n5;
// Order 5 or higher
return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat)
+ C6 * sin(6*lat) + C8 * sin(8*lat) + C10 * sin(10*lat));
}
};
}}} // namespace boost::geometry::formula
#endif // BOOST_GEOMETRY_FORMULAS_MERIDIAN_INVERSE_HPP

72
linx64/include/boost/geometry/formulas/meridian_segment.hpp Normal file
View File

@@ -0,0 +1,72 @@
// Boost.Geometry
// Copyright (c) 2017-2018 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_MERIDIAN_SEGMENT_HPP
#define BOOST_GEOMETRY_FORMULAS_MERIDIAN_SEGMENT_HPP
#include <boost/math/constants/constants.hpp>
#include <boost/geometry/core/radius.hpp>
#include <boost/geometry/util/condition.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
namespace boost { namespace geometry { namespace formula
{
/*!
\brief Test if a segment is meridian or not.
*/
class meridian_segment
{
public :
enum SegmentType {NonMeridian, MeridianCrossingPole, MeridianNotCrossingPole};
template <typename T>
static inline SegmentType is_meridian(T lon1, T lat1, T lon2, T lat2)
{
SegmentType res = NonMeridian;
T diff = geometry::math::longitude_distance_signed<geometry::radian>(lon1, lon2);
if ( meridian_not_crossing_pole(lat1, lat2, diff) )
{
res = MeridianNotCrossingPole;
}
else if ( meridian_crossing_pole(diff) )
{
res = MeridianCrossingPole;
}
return res;
}
template <typename T>
static bool meridian_not_crossing_pole(T lat1, T lat2, T diff)
{
T half_pi = math::half_pi<T>();
return math::equals(diff, T(0)) ||
(math::equals(lat2, half_pi) && math::equals(lat1, -half_pi));
}
template <typename T>
static bool meridian_crossing_pole(T diff)
{
return math::equals(math::abs(diff), math::pi<T>());
}
};
}}} // namespace boost::geometry::formula
#endif //BOOST_GEOMETRY_FORMULAS_MERIDIAN_SEGMENT_HPP

112
linx64/include/boost/geometry/formulas/quarter_meridian.hpp Normal file
View File

@@ -0,0 +1,112 @@
// Boost.Geometry
// Copyright (c) 2018 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_QUARTER_MERIDIAN_HPP
#define BOOST_GEOMETRY_FORMULAS_QUARTER_MERIDIAN_HPP
#include <boost/geometry/algorithms/not_implemented.hpp>
#include <boost/geometry/core/radius.hpp>
#include <boost/geometry/core/tag.hpp>
#include <boost/geometry/core/tags.hpp>
#include <boost/geometry/formulas/flattening.hpp>
#include <boost/geometry/util/math.hpp>
namespace boost { namespace geometry
{
#ifndef DOXYGEN_NO_DISPATCH
namespace formula_dispatch
{
template <typename ResultType, typename Geometry, typename Tag = typename tag<Geometry>::type>
struct quarter_meridian
: not_implemented<Tag>
{};
template <typename ResultType, typename Geometry>
struct quarter_meridian<ResultType, Geometry, srs_spheroid_tag>
{
//https://en.wikipedia.org/wiki/Meridian_arc#Generalized_series
//http://www.wolframalpha.com/input/?i=(sum(((2*j-3)!!%2F(2*j)!!)%5E2*n%5E(2*j),j,0,8))
static inline ResultType apply(Geometry const& geometry)
{
//order 8 expansion
ResultType const C[] =
{
1073741824,
268435456,
16777216,
4194304,
1638400,
802816,
451584,
278784,
184041
};
ResultType const c2 = 2;
ResultType const c4 = 4;
ResultType const f = formula::flattening<ResultType>(geometry);
ResultType const n = f / (c2 - f);
ResultType const ab4 = (get_radius<0>(geometry)
+ get_radius<2>(geometry)) / c4;
return geometry::math::pi<ResultType>() * ab4 *
horner_evaluate(n*n, C, C+8) / C[0];
}
private :
//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);
if (begin == end)
{
return result;
}
IteratorType it = end;
do
{
result = result * x + *--it;
}
while (it != begin);
return result;
}
};
} // namespace formula_dispatch
#endif // DOXYGEN_NO_DISPATCH
#ifndef DOXYGEN_NO_DETAIL
namespace formula
{
template <typename ResultType, typename Geometry>
ResultType quarter_meridian(Geometry const& geometry)
{
return formula_dispatch::quarter_meridian<ResultType, Geometry>::apply(geometry);
}
} // namespace formula
#endif // DOXYGEN_NO_DETAIL
}} // namespace boost::geometry
#endif // BOOST_GEOMETRY_FORMULAS_QUARTER_MERIDIAN_HPP

89
linx64/include/boost/geometry/formulas/sjoberg_intersection.hpp
View File

@@ -1,6 +1,6 @@
// Boost.Geometry
// Copyright (c) 2016 Oracle and/or its affiliates.
// Copyright (c) 2016-2017 Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
@@ -15,10 +15,10 @@
#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/util/normalize_spheroidal_coordinates.hpp>
#include <boost/geometry/formulas/flattening.hpp>
#include <boost/geometry/formulas/spherical.hpp>
@@ -97,10 +97,16 @@ private:
CT const sin_dlon1 = sin(dlon1);
CT const dlon2 = lon_b2 - lon2;
CT const sin_dlon2 = sin(dlon2);
CT const cos_dlon1 = cos(dlon1);
CT const cos_dlon2 = cos(dlon2);
CT const tan_alpha1_x = cos_lat1 * tan_lat_a2 - sin_lat1 * cos_dlon1;
CT const tan_alpha2_x = cos_lat2 * tan_lat_b2 - sin_lat2 * cos_dlon2;
CT const c0 = 0;
bool const is_vertical1 = math::equals(sin_dlon1, c0);
bool const is_vertical2 = math::equals(sin_dlon2, c0);
bool const is_vertical1 = math::equals(sin_dlon1, c0) || math::equals(tan_alpha1_x, c0);
bool const is_vertical2 = math::equals(sin_dlon2, c0) || math::equals(tan_alpha2_x, c0);
CT tan_alpha1 = 0;
CT tan_alpha2 = 0;
@@ -112,27 +118,18 @@ private:
}
else if (is_vertical1)
{
CT const cos_dlon2 = cos(dlon2);
CT const tan_alpha2_x = cos_lat2 * tan_lat_b2 - sin_lat2 * cos_dlon2;
tan_alpha2 = sin_dlon2 / tan_alpha2_x;
lon = lon1;
}
else if (is_vertical2)
{
CT const cos_dlon1 = cos(dlon1);
CT const tan_alpha1_x = cos_lat1 * tan_lat_a2 - sin_lat1 * cos_dlon1;
tan_alpha1 = sin_dlon1 / tan_alpha1_x;
lon = lon2;
}
else
{
CT const cos_dlon1 = cos(dlon1);
CT const cos_dlon2 = cos(dlon2);
CT const tan_alpha1_x = cos_lat1 * tan_lat_a2 - sin_lat1 * cos_dlon1;
CT const tan_alpha2_x = cos_lat2 * tan_lat_b2 - sin_lat2 * cos_dlon2;
tan_alpha1 = sin_dlon1 / tan_alpha1_x;
tan_alpha2 = sin_dlon2 / tan_alpha2_x;
@@ -761,16 +758,36 @@ public:
sjoberg_intersection_spherical_02<CT>::apply_alt(lon_a1, lat_a1, lon_a2, lat_a2,
lon_b1, lat_b1, lon_b2, lat_b2,
lon_sph, tan_lat_sph);
// Return for sphere
if (math::equals(f, c0))
{
lon = lon_sph;
lat = atan(tan_lat_sph);
return true;
}
t = one_minus_f * tan_lat_sph; // tan(beta)
}
}
// TODO: no need to calculate atan here if reduced latitudes were used
// instead of latitudes above, in sjoberg_intersection_spherical_02
CT const beta = atan(t);
if (enable_02 && newton_method(geod1, geod2, beta, t, lon1_minus_lon2, lon, lat))
if (enable_02 && newton_method(geod1, geod2, beta, t, lon1_minus_lon2, lon_sph, lon, lat))
{
return true;
// TODO: Newton's method may return wrong result in some specific cases
// Detected for sphere and nearly sphere, e.g. A=6371228, B=6371227
// and segments s1=(-121 -19,37 8) and s2=(-19 -15,-104 -58)
// It's unclear if this is a bug or a characteristic of this method
// so until this is investigated check if the resulting longitude is
// between endpoints of the segments. It should be since before calling
// this formula sides of endpoints WRT other segments are checked.
if ( is_result_longitude_ok(geod1, lon_a1, lon_a2, lon)
&& is_result_longitude_ok(geod2, lon_b1, lon_b2, lon) )
{
return true;
}
}
return converge_07(geod1, geod2, beta, t, lon1_minus_lon2, lon_sph, lon, lat);
@@ -778,7 +795,7 @@ public:
private:
static inline bool newton_method(geodesic_type const& geod1, geodesic_type const& geod2, // in
CT beta, CT t, CT const& lon1_minus_lon2, // in
CT beta, CT t, CT const& lon1_minus_lon2, CT const& lon_sph, // in
CT & lon, CT & lat) // out
{
CT const c0 = 0;
@@ -789,6 +806,13 @@ private:
CT lon1_diff = 0;
CT lon2_diff = 0;
// The segment is vertical and intersection point is behind the vertex
// this method is unable to calculate correct result
if (geod1.is_Cj_zero && math::abs(geod1.lonj - lon_sph) > math::half_pi<CT>())
return false;
if (geod2.is_Cj_zero && math::abs(geod2.lonj - lon_sph) > math::half_pi<CT>())
return false;
CT abs_dbeta_last = 0;
// [Sjoberg02] converges faster than solution in [Sjoberg07]
@@ -840,6 +864,8 @@ private:
// In this case sin_beta/sqrt_1_Cj_sqr and t/t0j is not in [-1, 1]
// 3. f1-f2 may be 0 which means that the intermediate point is on the vertex
// In this case it's not possible to check if this is the correct result
// 4. f1-f2 may also be 0 in other cases, e.g.
// geodesics are symetrical wrt equator and longitude directions are different
CT const dbeta_denom = f1 - f2;
//CT const dbeta_denom = math::abs(f1) + math::abs(f2);
@@ -882,6 +908,28 @@ private:
return true;
}
static inline bool is_result_longitude_ok(geodesic_type const& geod,
CT const& lon1, CT const& lon2, CT const& lon)
{
CT const c0 = 0;
if (geod.is_Cj_zero)
return true; // don't check vertical segment
CT dist1p = math::longitude_distance_signed<radian>(lon1, lon);
CT dist12 = math::longitude_distance_signed<radian>(lon1, lon2);
if (dist12 < c0)
{
dist1p = -dist1p;
dist12 = -dist12;
}
return (c0 <= dist1p && dist1p <= dist12)
|| math::equals(dist1p, c0)
|| math::equals(dist1p, dist12);
}
struct geodesics_type
{
geodesics_type(geodesic_type const& g1, geodesic_type const& g2)
@@ -1126,6 +1174,7 @@ private:
geodesic_type const& geod1, geodesic_type const& geod2,
CT & t1, CT& t2)
{
CT const c0 = 0;
CT const c1 = 1;
CT const c2 = 2;
@@ -1150,8 +1199,8 @@ private:
CT const K_t01_t02 = K * geod1.t0j * geod2.t0j;
CT const D1 = math::sqrt(A + B);
CT const D2 = math::sqrt(A - B);
CT const t_new1 = K_t01_t02 / D1;
CT const t_new2 = K_t01_t02 / D2;
CT const t_new1 = math::equals(D1, c0) ? c0 : K_t01_t02 / D1;
CT const t_new2 = math::equals(D2, c0) ? c0 : K_t01_t02 / D2;
CT const t_new3 = -t_new1;
CT const t_new4 = -t_new2;

67
linx64/include/boost/geometry/formulas/spherical.hpp
View File

@@ -1,6 +1,6 @@
// Boost.Geometry
// Copyright (c) 2016, Oracle and/or its affiliates.
// Copyright (c) 2016-2018, 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
@@ -13,8 +13,10 @@
#include <boost/geometry/core/coordinate_system.hpp>
#include <boost/geometry/core/coordinate_type.hpp>
#include <boost/geometry/core/cs.hpp>
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/core/radian_access.hpp>
#include <boost/geometry/core/radius.hpp>
//#include <boost/geometry/arithmetic/arithmetic.hpp>
#include <boost/geometry/arithmetic/cross_product.hpp>
@@ -24,6 +26,8 @@
#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
#include <boost/geometry/util/select_coordinate_type.hpp>
#include <boost/geometry/formulas/result_direct.hpp>
namespace boost { namespace geometry {
namespace formula {
@@ -95,7 +99,7 @@ static inline PointSph cart3d_to_sph(Point3d const& point_3d)
math::normalize_spheroidal_coordinates
<
typename coordinate_system<PointSph>::type::units,
typename detail::cs_angular_units<PointSph>::type,
coord_t
>(lon, lat);
@@ -180,6 +184,7 @@ inline T spherical_azimuth(T const& lon1, T const& lat1, T const& lon2, T const&
template <typename T>
inline int azimuth_side_value(T const& azi_a1_p, T const& azi_a1_a2)
{
T const c0 = 0;
T const pi = math::pi<T>();
T const two_pi = math::two_pi<T>();
@@ -210,13 +215,69 @@ inline int azimuth_side_value(T const& azi_a1_p, T const& azi_a1_a2)
// the difference to 0 as well
// positive azimuth is on the right side
return math::equals(a_diff, 0)
return math::equals(a_diff, c0)
|| math::equals(a_diff, pi)
|| math::equals(a_diff, -pi) ? 0
: a_diff > 0 ? -1 // right
: 1; // left
}
template
<
bool Coordinates,
bool ReverseAzimuth,
typename CT,
typename Sphere
>
inline result_direct<CT> spherical_direct(CT const& lon1,
CT const& lat1,
CT const& sig12,
CT const& alp1,
Sphere const& sphere)
{
result_direct<CT> result;
CT const sin_alp1 = sin(alp1);
CT const sin_lat1 = sin(lat1);
CT const cos_alp1 = cos(alp1);
CT const cos_lat1 = cos(lat1);
CT const norm = math::sqrt(cos_alp1 * cos_alp1 + sin_alp1 * sin_alp1
* sin_lat1 * sin_lat1);
CT const alp0 = atan2(sin_alp1 * cos_lat1, norm);
CT const sig1 = atan2(sin_lat1, cos_alp1 * cos_lat1);
CT const sig2 = sig1 + sig12 / get_radius<0>(sphere);
CT const cos_sig2 = cos(sig2);
CT const sin_alp0 = sin(alp0);
CT const cos_alp0 = cos(alp0);
if (Coordinates)
{
CT const sin_sig2 = sin(sig2);
CT const sin_sig1 = sin(sig1);
CT const cos_sig1 = cos(sig1);
CT const norm2 = math::sqrt(cos_alp0 * cos_alp0 * cos_sig2 * cos_sig2
+ sin_alp0 * sin_alp0);
CT const lat2 = atan2(cos_alp0 * sin_sig2, norm2);
CT const omg1 = atan2(sin_alp0 * sin_sig1, cos_sig1);
CT const lon2 = atan2(sin_alp0 * sin_sig2, cos_sig2);
result.lon2 = lon1 + lon2 - omg1;
result.lat2 = lat2;
}
if (ReverseAzimuth)
{
CT const alp2 = atan2(sin_alp0, cos_alp0 * cos_sig2);
result.reverse_azimuth = alp2;
}
return result;
}
} // namespace formula
}} // namespace boost::geometry

57
linx64/include/boost/geometry/formulas/thomas_direct.hpp
View File

@@ -1,7 +1,8 @@
// Boost.Geometry
// Copyright (c) 2016 Oracle and/or its affiliates.
// Copyright (c) 2016-2018 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,
@@ -15,7 +16,6 @@
#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>
@@ -31,16 +31,16 @@ 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.
Forsyth-Andoyer-Lambert type approximation with first/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 SecondOrder = true,
bool EnableCoordinates = true,
bool EnableReverseAzimuth = false,
bool EnableReducedLength = false,
@@ -67,13 +67,6 @@ public:
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;
@@ -92,7 +85,7 @@ public:
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));
@@ -109,9 +102,18 @@ public:
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;
CT D = 0;
CT P = 0;
if ( BOOST_GEOMETRY_CONDITION(SecondOrder) )
{
D = (c1 - C2) * (c1 - C2 - C1 * M);
P = C2 * (c1 + C1 * M / c2) / D;
}
else
{
D = c1 - c2 * C2 - C1 * M;
P = 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
@@ -133,9 +135,14 @@ public:
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 X = 0;
CT d_sigma = d - Y;
if ( BOOST_GEOMETRY_CONDITION(SecondOrder) )
{
X = math::sqr(C2) * sin_d * cos_d * (2 * math::sqr(V) - c1);
d_sigma += X;
}
CT const sin_d_sigma = sin(d_sigma);
CT const cos_d_sigma = cos(d_sigma);
@@ -152,11 +159,16 @@ public:
if (BOOST_GEOMETRY_CONDITION(CalcCoordinates))
{
CT const S_sigma = c2 * sigma1 - d_sigma;
CT const cos_S_sigma = cos(S_sigma);
CT cos_S_sigma = 0;
CT H = C1 * d_sigma;
if ( BOOST_GEOMETRY_CONDITION(SecondOrder) )
{
cos_S_sigma = cos(S_sigma);
H = H * (c1 - C2) - C1 * C2 * sin_d_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))
@@ -169,7 +181,8 @@ public:
{
CT const sigma2 = S_sigma - sigma1;
//theta2 = asin(cos(sigma2)) <=> sin_theta0 = 1
CT const tan_theta2 = cos(sigma2) / sin(sigma2);
// NOTE: cos(sigma2) defines the sign of tan_theta2
CT const tan_theta2 = cos(sigma2) / math::abs(sin(sigma2));
result.lat2 = atan(tan_theta2 / one_minus_f);
}
@@ -177,7 +190,7 @@ public:
{
result.lat2 = -result.lat2;
}
}
}
if (BOOST_GEOMETRY_CONDITION(CalcQuantities))
{

8
linx64/include/boost/geometry/formulas/thomas_inverse.hpp
View File

@@ -1,7 +1,8 @@
// Boost.Geometry
// Copyright (c) 2015-2016 Oracle and/or its affiliates.
// Copyright (c) 2015-2018 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,
@@ -15,7 +16,6 @@
#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>
@@ -178,7 +178,7 @@ public:
CT const pi = math::pi<CT>();
if (BOOST_GEOMETRY_CONDITION(EnableAzimuth))
if (BOOST_GEOMETRY_CONDITION(CalcFwdAzimuth))
{
CT alpha1 = v + u;
if (alpha1 > pi)
@@ -189,7 +189,7 @@ public:
result.azimuth = alpha1;
}
if (BOOST_GEOMETRY_CONDITION(EnableReverseAzimuth))
if (BOOST_GEOMETRY_CONDITION(CalcRevAzimuth))
{
CT alpha2 = pi - (v - u);
if (alpha2 > pi)

43
linx64/include/boost/geometry/formulas/unit_spheroid.hpp Normal file
View File

@@ -0,0 +1,43 @@
// Boost.Geometry
// Copyright (c) 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_UNIT_SPHEROID_HPP
#define BOOST_GEOMETRY_FORMULAS_UNIT_SPHEROID_HPP
#include <boost/geometry/core/radius.hpp>
namespace boost { namespace geometry
{
#ifndef DOXYGEN_NO_DETAIL
namespace formula
{
template <typename ResultType, typename Spheroid>
inline ResultType unit_spheroid_b(Spheroid const& spheroid)
{
return ResultType(get_radius<2>(spheroid))
/ ResultType(get_radius<0>(spheroid));
}
template <typename ResultSpheroid, typename Spheroid>
inline ResultSpheroid unit_spheroid(Spheroid const& spheroid)
{
typedef typename radius_type<ResultSpheroid>::type radius_t;
return ResultSpheroid(radius_t(1),
unit_spheroid_b<radius_t>(spheroid));
}
} // namespace formula
#endif // DOXYGEN_NO_DETAIL
}} // namespace boost::geometry
#endif // BOOST_GEOMETRY_FORMULAS_UNIT_SPHEROID_HPP

4
linx64/include/boost/geometry/formulas/vertex_latitude.hpp
View File

@@ -12,11 +12,13 @@
#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
{

354
linx64/include/boost/geometry/formulas/vertex_longitude.hpp Normal file
View File

@@ -0,0 +1,354 @@
// Boost.Geometry
// Copyright (c) 2016-2018 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_LONGITUDE_HPP
#define BOOST_GEOMETRY_FORMULAS_MAXIMUM_LONGITUDE_HPP
#include <boost/geometry/formulas/spherical.hpp>
#include <boost/geometry/formulas/flattening.hpp>
#include <boost/mpl/assert.hpp>
#include <boost/math/special_functions/hypot.hpp>
namespace boost { namespace geometry { namespace formula
{
/*!
\brief Algorithm to compute the vertex longitude of a geodesic segment. Vertex is
a point on the geodesic that maximizes (or minimizes) the latitude. The algorithm
is given the vertex latitude.
*/
//Classes for spesific CS
template <typename CT>
class vertex_longitude_on_sphere
{
public:
template <typename T>
static inline CT apply(T const& lat1, //segment point 1
T const& lat2, //segment point 2
T const& lat3, //vertex latitude
T const& sin_l12,
T const& cos_l12) //lon1 -lon2
{
//https://en.wikipedia.org/wiki/Great-circle_navigation#Finding_way-points
CT const A = sin(lat1) * cos(lat2) * cos(lat3) * sin_l12;
CT const B = sin(lat1) * cos(lat2) * cos(lat3) * cos_l12
- cos(lat1) * sin(lat2) * cos(lat3);
CT lon = atan2(B, A);
return lon + math::pi<CT>();
}
};
template <typename CT>
class vertex_longitude_on_spheroid
{
template<typename T>
static inline void normalize(T& x, T& y)
{
T h = boost::math::hypot(x, y);
x /= h;
y /= h;
}
public:
template <typename T, typename Spheroid>
static inline CT apply(T const& lat1, //segment point 1
T const& lat2, //segment point 2
T const& lat3, //vertex latitude
T& alp1,
Spheroid const& spheroid)
{
// We assume that segment points lay on different side w.r.t.
// the vertex
// Constants
CT const c0 = 0;
CT const c2 = 2;
CT const half_pi = math::pi<CT>() / c2;
if (math::equals(lat1, half_pi)
|| math::equals(lat2, half_pi)
|| math::equals(lat1, -half_pi)
|| math::equals(lat2, -half_pi))
{
// one segment point is the pole
return c0;
}
// More constants
CT const f = flattening<CT>(spheroid);
CT const pi = math::pi<CT>();
CT const c1 = 1;
CT const cminus1 = -1;
// First, compute longitude on auxiliary sphere
CT const one_minus_f = c1 - f;
CT const bet1 = atan(one_minus_f * tan(lat1));
CT const bet2 = atan(one_minus_f * tan(lat2));
CT const bet3 = atan(one_minus_f * tan(lat3));
CT cos_bet1 = cos(bet1);
CT cos_bet2 = cos(bet2);
CT const sin_bet1 = sin(bet1);
CT const sin_bet2 = sin(bet2);
CT const sin_bet3 = sin(bet3);
CT omg12 = 0;
if (bet1 < c0)
{
cos_bet1 *= cminus1;
omg12 += pi;
}
if (bet2 < c0)
{
cos_bet2 *= cminus1;
omg12 += pi;
}
CT const sin_alp1 = sin(alp1);
CT const cos_alp1 = math::sqrt(c1 - math::sqr(sin_alp1));
CT const norm = math::sqrt(math::sqr(cos_alp1) + math::sqr(sin_alp1 * sin_bet1));
CT const sin_alp0 = sin(atan2(sin_alp1 * cos_bet1, norm));
BOOST_ASSERT(cos_bet2 != c0);
CT const sin_alp2 = sin_alp1 * cos_bet1 / cos_bet2;
CT const cos_alp0 = math::sqrt(c1 - math::sqr(sin_alp0));
CT const cos_alp2 = math::sqrt(c1 - math::sqr(sin_alp2));
CT const sig1 = atan2(sin_bet1, cos_alp1 * cos_bet1);
CT const sig2 = atan2(sin_bet2, -cos_alp2 * cos_bet2); //lat3 is a vertex
CT const cos_sig1 = cos(sig1);
CT const sin_sig1 = math::sqrt(c1 - math::sqr(cos_sig1));
CT const cos_sig2 = cos(sig2);
CT const sin_sig2 = math::sqrt(c1 - math::sqr(cos_sig2));
CT const omg1 = atan2(sin_alp0 * sin_sig1, cos_sig1);
CT const omg2 = atan2(sin_alp0 * sin_sig2, cos_sig2);
omg12 += omg1 - omg2;
CT const sin_omg12 = sin(omg12);
CT const cos_omg12 = cos(omg12);
CT omg13 = geometry::formula::vertex_longitude_on_sphere<CT>
::apply(bet1, bet2, bet3, sin_omg12, cos_omg12);
if (lat1 * lat2 < c0)//different hemispheres
{
if ((lat2 - lat1) * lat3 > c0)// ascending segment
{
omg13 = pi - omg13;
}
}
// Second, compute the ellipsoidal longitude
CT const e2 = f * (c2 - f);
CT const ep = math::sqrt(e2 / (c1 - e2));
CT const k2 = math::sqr(ep * cos_alp0);
CT const sqrt_k2_plus_one = math::sqrt(c1 + k2);
CT const eps = (sqrt_k2_plus_one - c1) / (sqrt_k2_plus_one + c1);
CT const eps2 = eps * eps;
CT const n = f / (c2 - f);
// sig3 is the length from equator to the vertex
CT sig3;
if(sin_bet3 > c0)
{
sig3 = half_pi;
} else {
sig3 = -half_pi;
}
CT const cos_sig3 = 0;
CT const sin_sig3 = 1;
CT sig13 = sig3 - sig1;
if (sig13 > pi)
{
sig13 -= 2 * pi;
}
// Order 2 approximation
CT const c1over2 = 0.5;
CT const c1over4 = 0.25;
CT const c1over8 = 0.125;
CT const c1over16 = 0.0625;
CT const c4 = 4;
CT const c8 = 8;
CT const A3 = 1 - (c1over2 - c1over2 * n) * eps - c1over4 * eps2;
CT const C31 = (c1over4 - c1over4 * n) * eps + c1over8 * eps2;
CT const C32 = c1over16 * eps2;
CT const sin2_sig3 = c2 * cos_sig3 * sin_sig3;
CT const sin4_sig3 = sin_sig3 * (-c4 * cos_sig3
+ c8 * cos_sig3 * cos_sig3 * cos_sig3);
CT const sin2_sig1 = c2 * cos_sig1 * sin_sig1;
CT const sin4_sig1 = sin_sig1 * (-c4 * cos_sig1
+ c8 * cos_sig1 * cos_sig1 * cos_sig1);
CT const I3 = A3 * (sig13
+ C31 * (sin2_sig3 - sin2_sig1)
+ C32 * (sin4_sig3 - sin4_sig1));
CT const sign = bet3 >= c0
? c1
: cminus1;
CT const dlon_max = omg13 - sign * f * sin_alp0 * I3;
return dlon_max;
}
};
//CS_tag dispatching
template <typename CT, typename CS_Tag>
struct compute_vertex_lon
{
BOOST_MPL_ASSERT_MSG
(
false, NOT_IMPLEMENTED_FOR_THIS_COORDINATE_SYSTEM, (types<CS_Tag>)
);
};
template <typename CT>
struct compute_vertex_lon<CT, spherical_equatorial_tag>
{
template <typename Strategy>
static inline CT apply(CT const& lat1,
CT const& lat2,
CT const& vertex_lat,
CT const& sin_l12,
CT const& cos_l12,
CT,
Strategy)
{
return vertex_longitude_on_sphere<CT>
::apply(lat1,
lat2,
vertex_lat,
sin_l12,
cos_l12);
}
};
template <typename CT>
struct compute_vertex_lon<CT, geographic_tag>
{
template <typename Strategy>
static inline CT apply(CT const& lat1,
CT const& lat2,
CT const& vertex_lat,
CT,
CT,
CT& alp1,
Strategy const& azimuth_strategy)
{
return vertex_longitude_on_spheroid<CT>
::apply(lat1,
lat2,
vertex_lat,
alp1,
azimuth_strategy.model());
}
};
// Vertex longitude interface
// Assume that lon1 < lon2 and vertex_lat is the latitude of the vertex
template <typename CT, typename CS_Tag>
class vertex_longitude
{
public :
template <typename Strategy>
static inline CT apply(CT& lon1,
CT& lat1,
CT& lon2,
CT& lat2,
CT const& vertex_lat,
CT& alp1,
Strategy const& azimuth_strategy)
{
CT const c0 = 0;
CT pi = math::pi<CT>();
//Vertex is a segment's point
if (math::equals(vertex_lat, lat1))
{
return lon1;
}
if (math::equals(vertex_lat, lat2))
{
return lon2;
}
//Segment lay on meridian
if (math::equals(lon1, lon2))
{
return (std::max)(lat1, lat2);
}
BOOST_ASSERT(lon1 < lon2);
CT dlon = compute_vertex_lon<CT, CS_Tag>::apply(lat1, lat2,
vertex_lat,
sin(lon1 - lon2),
cos(lon1 - lon2),
alp1,
azimuth_strategy);
CT vertex_lon = std::fmod(lon1 + dlon, 2 * pi);
if (vertex_lat < c0)
{
vertex_lon -= pi;
}
if (std::abs(lon1 - lon2) > pi)
{
vertex_lon -= pi;
}
return vertex_lon;
}
};
template <typename CT>
class vertex_longitude<CT, cartesian_tag>
{
public :
template <typename Strategy>
static inline CT apply(CT& /*lon1*/,
CT& /*lat1*/,
CT& lon2,
CT& /*lat2*/,
CT const& /*vertex_lat*/,
CT& /*alp1*/,
Strategy const& /*azimuth_strategy*/)
{
return lon2;
}
};
}}} // namespace boost::geometry::formula
#endif // BOOST_GEOMETRY_FORMULAS_MAXIMUM_LONGITUDE_HPP

12
linx64/include/boost/geometry/formulas/vincenty_direct.hpp
View File

@@ -2,8 +2,8 @@
// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
// This file was modified by Oracle on 2014, 2016.
// Modifications copyright (c) 2014-2016 Oracle and/or its affiliates.
// This file was modified by Oracle on 2014, 2016, 2017.
// Modifications copyright (c) 2014-2017 Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
@@ -18,7 +18,6 @@
#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>
@@ -75,13 +74,6 @@ public:
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);

9
linx64/include/boost/geometry/formulas/vincenty_inverse.hpp
View File

@@ -1,6 +1,7 @@
// Boost.Geometry
// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
// Copyright (c) 2018 Adam Wulkiewicz, Lodz, Poland.
// This file was modified by Oracle on 2014, 2016, 2017.
// Modifications copyright (c) 2014-2017 Oracle and/or its affiliates.
@@ -18,7 +19,6 @@
#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>
@@ -79,6 +79,7 @@ public:
return result;
}
CT const c0 = 0;
CT const c1 = 1;
CT const c2 = 2;
CT const c3 = 3;
@@ -145,7 +146,7 @@ public:
CT cos_sigma = sin_U1 * sin_U2 + cos_U1 * cos_U2 * cos_lambda; // (15)
sin_alpha = cos_U1 * cos_U2 * sin_lambda / sin_sigma; // (17)
cos2_alpha = c1 - math::sqr(sin_alpha);
cos_2sigma_m = math::equals(cos2_alpha, 0) ? 0 : cos_sigma - c2 * sin_U1 * sin_U2 / cos2_alpha; // (18)
cos_2sigma_m = math::equals(cos2_alpha, c0) ? c0 : 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)
@@ -163,10 +164,6 @@ public:
{
// 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;