summaryrefslogtreecommitdiff
path: root/boost/geometry/formulas
diff options
context:
space:
mode:
Diffstat (limited to 'boost/geometry/formulas')
-rw-r--r--boost/geometry/formulas/area_formulas.hpp170
-rw-r--r--boost/geometry/formulas/differential_quantities.hpp27
-rw-r--r--boost/geometry/formulas/geographic.hpp16
-rw-r--r--boost/geometry/formulas/gnomonic_intersection.hpp11
-rw-r--r--boost/geometry/formulas/gnomonic_spheroid.hpp6
-rw-r--r--boost/geometry/formulas/interpolate_point_spherical.hpp5
-rw-r--r--boost/geometry/formulas/karney_direct.hpp49
-rw-r--r--boost/geometry/formulas/karney_inverse.hpp9
-rw-r--r--boost/geometry/formulas/meridian_direct.hpp11
-rw-r--r--boost/geometry/formulas/meridian_inverse.hpp13
-rw-r--r--boost/geometry/formulas/sjoberg_intersection.hpp43
-rw-r--r--boost/geometry/formulas/spherical.hpp2
-rw-r--r--boost/geometry/formulas/thomas_inverse.hpp6
-rw-r--r--boost/geometry/formulas/vertex_longitude.hpp2
-rw-r--r--boost/geometry/formulas/vincenty_direct.hpp2
-rw-r--r--boost/geometry/formulas/vincenty_inverse.hpp4
16 files changed, 199 insertions, 177 deletions
diff --git a/boost/geometry/formulas/area_formulas.hpp b/boost/geometry/formulas/area_formulas.hpp
index d459528c27..b8da8a6727 100644
--- a/boost/geometry/formulas/area_formulas.hpp
+++ b/boost/geometry/formulas/area_formulas.hpp
@@ -1,7 +1,8 @@
// Boost.Geometry
-// Copyright (c) 2015-2021 Oracle and/or its affiliates.
+// Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
+// Copyright (c) 2015-2022 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
@@ -342,6 +343,15 @@ public:
}
}
+ static inline CT trapezoidal_formula(CT lat1r, CT lat2r, CT lon21r)
+ {
+ CT const c1 = CT(1);
+ CT const c2 = CT(2);
+ CT const tan_lat1 = tan(lat1r / c2);
+ CT const tan_lat2 = tan(lat2r / c2);
+
+ return c2 * atan(((tan_lat1 + tan_lat2) / (c1 + tan_lat1 * tan_lat2))* tan(lon21r / c2));
+ }
/*
Compute the spherical excess of a geodesic (or shperical) segment
@@ -354,43 +364,44 @@ public:
static inline CT spherical(PointOfSegment const& p1,
PointOfSegment const& p2)
{
+ CT const pi = math::pi<CT>();
+
CT excess;
- if (LongSegment) // not for segments parallel to equator
+ CT const lon1r = get_as_radian<0>(p1);
+ CT const lat1r = get_as_radian<1>(p1);
+ CT const lon2r = get_as_radian<0>(p2);
+ CT const lat2r = get_as_radian<1>(p2);
+
+ CT lon12r = lon2r - lon1r;
+ math::normalize_longitude<radian, CT>(lon12r);
+
+ if (lon12r == pi || lon12r == -pi)
{
- 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));
+ return pi;
+ }
- CT omg12 = geometry::get_as_radian<0>(p1)
- - geometry::get_as_radian<0>(p2);
- CT comg12 = cos(omg12);
- CT somg12 = sin(omg12);
+ if (BOOST_GEOMETRY_CONDITION(LongSegment) && lat1r != lat2r) // not for segments parallel to equator
+ {
+ CT const cbet1 = cos(lat1r);
+ CT const sbet1 = sin(lat1r);
+ CT const cbet2 = cos(lat2r);
+ CT const sbet2 = sin(lat2r);
- CT alp1 = atan2(cbet1 * sbet2
- - sbet1 * cbet2 * comg12,
- cbet2 * somg12);
+ CT const omg12 = lon2r - lon1r;
+ CT const comg12 = cos(omg12);
+ CT const somg12 = sin(omg12);
- CT alp2 = atan2(cbet1 * sbet2 * comg12
- - sbet1 * cbet2,
- cbet1 * somg12);
+ CT const cbet1_sbet2 = cbet1 * sbet2;
+ CT const sbet1_cbet2 = sbet1 * cbet2;
+ CT const alp1 = atan2(cbet1_sbet2 - sbet1_cbet2 * comg12, cbet2 * somg12);
+ CT const 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));
+ excess = trapezoidal_formula(lat1r, lat2r, lon12r);
}
return excess;
@@ -437,8 +448,13 @@ public:
// Constants
+ CT const c0 = CT(0);
+ CT const c1 = CT(1);
+ CT const c2 = CT(2);
+ CT const pi = math::pi<CT>();
+ CT const half_pi = pi / c2;
CT const ep = spheroid_const.m_ep;
- CT const one_minus_f = CT(1) - spheroid_const.m_f;
+ CT const one_minus_f = c1 - spheroid_const.m_f;
// Basic trigonometric computations
// the compiler could optimize here using sincos function
@@ -472,43 +488,47 @@ public:
CT excess;
- auto const half_pi = math::pi<CT>() / 2;
- bool meridian = lon2r - lon1r == CT(0)
- || lat1r == half_pi || lat1r == -half_pi
- || lat2r == half_pi || lat2r == -half_pi;
+ CT lon12r = lon2r - lon1r;
+ math::normalize_longitude<radian, CT>(lon12r);
- if (!meridian && (i_res.distance)
- < mean_radius<CT>(spheroid_const.m_spheroid) / CT(638)) // short segment
+ // Comparing with "==" works with all test cases here, but could potential create numerical issues
+ if (lon12r == pi || lon12r == -pi)
{
- CT tan_lat1 = tan(lat1r / 2.0);
- CT tan_lat2 = tan(lat2r / 2.0);
-
- excess = CT(2.0)
- * atan(((tan_lat1 + tan_lat2) / (CT(1) + tan_lat1 * tan_lat2))
- * tan((lon2r - lon1r) / 2));
+ result.spherical_term = pi;
}
else
{
- /* in some cases this formula gives more accurate results
- *
- * 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));
- */
-
- excess = alp2 - alp1;
+ bool const meridian = lon12r == c0
+ || lat1r == half_pi || lat1r == -half_pi
+ || lat2r == half_pi || lat2r == -half_pi;
+
+ if (!meridian && (i_res.distance)
+ < mean_radius<CT>(spheroid_const.m_spheroid) / CT(638)) // short segment
+ {
+ excess = trapezoidal_formula(lat1r, lat2r, lon12r);
+ }
+ else
+ {
+ /* in some cases this formula gives more accurate results
+ 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));
+ */
+
+ excess = alp2 - alp1;
+ }
+
+ result.spherical_term = excess;
}
- result.spherical_term = excess;
-
// Ellipsoidal term computation (uses integral approximation)
- CT const cos_alp0 = math::sqrt(CT(1) - math::sqr(sin_alp0));
+ CT const cos_alp0 = math::sqrt(c1 - math::sqr(sin_alp0));
//CT const cos_alp0 = hypot(cos_alp1, sin_alp1 * sin_bet1);
CT cos_sig1 = cos_alp1 * cos_bet1;
CT cos_sig2 = cos_alp2 * cos_bet2;
@@ -523,8 +543,8 @@ public:
if (ExpandEpsN) // expand by eps and n
{
CT const k2 = math::sqr(ep * cos_alp0);
- CT const sqrt_k2_plus_one = math::sqrt(CT(1) + k2);
- CT const eps = (sqrt_k2_plus_one - CT(1)) / (sqrt_k2_plus_one + CT(1));
+ CT const sqrt_k2_plus_one = math::sqrt(c1 + k2);
+ CT const eps = (sqrt_k2_plus_one - c1) / (sqrt_k2_plus_one + c1);
// Generate and evaluate the polynomials on eps (i.e. var2 = eps)
// to get the final series coefficients
@@ -563,20 +583,26 @@ public:
static inline bool crosses_prime_meridian(PointOfSegment const& p1,
PointOfSegment const& p2)
{
- 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);
+ CT const pi = geometry::math::pi<CT>();
+ CT const two_pi = geometry::math::two_pi<CT>();
+
+ CT const lon1r = get_as_radian<0>(p1);
+ CT const lon2r = get_as_radian<0>(p2);
+
+ CT lon12 = lon2r - lon1r;
+ math::normalize_longitude<radian, CT>(lon12);
+
+ // Comparing with "==" works with all test cases here, but could potential create numerical issues
+ if (lon12 == pi || lon12 == -pi)
+ {
+ return true;
+ }
+
+ CT const p1_lon = lon1r - ( std::floor( lon1r / two_pi ) * two_pi );
+ CT const p2_lon = lon2r - ( std::floor( lon2r / two_pi ) * two_pi );
+
+ CT const max_lon = (std::max)(p1_lon, p2_lon);
+ CT const min_lon = (std::min)(p1_lon, p2_lon);
return max_lon > pi && min_lon < pi && max_lon - min_lon > pi;
}
diff --git a/boost/geometry/formulas/differential_quantities.hpp b/boost/geometry/formulas/differential_quantities.hpp
index 789275f145..1df9adb8e6 100644
--- a/boost/geometry/formulas/differential_quantities.hpp
+++ b/boost/geometry/formulas/differential_quantities.hpp
@@ -1,7 +1,8 @@
// Boost.Geometry
-// Copyright (c) 2016-2019 Oracle and/or its affiliates.
+// Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
+// Copyright (c) 2016-2019 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,
@@ -67,7 +68,7 @@ public:
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))
{
@@ -80,7 +81,7 @@ public:
CT m12 = azi_sign * sin(sig_12) * b;
reduced_length = m12;
}
-
+
if (BOOST_GEOMETRY_CONDITION(EnableGeodesicScale))
{
CT M12 = cos(sig_12);
@@ -162,7 +163,7 @@ private:
CT const& sin_sig2, CT const& cos_sig2,
CT const& cos_alp0_sqr, CT const& f)
{
- if (Order == 0)
+ if (BOOST_GEOMETRY_CONDITION(Order == 0))
{
return 0;
}
@@ -176,15 +177,15 @@ private:
CT const sin_2sig_12 = sin_2sig2 - sin_2sig1;
CT const L1 = sig_12 - sin_2sig_12 / c2;
- if (Order == 1)
+ if (BOOST_GEOMETRY_CONDITION(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;
@@ -195,7 +196,7 @@ private:
+ (c12 * cos_alp0_sqr - c24) * sig_12)
/ c16;
- if (Order == 2)
+ if (BOOST_GEOMETRY_CONDITION(Order == 2))
{
return cos_alp0_sqr * f * (L1 + f * L2);
}
@@ -241,7 +242,7 @@ private:
CT const& sin_sig2, CT const& cos_sig2,
CT const& cos_alp0_sqr, CT const& ep_sqr)
{
- if (Order == 0)
+ if (BOOST_GEOMETRY_CONDITION(Order == 0))
{
return 0;
}
@@ -259,21 +260,21 @@ private:
CT const L1 = (c2 * sig_12 - sin_2sig_12) / c4;
- if (Order == 1)
+ if (BOOST_GEOMETRY_CONDITION(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)
+ if (BOOST_GEOMETRY_CONDITION(Order == 2))
{
return c2a0ep2 * (L1 + c2a0ep2 * L2);
}
diff --git a/boost/geometry/formulas/geographic.hpp b/boost/geometry/formulas/geographic.hpp
index 9501526ba8..cc726876aa 100644
--- a/boost/geometry/formulas/geographic.hpp
+++ b/boost/geometry/formulas/geographic.hpp
@@ -29,7 +29,7 @@
#include <boost/geometry/util/select_coordinate_type.hpp>
namespace boost { namespace geometry {
-
+
namespace formula {
template <typename Point3d, typename PointGeo, typename Spheroid>
@@ -107,11 +107,11 @@ inline PointGeo cart3d_to_geo(Point3d const& point_3d, Spheroid const& spheroid)
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
@@ -141,8 +141,8 @@ inline PointGeo cart3d_to_geo(Point3d const& point_3d, Spheroid const& spheroid)
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;
-
+ typedef typename coordinate_type<Point3d>::type coord_t;
+
// len_xy = sqrt(x^2 + y^2)
// r = len_xy - |z / tan(lat)|
// assuming h = 0
@@ -151,7 +151,7 @@ inline Point3d projected_to_xy(Point3d const& point_3d, Spheroid const& spheroid
// 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);
@@ -178,7 +178,7 @@ inline Point3d projected_to_surface(Point3d const& direction, Spheroid const& sp
//(x*x+y*y)/(a*a) + z*z/(b*b) = 1
// x = d.x * t
// y = d.y * t
- // z = d.z * 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);
@@ -217,7 +217,7 @@ inline bool projected_to_surface(Point3d const& origin, Point3d const& direction
//(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
+ // 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);
diff --git a/boost/geometry/formulas/gnomonic_intersection.hpp b/boost/geometry/formulas/gnomonic_intersection.hpp
index 7e1b7bcdab..64795f7625 100644
--- a/boost/geometry/formulas/gnomonic_intersection.hpp
+++ b/boost/geometry/formulas/gnomonic_intersection.hpp
@@ -1,7 +1,8 @@
// Boost.Geometry
-// Copyright (c) 2016 Oracle and/or its affiliates.
+// Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
+// 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,
@@ -60,7 +61,7 @@ public:
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,
@@ -80,7 +81,7 @@ public:
CT xa1, ya1, xa2, ya2;
CT xb1, yb1, xb2, yb2;
CT x, y;
- double lat1, lon1;
+ CT lon1, lat1;
bool ok = gnom_t::forward(lon, lat, lona1, lata1, xa1, ya1, spheroid)
&& gnom_t::forward(lon, lat, lona2, lata2, xa2, ya2, spheroid)
@@ -93,7 +94,7 @@ public:
{
return false;
}
-
+
if (math::equals(lat1, lat) && math::equals(lon1, lon))
{
break;
@@ -137,7 +138,7 @@ private:
x = get<0>(p) / z;
y = get<1>(p) / z;
-
+
return true;
}
};
diff --git a/boost/geometry/formulas/gnomonic_spheroid.hpp b/boost/geometry/formulas/gnomonic_spheroid.hpp
index 4710c6c063..75544fc7db 100644
--- a/boost/geometry/formulas/gnomonic_spheroid.hpp
+++ b/boost/geometry/formulas/gnomonic_spheroid.hpp
@@ -66,7 +66,7 @@ public:
CT rho = m / M;
x = sin(i_res.azimuth) * rho;
y = cos(i_res.azimuth) * rho;
-
+
return true;
}
@@ -78,7 +78,7 @@ public:
{
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);
@@ -95,7 +95,7 @@ public:
// 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;
diff --git a/boost/geometry/formulas/interpolate_point_spherical.hpp b/boost/geometry/formulas/interpolate_point_spherical.hpp
index 2e712d44ce..ad0b746997 100644
--- a/boost/geometry/formulas/interpolate_point_spherical.hpp
+++ b/boost/geometry/formulas/interpolate_point_spherical.hpp
@@ -1,6 +1,6 @@
// Boost.Geometry
-// Copyright (c) 2019-2021 Oracle and/or its affiliates.
+// Copyright (c) 2019-2023 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
@@ -12,6 +12,7 @@
#ifndef BOOST_GEOMETRY_FORMULAS_INTERPOLATE_POINT_SPHERICAL_HPP
#define BOOST_GEOMETRY_FORMULAS_INTERPOLATE_POINT_SPHERICAL_HPP
+#include <boost/geometry/arithmetic/normalize.hpp>
#include <boost/geometry/formulas/spherical.hpp>
#include <boost/geometry/geometries/point.hpp>
@@ -29,7 +30,7 @@ public :
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);
diff --git a/boost/geometry/formulas/karney_direct.hpp b/boost/geometry/formulas/karney_direct.hpp
index 463697e4bc..bfe35ddc49 100644
--- a/boost/geometry/formulas/karney_direct.hpp
+++ b/boost/geometry/formulas/karney_direct.hpp
@@ -1,12 +1,14 @@
// Boost.Geometry
// Copyright (c) 2018 Adeel Ahmad, Islamabad, Pakistan.
+// Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
// Contributed and/or modified by Adeel Ahmad,
// as part of Google Summer of Code 2018 program.
-// This file was modified by Oracle on 2018-2020.
-// Modifications copyright (c) 2018-2020 Oracle and/or its affiliates.
+// This file was modified by Oracle on 2018-2022.
+// Modifications copyright (c) 2018-2022 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,
@@ -30,8 +32,6 @@
#define BOOST_GEOMETRY_FORMULAS_KARNEY_DIRECT_HPP
-#include <boost/array.hpp>
-
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/hypot.hpp>
@@ -82,10 +82,10 @@ public:
{
result_type result;
- CT lon1 = lo1;
- CT const lat1 = la1;
+ CT lon1 = lo1 * math::r2d<CT>();
+ CT const lat1 = la1 * math::r2d<CT>();
- Azi azi12 = azimuth12;
+ Azi azi12 = azimuth12 * math::r2d<CT>();
math::normalize_azimuth<degree, Azi>(azi12);
CT const c0 = 0;
@@ -151,10 +151,9 @@ public:
// 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 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);
@@ -169,9 +168,6 @@ public:
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))
@@ -182,9 +178,6 @@ public:
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;
@@ -201,15 +194,11 @@ public:
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));
+ CT const sin_cos_res = se::sin_cos_series(sin_sigma2, cos_sigma2, coeffs_C3);
+ CT const lam12 = omega12 + A3c * (sigma12 + (sin_cos_res - B31));
- // Convert to radians to get the
- // longitudinal difference.
- CT lon12 = lam12 / math::d2r<CT>();
+ // Convert to degrees to get the longitudinal difference.
+ CT lon12 = lam12 * math::r2d<CT>();
// Add the longitude at first point to the longitudinal
// difference and normalize the result.
@@ -224,6 +213,8 @@ public:
// otherwise differential quantities are calculated incorrectly.
// But here it's ok since result.lon2 is not used after this point.
math::normalize_longitude<degree, CT>(result.lon2);
+
+ result.lon2 *= math::d2r<CT>();
}
if (BOOST_GEOMETRY_CONDITION(CalcQuantities))
@@ -252,12 +243,10 @@ public:
cos_sigma1 * cos_sigma2 * J12);
// Find the geodesic scale.
- CT const t = k2 * (sin_sigma2 - sin_sigma1) *
- (sin_sigma2 + sin_sigma1) / (dn1 + dn2);
+ 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;
+ result.geodesic_scale = cos_sigma12 + (t * sin_sigma2 - cos_sigma2 * J12) *
+ sin_sigma1 / dn1;
}
return result;
diff --git a/boost/geometry/formulas/karney_inverse.hpp b/boost/geometry/formulas/karney_inverse.hpp
index 936fbe22c4..47efc9ff5e 100644
--- a/boost/geometry/formulas/karney_inverse.hpp
+++ b/boost/geometry/formulas/karney_inverse.hpp
@@ -1,6 +1,7 @@
// Boost.Geometry
// Copyright (c) 2018 Adeel Ahmad, Islamabad, Pakistan.
+// Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
// Contributed and/or modified by Adeel Ahmad, as part of Google Summer of Code 2018 program.
@@ -267,8 +268,8 @@ public:
CT sin_sigma2 = sin_beta2;
CT cos_sigma2 = cos_alpha2 * cos_beta2;
- CT sigma12 = std::atan2((std::max)(c0, cos_sigma1 * sin_sigma2 - sin_sigma1 * cos_sigma2),
- cos_sigma1 * cos_sigma2 + sin_sigma1 * sin_sigma2);
+ sigma12 = std::atan2((std::max)(c0, cos_sigma1 * sin_sigma2 - sin_sigma1 * cos_sigma2),
+ cos_sigma1 * cos_sigma2 + sin_sigma1 * sin_sigma2);
CT dummy;
meridian_length(n, ep2, sigma12, sin_sigma1, cos_sigma1, dn1,
@@ -281,7 +282,7 @@ public:
{
if (sigma12 < c3 * tiny)
{
- sigma12 = m12x = s12x = c0;
+ sigma12 = m12x = s12x = c0;
}
m12x *= b;
@@ -381,7 +382,7 @@ public:
cos_alpha1 / sin_alpha1 > cos_alpha1b / sin_alpha1b))
{
sin_alpha1b = sin_alpha1;
- cos_alpha1b = cos_alpha1;
+ cos_alpha1b = cos_alpha1;
}
else if (v < c0 && (iteration > max_iterations ||
cos_alpha1 / sin_alpha1 < cos_alpha1a / sin_alpha1a))
diff --git a/boost/geometry/formulas/meridian_direct.hpp b/boost/geometry/formulas/meridian_direct.hpp
index f8e73f57c6..330af6703f 100644
--- a/boost/geometry/formulas/meridian_direct.hpp
+++ b/boost/geometry/formulas/meridian_direct.hpp
@@ -1,7 +1,8 @@
// Boost.Geometry
-// Copyright (c) 2018 Oracle and/or its affiliates.
+// Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
+// 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
@@ -127,14 +128,14 @@ public:
CT mp = formula::quarter_meridian<CT>(spheroid);
CT mu = geometry::math::pi<CT>()/CT(2) * m / mp;
- if (Order == 0)
+ if (BOOST_GEOMETRY_CONDITION(Order == 0))
{
return mu;
}
CT H2 = 1.5 * n;
- if (Order == 1)
+ if (BOOST_GEOMETRY_CONDITION(Order == 1))
{
return mu + H2 * sin(2*mu);
}
@@ -142,7 +143,7 @@ public:
CT n2 = n * n;
CT H4 = 1.3125 * n2;
- if (Order == 2)
+ if (BOOST_GEOMETRY_CONDITION(Order == 2))
{
return mu + H2 * sin(2*mu) + H4 * sin(4*mu);
}
@@ -151,7 +152,7 @@ public:
H2 -= 0.84375 * n3;
CT H6 = 1.572916667 * n3;
- if (Order == 3)
+ if (BOOST_GEOMETRY_CONDITION(Order == 3))
{
return mu + H2 * sin(2*mu) + H4 * sin(4*mu) + H6 * sin(6*mu);
}
diff --git a/boost/geometry/formulas/meridian_inverse.hpp b/boost/geometry/formulas/meridian_inverse.hpp
index 43bec19970..c21196f5c9 100644
--- a/boost/geometry/formulas/meridian_inverse.hpp
+++ b/boost/geometry/formulas/meridian_inverse.hpp
@@ -1,7 +1,8 @@
// Boost.Geometry
-// Copyright (c) 2017-2018 Oracle and/or its affiliates.
+// Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
+// 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,
@@ -118,14 +119,14 @@ public :
CT M = a/(1+n);
CT C0 = 1;
- if (Order == 0)
+ if (BOOST_GEOMETRY_CONDITION(Order == 0))
{
return M * C0 * lat;
}
CT C2 = -1.5 * n;
- if (Order == 1)
+ if (BOOST_GEOMETRY_CONDITION(Order == 1))
{
return M * (C0 * lat + C2 * sin(2*lat));
}
@@ -134,7 +135,7 @@ public :
C0 += .25 * n2;
CT C4 = 0.9375 * n2;
- if (Order == 2)
+ if (BOOST_GEOMETRY_CONDITION(Order == 2))
{
return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat));
}
@@ -143,7 +144,7 @@ public :
C2 += 0.1875 * n3;
CT C6 = -0.729166667 * n3;
- if (Order == 3)
+ if (BOOST_GEOMETRY_CONDITION(Order == 3))
{
return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat)
+ C6 * sin(6*lat));
@@ -153,7 +154,7 @@ public :
C4 -= 0.234375 * n4;
CT C8 = 0.615234375 * n4;
- if (Order == 4)
+ if (BOOST_GEOMETRY_CONDITION(Order == 4))
{
return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat)
+ C6 * sin(6*lat) + C8 * sin(8*lat));
diff --git a/boost/geometry/formulas/sjoberg_intersection.hpp b/boost/geometry/formulas/sjoberg_intersection.hpp
index 12abb76757..af56f285f4 100644
--- a/boost/geometry/formulas/sjoberg_intersection.hpp
+++ b/boost/geometry/formulas/sjoberg_intersection.hpp
@@ -1,7 +1,8 @@
// Boost.Geometry
-// Copyright (c) 2016-2019 Oracle and/or its affiliates.
+// Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
+// Copyright (c) 2016-2019 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,
@@ -71,7 +72,7 @@ struct sjoberg_intersection_spherical_02
CT const tan_lat_a2 = tan(lat_a2);
CT const tan_lat_b2 = tan(lat_b2);
-
+
return apply(lon1, lon_a2, lon2, lon_b2,
sin_lon1, cos_lon1, sin_lat1, cos_lat1,
sin_lon2, cos_lon2, sin_lat2, cos_lat2,
@@ -103,7 +104,7 @@ private:
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) || math::equals(tan_alpha1_x, c0);
bool const is_vertical2 = math::equals(sin_dlon2, c0) || math::equals(tan_alpha2_x, c0);
@@ -132,7 +133,7 @@ private:
{
tan_alpha1 = sin_dlon1 / tan_alpha1_x;
tan_alpha2 = sin_dlon2 / tan_alpha2_x;
-
+
CT const T1 = tan_alpha1 * cos_lat1;
CT const T2 = tan_alpha2 * cos_lat2;
CT const T1T2 = T1*T2;
@@ -152,12 +153,12 @@ private:
CT const lon_dist2 = (std::max)((std::min)(math::longitude_difference<radian>(lon1, lon_2),
math::longitude_difference<radian>(lon_a2, lon_2)),
(std::min)(math::longitude_difference<radian>(lon2, lon_2),
- math::longitude_difference<radian>(lon_b2, lon_2)));
+ math::longitude_difference<radian>(lon_b2, lon_2)));
if (lon_dist2 < lon_dist1)
{
lon = lon_2;
}
-
+
CT const sin_lon = sin(lon);
CT const cos_lon = cos(lon);
@@ -177,7 +178,7 @@ private:
CT const lat_x_2 = tan_alpha2 * cos_lat2;
tan_lat = lat_y_2 / lat_x_2;
}
-
+
return true;
}
};
@@ -206,19 +207,19 @@ inline CT sjoberg_d_lambda_e_sqr(CT const& sin_betaj, CT const& sin_beta,
{
using math::detail::bounded;
- if (Order == 0)
+ if (BOOST_GEOMETRY_CONDITION(Order == 0))
{
return 0;
}
CT const c1 = 1;
CT const c2 = 2;
-
+
CT const asin_B = asin(bounded(sin_beta / sqrt_1_Cj_sqr, -c1, c1));
CT const asin_Bj = asin(sin_betaj / sqrt_1_Cj_sqr);
CT const L0 = (asin_B - asin_Bj) / c2;
- if (Order == 1)
+ if (BOOST_GEOMETRY_CONDITION(Order == 1))
{
return -Cj * e_sqr * L0;
}
@@ -237,7 +238,7 @@ inline CT sjoberg_d_lambda_e_sqr(CT const& sin_betaj, CT const& sin_beta,
CT const sqrt_Yj = math::sqrt(-Xj_sqr + one_minus_Cj_sqr);
CT const L1 = (Cj_sqr_plus_one * (asin_B - asin_Bj) + X * sqrt_Y - Xj * sqrt_Yj) / c16;
- if (Order == 2)
+ if (BOOST_GEOMETRY_CONDITION(Order == 2))
{
return -Cj * e_sqr * (L0 + e_sqr * L1);
}
@@ -251,7 +252,7 @@ inline CT sjoberg_d_lambda_e_sqr(CT const& sin_betaj, CT const& sin_beta,
CT const Fj = Xj * (-c2 * Xj_sqr + c3 * Cj_sqr + c5);
CT const L2 = (E * (asin_B - asin_Bj) + F * sqrt_Y - Fj * sqrt_Yj) / c128;
- if (Order == 3)
+ if (BOOST_GEOMETRY_CONDITION(Order == 3))
{
return -Cj * e_sqr * (L0 + e_sqr * (L1 + e_sqr * L2));
}
@@ -405,7 +406,7 @@ public:
{
return false;
}
-
+
// NOTE: beta may be slightly out of bounds here but d_lambda handles that
CT const dLj = d_lambda(sin_beta);
delta_k = sign_lon_diff * (/*asin_t_t0j*/ - asin_tj_t0j + dLj);
@@ -451,7 +452,7 @@ public:
delta_k_behind = sign_lon_diff * pi;
return true;
}
-
+
// beta out of bounds and not close
if (check_sin_beta
&& ! (is_sin_beta_ok(sin_beta)
@@ -459,7 +460,7 @@ public:
{
return false;
}
-
+
// NOTE: beta may be slightly out of bounds here but d_lambda handles that
CT const dLj = d_lambda(sin_beta);
delta_k_before = sign_lon_diff * (/*asin_t_t0j*/ - asin_tj_t0j + dLj);
@@ -768,7 +769,7 @@ public:
}
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
@@ -802,7 +803,7 @@ private:
CT const c1 = 1;
CT const e_sqr = geod1.e_sqr;
-
+
CT lon1_diff = 0;
CT lon2_diff = 0;
@@ -892,7 +893,7 @@ private:
if (math::equals(dbeta_denom, c0))
{
return false;
- }
+ }
// The sign of dbeta is changed WRT [Sjoberg02]
CT const dbeta = (lon1_minus_lon2 + lon1_diff - lon2_diff) / dbeta_denom;
@@ -1027,7 +1028,7 @@ private:
{
break;
}
-
+
bool try_t2 = false;
converge_07_result result_curr;
if (converge_07_step_one(CT(sin(beta1)), t1, lon1_minus_lon2, geodesics, lon_sph, result_curr))
@@ -1063,7 +1064,7 @@ private:
}
}
-
+
if (try_t2)
{
CT t2 = t;
@@ -1156,7 +1157,7 @@ private:
{
using math::detail::bounded;
CT const c1 = 1;
-
+
CT k_diff_before = 0;
CT k_diff_behind = 0;
diff --git a/boost/geometry/formulas/spherical.hpp b/boost/geometry/formulas/spherical.hpp
index c73b8b1555..3a388bfbce 100644
--- a/boost/geometry/formulas/spherical.hpp
+++ b/boost/geometry/formulas/spherical.hpp
@@ -29,7 +29,7 @@
#include <boost/geometry/formulas/result_direct.hpp>
namespace boost { namespace geometry {
-
+
namespace formula {
template <typename T>
diff --git a/boost/geometry/formulas/thomas_inverse.hpp b/boost/geometry/formulas/thomas_inverse.hpp
index cf69c9df1d..c93aaf838c 100644
--- a/boost/geometry/formulas/thomas_inverse.hpp
+++ b/boost/geometry/formulas/thomas_inverse.hpp
@@ -136,7 +136,7 @@ public:
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);
@@ -151,7 +151,7 @@ public:
//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
@@ -162,7 +162,7 @@ public:
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.
diff --git a/boost/geometry/formulas/vertex_longitude.hpp b/boost/geometry/formulas/vertex_longitude.hpp
index e8c5e8775b..bef8f5a3f6 100644
--- a/boost/geometry/formulas/vertex_longitude.hpp
+++ b/boost/geometry/formulas/vertex_longitude.hpp
@@ -212,7 +212,7 @@ public:
CT const sign = bet3 >= c0
? c1
: cminus1;
-
+
CT const dlon_max = omg13 - sign * f * sin_alp0 * I3;
return dlon_max;
diff --git a/boost/geometry/formulas/vincenty_direct.hpp b/boost/geometry/formulas/vincenty_direct.hpp
index bc150e6d43..274af3f2f2 100644
--- a/boost/geometry/formulas/vincenty_direct.hpp
+++ b/boost/geometry/formulas/vincenty_direct.hpp
@@ -141,7 +141,7 @@ public:
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)
diff --git a/boost/geometry/formulas/vincenty_inverse.hpp b/boost/geometry/formulas/vincenty_inverse.hpp
index a25e5a0843..24f285d9bc 100644
--- a/boost/geometry/formulas/vincenty_inverse.hpp
+++ b/boost/geometry/formulas/vincenty_inverse.hpp
@@ -159,7 +159,7 @@ public:
} 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) )
{
// Some types cannot divide by doubles
@@ -187,7 +187,7 @@ public:
result.distance = radius_b * A * (sigma - delta_sigma); // (19)
}
-
+
if ( BOOST_GEOMETRY_CONDITION(CalcAzimuths) )
{
if (BOOST_GEOMETRY_CONDITION(CalcFwdAzimuth))
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-24 20:50:17 +0000