stabilize build system: depends, installer, boost/bdb fixes, cross targets groundwork
This commit is contained in:
@@ -0,0 +1,260 @@
|
||||
// Boost.Geometry
|
||||
|
||||
// Copyright (c) 2015-2017 Oracle and/or its affiliates.
|
||||
|
||||
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
|
||||
|
||||
// Use, modification and distribution is subject to the Boost Software License,
|
||||
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
|
||||
// http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#ifndef BOOST_GEOMETRY_FORMULAS_ANDOYER_INVERSE_HPP
|
||||
#define BOOST_GEOMETRY_FORMULAS_ANDOYER_INVERSE_HPP
|
||||
|
||||
|
||||
#include <boost/math/constants/constants.hpp>
|
||||
|
||||
#include <boost/geometry/core/radius.hpp>
|
||||
#include <boost/geometry/core/srs.hpp>
|
||||
|
||||
#include <boost/geometry/util/condition.hpp>
|
||||
#include <boost/geometry/util/math.hpp>
|
||||
|
||||
#include <boost/geometry/formulas/differential_quantities.hpp>
|
||||
#include <boost/geometry/formulas/flattening.hpp>
|
||||
#include <boost/geometry/formulas/result_inverse.hpp>
|
||||
|
||||
|
||||
namespace boost { namespace geometry { namespace formula
|
||||
{
|
||||
|
||||
/*!
|
||||
\brief The solution of the inverse problem of geodesics on latlong coordinates,
|
||||
Forsyth-Andoyer-Lambert type approximation with first order terms.
|
||||
\author See
|
||||
- Technical Report: PAUL D. THOMAS, MATHEMATICAL MODELS FOR NAVIGATION SYSTEMS, 1965
|
||||
http://www.dtic.mil/docs/citations/AD0627893
|
||||
- Technical Report: PAUL D. THOMAS, SPHEROIDAL GEODESICS, REFERENCE SYSTEMS, AND LOCAL GEOMETRY, 1970
|
||||
http://www.dtic.mil/docs/citations/AD703541
|
||||
*/
|
||||
template <
|
||||
typename CT,
|
||||
bool EnableDistance,
|
||||
bool EnableAzimuth,
|
||||
bool EnableReverseAzimuth = false,
|
||||
bool EnableReducedLength = false,
|
||||
bool EnableGeodesicScale = false
|
||||
>
|
||||
class andoyer_inverse
|
||||
{
|
||||
static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale;
|
||||
static const bool CalcAzimuths = EnableAzimuth || EnableReverseAzimuth || CalcQuantities;
|
||||
static const bool CalcFwdAzimuth = EnableAzimuth || CalcQuantities;
|
||||
static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcQuantities;
|
||||
|
||||
public:
|
||||
typedef result_inverse<CT> result_type;
|
||||
|
||||
template <typename T1, typename T2, typename Spheroid>
|
||||
static inline result_type apply(T1 const& lon1,
|
||||
T1 const& lat1,
|
||||
T2 const& lon2,
|
||||
T2 const& lat2,
|
||||
Spheroid const& spheroid)
|
||||
{
|
||||
result_type result;
|
||||
|
||||
// coordinates in radians
|
||||
|
||||
if ( math::equals(lon1, lon2) && math::equals(lat1, lat2) )
|
||||
{
|
||||
return result;
|
||||
}
|
||||
|
||||
CT const c0 = CT(0);
|
||||
CT const c1 = CT(1);
|
||||
CT const pi = math::pi<CT>();
|
||||
CT const f = formula::flattening<CT>(spheroid);
|
||||
|
||||
CT const dlon = lon2 - lon1;
|
||||
CT const sin_dlon = sin(dlon);
|
||||
CT const cos_dlon = cos(dlon);
|
||||
CT const sin_lat1 = sin(lat1);
|
||||
CT const cos_lat1 = cos(lat1);
|
||||
CT const sin_lat2 = sin(lat2);
|
||||
CT const cos_lat2 = cos(lat2);
|
||||
|
||||
// H,G,T = infinity if cos_d = 1 or cos_d = -1
|
||||
// lat1 == +-90 && lat2 == +-90
|
||||
// lat1 == lat2 && lon1 == lon2
|
||||
CT cos_d = sin_lat1*sin_lat2 + cos_lat1*cos_lat2*cos_dlon;
|
||||
// on some platforms cos_d may be outside valid range
|
||||
if (cos_d < -c1)
|
||||
cos_d = -c1;
|
||||
else if (cos_d > c1)
|
||||
cos_d = c1;
|
||||
|
||||
CT const d = acos(cos_d); // [0, pi]
|
||||
CT const sin_d = sin(d); // [-1, 1]
|
||||
|
||||
if ( BOOST_GEOMETRY_CONDITION(EnableDistance) )
|
||||
{
|
||||
CT const K = math::sqr(sin_lat1-sin_lat2);
|
||||
CT const L = math::sqr(sin_lat1+sin_lat2);
|
||||
CT const three_sin_d = CT(3) * sin_d;
|
||||
|
||||
CT const one_minus_cos_d = c1 - cos_d;
|
||||
CT const one_plus_cos_d = c1 + cos_d;
|
||||
// cos_d = 1 or cos_d = -1 means that the points are antipodal
|
||||
|
||||
CT const H = math::equals(one_minus_cos_d, c0) ?
|
||||
c0 :
|
||||
(d + three_sin_d) / one_minus_cos_d;
|
||||
CT const G = math::equals(one_plus_cos_d, c0) ?
|
||||
c0 :
|
||||
(d - three_sin_d) / one_plus_cos_d;
|
||||
|
||||
CT const dd = -(f/CT(4))*(H*K+G*L);
|
||||
|
||||
CT const a = get_radius<0>(spheroid);
|
||||
|
||||
result.distance = a * (d + dd);
|
||||
}
|
||||
|
||||
if ( BOOST_GEOMETRY_CONDITION(CalcAzimuths) )
|
||||
{
|
||||
// sin_d = 0 <=> antipodal points
|
||||
if (math::equals(sin_d, c0))
|
||||
{
|
||||
// T = inf
|
||||
// dA = inf
|
||||
// azimuth = -inf
|
||||
result.azimuth = lat1 <= lat2 ? c0 : pi;
|
||||
}
|
||||
else
|
||||
{
|
||||
CT const c2 = CT(2);
|
||||
|
||||
CT A = c0;
|
||||
CT U = c0;
|
||||
if (math::equals(cos_lat2, c0))
|
||||
{
|
||||
if (sin_lat2 < c0)
|
||||
{
|
||||
A = pi;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
CT const tan_lat2 = sin_lat2/cos_lat2;
|
||||
CT const M = cos_lat1*tan_lat2-sin_lat1*cos_dlon;
|
||||
A = atan2(sin_dlon, M);
|
||||
CT const sin_2A = sin(c2*A);
|
||||
U = (f/ c2)*math::sqr(cos_lat1)*sin_2A;
|
||||
}
|
||||
|
||||
CT B = c0;
|
||||
CT V = c0;
|
||||
if (math::equals(cos_lat1, c0))
|
||||
{
|
||||
if (sin_lat1 < c0)
|
||||
{
|
||||
B = pi;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
CT const tan_lat1 = sin_lat1/cos_lat1;
|
||||
CT const N = cos_lat2*tan_lat1-sin_lat2*cos_dlon;
|
||||
B = atan2(sin_dlon, N);
|
||||
CT const sin_2B = sin(c2*B);
|
||||
V = (f/ c2)*math::sqr(cos_lat2)*sin_2B;
|
||||
}
|
||||
|
||||
CT const T = d / sin_d;
|
||||
|
||||
// even with sin_d == 0 checked above if the second point
|
||||
// is somewhere in the antipodal area T may still be great
|
||||
// therefore dA and dB may be great and the resulting azimuths
|
||||
// may be some more or less arbitrary angles
|
||||
|
||||
if (BOOST_GEOMETRY_CONDITION(CalcFwdAzimuth))
|
||||
{
|
||||
CT const dA = V*T - U;
|
||||
result.azimuth = A - dA;
|
||||
normalize_azimuth(result.azimuth, A, dA);
|
||||
}
|
||||
|
||||
if (BOOST_GEOMETRY_CONDITION(CalcRevAzimuth))
|
||||
{
|
||||
CT const dB = -U*T + V;
|
||||
result.reverse_azimuth = pi - B - dB;
|
||||
if (result.reverse_azimuth > pi)
|
||||
{
|
||||
result.reverse_azimuth -= 2 * pi;
|
||||
}
|
||||
normalize_azimuth(result.reverse_azimuth, B, dB);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (BOOST_GEOMETRY_CONDITION(CalcQuantities))
|
||||
{
|
||||
typedef differential_quantities<CT, EnableReducedLength, EnableGeodesicScale, 1> quantities;
|
||||
quantities::apply(dlon, sin_lat1, cos_lat1, sin_lat2, cos_lat2,
|
||||
result.azimuth, result.reverse_azimuth,
|
||||
get_radius<2>(spheroid), f,
|
||||
result.reduced_length, result.geodesic_scale);
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
private:
|
||||
static inline void normalize_azimuth(CT & azimuth, CT const& A, CT const& dA)
|
||||
{
|
||||
CT const c0 = 0;
|
||||
|
||||
if (A >= c0) // A indicates Eastern hemisphere
|
||||
{
|
||||
if (dA >= c0) // A altered towards 0
|
||||
{
|
||||
if (azimuth < c0)
|
||||
{
|
||||
azimuth = c0;
|
||||
}
|
||||
}
|
||||
else // dA < 0, A altered towards pi
|
||||
{
|
||||
CT const pi = math::pi<CT>();
|
||||
if (azimuth > pi)
|
||||
{
|
||||
azimuth = pi;
|
||||
}
|
||||
}
|
||||
}
|
||||
else // A indicates Western hemisphere
|
||||
{
|
||||
if (dA <= c0) // A altered towards 0
|
||||
{
|
||||
if (azimuth > c0)
|
||||
{
|
||||
azimuth = c0;
|
||||
}
|
||||
}
|
||||
else // dA > 0, A altered towards -pi
|
||||
{
|
||||
CT const minus_pi = -math::pi<CT>();
|
||||
if (azimuth < minus_pi)
|
||||
{
|
||||
azimuth = minus_pi;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
}}} // namespace boost::geometry::formula
|
||||
|
||||
|
||||
#endif // BOOST_GEOMETRY_FORMULAS_ANDOYER_INVERSE_HPP
|
||||
@@ -0,0 +1,578 @@
|
||||
// Boost.Geometry
|
||||
|
||||
// Copyright (c) 2015-2016 Oracle and/or its affiliates.
|
||||
|
||||
// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
|
||||
|
||||
// Use, modification and distribution is subject to the Boost Software License,
|
||||
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
|
||||
// http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#ifndef BOOST_GEOMETRY_FORMULAS_AREA_FORMULAS_HPP
|
||||
#define BOOST_GEOMETRY_FORMULAS_AREA_FORMULAS_HPP
|
||||
|
||||
#include <boost/geometry/formulas/flattening.hpp>
|
||||
#include <boost/math/special_functions/hypot.hpp>
|
||||
|
||||
namespace boost { namespace geometry { namespace formula
|
||||
{
|
||||
|
||||
/*!
|
||||
\brief Formulas for computing spherical and ellipsoidal polygon area.
|
||||
The current class computes the area of the trapezoid defined by a segment
|
||||
the two meridians passing by the endpoints and the equator.
|
||||
\author See
|
||||
- Danielsen JS, The area under the geodesic. Surv Rev 30(232):
|
||||
61–66, 1989
|
||||
- Charles F.F Karney, Algorithms for geodesics, 2011
|
||||
https://arxiv.org/pdf/1109.4448.pdf
|
||||
*/
|
||||
|
||||
template <
|
||||
typename CT,
|
||||
std::size_t SeriesOrder = 2,
|
||||
bool ExpandEpsN = true
|
||||
>
|
||||
class area_formulas
|
||||
{
|
||||
|
||||
public:
|
||||
|
||||
//TODO: move the following to a more general space to be used by other
|
||||
// classes as well
|
||||
/*
|
||||
Evaluate the polynomial in x using Horner's method.
|
||||
*/
|
||||
template <typename NT, typename IteratorType>
|
||||
static inline NT horner_evaluate(NT x,
|
||||
IteratorType begin,
|
||||
IteratorType end)
|
||||
{
|
||||
NT result(0);
|
||||
IteratorType it = end;
|
||||
do
|
||||
{
|
||||
result = result * x + *--it;
|
||||
}
|
||||
while (it != begin);
|
||||
return result;
|
||||
}
|
||||
|
||||
/*
|
||||
Clenshaw algorithm for summing trigonometric series
|
||||
https://en.wikipedia.org/wiki/Clenshaw_algorithm
|
||||
*/
|
||||
template <typename NT, typename IteratorType>
|
||||
static inline NT clenshaw_sum(NT cosx,
|
||||
IteratorType begin,
|
||||
IteratorType end)
|
||||
{
|
||||
IteratorType it = end;
|
||||
bool odd = true;
|
||||
CT b_k, b_k1(0), b_k2(0);
|
||||
do
|
||||
{
|
||||
CT c_k = odd ? *--it : NT(0);
|
||||
b_k = c_k + NT(2) * cosx * b_k1 - b_k2;
|
||||
b_k2 = b_k1;
|
||||
b_k1 = b_k;
|
||||
odd = !odd;
|
||||
}
|
||||
while (it != begin);
|
||||
|
||||
return *begin + b_k1 * cosx - b_k2;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
static inline void normalize(T& x, T& y)
|
||||
{
|
||||
T h = boost::math::hypot(x, y);
|
||||
x /= h;
|
||||
y /= h;
|
||||
}
|
||||
|
||||
/*
|
||||
Generate and evaluate the series expansion of the following integral
|
||||
|
||||
I4 = -integrate( (t(ep2) - t(k2*sin(sigma1)^2)) / (ep2 - k2*sin(sigma1)^2)
|
||||
* sin(sigma1)/2, sigma1, pi/2, sigma )
|
||||
where
|
||||
|
||||
t(x) = sqrt(1+1/x)*asinh(sqrt(x)) + x
|
||||
|
||||
valid for ep2 and k2 small. We substitute k2 = 4 * eps / (1 - eps)^2
|
||||
and ep2 = 4 * n / (1 - n)^2 and expand in eps and n.
|
||||
|
||||
The resulting sum of the series is of the form
|
||||
|
||||
sum(C4[l] * cos((2*l+1)*sigma), l, 0, maxpow-1) )
|
||||
|
||||
The above expansion is performed in Computer Algebra System Maxima.
|
||||
The C++ code (that yields the function evaluate_coeffs_n below) is generated
|
||||
by the following Maxima script and is based on script:
|
||||
http://geographiclib.sourceforge.net/html/geod.mac
|
||||
|
||||
// Maxima script begin
|
||||
taylordepth:5$
|
||||
ataylor(expr,var,ord):=expand(ratdisrep(taylor(expr,var,0,ord)))$
|
||||
jtaylor(expr,var1,var2,ord):=block([zz],expand(subst([zz=1],
|
||||
ratdisrep(taylor(subst([var1=zz*var1,var2=zz*var2],expr),zz,0,ord)))))$
|
||||
|
||||
compute(maxpow):=block([int,t,intexp,area, x,ep2,k2],
|
||||
maxpow:maxpow-1,
|
||||
t : sqrt(1+1/x) * asinh(sqrt(x)) + x,
|
||||
int:-(tf(ep2) - tf(k2*sin(sigma)^2)) / (ep2 - k2*sin(sigma)^2)
|
||||
* sin(sigma)/2,
|
||||
int:subst([tf(ep2)=subst([x=ep2],t),
|
||||
tf(k2*sin(sigma)^2)=subst([x=k2*sin(sigma)^2],t)],
|
||||
int),
|
||||
int:subst([abs(sin(sigma))=sin(sigma)],int),
|
||||
int:subst([k2=4*eps/(1-eps)^2,ep2=4*n/(1-n)^2],int),
|
||||
intexp:jtaylor(int,n,eps,maxpow),
|
||||
area:trigreduce(integrate(intexp,sigma)),
|
||||
area:expand(area-subst(sigma=%pi/2,area)),
|
||||
for i:0 thru maxpow do C4[i]:coeff(area,cos((2*i+1)*sigma)),
|
||||
if expand(area-sum(C4[i]*cos((2*i+1)*sigma),i,0,maxpow)) # 0
|
||||
then error("left over terms in I4"),
|
||||
'done)$
|
||||
|
||||
printcode(maxpow):=
|
||||
block([tab2:" ",tab3:" "],
|
||||
print(" switch (SeriesOrder) {"),
|
||||
for nn:1 thru maxpow do block([c],
|
||||
print(concat(tab2,"case ",string(nn-1),":")),
|
||||
c:0,
|
||||
for m:0 thru nn-1 do block(
|
||||
[q:jtaylor(subst([n=n],C4[m]),n,eps,nn-1),
|
||||
linel:1200],
|
||||
for j:m thru nn-1 do (
|
||||
print(concat(tab3,"coeffs_n[",c,"] = ",
|
||||
string(horner(coeff(q,eps,j))),";")),
|
||||
c:c+1)
|
||||
),
|
||||
print(concat(tab3,"break;"))),
|
||||
print(" }"),
|
||||
'done)$
|
||||
|
||||
maxpow:6$
|
||||
compute(maxpow)$
|
||||
printcode(maxpow)$
|
||||
// Maxima script end
|
||||
|
||||
In the resulting code we should replace each number x by CT(x)
|
||||
e.g. using the following scirpt:
|
||||
sed -e 's/[0-9]\+/CT(&)/g; s/\[CT(/\[/g; s/)\]/\]/g;
|
||||
s/case\sCT(/case /g; s/):/:/g'
|
||||
*/
|
||||
|
||||
static inline void evaluate_coeffs_n(CT n, CT coeffs_n[])
|
||||
{
|
||||
|
||||
switch (SeriesOrder) {
|
||||
case 0:
|
||||
coeffs_n[0] = CT(2)/CT(3);
|
||||
break;
|
||||
case 1:
|
||||
coeffs_n[0] = (CT(10)-CT(4)*n)/CT(15);
|
||||
coeffs_n[1] = -CT(1)/CT(5);
|
||||
coeffs_n[2] = CT(1)/CT(45);
|
||||
break;
|
||||
case 2:
|
||||
coeffs_n[0] = (n*(CT(8)*n-CT(28))+CT(70))/CT(105);
|
||||
coeffs_n[1] = (CT(16)*n-CT(7))/CT(35);
|
||||
coeffs_n[2] = -CT(2)/CT(105);
|
||||
coeffs_n[3] = (CT(7)-CT(16)*n)/CT(315);
|
||||
coeffs_n[4] = -CT(2)/CT(105);
|
||||
coeffs_n[5] = CT(4)/CT(525);
|
||||
break;
|
||||
case 3:
|
||||
coeffs_n[0] = (n*(n*(CT(4)*n+CT(24))-CT(84))+CT(210))/CT(315);
|
||||
coeffs_n[1] = ((CT(48)-CT(32)*n)*n-CT(21))/CT(105);
|
||||
coeffs_n[2] = (-CT(32)*n-CT(6))/CT(315);
|
||||
coeffs_n[3] = CT(11)/CT(315);
|
||||
coeffs_n[4] = (n*(CT(32)*n-CT(48))+CT(21))/CT(945);
|
||||
coeffs_n[5] = (CT(64)*n-CT(18))/CT(945);
|
||||
coeffs_n[6] = -CT(1)/CT(105);
|
||||
coeffs_n[7] = (CT(12)-CT(32)*n)/CT(1575);
|
||||
coeffs_n[8] = -CT(8)/CT(1575);
|
||||
coeffs_n[9] = CT(8)/CT(2205);
|
||||
break;
|
||||
case 4:
|
||||
coeffs_n[0] = (n*(n*(n*(CT(16)*n+CT(44))+CT(264))-CT(924))+CT(2310))/CT(3465);
|
||||
coeffs_n[1] = (n*(n*(CT(48)*n-CT(352))+CT(528))-CT(231))/CT(1155);
|
||||
coeffs_n[2] = (n*(CT(1088)*n-CT(352))-CT(66))/CT(3465);
|
||||
coeffs_n[3] = (CT(121)-CT(368)*n)/CT(3465);
|
||||
coeffs_n[4] = CT(4)/CT(1155);
|
||||
coeffs_n[5] = (n*((CT(352)-CT(48)*n)*n-CT(528))+CT(231))/CT(10395);
|
||||
coeffs_n[6] = ((CT(704)-CT(896)*n)*n-CT(198))/CT(10395);
|
||||
coeffs_n[7] = (CT(80)*n-CT(99))/CT(10395);
|
||||
coeffs_n[8] = CT(4)/CT(1155);
|
||||
coeffs_n[9] = (n*(CT(320)*n-CT(352))+CT(132))/CT(17325);
|
||||
coeffs_n[10] = (CT(384)*n-CT(88))/CT(17325);
|
||||
coeffs_n[11] = -CT(8)/CT(1925);
|
||||
coeffs_n[12] = (CT(88)-CT(256)*n)/CT(24255);
|
||||
coeffs_n[13] = -CT(16)/CT(8085);
|
||||
coeffs_n[14] = CT(64)/CT(31185);
|
||||
break;
|
||||
case 5:
|
||||
coeffs_n[0] = (n*(n*(n*(n*(CT(100)*n+CT(208))+CT(572))+CT(3432))-CT(12012))+CT(30030))
|
||||
/CT(45045);
|
||||
coeffs_n[1] = (n*(n*(n*(CT(64)*n+CT(624))-CT(4576))+CT(6864))-CT(3003))/CT(15015);
|
||||
coeffs_n[2] = (n*((CT(14144)-CT(10656)*n)*n-CT(4576))-CT(858))/CT(45045);
|
||||
coeffs_n[3] = ((-CT(224)*n-CT(4784))*n+CT(1573))/CT(45045);
|
||||
coeffs_n[4] = (CT(1088)*n+CT(156))/CT(45045);
|
||||
coeffs_n[5] = CT(97)/CT(15015);
|
||||
coeffs_n[6] = (n*(n*((-CT(64)*n-CT(624))*n+CT(4576))-CT(6864))+CT(3003))/CT(135135);
|
||||
coeffs_n[7] = (n*(n*(CT(5952)*n-CT(11648))+CT(9152))-CT(2574))/CT(135135);
|
||||
coeffs_n[8] = (n*(CT(5792)*n+CT(1040))-CT(1287))/CT(135135);
|
||||
coeffs_n[9] = (CT(468)-CT(2944)*n)/CT(135135);
|
||||
coeffs_n[10] = CT(1)/CT(9009);
|
||||
coeffs_n[11] = (n*((CT(4160)-CT(1440)*n)*n-CT(4576))+CT(1716))/CT(225225);
|
||||
coeffs_n[12] = ((CT(4992)-CT(8448)*n)*n-CT(1144))/CT(225225);
|
||||
coeffs_n[13] = (CT(1856)*n-CT(936))/CT(225225);
|
||||
coeffs_n[14] = CT(8)/CT(10725);
|
||||
coeffs_n[15] = (n*(CT(3584)*n-CT(3328))+CT(1144))/CT(315315);
|
||||
coeffs_n[16] = (CT(1024)*n-CT(208))/CT(105105);
|
||||
coeffs_n[17] = -CT(136)/CT(63063);
|
||||
coeffs_n[18] = (CT(832)-CT(2560)*n)/CT(405405);
|
||||
coeffs_n[19] = -CT(128)/CT(135135);
|
||||
coeffs_n[20] = CT(128)/CT(99099);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
/*
|
||||
Expand in k2 and ep2.
|
||||
*/
|
||||
static inline void evaluate_coeffs_ep(CT ep, CT coeffs_n[])
|
||||
{
|
||||
switch (SeriesOrder) {
|
||||
case 0:
|
||||
coeffs_n[0] = CT(2)/CT(3);
|
||||
break;
|
||||
case 1:
|
||||
coeffs_n[0] = (CT(10)-ep)/CT(15);
|
||||
coeffs_n[1] = -CT(1)/CT(20);
|
||||
coeffs_n[2] = CT(1)/CT(180);
|
||||
break;
|
||||
case 2:
|
||||
coeffs_n[0] = (ep*(CT(4)*ep-CT(7))+CT(70))/CT(105);
|
||||
coeffs_n[1] = (CT(4)*ep-CT(7))/CT(140);
|
||||
coeffs_n[2] = CT(1)/CT(42);
|
||||
coeffs_n[3] = (CT(7)-CT(4)*ep)/CT(1260);
|
||||
coeffs_n[4] = -CT(1)/CT(252);
|
||||
coeffs_n[5] = CT(1)/CT(2100);
|
||||
break;
|
||||
case 3:
|
||||
coeffs_n[0] = (ep*((CT(12)-CT(8)*ep)*ep-CT(21))+CT(210))/CT(315);
|
||||
coeffs_n[1] = ((CT(12)-CT(8)*ep)*ep-CT(21))/CT(420);
|
||||
coeffs_n[2] = (CT(3)-CT(2)*ep)/CT(126);
|
||||
coeffs_n[3] = -CT(1)/CT(72);
|
||||
coeffs_n[4] = (ep*(CT(8)*ep-CT(12))+CT(21))/CT(3780);
|
||||
coeffs_n[5] = (CT(2)*ep-CT(3))/CT(756);
|
||||
coeffs_n[6] = CT(1)/CT(360);
|
||||
coeffs_n[7] = (CT(3)-CT(2)*ep)/CT(6300);
|
||||
coeffs_n[8] = -CT(1)/CT(1800);
|
||||
coeffs_n[9] = CT(1)/CT(17640);
|
||||
break;
|
||||
case 4:
|
||||
coeffs_n[0] = (ep*(ep*(ep*(CT(64)*ep-CT(88))+CT(132))-CT(231))+CT(2310))/CT(3465);
|
||||
coeffs_n[1] = (ep*(ep*(CT(64)*ep-CT(88))+CT(132))-CT(231))/CT(4620);
|
||||
coeffs_n[2] = (ep*(CT(16)*ep-CT(22))+CT(33))/CT(1386);
|
||||
coeffs_n[3] = (CT(8)*ep-CT(11))/CT(792);
|
||||
coeffs_n[4] = CT(1)/CT(110);
|
||||
coeffs_n[5] = (ep*((CT(88)-CT(64)*ep)*ep-CT(132))+CT(231))/CT(41580);
|
||||
coeffs_n[6] = ((CT(22)-CT(16)*ep)*ep-CT(33))/CT(8316);
|
||||
coeffs_n[7] = (CT(11)-CT(8)*ep)/CT(3960);
|
||||
coeffs_n[8] = -CT(1)/CT(495);
|
||||
coeffs_n[9] = (ep*(CT(16)*ep-CT(22))+CT(33))/CT(69300);
|
||||
coeffs_n[10] = (CT(8)*ep-CT(11))/CT(19800);
|
||||
coeffs_n[11] = CT(1)/CT(1925);
|
||||
coeffs_n[12] = (CT(11)-CT(8)*ep)/CT(194040);
|
||||
coeffs_n[13] = -CT(1)/CT(10780);
|
||||
coeffs_n[14] = CT(1)/CT(124740);
|
||||
break;
|
||||
case 5:
|
||||
coeffs_n[0] = (ep*(ep*(ep*((CT(832)-CT(640)*ep)*ep-CT(1144))+CT(1716))-CT(3003))+CT(30030))/CT(45045);
|
||||
coeffs_n[1] = (ep*(ep*((CT(832)-CT(640)*ep)*ep-CT(1144))+CT(1716))-CT(3003))/CT(60060);
|
||||
coeffs_n[2] = (ep*((CT(208)-CT(160)*ep)*ep-CT(286))+CT(429))/CT(18018);
|
||||
coeffs_n[3] = ((CT(104)-CT(80)*ep)*ep-CT(143))/CT(10296);
|
||||
coeffs_n[4] = (CT(13)-CT(10)*ep)/CT(1430);
|
||||
coeffs_n[5] = -CT(1)/CT(156);
|
||||
coeffs_n[6] = (ep*(ep*(ep*(CT(640)*ep-CT(832))+CT(1144))-CT(1716))+CT(3003))/CT(540540);
|
||||
coeffs_n[7] = (ep*(ep*(CT(160)*ep-CT(208))+CT(286))-CT(429))/CT(108108);
|
||||
coeffs_n[8] = (ep*(CT(80)*ep-CT(104))+CT(143))/CT(51480);
|
||||
coeffs_n[9] = (CT(10)*ep-CT(13))/CT(6435);
|
||||
coeffs_n[10] = CT(5)/CT(3276);
|
||||
coeffs_n[11] = (ep*((CT(208)-CT(160)*ep)*ep-CT(286))+CT(429))/CT(900900);
|
||||
coeffs_n[12] = ((CT(104)-CT(80)*ep)*ep-CT(143))/CT(257400);
|
||||
coeffs_n[13] = (CT(13)-CT(10)*ep)/CT(25025);
|
||||
coeffs_n[14] = -CT(1)/CT(2184);
|
||||
coeffs_n[15] = (ep*(CT(80)*ep-CT(104))+CT(143))/CT(2522520);
|
||||
coeffs_n[16] = (CT(10)*ep-CT(13))/CT(140140);
|
||||
coeffs_n[17] = CT(5)/CT(45864);
|
||||
coeffs_n[18] = (CT(13)-CT(10)*ep)/CT(1621620);
|
||||
coeffs_n[19] = -CT(1)/CT(58968);
|
||||
coeffs_n[20] = CT(1)/CT(792792);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
/*
|
||||
Given the set of coefficients coeffs1[] evaluate on var2 and return
|
||||
the set of coefficients coeffs2[]
|
||||
*/
|
||||
static inline void evaluate_coeffs_var2(CT var2,
|
||||
CT coeffs1[],
|
||||
CT coeffs2[]){
|
||||
std::size_t begin(0), end(0);
|
||||
for(std::size_t i = 0; i <= SeriesOrder; i++){
|
||||
end = begin + SeriesOrder + 1 - i;
|
||||
coeffs2[i] = ((i==0) ? CT(1) : pow(var2,int(i)))
|
||||
* horner_evaluate(var2, coeffs1 + begin, coeffs1 + end);
|
||||
begin = end;
|
||||
}
|
||||
}
|
||||
|
||||
/*
|
||||
Compute the spherical excess of a geodesic (or shperical) segment
|
||||
*/
|
||||
template <
|
||||
bool LongSegment,
|
||||
typename PointOfSegment
|
||||
>
|
||||
static inline CT spherical(PointOfSegment const& p1,
|
||||
PointOfSegment const& p2)
|
||||
{
|
||||
CT excess;
|
||||
|
||||
if(LongSegment) // not for segments parallel to equator
|
||||
{
|
||||
CT cbet1 = cos(geometry::get_as_radian<1>(p1));
|
||||
CT sbet1 = sin(geometry::get_as_radian<1>(p1));
|
||||
CT cbet2 = cos(geometry::get_as_radian<1>(p2));
|
||||
CT sbet2 = sin(geometry::get_as_radian<1>(p2));
|
||||
|
||||
CT omg12 = geometry::get_as_radian<0>(p1)
|
||||
- geometry::get_as_radian<0>(p2);
|
||||
CT comg12 = cos(omg12);
|
||||
CT somg12 = sin(omg12);
|
||||
|
||||
CT alp1 = atan2(cbet1 * sbet2
|
||||
- sbet1 * cbet2 * comg12,
|
||||
cbet2 * somg12);
|
||||
|
||||
CT alp2 = atan2(cbet1 * sbet2 * comg12
|
||||
- sbet1 * cbet2,
|
||||
cbet1 * somg12);
|
||||
|
||||
excess = alp2 - alp1;
|
||||
|
||||
} else {
|
||||
|
||||
// Trapezoidal formula
|
||||
|
||||
CT tan_lat1 =
|
||||
tan(geometry::get_as_radian<1>(p1) / 2.0);
|
||||
CT tan_lat2 =
|
||||
tan(geometry::get_as_radian<1>(p2) / 2.0);
|
||||
|
||||
excess = CT(2.0)
|
||||
* atan(((tan_lat1 + tan_lat2) / (CT(1) + tan_lat1 * tan_lat2))
|
||||
* tan((geometry::get_as_radian<0>(p2)
|
||||
- geometry::get_as_radian<0>(p1)) / 2));
|
||||
}
|
||||
|
||||
return excess;
|
||||
}
|
||||
|
||||
struct return_type_ellipsoidal
|
||||
{
|
||||
return_type_ellipsoidal()
|
||||
: spherical_term(0),
|
||||
ellipsoidal_term(0)
|
||||
{}
|
||||
|
||||
CT spherical_term;
|
||||
CT ellipsoidal_term;
|
||||
};
|
||||
|
||||
/*
|
||||
Compute the ellipsoidal correction of a geodesic (or shperical) segment
|
||||
*/
|
||||
template <
|
||||
template <typename, bool, bool, bool, bool, bool> class Inverse,
|
||||
//typename AzimuthStrategy,
|
||||
typename PointOfSegment,
|
||||
typename SpheroidConst
|
||||
>
|
||||
static inline return_type_ellipsoidal ellipsoidal(PointOfSegment const& p1,
|
||||
PointOfSegment const& p2,
|
||||
SpheroidConst spheroid_const)
|
||||
{
|
||||
return_type_ellipsoidal result;
|
||||
|
||||
// Azimuth Approximation
|
||||
|
||||
typedef Inverse<CT, false, true, true, false, false> inverse_type;
|
||||
typedef typename inverse_type::result_type inverse_result;
|
||||
|
||||
inverse_result i_res = inverse_type::apply(get_as_radian<0>(p1),
|
||||
get_as_radian<1>(p1),
|
||||
get_as_radian<0>(p2),
|
||||
get_as_radian<1>(p2),
|
||||
spheroid_const.m_spheroid);
|
||||
|
||||
CT alp1 = i_res.azimuth;
|
||||
CT alp2 = i_res.reverse_azimuth;
|
||||
|
||||
// Constants
|
||||
|
||||
CT const ep = spheroid_const.m_ep;
|
||||
CT const f = formula::flattening<CT>(spheroid_const.m_spheroid);
|
||||
CT const one_minus_f = CT(1) - f;
|
||||
std::size_t const series_order_plus_one = SeriesOrder + 1;
|
||||
std::size_t const series_order_plus_two = SeriesOrder + 2;
|
||||
|
||||
// Basic trigonometric computations
|
||||
|
||||
CT tan_bet1 = tan(get_as_radian<1>(p1)) * one_minus_f;
|
||||
CT tan_bet2 = tan(get_as_radian<1>(p2)) * one_minus_f;
|
||||
CT cos_bet1 = cos(atan(tan_bet1));
|
||||
CT cos_bet2 = cos(atan(tan_bet2));
|
||||
CT sin_bet1 = tan_bet1 * cos_bet1;
|
||||
CT sin_bet2 = tan_bet2 * cos_bet2;
|
||||
CT sin_alp1 = sin(alp1);
|
||||
CT cos_alp1 = cos(alp1);
|
||||
CT cos_alp2 = cos(alp2);
|
||||
CT sin_alp0 = sin_alp1 * cos_bet1;
|
||||
|
||||
// Spherical term computation
|
||||
|
||||
CT sin_omg1 = sin_alp0 * sin_bet1;
|
||||
CT cos_omg1 = cos_alp1 * cos_bet1;
|
||||
CT sin_omg2 = sin_alp0 * sin_bet2;
|
||||
CT cos_omg2 = cos_alp2 * cos_bet2;
|
||||
CT cos_omg12 = cos_omg1 * cos_omg2 + sin_omg1 * sin_omg2;
|
||||
CT excess;
|
||||
|
||||
bool meridian = get<0>(p2) - get<0>(p1) == CT(0)
|
||||
|| get<1>(p1) == CT(90) || get<1>(p1) == -CT(90)
|
||||
|| get<1>(p2) == CT(90) || get<1>(p2) == -CT(90);
|
||||
|
||||
if (!meridian && cos_omg12 > -CT(0.7)
|
||||
&& sin_bet2 - sin_bet1 < CT(1.75)) // short segment
|
||||
{
|
||||
CT sin_omg12 = cos_omg1 * sin_omg2 - sin_omg1 * cos_omg2;
|
||||
normalize(sin_omg12, cos_omg12);
|
||||
|
||||
CT cos_omg12p1 = CT(1) + cos_omg12;
|
||||
CT cos_bet1p1 = CT(1) + cos_bet1;
|
||||
CT cos_bet2p1 = CT(1) + cos_bet2;
|
||||
excess = CT(2) * atan2(sin_omg12 * (sin_bet1 * cos_bet2p1 + sin_bet2 * cos_bet1p1),
|
||||
cos_omg12p1 * (sin_bet1 * sin_bet2 + cos_bet1p1 * cos_bet2p1));
|
||||
}
|
||||
else
|
||||
{
|
||||
/*
|
||||
CT sin_alp2 = sin(alp2);
|
||||
CT sin_alp12 = sin_alp2 * cos_alp1 - cos_alp2 * sin_alp1;
|
||||
CT cos_alp12 = cos_alp2 * cos_alp1 + sin_alp2 * sin_alp1;
|
||||
excess = atan2(sin_alp12, cos_alp12);
|
||||
*/
|
||||
excess = alp2 - alp1;
|
||||
}
|
||||
|
||||
result.spherical_term = excess;
|
||||
|
||||
// Ellipsoidal term computation (uses integral approximation)
|
||||
|
||||
CT cos_alp0 = math::sqrt(CT(1) - math::sqr(sin_alp0));
|
||||
CT cos_sig1 = cos_alp1 * cos_bet1;
|
||||
CT cos_sig2 = cos_alp2 * cos_bet2;
|
||||
CT sin_sig1 = sin_bet1;
|
||||
CT sin_sig2 = sin_bet2;
|
||||
|
||||
normalize(sin_sig1, cos_sig1);
|
||||
normalize(sin_sig2, cos_sig2);
|
||||
|
||||
CT coeffs[SeriesOrder + 1];
|
||||
const std::size_t coeffs_var_size = (series_order_plus_two
|
||||
* series_order_plus_one) / 2;
|
||||
CT coeffs_var[coeffs_var_size];
|
||||
|
||||
if(ExpandEpsN){ // expand by eps and n
|
||||
|
||||
CT k2 = math::sqr(ep * cos_alp0);
|
||||
CT sqrt_k2_plus_one = math::sqrt(CT(1) + k2);
|
||||
CT eps = (sqrt_k2_plus_one - CT(1)) / (sqrt_k2_plus_one + CT(1));
|
||||
CT n = f / (CT(2) - f);
|
||||
|
||||
// Generate and evaluate the polynomials on n
|
||||
// to get the series coefficients (that depend on eps)
|
||||
evaluate_coeffs_n(n, coeffs_var);
|
||||
|
||||
// Generate and evaluate the polynomials on eps (i.e. var2 = eps)
|
||||
// to get the final series coefficients
|
||||
evaluate_coeffs_var2(eps, coeffs_var, coeffs);
|
||||
|
||||
}else{ // expand by k2 and ep
|
||||
|
||||
CT k2 = math::sqr(ep * cos_alp0);
|
||||
CT ep2 = math::sqr(ep);
|
||||
|
||||
// Generate and evaluate the polynomials on ep2
|
||||
evaluate_coeffs_ep(ep2, coeffs_var);
|
||||
|
||||
// Generate and evaluate the polynomials on k2 (i.e. var2 = k2)
|
||||
evaluate_coeffs_var2(k2, coeffs_var, coeffs);
|
||||
|
||||
}
|
||||
|
||||
// Evaluate the trigonometric sum
|
||||
CT I12 = clenshaw_sum(cos_sig2, coeffs, coeffs + series_order_plus_one)
|
||||
- clenshaw_sum(cos_sig1, coeffs, coeffs + series_order_plus_one);
|
||||
|
||||
// The part of the ellipsodal correction that depends on
|
||||
// point coordinates
|
||||
result.ellipsoidal_term = cos_alp0 * sin_alp0 * I12;
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
// Keep track whenever a segment crosses the prime meridian
|
||||
// First normalize to [0,360)
|
||||
template <typename PointOfSegment, typename StateType>
|
||||
static inline int crosses_prime_meridian(PointOfSegment const& p1,
|
||||
PointOfSegment const& p2,
|
||||
StateType& state)
|
||||
{
|
||||
CT const pi
|
||||
= geometry::math::pi<CT>();
|
||||
CT const two_pi
|
||||
= geometry::math::two_pi<CT>();
|
||||
|
||||
CT p1_lon = get_as_radian<0>(p1)
|
||||
- ( floor( get_as_radian<0>(p1) / two_pi )
|
||||
* two_pi );
|
||||
CT p2_lon = get_as_radian<0>(p2)
|
||||
- ( floor( get_as_radian<0>(p2) / two_pi )
|
||||
* two_pi );
|
||||
|
||||
CT max_lon = (std::max)(p1_lon, p2_lon);
|
||||
CT min_lon = (std::min)(p1_lon, p2_lon);
|
||||
|
||||
if(max_lon > pi && min_lon < pi && max_lon - min_lon > pi)
|
||||
{
|
||||
state.m_crosses_prime_meridian++;
|
||||
}
|
||||
|
||||
return state.m_crosses_prime_meridian;
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
}}} // namespace boost::geometry::formula
|
||||
|
||||
|
||||
#endif // BOOST_GEOMETRY_FORMULAS_AREA_FORMULAS_HPP
|
||||
+303
@@ -0,0 +1,303 @@
|
||||
// Boost.Geometry
|
||||
|
||||
// Copyright (c) 2016-2017 Oracle and/or its affiliates.
|
||||
|
||||
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
|
||||
|
||||
// Use, modification and distribution is subject to the Boost Software License,
|
||||
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
|
||||
// http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#ifndef BOOST_GEOMETRY_FORMULAS_INVERSE_DIFFERENTIAL_QUANTITIES_HPP
|
||||
#define BOOST_GEOMETRY_FORMULAS_INVERSE_DIFFERENTIAL_QUANTITIES_HPP
|
||||
|
||||
|
||||
#include <boost/geometry/util/condition.hpp>
|
||||
#include <boost/geometry/util/math.hpp>
|
||||
|
||||
|
||||
namespace boost { namespace geometry { namespace formula
|
||||
{
|
||||
|
||||
/*!
|
||||
\brief The solution of a part of the inverse problem - differential quantities.
|
||||
\author See
|
||||
- Charles F.F Karney, Algorithms for geodesics, 2011
|
||||
https://arxiv.org/pdf/1109.4448.pdf
|
||||
*/
|
||||
template <
|
||||
typename CT,
|
||||
bool EnableReducedLength,
|
||||
bool EnableGeodesicScale,
|
||||
unsigned int Order = 2,
|
||||
bool ApproxF = true
|
||||
>
|
||||
class differential_quantities
|
||||
{
|
||||
public:
|
||||
static inline void apply(CT const& lon1, CT const& lat1,
|
||||
CT const& lon2, CT const& lat2,
|
||||
CT const& azimuth, CT const& reverse_azimuth,
|
||||
CT const& b, CT const& f,
|
||||
CT & reduced_length, CT & geodesic_scale)
|
||||
{
|
||||
CT const dlon = lon2 - lon1;
|
||||
CT const sin_lat1 = sin(lat1);
|
||||
CT const cos_lat1 = cos(lat1);
|
||||
CT const sin_lat2 = sin(lat2);
|
||||
CT const cos_lat2 = cos(lat2);
|
||||
|
||||
apply(dlon, sin_lat1, cos_lat1, sin_lat2, cos_lat2,
|
||||
azimuth, reverse_azimuth,
|
||||
b, f,
|
||||
reduced_length, geodesic_scale);
|
||||
}
|
||||
|
||||
static inline void apply(CT const& dlon,
|
||||
CT const& sin_lat1, CT const& cos_lat1,
|
||||
CT const& sin_lat2, CT const& cos_lat2,
|
||||
CT const& azimuth, CT const& reverse_azimuth,
|
||||
CT const& b, CT const& f,
|
||||
CT & reduced_length, CT & geodesic_scale)
|
||||
{
|
||||
CT const c0 = 0;
|
||||
CT const c1 = 1;
|
||||
CT const one_minus_f = c1 - f;
|
||||
|
||||
CT sin_bet1 = one_minus_f * sin_lat1;
|
||||
CT sin_bet2 = one_minus_f * sin_lat2;
|
||||
|
||||
// equator
|
||||
if (math::equals(sin_bet1, c0) && math::equals(sin_bet2, c0))
|
||||
{
|
||||
CT const sig_12 = math::abs(dlon) / one_minus_f;
|
||||
if (BOOST_GEOMETRY_CONDITION(EnableReducedLength))
|
||||
{
|
||||
CT m12 = sin(sig_12) * b;
|
||||
reduced_length = m12;
|
||||
}
|
||||
|
||||
if (BOOST_GEOMETRY_CONDITION(EnableGeodesicScale))
|
||||
{
|
||||
CT M12 = cos(sig_12);
|
||||
geodesic_scale = M12;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
CT const c2 = 2;
|
||||
CT const e2 = f * (c2 - f);
|
||||
CT const ep2 = e2 / math::sqr(one_minus_f);
|
||||
|
||||
CT const sin_alp1 = sin(azimuth);
|
||||
CT const cos_alp1 = cos(azimuth);
|
||||
//CT const sin_alp2 = sin(reverse_azimuth);
|
||||
CT const cos_alp2 = cos(reverse_azimuth);
|
||||
|
||||
CT cos_bet1 = cos_lat1;
|
||||
CT cos_bet2 = cos_lat2;
|
||||
|
||||
normalize(sin_bet1, cos_bet1);
|
||||
normalize(sin_bet2, cos_bet2);
|
||||
|
||||
CT sin_sig1 = sin_bet1;
|
||||
CT cos_sig1 = cos_alp1 * cos_bet1;
|
||||
CT sin_sig2 = sin_bet2;
|
||||
CT cos_sig2 = cos_alp2 * cos_bet2;
|
||||
|
||||
normalize(sin_sig1, cos_sig1);
|
||||
normalize(sin_sig2, cos_sig2);
|
||||
|
||||
CT const sin_alp0 = sin_alp1 * cos_bet1;
|
||||
CT const cos_alp0_sqr = c1 - math::sqr(sin_alp0);
|
||||
|
||||
CT const J12 = BOOST_GEOMETRY_CONDITION(ApproxF) ?
|
||||
J12_f(sin_sig1, cos_sig1, sin_sig2, cos_sig2, cos_alp0_sqr, f) :
|
||||
J12_ep_sqr(sin_sig1, cos_sig1, sin_sig2, cos_sig2, cos_alp0_sqr, ep2) ;
|
||||
|
||||
CT const dn1 = math::sqrt(c1 + ep2 * math::sqr(sin_bet1));
|
||||
CT const dn2 = math::sqrt(c1 + ep2 * math::sqr(sin_bet2));
|
||||
|
||||
if (BOOST_GEOMETRY_CONDITION(EnableReducedLength))
|
||||
{
|
||||
CT const m12_b = dn2 * (cos_sig1 * sin_sig2)
|
||||
- dn1 * (sin_sig1 * cos_sig2)
|
||||
- cos_sig1 * cos_sig2 * J12;
|
||||
CT const m12 = m12_b * b;
|
||||
|
||||
reduced_length = m12;
|
||||
}
|
||||
|
||||
if (BOOST_GEOMETRY_CONDITION(EnableGeodesicScale))
|
||||
{
|
||||
CT const cos_sig12 = cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2;
|
||||
CT const t = ep2 * (cos_bet1 - cos_bet2) * (cos_bet1 + cos_bet2) / (dn1 + dn2);
|
||||
CT const M12 = cos_sig12 + (t * sin_sig2 - cos_sig2 * J12) * sin_sig1 / dn1;
|
||||
|
||||
geodesic_scale = M12;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
private:
|
||||
/*! Approximation of J12, expanded into taylor series in f
|
||||
Maxima script:
|
||||
ep2: f * (2-f) / ((1-f)^2);
|
||||
k2: ca02 * ep2;
|
||||
assume(f < 1);
|
||||
assume(sig > 0);
|
||||
I1(sig):= integrate(sqrt(1 + k2 * sin(s)^2), s, 0, sig);
|
||||
I2(sig):= integrate(1/sqrt(1 + k2 * sin(s)^2), s, 0, sig);
|
||||
J(sig):= I1(sig) - I2(sig);
|
||||
S: taylor(J(sig), f, 0, 3);
|
||||
S1: factor( 2*integrate(sin(s)^2,s,0,sig)*ca02*f );
|
||||
S2: factor( ((integrate(-6*ca02^2*sin(s)^4+6*ca02*sin(s)^2,s,0,sig)+integrate(-2*ca02^2*sin(s)^4+6*ca02*sin(s)^2,s,0,sig))*f^2)/4 );
|
||||
S3: factor( ((integrate(30*ca02^3*sin(s)^6-54*ca02^2*sin(s)^4+24*ca02*sin(s)^2,s,0,sig)+integrate(6*ca02^3*sin(s)^6-18*ca02^2*sin(s)^4+24*ca02*sin(s)^2,s,0,sig))*f^3)/12 );
|
||||
*/
|
||||
static inline CT J12_f(CT const& sin_sig1, CT const& cos_sig1,
|
||||
CT const& sin_sig2, CT const& cos_sig2,
|
||||
CT const& cos_alp0_sqr, CT const& f)
|
||||
{
|
||||
if (Order == 0)
|
||||
{
|
||||
return 0;
|
||||
}
|
||||
|
||||
CT const c2 = 2;
|
||||
|
||||
CT const sig_12 = atan2(cos_sig1 * sin_sig2 - sin_sig1 * cos_sig2,
|
||||
cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2);
|
||||
CT const sin_2sig1 = c2 * cos_sig1 * sin_sig1; // sin(2sig1)
|
||||
CT const sin_2sig2 = c2 * cos_sig2 * sin_sig2; // sin(2sig2)
|
||||
CT const sin_2sig_12 = sin_2sig2 - sin_2sig1;
|
||||
CT const L1 = sig_12 - sin_2sig_12 / c2;
|
||||
|
||||
if (Order == 1)
|
||||
{
|
||||
return cos_alp0_sqr * f * L1;
|
||||
}
|
||||
|
||||
CT const sin_4sig1 = c2 * sin_2sig1 * (math::sqr(cos_sig1) - math::sqr(sin_sig1)); // sin(4sig1)
|
||||
CT const sin_4sig2 = c2 * sin_2sig2 * (math::sqr(cos_sig2) - math::sqr(sin_sig2)); // sin(4sig2)
|
||||
CT const sin_4sig_12 = sin_4sig2 - sin_4sig1;
|
||||
|
||||
CT const c8 = 8;
|
||||
CT const c12 = 12;
|
||||
CT const c16 = 16;
|
||||
CT const c24 = 24;
|
||||
|
||||
CT const L2 = -( cos_alp0_sqr * sin_4sig_12
|
||||
+ (-c8 * cos_alp0_sqr + c12) * sin_2sig_12
|
||||
+ (c12 * cos_alp0_sqr - c24) * sig_12)
|
||||
/ c16;
|
||||
|
||||
if (Order == 2)
|
||||
{
|
||||
return cos_alp0_sqr * f * (L1 + f * L2);
|
||||
}
|
||||
|
||||
CT const c4 = 4;
|
||||
CT const c9 = 9;
|
||||
CT const c48 = 48;
|
||||
CT const c60 = 60;
|
||||
CT const c64 = 64;
|
||||
CT const c96 = 96;
|
||||
CT const c128 = 128;
|
||||
CT const c144 = 144;
|
||||
|
||||
CT const cos_alp0_quad = math::sqr(cos_alp0_sqr);
|
||||
CT const sin3_2sig1 = math::sqr(sin_2sig1) * sin_2sig1;
|
||||
CT const sin3_2sig2 = math::sqr(sin_2sig2) * sin_2sig2;
|
||||
CT const sin3_2sig_12 = sin3_2sig2 - sin3_2sig1;
|
||||
|
||||
CT const A = (c9 * cos_alp0_quad - c12 * cos_alp0_sqr) * sin_4sig_12;
|
||||
CT const B = c4 * cos_alp0_quad * sin3_2sig_12;
|
||||
CT const C = (-c48 * cos_alp0_quad + c96 * cos_alp0_sqr - c64) * sin_2sig_12;
|
||||
CT const D = (c60 * cos_alp0_quad - c144 * cos_alp0_sqr + c128) * sig_12;
|
||||
|
||||
CT const L3 = (A + B + C + D) / c64;
|
||||
|
||||
// Order 3 and higher
|
||||
return cos_alp0_sqr * f * (L1 + f * (L2 + f * L3));
|
||||
}
|
||||
|
||||
/*! Approximation of J12, expanded into taylor series in e'^2
|
||||
Maxima script:
|
||||
k2: ca02 * ep2;
|
||||
assume(sig > 0);
|
||||
I1(sig):= integrate(sqrt(1 + k2 * sin(s)^2), s, 0, sig);
|
||||
I2(sig):= integrate(1/sqrt(1 + k2 * sin(s)^2), s, 0, sig);
|
||||
J(sig):= I1(sig) - I2(sig);
|
||||
S: taylor(J(sig), ep2, 0, 3);
|
||||
S1: factor( integrate(sin(s)^2,s,0,sig)*ca02*ep2 );
|
||||
S2: factor( (integrate(sin(s)^4,s,0,sig)*ca02^2*ep2^2)/2 );
|
||||
S3: factor( (3*integrate(sin(s)^6,s,0,sig)*ca02^3*ep2^3)/8 );
|
||||
*/
|
||||
static inline CT J12_ep_sqr(CT const& sin_sig1, CT const& cos_sig1,
|
||||
CT const& sin_sig2, CT const& cos_sig2,
|
||||
CT const& cos_alp0_sqr, CT const& ep_sqr)
|
||||
{
|
||||
if (Order == 0)
|
||||
{
|
||||
return 0;
|
||||
}
|
||||
|
||||
CT const c2 = 2;
|
||||
CT const c4 = 4;
|
||||
|
||||
CT const c2a0ep2 = cos_alp0_sqr * ep_sqr;
|
||||
|
||||
CT const sig_12 = atan2(cos_sig1 * sin_sig2 - sin_sig1 * cos_sig2,
|
||||
cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2); // sig2 - sig1
|
||||
CT const sin_2sig1 = c2 * cos_sig1 * sin_sig1; // sin(2sig1)
|
||||
CT const sin_2sig2 = c2 * cos_sig2 * sin_sig2; // sin(2sig2)
|
||||
CT const sin_2sig_12 = sin_2sig2 - sin_2sig1;
|
||||
|
||||
CT const L1 = (c2 * sig_12 - sin_2sig_12) / c4;
|
||||
|
||||
if (Order == 1)
|
||||
{
|
||||
return c2a0ep2 * L1;
|
||||
}
|
||||
|
||||
CT const c8 = 8;
|
||||
CT const c64 = 64;
|
||||
|
||||
CT const sin_4sig1 = c2 * sin_2sig1 * (math::sqr(cos_sig1) - math::sqr(sin_sig1)); // sin(4sig1)
|
||||
CT const sin_4sig2 = c2 * sin_2sig2 * (math::sqr(cos_sig2) - math::sqr(sin_sig2)); // sin(4sig2)
|
||||
CT const sin_4sig_12 = sin_4sig2 - sin_4sig1;
|
||||
|
||||
CT const L2 = (sin_4sig_12 - c8 * sin_2sig_12 + 12 * sig_12) / c64;
|
||||
|
||||
if (Order == 2)
|
||||
{
|
||||
return c2a0ep2 * (L1 + c2a0ep2 * L2);
|
||||
}
|
||||
|
||||
CT const sin3_2sig1 = math::sqr(sin_2sig1) * sin_2sig1;
|
||||
CT const sin3_2sig2 = math::sqr(sin_2sig2) * sin_2sig2;
|
||||
CT const sin3_2sig_12 = sin3_2sig2 - sin3_2sig1;
|
||||
|
||||
CT const c9 = 9;
|
||||
CT const c48 = 48;
|
||||
CT const c60 = 60;
|
||||
CT const c512 = 512;
|
||||
|
||||
CT const L3 = (c9 * sin_4sig_12 + c4 * sin3_2sig_12 - c48 * sin_2sig_12 + c60 * sig_12) / c512;
|
||||
|
||||
// Order 3 and higher
|
||||
return c2a0ep2 * (L1 + c2a0ep2 * (L2 + c2a0ep2 * L3));
|
||||
}
|
||||
|
||||
static inline void normalize(CT & x, CT & y)
|
||||
{
|
||||
CT const len = math::sqrt(math::sqr(x) + math::sqr(y));
|
||||
x /= len;
|
||||
y /= len;
|
||||
}
|
||||
};
|
||||
|
||||
}}} // namespace boost::geometry::formula
|
||||
|
||||
|
||||
#endif // BOOST_GEOMETRY_FORMULAS_INVERSE_DIFFERENTIAL_QUANTITIES_HPP
|
||||
@@ -0,0 +1,70 @@
|
||||
// Boost.Geometry
|
||||
|
||||
// Copyright (c) 2016 Oracle and/or its affiliates.
|
||||
|
||||
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
|
||||
|
||||
// Use, modification and distribution is subject to the Boost Software License,
|
||||
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
|
||||
// http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#ifndef BOOST_GEOMETRY_FORMULAS_ECCENCRICITY_SQR_HPP
|
||||
#define BOOST_GEOMETRY_FORMULAS_ECCENCRICITY_SQR_HPP
|
||||
|
||||
#include <boost/geometry/core/radius.hpp>
|
||||
#include <boost/geometry/core/tag.hpp>
|
||||
#include <boost/geometry/core/tags.hpp>
|
||||
|
||||
#include <boost/geometry/algorithms/not_implemented.hpp>
|
||||
|
||||
namespace boost { namespace geometry
|
||||
{
|
||||
|
||||
#ifndef DOXYGEN_NO_DISPATCH
|
||||
namespace formula_dispatch
|
||||
{
|
||||
|
||||
template <typename ResultType, typename Geometry, typename Tag = typename tag<Geometry>::type>
|
||||
struct eccentricity_sqr
|
||||
: not_implemented<Tag>
|
||||
{};
|
||||
|
||||
template <typename ResultType, typename Geometry>
|
||||
struct eccentricity_sqr<ResultType, Geometry, srs_sphere_tag>
|
||||
{
|
||||
static inline ResultType apply(Geometry const& /*geometry*/)
|
||||
{
|
||||
return ResultType(0);
|
||||
}
|
||||
};
|
||||
|
||||
template <typename ResultType, typename Geometry>
|
||||
struct eccentricity_sqr<ResultType, Geometry, srs_spheroid_tag>
|
||||
{
|
||||
static inline ResultType apply(Geometry const& geometry)
|
||||
{
|
||||
// 1 - (b / a)^2
|
||||
return ResultType(1) - math::sqr(ResultType(get_radius<2>(geometry))
|
||||
/ ResultType(get_radius<0>(geometry)));
|
||||
}
|
||||
};
|
||||
|
||||
} // namespace formula_dispatch
|
||||
#endif // DOXYGEN_NO_DISPATCH
|
||||
|
||||
#ifndef DOXYGEN_NO_DETAIL
|
||||
namespace formula
|
||||
{
|
||||
|
||||
template <typename ResultType, typename Geometry>
|
||||
ResultType eccentricity_sqr(Geometry const& geometry)
|
||||
{
|
||||
return formula_dispatch::eccentricity_sqr<ResultType, Geometry>::apply(geometry);
|
||||
}
|
||||
|
||||
} // namespace formula
|
||||
#endif // DOXYGEN_NO_DETAIL
|
||||
|
||||
}} // namespace boost::geometry
|
||||
|
||||
#endif // BOOST_GEOMETRY_FORMULAS_ECCENCRICITY_SQR_HPP
|
||||
@@ -0,0 +1,70 @@
|
||||
// Boost.Geometry
|
||||
|
||||
// Copyright (c) 2014-2016 Oracle and/or its affiliates.
|
||||
|
||||
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
|
||||
|
||||
// Use, modification and distribution is subject to the Boost Software License,
|
||||
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
|
||||
// http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#ifndef BOOST_GEOMETRY_FORMULAS_FLATTENING_HPP
|
||||
#define BOOST_GEOMETRY_FORMULAS_FLATTENING_HPP
|
||||
|
||||
#include <boost/geometry/core/radius.hpp>
|
||||
#include <boost/geometry/core/tag.hpp>
|
||||
#include <boost/geometry/core/tags.hpp>
|
||||
|
||||
#include <boost/geometry/algorithms/not_implemented.hpp>
|
||||
|
||||
namespace boost { namespace geometry
|
||||
{
|
||||
|
||||
#ifndef DOXYGEN_NO_DISPATCH
|
||||
namespace formula_dispatch
|
||||
{
|
||||
|
||||
template <typename ResultType, typename Geometry, typename Tag = typename tag<Geometry>::type>
|
||||
struct flattening
|
||||
: not_implemented<Tag>
|
||||
{};
|
||||
|
||||
template <typename ResultType, typename Geometry>
|
||||
struct flattening<ResultType, Geometry, srs_sphere_tag>
|
||||
{
|
||||
static inline ResultType apply(Geometry const& /*geometry*/)
|
||||
{
|
||||
return ResultType(0);
|
||||
}
|
||||
};
|
||||
|
||||
template <typename ResultType, typename Geometry>
|
||||
struct flattening<ResultType, Geometry, srs_spheroid_tag>
|
||||
{
|
||||
static inline ResultType apply(Geometry const& geometry)
|
||||
{
|
||||
// (a - b) / a
|
||||
return ResultType(get_radius<0>(geometry) - get_radius<2>(geometry))
|
||||
/ ResultType(get_radius<0>(geometry));
|
||||
}
|
||||
};
|
||||
|
||||
} // namespace formula_dispatch
|
||||
#endif // DOXYGEN_NO_DISPATCH
|
||||
|
||||
#ifndef DOXYGEN_NO_DETAIL
|
||||
namespace formula
|
||||
{
|
||||
|
||||
template <typename ResultType, typename Geometry>
|
||||
ResultType flattening(Geometry const& geometry)
|
||||
{
|
||||
return formula_dispatch::flattening<ResultType, Geometry>::apply(geometry);
|
||||
}
|
||||
|
||||
} // namespace formula
|
||||
#endif // DOXYGEN_NO_DETAIL
|
||||
|
||||
}} // namespace boost::geometry
|
||||
|
||||
#endif // BOOST_GEOMETRY_FORMULAS_FLATTENING_HPP
|
||||
@@ -0,0 +1,457 @@
|
||||
// Boost.Geometry
|
||||
|
||||
// Copyright (c) 2016, Oracle and/or its affiliates.
|
||||
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
|
||||
|
||||
// Use, modification and distribution is subject to the Boost Software License,
|
||||
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
|
||||
// http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#ifndef BOOST_GEOMETRY_FORMULAS_GEOGRAPHIC_HPP
|
||||
#define BOOST_GEOMETRY_FORMULAS_GEOGRAPHIC_HPP
|
||||
|
||||
#include <boost/geometry/core/coordinate_system.hpp>
|
||||
#include <boost/geometry/core/coordinate_type.hpp>
|
||||
#include <boost/geometry/core/access.hpp>
|
||||
#include <boost/geometry/core/radian_access.hpp>
|
||||
|
||||
#include <boost/geometry/arithmetic/arithmetic.hpp>
|
||||
#include <boost/geometry/arithmetic/cross_product.hpp>
|
||||
#include <boost/geometry/arithmetic/dot_product.hpp>
|
||||
#include <boost/geometry/arithmetic/normalize.hpp>
|
||||
|
||||
#include <boost/geometry/formulas/eccentricity_sqr.hpp>
|
||||
#include <boost/geometry/formulas/flattening.hpp>
|
||||
|
||||
#include <boost/geometry/util/math.hpp>
|
||||
#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
|
||||
#include <boost/geometry/util/select_coordinate_type.hpp>
|
||||
|
||||
namespace boost { namespace geometry {
|
||||
|
||||
namespace formula {
|
||||
|
||||
template <typename Point3d, typename PointGeo, typename Spheroid>
|
||||
inline Point3d geo_to_cart3d(PointGeo const& point_geo, Spheroid const& spheroid)
|
||||
{
|
||||
typedef typename coordinate_type<Point3d>::type calc_t;
|
||||
|
||||
calc_t const c1 = 1;
|
||||
calc_t const e_sqr = eccentricity_sqr<calc_t>(spheroid);
|
||||
|
||||
calc_t const lon = get_as_radian<0>(point_geo);
|
||||
calc_t const lat = get_as_radian<1>(point_geo);
|
||||
|
||||
Point3d res;
|
||||
|
||||
calc_t const sin_lat = sin(lat);
|
||||
|
||||
// "unit" spheroid, a = 1
|
||||
calc_t const N = c1 / math::sqrt(c1 - e_sqr * math::sqr(sin_lat));
|
||||
calc_t const N_cos_lat = N * cos(lat);
|
||||
|
||||
set<0>(res, N_cos_lat * cos(lon));
|
||||
set<1>(res, N_cos_lat * sin(lon));
|
||||
set<2>(res, N * (c1 - e_sqr) * sin_lat);
|
||||
|
||||
return res;
|
||||
}
|
||||
|
||||
template <typename PointGeo, typename Spheroid, typename Point3d>
|
||||
inline void geo_to_cart3d(PointGeo const& point_geo, Point3d & result, Point3d & north, Point3d & east, Spheroid const& spheroid)
|
||||
{
|
||||
typedef typename coordinate_type<Point3d>::type calc_t;
|
||||
|
||||
calc_t const c1 = 1;
|
||||
calc_t const e_sqr = eccentricity_sqr<calc_t>(spheroid);
|
||||
|
||||
calc_t const lon = get_as_radian<0>(point_geo);
|
||||
calc_t const lat = get_as_radian<1>(point_geo);
|
||||
|
||||
calc_t const sin_lon = sin(lon);
|
||||
calc_t const cos_lon = cos(lon);
|
||||
calc_t const sin_lat = sin(lat);
|
||||
calc_t const cos_lat = cos(lat);
|
||||
|
||||
// "unit" spheroid, a = 1
|
||||
calc_t const N = c1 / math::sqrt(c1 - e_sqr * math::sqr(sin_lat));
|
||||
calc_t const N_cos_lat = N * cos_lat;
|
||||
|
||||
set<0>(result, N_cos_lat * cos_lon);
|
||||
set<1>(result, N_cos_lat * sin_lon);
|
||||
set<2>(result, N * (c1 - e_sqr) * sin_lat);
|
||||
|
||||
set<0>(east, -sin_lon);
|
||||
set<1>(east, cos_lon);
|
||||
set<2>(east, 0);
|
||||
|
||||
set<0>(north, -sin_lat * cos_lon);
|
||||
set<1>(north, -sin_lat * sin_lon);
|
||||
set<2>(north, cos_lat);
|
||||
}
|
||||
|
||||
template <typename PointGeo, typename Point3d, typename Spheroid>
|
||||
inline PointGeo cart3d_to_geo(Point3d const& point_3d, Spheroid const& spheroid)
|
||||
{
|
||||
typedef typename coordinate_type<PointGeo>::type coord_t;
|
||||
typedef typename coordinate_type<Point3d>::type calc_t;
|
||||
|
||||
calc_t const c1 = 1;
|
||||
//calc_t const c2 = 2;
|
||||
calc_t const e_sqr = eccentricity_sqr<calc_t>(spheroid);
|
||||
|
||||
calc_t const x = get<0>(point_3d);
|
||||
calc_t const y = get<1>(point_3d);
|
||||
calc_t const z = get<2>(point_3d);
|
||||
calc_t const xy_l = math::sqrt(math::sqr(x) + math::sqr(y));
|
||||
|
||||
calc_t const lonr = atan2(y, x);
|
||||
|
||||
// NOTE: Alternative version
|
||||
// http://www.iag-aig.org/attach/989c8e501d9c5b5e2736955baf2632f5/V60N2_5FT.pdf
|
||||
// calc_t const lonr = c2 * atan2(y, x + xy_l);
|
||||
|
||||
calc_t const latr = atan2(z, (c1 - e_sqr) * xy_l);
|
||||
|
||||
// NOTE: If h is equal to 0 then there is no need to improve value of latitude
|
||||
// because then N_i / (N_i + h_i) = 1
|
||||
// http://www.navipedia.net/index.php/Ellipsoidal_and_Cartesian_Coordinates_Conversion
|
||||
|
||||
PointGeo res;
|
||||
|
||||
set_from_radian<0>(res, lonr);
|
||||
set_from_radian<1>(res, latr);
|
||||
|
||||
coord_t lon = get<0>(res);
|
||||
coord_t lat = get<1>(res);
|
||||
|
||||
math::normalize_spheroidal_coordinates
|
||||
<
|
||||
typename coordinate_system<PointGeo>::type::units,
|
||||
coord_t
|
||||
>(lon, lat);
|
||||
|
||||
set<0>(res, lon);
|
||||
set<1>(res, lat);
|
||||
|
||||
return res;
|
||||
}
|
||||
|
||||
template <typename Point3d, typename Spheroid>
|
||||
inline Point3d projected_to_xy(Point3d const& point_3d, Spheroid const& spheroid)
|
||||
{
|
||||
typedef typename coordinate_type<Point3d>::type coord_t;
|
||||
|
||||
// len_xy = sqrt(x^2 + y^2)
|
||||
// r = len_xy - |z / tan(lat)|
|
||||
// assuming h = 0
|
||||
// lat = atan2(z, (1 - e^2) * len_xy);
|
||||
// |z / tan(lat)| = (1 - e^2) * len_xy
|
||||
// r = e^2 * len_xy
|
||||
// x_res = r * cos(lon) = e^2 * len_xy * x / len_xy = e^2 * x
|
||||
// y_res = r * sin(lon) = e^2 * len_xy * y / len_xy = e^2 * y
|
||||
|
||||
coord_t const c0 = 0;
|
||||
coord_t const e_sqr = formula::eccentricity_sqr<coord_t>(spheroid);
|
||||
|
||||
Point3d res;
|
||||
|
||||
set<0>(res, e_sqr * get<0>(point_3d));
|
||||
set<1>(res, e_sqr * get<1>(point_3d));
|
||||
set<2>(res, c0);
|
||||
|
||||
return res;
|
||||
}
|
||||
|
||||
template <typename Point3d, typename Spheroid>
|
||||
inline Point3d projected_to_surface(Point3d const& direction, Spheroid const& spheroid)
|
||||
{
|
||||
typedef typename coordinate_type<Point3d>::type coord_t;
|
||||
|
||||
//coord_t const c0 = 0;
|
||||
coord_t const c2 = 2;
|
||||
coord_t const c4 = 4;
|
||||
|
||||
// calculate the point of intersection of a ray and spheroid's surface
|
||||
// the origin is the origin of the coordinate system
|
||||
//(x*x+y*y)/(a*a) + z*z/(b*b) = 1
|
||||
// x = d.x * t
|
||||
// y = d.y * t
|
||||
// z = d.z * t
|
||||
coord_t const dx = get<0>(direction);
|
||||
coord_t const dy = get<1>(direction);
|
||||
coord_t const dz = get<2>(direction);
|
||||
|
||||
//coord_t const a_sqr = math::sqr(get_radius<0>(spheroid));
|
||||
//coord_t const b_sqr = math::sqr(get_radius<2>(spheroid));
|
||||
// "unit" spheroid, a = 1
|
||||
coord_t const a_sqr = 1;
|
||||
coord_t const b_sqr = math::sqr(get_radius<2>(spheroid) / get_radius<0>(spheroid));
|
||||
|
||||
coord_t const param_a = (dx*dx + dy*dy) / a_sqr + dz*dz / b_sqr;
|
||||
coord_t const delta = c4 * param_a;
|
||||
// delta >= 0
|
||||
coord_t const t = math::sqrt(delta) / (c2 * param_a);
|
||||
|
||||
// result = direction * t
|
||||
Point3d result = direction;
|
||||
multiply_value(result, t);
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
template <typename Point3d, typename Spheroid>
|
||||
inline bool projected_to_surface(Point3d const& origin, Point3d const& direction, Point3d & result1, Point3d & result2, Spheroid const& spheroid)
|
||||
{
|
||||
typedef typename coordinate_type<Point3d>::type coord_t;
|
||||
|
||||
coord_t const c0 = 0;
|
||||
coord_t const c1 = 1;
|
||||
coord_t const c2 = 2;
|
||||
coord_t const c4 = 4;
|
||||
|
||||
// calculate the point of intersection of a ray and spheroid's surface
|
||||
//(x*x+y*y)/(a*a) + z*z/(b*b) = 1
|
||||
// x = o.x + d.x * t
|
||||
// y = o.y + d.y * t
|
||||
// z = o.z + d.z * t
|
||||
coord_t const ox = get<0>(origin);
|
||||
coord_t const oy = get<1>(origin);
|
||||
coord_t const oz = get<2>(origin);
|
||||
coord_t const dx = get<0>(direction);
|
||||
coord_t const dy = get<1>(direction);
|
||||
coord_t const dz = get<2>(direction);
|
||||
|
||||
//coord_t const a_sqr = math::sqr(get_radius<0>(spheroid));
|
||||
//coord_t const b_sqr = math::sqr(get_radius<2>(spheroid));
|
||||
// "unit" spheroid, a = 1
|
||||
coord_t const a_sqr = 1;
|
||||
coord_t const b_sqr = math::sqr(get_radius<2>(spheroid) / get_radius<0>(spheroid));
|
||||
|
||||
coord_t const param_a = (dx*dx + dy*dy) / a_sqr + dz*dz / b_sqr;
|
||||
coord_t const param_b = c2 * ((ox*dx + oy*dy) / a_sqr + oz*dz / b_sqr);
|
||||
coord_t const param_c = (ox*ox + oy*oy) / a_sqr + oz*oz / b_sqr - c1;
|
||||
|
||||
coord_t const delta = math::sqr(param_b) - c4 * param_a*param_c;
|
||||
|
||||
// equals() ?
|
||||
if (delta < c0 || param_a == 0)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
// result = origin + direction * t
|
||||
|
||||
coord_t const sqrt_delta = math::sqrt(delta);
|
||||
coord_t const two_a = c2 * param_a;
|
||||
|
||||
coord_t const t1 = (-param_b + sqrt_delta) / two_a;
|
||||
result1 = direction;
|
||||
multiply_value(result1, t1);
|
||||
add_point(result1, origin);
|
||||
|
||||
coord_t const t2 = (-param_b - sqrt_delta) / two_a;
|
||||
result2 = direction;
|
||||
multiply_value(result2, t2);
|
||||
add_point(result2, origin);
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
template <typename Point3d, typename Spheroid>
|
||||
inline bool great_elliptic_intersection(Point3d const& a1, Point3d const& a2,
|
||||
Point3d const& b1, Point3d const& b2,
|
||||
Point3d & result,
|
||||
Spheroid const& spheroid)
|
||||
{
|
||||
typedef typename coordinate_type<Point3d>::type coord_t;
|
||||
|
||||
coord_t c0 = 0;
|
||||
coord_t c1 = 1;
|
||||
|
||||
Point3d n1 = cross_product(a1, a2);
|
||||
Point3d n2 = cross_product(b1, b2);
|
||||
|
||||
// intersection direction
|
||||
Point3d id = cross_product(n1, n2);
|
||||
coord_t id_len_sqr = dot_product(id, id);
|
||||
|
||||
if (math::equals(id_len_sqr, c0))
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
// no need to normalize a1 and a2 because the intersection point on
|
||||
// the opposite side of the globe is at the same distance from the origin
|
||||
coord_t cos_a1i = dot_product(a1, id);
|
||||
coord_t cos_a2i = dot_product(a2, id);
|
||||
coord_t gri = math::detail::greatest(cos_a1i, cos_a2i);
|
||||
Point3d neg_id = id;
|
||||
multiply_value(neg_id, -c1);
|
||||
coord_t cos_a1ni = dot_product(a1, neg_id);
|
||||
coord_t cos_a2ni = dot_product(a2, neg_id);
|
||||
coord_t grni = math::detail::greatest(cos_a1ni, cos_a2ni);
|
||||
|
||||
if (gri >= grni)
|
||||
{
|
||||
result = projected_to_surface(id, spheroid);
|
||||
}
|
||||
else
|
||||
{
|
||||
result = projected_to_surface(neg_id, spheroid);
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
template <typename Point3d1, typename Point3d2>
|
||||
static inline int elliptic_side_value(Point3d1 const& origin, Point3d1 const& norm, Point3d2 const& pt)
|
||||
{
|
||||
typedef typename coordinate_type<Point3d1>::type calc_t;
|
||||
calc_t c0 = 0;
|
||||
|
||||
// vector oposite to pt - origin
|
||||
// only for the purpose of assigning origin
|
||||
Point3d1 vec = origin;
|
||||
subtract_point(vec, pt);
|
||||
|
||||
calc_t d = dot_product(norm, vec);
|
||||
|
||||
// since the vector is opposite the signs are opposite
|
||||
return math::equals(d, c0) ? 0
|
||||
: d < c0 ? 1
|
||||
: -1; // d > 0
|
||||
}
|
||||
|
||||
template <typename Point3d, typename Spheroid>
|
||||
inline bool planes_spheroid_intersection(Point3d const& o1, Point3d const& n1,
|
||||
Point3d const& o2, Point3d const& n2,
|
||||
Point3d & ip1, Point3d & ip2,
|
||||
Spheroid const& spheroid)
|
||||
{
|
||||
typedef typename coordinate_type<Point3d>::type coord_t;
|
||||
|
||||
coord_t c0 = 0;
|
||||
coord_t c1 = 1;
|
||||
|
||||
// Below
|
||||
// n . (p - o) = 0
|
||||
// n . p - n . o = 0
|
||||
// n . p + d = 0
|
||||
// n . p = h
|
||||
|
||||
// intersection direction
|
||||
Point3d id = cross_product(n1, n2);
|
||||
|
||||
if (math::equals(dot_product(id, id), c0))
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
coord_t dot_n1_n2 = dot_product(n1, n2);
|
||||
coord_t dot_n1_n2_sqr = math::sqr(dot_n1_n2);
|
||||
|
||||
coord_t h1 = dot_product(n1, o1);
|
||||
coord_t h2 = dot_product(n2, o2);
|
||||
|
||||
coord_t denom = c1 - dot_n1_n2_sqr;
|
||||
coord_t C1 = (h1 - h2 * dot_n1_n2) / denom;
|
||||
coord_t C2 = (h2 - h1 * dot_n1_n2) / denom;
|
||||
|
||||
// C1 * n1 + C2 * n2
|
||||
Point3d C1_n1 = n1;
|
||||
multiply_value(C1_n1, C1);
|
||||
Point3d C2_n2 = n2;
|
||||
multiply_value(C2_n2, C2);
|
||||
Point3d io = C1_n1;
|
||||
add_point(io, C2_n2);
|
||||
|
||||
if (! projected_to_surface(io, id, ip1, ip2, spheroid))
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
template <typename Point3d, typename Spheroid>
|
||||
inline void experimental_elliptic_plane(Point3d const& p1, Point3d const& p2,
|
||||
Point3d & v1, Point3d & v2,
|
||||
Point3d & origin, Point3d & normal,
|
||||
Spheroid const& spheroid)
|
||||
{
|
||||
typedef typename coordinate_type<Point3d>::type coord_t;
|
||||
|
||||
Point3d xy1 = projected_to_xy(p1, spheroid);
|
||||
Point3d xy2 = projected_to_xy(p2, spheroid);
|
||||
|
||||
// origin = (xy1 + xy2) / 2
|
||||
origin = xy1;
|
||||
add_point(origin, xy2);
|
||||
multiply_value(origin, coord_t(0.5));
|
||||
|
||||
// v1 = p1 - origin
|
||||
v1 = p1;
|
||||
subtract_point(v1, origin);
|
||||
// v2 = p2 - origin
|
||||
v2 = p2;
|
||||
subtract_point(v2, origin);
|
||||
|
||||
normal = cross_product(v1, v2);
|
||||
}
|
||||
|
||||
template <typename Point3d, typename Spheroid>
|
||||
inline void experimental_elliptic_plane(Point3d const& p1, Point3d const& p2,
|
||||
Point3d & origin, Point3d & normal,
|
||||
Spheroid const& spheroid)
|
||||
{
|
||||
Point3d v1, v2;
|
||||
experimental_elliptic_plane(p1, p2, v1, v2, origin, normal, spheroid);
|
||||
}
|
||||
|
||||
template <typename Point3d, typename Spheroid>
|
||||
inline bool experimental_elliptic_intersection(Point3d const& a1, Point3d const& a2,
|
||||
Point3d const& b1, Point3d const& b2,
|
||||
Point3d & result,
|
||||
Spheroid const& spheroid)
|
||||
{
|
||||
typedef typename coordinate_type<Point3d>::type coord_t;
|
||||
|
||||
coord_t c0 = 0;
|
||||
coord_t c1 = 1;
|
||||
|
||||
Point3d a1v, a2v, o1, n1;
|
||||
experimental_elliptic_plane(a1, a2, a1v, a2v, o1, n1, spheroid);
|
||||
Point3d b1v, b2v, o2, n2;
|
||||
experimental_elliptic_plane(b1, b2, b1v, b2v, o2, n2, spheroid);
|
||||
|
||||
if (! detail::vec_normalize(n1) || ! detail::vec_normalize(n2))
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
Point3d ip1_s, ip2_s;
|
||||
if (! planes_spheroid_intersection(o1, n1, o2, n2, ip1_s, ip2_s, spheroid))
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
// NOTE: simplified test, may not work in all cases
|
||||
coord_t dot_a1i1 = dot_product(a1, ip1_s);
|
||||
coord_t dot_a2i1 = dot_product(a2, ip1_s);
|
||||
coord_t gri1 = math::detail::greatest(dot_a1i1, dot_a2i1);
|
||||
coord_t dot_a1i2 = dot_product(a1, ip2_s);
|
||||
coord_t dot_a2i2 = dot_product(a2, ip2_s);
|
||||
coord_t gri2 = math::detail::greatest(dot_a1i2, dot_a2i2);
|
||||
|
||||
result = gri1 >= gri2 ? ip1_s : ip2_s;
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
} // namespace formula
|
||||
|
||||
}} // namespace boost::geometry
|
||||
|
||||
#endif // BOOST_GEOMETRY_FORMULAS_GEOGRAPHIC_HPP
|
||||
+148
@@ -0,0 +1,148 @@
|
||||
// Boost.Geometry
|
||||
|
||||
// Copyright (c) 2016 Oracle and/or its affiliates.
|
||||
|
||||
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
|
||||
|
||||
// Use, modification and distribution is subject to the Boost Software License,
|
||||
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
|
||||
// http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#ifndef BOOST_GEOMETRY_FORMULAS_GNOMONIC_INTERSECTION_HPP
|
||||
#define BOOST_GEOMETRY_FORMULAS_GNOMONIC_INTERSECTION_HPP
|
||||
|
||||
#include <boost/geometry/core/access.hpp>
|
||||
#include <boost/geometry/core/cs.hpp>
|
||||
|
||||
#include <boost/geometry/arithmetic/cross_product.hpp>
|
||||
#include <boost/geometry/formulas/gnomonic_spheroid.hpp>
|
||||
#include <boost/geometry/geometries/point.hpp>
|
||||
#include <boost/geometry/util/math.hpp>
|
||||
|
||||
|
||||
namespace boost { namespace geometry { namespace formula
|
||||
{
|
||||
|
||||
/*!
|
||||
\brief The intersection of two geodesics using spheroidal gnomonic projection
|
||||
as proposed by Karney.
|
||||
\author See
|
||||
- Charles F.F Karney, Algorithms for geodesics, 2011
|
||||
https://arxiv.org/pdf/1109.4448.pdf
|
||||
- GeographicLib forum thread: Intersection between two geodesic lines
|
||||
https://sourceforge.net/p/geographiclib/discussion/1026621/thread/21aaff9f/
|
||||
*/
|
||||
template
|
||||
<
|
||||
typename CT,
|
||||
template <typename, bool, bool, bool, bool, bool> class Inverse,
|
||||
template <typename, bool, bool, bool, bool> class Direct
|
||||
>
|
||||
class gnomonic_intersection
|
||||
{
|
||||
public:
|
||||
template <typename T1, typename T2, typename Spheroid>
|
||||
static inline bool apply(T1 const& lona1, T1 const& lata1,
|
||||
T1 const& lona2, T1 const& lata2,
|
||||
T2 const& lonb1, T2 const& latb1,
|
||||
T2 const& lonb2, T2 const& latb2,
|
||||
CT & lon, CT & lat,
|
||||
Spheroid const& spheroid)
|
||||
{
|
||||
CT const lon_a1 = lona1;
|
||||
CT const lat_a1 = lata1;
|
||||
CT const lon_a2 = lona2;
|
||||
CT const lat_a2 = lata2;
|
||||
CT const lon_b1 = lonb1;
|
||||
CT const lat_b1 = latb1;
|
||||
CT const lon_b2 = lonb2;
|
||||
CT const lat_b2 = latb2;
|
||||
|
||||
return apply(lon_a1, lat_a1, lon_a2, lat_a2, lon_b1, lat_b1, lon_b2, lat_b2, lon, lat, spheroid);
|
||||
}
|
||||
|
||||
template <typename Spheroid>
|
||||
static inline bool apply(CT const& lona1, CT const& lata1,
|
||||
CT const& lona2, CT const& lata2,
|
||||
CT const& lonb1, CT const& latb1,
|
||||
CT const& lonb2, CT const& latb2,
|
||||
CT & lon, CT & lat,
|
||||
Spheroid const& spheroid)
|
||||
{
|
||||
typedef gnomonic_spheroid<CT, Inverse, Direct> gnom_t;
|
||||
|
||||
lon = (lona1 + lona2 + lonb1 + lonb2) / 4;
|
||||
lat = (lata1 + lata2 + latb1 + latb2) / 4;
|
||||
// TODO: consider normalizing lon
|
||||
|
||||
for (int i = 0; i < 10; ++i)
|
||||
{
|
||||
CT xa1, ya1, xa2, ya2;
|
||||
CT xb1, yb1, xb2, yb2;
|
||||
CT x, y;
|
||||
double lat1, lon1;
|
||||
|
||||
bool ok = gnom_t::forward(lon, lat, lona1, lata1, xa1, ya1, spheroid)
|
||||
&& gnom_t::forward(lon, lat, lona2, lata2, xa2, ya2, spheroid)
|
||||
&& gnom_t::forward(lon, lat, lonb1, latb1, xb1, yb1, spheroid)
|
||||
&& gnom_t::forward(lon, lat, lonb2, latb2, xb2, yb2, spheroid)
|
||||
&& intersect(xa1, ya1, xa2, ya2, xb1, yb1, xb2, yb2, x, y)
|
||||
&& gnom_t::inverse(lon, lat, x, y, lon1, lat1, spheroid);
|
||||
|
||||
if (! ok)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
if (math::equals(lat1, lat) && math::equals(lon1, lon))
|
||||
{
|
||||
break;
|
||||
}
|
||||
|
||||
lat = lat1;
|
||||
lon = lon1;
|
||||
}
|
||||
|
||||
// NOTE: true is also returned if the number of iterations is too great
|
||||
// which means that the accuracy of the result is low
|
||||
return true;
|
||||
}
|
||||
|
||||
private:
|
||||
static inline bool intersect(CT const& xa1, CT const& ya1, CT const& xa2, CT const& ya2,
|
||||
CT const& xb1, CT const& yb1, CT const& xb2, CT const& yb2,
|
||||
CT & x, CT & y)
|
||||
{
|
||||
typedef model::point<CT, 3, cs::cartesian> v3d_t;
|
||||
|
||||
CT const c0 = 0;
|
||||
CT const c1 = 1;
|
||||
|
||||
v3d_t const va1(xa1, ya1, c1);
|
||||
v3d_t const va2(xa2, ya2, c1);
|
||||
v3d_t const vb1(xb1, yb1, c1);
|
||||
v3d_t const vb2(xb2, yb2, c1);
|
||||
|
||||
v3d_t const la = cross_product(va1, va2);
|
||||
v3d_t const lb = cross_product(vb1, vb2);
|
||||
v3d_t const p = cross_product(la, lb);
|
||||
|
||||
CT const z = get<2>(p);
|
||||
|
||||
if (math::equals(z, c0))
|
||||
{
|
||||
// degenerated or collinear segments
|
||||
return false;
|
||||
}
|
||||
|
||||
x = get<0>(p) / z;
|
||||
y = get<1>(p) / z;
|
||||
|
||||
return true;
|
||||
}
|
||||
};
|
||||
|
||||
}}} // namespace boost::geometry::formula
|
||||
|
||||
|
||||
#endif // BOOST_GEOMETRY_FORMULAS_GNOMONIC_INTERSECTION_HPP
|
||||
@@ -0,0 +1,125 @@
|
||||
// Boost.Geometry
|
||||
|
||||
// Copyright (c) 2016 Oracle and/or its affiliates.
|
||||
|
||||
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
|
||||
|
||||
// Use, modification and distribution is subject to the Boost Software License,
|
||||
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
|
||||
// http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#ifndef BOOST_GEOMETRY_FORMULAS_GNOMONIC_SPHEROID_HPP
|
||||
#define BOOST_GEOMETRY_FORMULAS_GNOMONIC_SPHEROID_HPP
|
||||
|
||||
|
||||
#include <boost/geometry/core/radius.hpp>
|
||||
|
||||
#include <boost/geometry/util/condition.hpp>
|
||||
#include <boost/geometry/util/math.hpp>
|
||||
|
||||
#include <boost/geometry/formulas/andoyer_inverse.hpp>
|
||||
#include <boost/geometry/formulas/flattening.hpp>
|
||||
#include <boost/geometry/formulas/thomas_inverse.hpp>
|
||||
#include <boost/geometry/formulas/vincenty_direct.hpp>
|
||||
#include <boost/geometry/formulas/vincenty_inverse.hpp>
|
||||
|
||||
|
||||
namespace boost { namespace geometry { namespace formula
|
||||
{
|
||||
|
||||
/*!
|
||||
\brief Gnomonic projection on spheroid (ellipsoid of revolution).
|
||||
\author See
|
||||
- Charles F.F Karney, Algorithms for geodesics, 2011
|
||||
https://arxiv.org/pdf/1109.4448.pdf
|
||||
*/
|
||||
template <
|
||||
typename CT,
|
||||
template <typename, bool, bool, bool, bool ,bool> class Inverse,
|
||||
template <typename, bool, bool, bool, bool> class Direct
|
||||
>
|
||||
class gnomonic_spheroid
|
||||
{
|
||||
typedef Inverse<CT, false, true, true, true, true> inverse_type;
|
||||
typedef typename inverse_type::result_type inverse_result;
|
||||
|
||||
typedef Direct<CT, false, false, true, true> direct_quantities_type;
|
||||
typedef Direct<CT, true, false, false, false> direct_coordinates_type;
|
||||
typedef typename direct_coordinates_type::result_type direct_result;
|
||||
|
||||
public:
|
||||
template <typename Spheroid>
|
||||
static inline bool forward(CT const& lon0, CT const& lat0,
|
||||
CT const& lon, CT const& lat,
|
||||
CT & x, CT & y,
|
||||
Spheroid const& spheroid)
|
||||
{
|
||||
inverse_result i_res = inverse_type::apply(lon0, lat0, lon, lat, spheroid);
|
||||
CT const& m = i_res.reduced_length;
|
||||
CT const& M = i_res.geodesic_scale;
|
||||
|
||||
if (math::smaller_or_equals(M, CT(0)))
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
CT rho = m / M;
|
||||
x = sin(i_res.azimuth) * rho;
|
||||
y = cos(i_res.azimuth) * rho;
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
template <typename Spheroid>
|
||||
static inline bool inverse(CT const& lon0, CT const& lat0,
|
||||
CT const& x, CT const& y,
|
||||
CT & lon, CT & lat,
|
||||
Spheroid const& spheroid)
|
||||
{
|
||||
CT const a = get_radius<0>(spheroid);
|
||||
CT const ds_threshold = a * std::numeric_limits<CT>::epsilon(); // TODO: 0 for non-fundamental type
|
||||
|
||||
CT const azimuth = atan2(x, y);
|
||||
CT const rho = math::sqrt(math::sqr(x) + math::sqr(y)); // use hypot?
|
||||
CT distance = a * atan(rho / a);
|
||||
|
||||
bool found = false;
|
||||
for (int i = 0 ; i < 10 ; ++i)
|
||||
{
|
||||
direct_result d_res = direct_quantities_type::apply(lon0, lat0, distance, azimuth, spheroid);
|
||||
CT const& m = d_res.reduced_length;
|
||||
CT const& M = d_res.geodesic_scale;
|
||||
|
||||
if (math::smaller_or_equals(M, CT(0)))
|
||||
{
|
||||
// found = false;
|
||||
return found;
|
||||
}
|
||||
|
||||
CT const drho = m / M - rho; // rho = m / M
|
||||
CT const ds = drho * math::sqr(M); // drho/ds = 1/M^2
|
||||
distance -= ds;
|
||||
|
||||
// ds_threshold may be 0
|
||||
if (math::abs(ds) <= ds_threshold)
|
||||
{
|
||||
found = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if (found)
|
||||
{
|
||||
direct_result d_res = direct_coordinates_type::apply(lon0, lat0, distance, azimuth, spheroid);
|
||||
lon = d_res.lon2;
|
||||
lat = d_res.lat2;
|
||||
}
|
||||
|
||||
return found;
|
||||
}
|
||||
};
|
||||
|
||||
}}} // namespace boost::geometry::formula
|
||||
|
||||
|
||||
#endif // BOOST_GEOMETRY_FORMULAS_GNOMONIC_SPHEROID_HPP
|
||||
@@ -0,0 +1,39 @@
|
||||
// Boost.Geometry
|
||||
|
||||
// Copyright (c) 2016 Oracle and/or its affiliates.
|
||||
|
||||
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
|
||||
|
||||
// Use, modification and distribution is subject to the Boost Software License,
|
||||
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
|
||||
// http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#ifndef BOOST_GEOMETRY_FORMULAS_RESULT_DIRECT_HPP
|
||||
#define BOOST_GEOMETRY_FORMULAS_RESULT_DIRECT_HPP
|
||||
|
||||
|
||||
namespace boost { namespace geometry { namespace formula
|
||||
{
|
||||
|
||||
template <typename T>
|
||||
struct result_direct
|
||||
{
|
||||
result_direct()
|
||||
: lon2(0)
|
||||
, lat2(0)
|
||||
, reverse_azimuth(0)
|
||||
, reduced_length(0)
|
||||
, geodesic_scale(1)
|
||||
{}
|
||||
|
||||
T lon2;
|
||||
T lat2;
|
||||
T reverse_azimuth;
|
||||
T reduced_length;
|
||||
T geodesic_scale;
|
||||
};
|
||||
|
||||
}}} // namespace boost::geometry::formula
|
||||
|
||||
|
||||
#endif // BOOST_GEOMETRY_FORMULAS_RESULT_DIRECT_HPP
|
||||
@@ -0,0 +1,39 @@
|
||||
// Boost.Geometry
|
||||
|
||||
// Copyright (c) 2015-2016 Oracle and/or its affiliates.
|
||||
|
||||
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
|
||||
|
||||
// Use, modification and distribution is subject to the Boost Software License,
|
||||
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
|
||||
// http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#ifndef BOOST_GEOMETRY_FORMULAS_RESULT_INVERSE_HPP
|
||||
#define BOOST_GEOMETRY_FORMULAS_RESULT_INVERSE_HPP
|
||||
|
||||
|
||||
namespace boost { namespace geometry { namespace formula
|
||||
{
|
||||
|
||||
template <typename T>
|
||||
struct result_inverse
|
||||
{
|
||||
result_inverse()
|
||||
: distance(0)
|
||||
, azimuth(0)
|
||||
, reverse_azimuth(0)
|
||||
, reduced_length(0)
|
||||
, geodesic_scale(1)
|
||||
{}
|
||||
|
||||
T distance;
|
||||
T azimuth;
|
||||
T reverse_azimuth;
|
||||
T reduced_length;
|
||||
T geodesic_scale;
|
||||
};
|
||||
|
||||
}}} // namespace boost::geometry::formula
|
||||
|
||||
|
||||
#endif // BOOST_GEOMETRY_FORMULAS_RESULT_INVERSE_HPP
|
||||
+1219
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,224 @@
|
||||
// Boost.Geometry
|
||||
|
||||
// Copyright (c) 2016, Oracle and/or its affiliates.
|
||||
// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
|
||||
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
|
||||
|
||||
// Use, modification and distribution is subject to the Boost Software License,
|
||||
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
|
||||
// http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#ifndef BOOST_GEOMETRY_FORMULAS_SPHERICAL_HPP
|
||||
#define BOOST_GEOMETRY_FORMULAS_SPHERICAL_HPP
|
||||
|
||||
#include <boost/geometry/core/coordinate_system.hpp>
|
||||
#include <boost/geometry/core/coordinate_type.hpp>
|
||||
#include <boost/geometry/core/access.hpp>
|
||||
#include <boost/geometry/core/radian_access.hpp>
|
||||
|
||||
//#include <boost/geometry/arithmetic/arithmetic.hpp>
|
||||
#include <boost/geometry/arithmetic/cross_product.hpp>
|
||||
#include <boost/geometry/arithmetic/dot_product.hpp>
|
||||
|
||||
#include <boost/geometry/util/math.hpp>
|
||||
#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
|
||||
#include <boost/geometry/util/select_coordinate_type.hpp>
|
||||
|
||||
namespace boost { namespace geometry {
|
||||
|
||||
namespace formula {
|
||||
|
||||
template <typename T>
|
||||
struct result_spherical
|
||||
{
|
||||
result_spherical()
|
||||
: azimuth(0)
|
||||
, reverse_azimuth(0)
|
||||
{}
|
||||
|
||||
T azimuth;
|
||||
T reverse_azimuth;
|
||||
};
|
||||
|
||||
template <typename T>
|
||||
static inline void sph_to_cart3d(T const& lon, T const& lat, T & x, T & y, T & z)
|
||||
{
|
||||
T const cos_lat = cos(lat);
|
||||
x = cos_lat * cos(lon);
|
||||
y = cos_lat * sin(lon);
|
||||
z = sin(lat);
|
||||
}
|
||||
|
||||
template <typename Point3d, typename PointSph>
|
||||
static inline Point3d sph_to_cart3d(PointSph const& point_sph)
|
||||
{
|
||||
typedef typename coordinate_type<Point3d>::type calc_t;
|
||||
|
||||
calc_t const lon = get_as_radian<0>(point_sph);
|
||||
calc_t const lat = get_as_radian<1>(point_sph);
|
||||
calc_t x, y, z;
|
||||
sph_to_cart3d(lon, lat, x, y, z);
|
||||
|
||||
Point3d res;
|
||||
set<0>(res, x);
|
||||
set<1>(res, y);
|
||||
set<2>(res, z);
|
||||
|
||||
return res;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
static inline void cart3d_to_sph(T const& x, T const& y, T const& z, T & lon, T & lat)
|
||||
{
|
||||
lon = atan2(y, x);
|
||||
lat = asin(z);
|
||||
}
|
||||
|
||||
template <typename PointSph, typename Point3d>
|
||||
static inline PointSph cart3d_to_sph(Point3d const& point_3d)
|
||||
{
|
||||
typedef typename coordinate_type<PointSph>::type coord_t;
|
||||
typedef typename coordinate_type<Point3d>::type calc_t;
|
||||
|
||||
calc_t const x = get<0>(point_3d);
|
||||
calc_t const y = get<1>(point_3d);
|
||||
calc_t const z = get<2>(point_3d);
|
||||
calc_t lonr, latr;
|
||||
cart3d_to_sph(x, y, z, lonr, latr);
|
||||
|
||||
PointSph res;
|
||||
set_from_radian<0>(res, lonr);
|
||||
set_from_radian<1>(res, latr);
|
||||
|
||||
coord_t lon = get<0>(res);
|
||||
coord_t lat = get<1>(res);
|
||||
|
||||
math::normalize_spheroidal_coordinates
|
||||
<
|
||||
typename coordinate_system<PointSph>::type::units,
|
||||
coord_t
|
||||
>(lon, lat);
|
||||
|
||||
set<0>(res, lon);
|
||||
set<1>(res, lat);
|
||||
|
||||
return res;
|
||||
}
|
||||
|
||||
// -1 right
|
||||
// 1 left
|
||||
// 0 on
|
||||
template <typename Point3d1, typename Point3d2>
|
||||
static inline int sph_side_value(Point3d1 const& norm, Point3d2 const& pt)
|
||||
{
|
||||
typedef typename select_coordinate_type<Point3d1, Point3d2>::type calc_t;
|
||||
calc_t c0 = 0;
|
||||
calc_t d = dot_product(norm, pt);
|
||||
return math::equals(d, c0) ? 0
|
||||
: d > c0 ? 1
|
||||
: -1; // d < 0
|
||||
}
|
||||
|
||||
template <typename CT, bool ReverseAzimuth, typename T1, typename T2>
|
||||
static inline result_spherical<CT> spherical_azimuth(T1 const& lon1,
|
||||
T1 const& lat1,
|
||||
T2 const& lon2,
|
||||
T2 const& lat2)
|
||||
{
|
||||
typedef result_spherical<CT> result_type;
|
||||
result_type result;
|
||||
|
||||
// http://williams.best.vwh.net/avform.htm#Crs
|
||||
// https://en.wikipedia.org/wiki/Great-circle_navigation
|
||||
CT dlon = lon2 - lon1;
|
||||
|
||||
// An optimization which should kick in often for Boxes
|
||||
//if ( math::equals(dlon, ReturnType(0)) )
|
||||
//if ( get<0>(p1) == get<0>(p2) )
|
||||
//{
|
||||
// return - sin(get_as_radian<1>(p1)) * cos_p2lat);
|
||||
//}
|
||||
|
||||
CT const cos_dlon = cos(dlon);
|
||||
CT const sin_dlon = sin(dlon);
|
||||
CT const cos_lat1 = cos(lat1);
|
||||
CT const cos_lat2 = cos(lat2);
|
||||
CT const sin_lat1 = sin(lat1);
|
||||
CT const sin_lat2 = sin(lat2);
|
||||
|
||||
{
|
||||
// "An alternative formula, not requiring the pre-computation of d"
|
||||
// In the formula below dlon is used as "d"
|
||||
CT const y = sin_dlon * cos_lat2;
|
||||
CT const x = cos_lat1 * sin_lat2 - sin_lat1 * cos_lat2 * cos_dlon;
|
||||
result.azimuth = atan2(y, x);
|
||||
}
|
||||
|
||||
if (ReverseAzimuth)
|
||||
{
|
||||
CT const y = sin_dlon * cos_lat1;
|
||||
CT const x = sin_lat2 * cos_lat1 * cos_dlon - cos_lat2 * sin_lat1;
|
||||
result.reverse_azimuth = atan2(y, x);
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
template <typename ReturnType, typename T1, typename T2>
|
||||
inline ReturnType spherical_azimuth(T1 const& lon1, T1 const& lat1,
|
||||
T2 const& lon2, T2 const& lat2)
|
||||
{
|
||||
return spherical_azimuth<ReturnType, false>(lon1, lat1, lon2, lat2).azimuth;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
inline T spherical_azimuth(T const& lon1, T const& lat1, T const& lon2, T const& lat2)
|
||||
{
|
||||
return spherical_azimuth<T, false>(lon1, lat1, lon2, lat2).azimuth;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
inline int azimuth_side_value(T const& azi_a1_p, T const& azi_a1_a2)
|
||||
{
|
||||
T const pi = math::pi<T>();
|
||||
T const two_pi = math::two_pi<T>();
|
||||
|
||||
// instead of the formula from XTD
|
||||
//calc_t a_diff = asin(sin(azi_a1_p - azi_a1_a2));
|
||||
|
||||
T a_diff = azi_a1_p - azi_a1_a2;
|
||||
// normalize, angle in [-pi, pi]
|
||||
while (a_diff > pi)
|
||||
a_diff -= two_pi;
|
||||
while (a_diff < -pi)
|
||||
a_diff += two_pi;
|
||||
|
||||
// NOTE: in general it shouldn't be required to support the pi/-pi case
|
||||
// because in non-cartesian systems it makes sense to check the side
|
||||
// only "between" the endpoints.
|
||||
// However currently the winding strategy calls the side strategy
|
||||
// for vertical segments to check if the point is "between the endpoints.
|
||||
// This could be avoided since the side strategy is not required for that
|
||||
// because meridian is the shortest path. So a difference of
|
||||
// longitudes would be sufficient (of course normalized to [-pi, pi]).
|
||||
|
||||
// NOTE: with the above said, the pi/-pi check is temporary
|
||||
// however in case if this was required
|
||||
// the geodesics on ellipsoid aren't "symmetrical"
|
||||
// therefore instead of comparing a_diff to pi and -pi
|
||||
// one should probably use inverse azimuths and compare
|
||||
// the difference to 0 as well
|
||||
|
||||
// positive azimuth is on the right side
|
||||
return math::equals(a_diff, 0)
|
||||
|| math::equals(a_diff, pi)
|
||||
|| math::equals(a_diff, -pi) ? 0
|
||||
: a_diff > 0 ? -1 // right
|
||||
: 1; // left
|
||||
}
|
||||
|
||||
} // namespace formula
|
||||
|
||||
}} // namespace boost::geometry
|
||||
|
||||
#endif // BOOST_GEOMETRY_FORMULAS_SPHERICAL_HPP
|
||||
@@ -0,0 +1,248 @@
|
||||
// Boost.Geometry
|
||||
|
||||
// Copyright (c) 2016 Oracle and/or its affiliates.
|
||||
|
||||
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
|
||||
|
||||
// Use, modification and distribution is subject to the Boost Software License,
|
||||
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
|
||||
// http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#ifndef BOOST_GEOMETRY_FORMULAS_THOMAS_DIRECT_HPP
|
||||
#define BOOST_GEOMETRY_FORMULAS_THOMAS_DIRECT_HPP
|
||||
|
||||
|
||||
#include <boost/math/constants/constants.hpp>
|
||||
|
||||
#include <boost/geometry/core/radius.hpp>
|
||||
#include <boost/geometry/core/srs.hpp>
|
||||
|
||||
#include <boost/geometry/util/condition.hpp>
|
||||
#include <boost/geometry/util/math.hpp>
|
||||
|
||||
#include <boost/geometry/formulas/differential_quantities.hpp>
|
||||
#include <boost/geometry/formulas/flattening.hpp>
|
||||
#include <boost/geometry/formulas/result_direct.hpp>
|
||||
|
||||
|
||||
namespace boost { namespace geometry { namespace formula
|
||||
{
|
||||
|
||||
|
||||
/*!
|
||||
\brief The solution of the direct problem of geodesics on latlong coordinates,
|
||||
Forsyth-Andoyer-Lambert type approximation with second order terms.
|
||||
\author See
|
||||
- Technical Report: PAUL D. THOMAS, MATHEMATICAL MODELS FOR NAVIGATION SYSTEMS, 1965
|
||||
http://www.dtic.mil/docs/citations/AD0627893
|
||||
- Technical Report: PAUL D. THOMAS, SPHEROIDAL GEODESICS, REFERENCE SYSTEMS, AND LOCAL GEOMETRY, 1970
|
||||
http://www.dtic.mil/docs/citations/AD0703541
|
||||
|
||||
*/
|
||||
template <
|
||||
typename CT,
|
||||
bool EnableCoordinates = true,
|
||||
bool EnableReverseAzimuth = false,
|
||||
bool EnableReducedLength = false,
|
||||
bool EnableGeodesicScale = false
|
||||
>
|
||||
class thomas_direct
|
||||
{
|
||||
static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale;
|
||||
static const bool CalcCoordinates = EnableCoordinates || CalcQuantities;
|
||||
static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcCoordinates || CalcQuantities;
|
||||
|
||||
public:
|
||||
typedef result_direct<CT> result_type;
|
||||
|
||||
template <typename T, typename Dist, typename Azi, typename Spheroid>
|
||||
static inline result_type apply(T const& lo1,
|
||||
T const& la1,
|
||||
Dist const& distance,
|
||||
Azi const& azimuth12,
|
||||
Spheroid const& spheroid)
|
||||
{
|
||||
result_type result;
|
||||
|
||||
CT const lon1 = lo1;
|
||||
CT const lat1 = la1;
|
||||
|
||||
if ( math::equals(distance, Dist(0)) || distance < Dist(0) )
|
||||
{
|
||||
result.lon2 = lon1;
|
||||
result.lat2 = lat1;
|
||||
return result;
|
||||
}
|
||||
|
||||
CT const c0 = 0;
|
||||
CT const c1 = 1;
|
||||
CT const c2 = 2;
|
||||
CT const c4 = 4;
|
||||
|
||||
CT const a = CT(get_radius<0>(spheroid));
|
||||
CT const b = CT(get_radius<2>(spheroid));
|
||||
CT const f = formula::flattening<CT>(spheroid);
|
||||
CT const one_minus_f = c1 - f;
|
||||
|
||||
CT const pi = math::pi<CT>();
|
||||
CT const pi_half = pi / c2;
|
||||
|
||||
// keep azimuth small - experiments show low accuracy
|
||||
// if the azimuth is closer to (+-)180 deg.
|
||||
CT azi12_alt = azimuth12;
|
||||
CT lat1_alt = lat1;
|
||||
bool alter_result = vflip_if_south(lat1, azimuth12, lat1_alt, azi12_alt);
|
||||
|
||||
CT const theta1 = math::equals(lat1_alt, pi_half) ? lat1_alt :
|
||||
math::equals(lat1_alt, -pi_half) ? lat1_alt :
|
||||
atan(one_minus_f * tan(lat1_alt));
|
||||
CT const sin_theta1 = sin(theta1);
|
||||
CT const cos_theta1 = cos(theta1);
|
||||
|
||||
CT const sin_a12 = sin(azi12_alt);
|
||||
CT const cos_a12 = cos(azi12_alt);
|
||||
|
||||
CT const M = cos_theta1 * sin_a12; // cos_theta0
|
||||
CT const theta0 = acos(M);
|
||||
CT const sin_theta0 = sin(theta0);
|
||||
|
||||
CT const N = cos_theta1 * cos_a12;
|
||||
CT const C1 = f * M; // lower-case c1 in the technical report
|
||||
CT const C2 = f * (c1 - math::sqr(M)) / c4; // lower-case c2 in the technical report
|
||||
CT const D = (c1 - C2) * (c1 - C2 - C1 * M);
|
||||
CT const P = C2 * (c1 + C1 * M / c2) / D;
|
||||
|
||||
// special case for equator:
|
||||
// sin_theta0 = 0 <=> lat1 = 0 ^ |azimuth12| = pi/2
|
||||
// NOTE: in this case it doesn't matter what's the value of cos_sigma1 because
|
||||
// theta1=0, theta0=0, M=1|-1, C2=0 so X=0 and Y=0 so d_sigma=d
|
||||
// cos_a12=0 so N=0, therefore
|
||||
// lat2=0, azi21=pi/2|-pi/2
|
||||
// d_eta = atan2(sin_d_sigma, cos_d_sigma)
|
||||
// H = C1 * d_sigma
|
||||
CT const cos_sigma1 = math::equals(sin_theta0, c0)
|
||||
? c1
|
||||
: normalized1_1(sin_theta1 / sin_theta0);
|
||||
CT const sigma1 = acos(cos_sigma1);
|
||||
CT const d = distance / (a * D);
|
||||
CT const u = 2 * (sigma1 - d);
|
||||
CT const cos_d = cos(d);
|
||||
CT const sin_d = sin(d);
|
||||
CT const cos_u = cos(u);
|
||||
CT const sin_u = sin(u);
|
||||
|
||||
CT const W = c1 - c2 * P * cos_u;
|
||||
CT const V = cos_u * cos_d - sin_u * sin_d;
|
||||
CT const X = math::sqr(C2) * sin_d * cos_d * (2 * math::sqr(V) - c1);
|
||||
CT const Y = c2 * P * V * W * sin_d;
|
||||
CT const d_sigma = d + X - Y;
|
||||
CT const sin_d_sigma = sin(d_sigma);
|
||||
CT const cos_d_sigma = cos(d_sigma);
|
||||
|
||||
if (BOOST_GEOMETRY_CONDITION(CalcRevAzimuth))
|
||||
{
|
||||
result.reverse_azimuth = atan2(M, N * cos_d_sigma - sin_theta1 * sin_d_sigma);
|
||||
|
||||
if (alter_result)
|
||||
{
|
||||
vflip_rev_azi(result.reverse_azimuth, azimuth12);
|
||||
}
|
||||
}
|
||||
|
||||
if (BOOST_GEOMETRY_CONDITION(CalcCoordinates))
|
||||
{
|
||||
CT const S_sigma = c2 * sigma1 - d_sigma;
|
||||
CT const cos_S_sigma = cos(S_sigma);
|
||||
CT const d_eta = atan2(sin_d_sigma * sin_a12, cos_theta1 * cos_d_sigma - sin_theta1 * sin_d_sigma * cos_a12);
|
||||
CT const H = C1 * (c1 - C2) * d_sigma - C1 * C2 * sin_d_sigma * cos_S_sigma;
|
||||
CT const d_lambda = d_eta - H;
|
||||
|
||||
result.lon2 = lon1 + d_lambda;
|
||||
|
||||
if (! math::equals(M, c0))
|
||||
{
|
||||
CT const sin_a21 = sin(result.reverse_azimuth);
|
||||
CT const tan_theta2 = (sin_theta1 * cos_d_sigma + N * sin_d_sigma) * sin_a21 / M;
|
||||
result.lat2 = atan(tan_theta2 / one_minus_f);
|
||||
}
|
||||
else
|
||||
{
|
||||
CT const sigma2 = S_sigma - sigma1;
|
||||
//theta2 = asin(cos(sigma2)) <=> sin_theta0 = 1
|
||||
CT const tan_theta2 = cos(sigma2) / sin(sigma2);
|
||||
result.lat2 = atan(tan_theta2 / one_minus_f);
|
||||
}
|
||||
|
||||
if (alter_result)
|
||||
{
|
||||
result.lat2 = -result.lat2;
|
||||
}
|
||||
}
|
||||
|
||||
if (BOOST_GEOMETRY_CONDITION(CalcQuantities))
|
||||
{
|
||||
typedef differential_quantities<CT, EnableReducedLength, EnableGeodesicScale, 2> quantities;
|
||||
quantities::apply(lon1, lat1, result.lon2, result.lat2,
|
||||
azimuth12, result.reverse_azimuth,
|
||||
b, f,
|
||||
result.reduced_length, result.geodesic_scale);
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
private:
|
||||
static inline bool vflip_if_south(CT const& lat1, CT const& azi12, CT & lat1_alt, CT & azi12_alt)
|
||||
{
|
||||
CT const c2 = 2;
|
||||
CT const pi = math::pi<CT>();
|
||||
CT const pi_half = pi / c2;
|
||||
|
||||
if (azi12 > pi_half)
|
||||
{
|
||||
azi12_alt = pi - azi12;
|
||||
lat1_alt = -lat1;
|
||||
return true;
|
||||
}
|
||||
else if (azi12 < -pi_half)
|
||||
{
|
||||
azi12_alt = -pi - azi12;
|
||||
lat1_alt = -lat1;
|
||||
return true;
|
||||
}
|
||||
|
||||
return false;
|
||||
}
|
||||
|
||||
static inline void vflip_rev_azi(CT & rev_azi, CT const& azimuth12)
|
||||
{
|
||||
CT const c0 = 0;
|
||||
CT const pi = math::pi<CT>();
|
||||
|
||||
if (rev_azi == c0)
|
||||
{
|
||||
rev_azi = azimuth12 >= 0 ? pi : -pi;
|
||||
}
|
||||
else if (rev_azi > c0)
|
||||
{
|
||||
rev_azi = pi - rev_azi;
|
||||
}
|
||||
else
|
||||
{
|
||||
rev_azi = -pi - rev_azi;
|
||||
}
|
||||
}
|
||||
|
||||
static inline CT normalized1_1(CT const& value)
|
||||
{
|
||||
CT const c1 = 1;
|
||||
return value > c1 ? c1 :
|
||||
value < -c1 ? -c1 :
|
||||
value;
|
||||
}
|
||||
};
|
||||
|
||||
}}} // namespace boost::geometry::formula
|
||||
|
||||
|
||||
#endif // BOOST_GEOMETRY_FORMULAS_THOMAS_DIRECT_HPP
|
||||
@@ -0,0 +1,220 @@
|
||||
// Boost.Geometry
|
||||
|
||||
// Copyright (c) 2015-2016 Oracle and/or its affiliates.
|
||||
|
||||
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
|
||||
|
||||
// Use, modification and distribution is subject to the Boost Software License,
|
||||
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
|
||||
// http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#ifndef BOOST_GEOMETRY_FORMULAS_THOMAS_INVERSE_HPP
|
||||
#define BOOST_GEOMETRY_FORMULAS_THOMAS_INVERSE_HPP
|
||||
|
||||
|
||||
#include <boost/math/constants/constants.hpp>
|
||||
|
||||
#include <boost/geometry/core/radius.hpp>
|
||||
#include <boost/geometry/core/srs.hpp>
|
||||
|
||||
#include <boost/geometry/util/condition.hpp>
|
||||
#include <boost/geometry/util/math.hpp>
|
||||
|
||||
#include <boost/geometry/formulas/differential_quantities.hpp>
|
||||
#include <boost/geometry/formulas/flattening.hpp>
|
||||
#include <boost/geometry/formulas/result_inverse.hpp>
|
||||
|
||||
|
||||
namespace boost { namespace geometry { namespace formula
|
||||
{
|
||||
|
||||
/*!
|
||||
\brief The solution of the inverse problem of geodesics on latlong coordinates,
|
||||
Forsyth-Andoyer-Lambert type approximation with second order terms.
|
||||
\author See
|
||||
- Technical Report: PAUL D. THOMAS, MATHEMATICAL MODELS FOR NAVIGATION SYSTEMS, 1965
|
||||
http://www.dtic.mil/docs/citations/AD0627893
|
||||
- Technical Report: PAUL D. THOMAS, SPHEROIDAL GEODESICS, REFERENCE SYSTEMS, AND LOCAL GEOMETRY, 1970
|
||||
http://www.dtic.mil/docs/citations/AD0703541
|
||||
*/
|
||||
template <
|
||||
typename CT,
|
||||
bool EnableDistance,
|
||||
bool EnableAzimuth,
|
||||
bool EnableReverseAzimuth = false,
|
||||
bool EnableReducedLength = false,
|
||||
bool EnableGeodesicScale = false
|
||||
>
|
||||
class thomas_inverse
|
||||
{
|
||||
static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale;
|
||||
static const bool CalcAzimuths = EnableAzimuth || EnableReverseAzimuth || CalcQuantities;
|
||||
static const bool CalcFwdAzimuth = EnableAzimuth || CalcQuantities;
|
||||
static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcQuantities;
|
||||
|
||||
public:
|
||||
typedef result_inverse<CT> result_type;
|
||||
|
||||
template <typename T1, typename T2, typename Spheroid>
|
||||
static inline result_type apply(T1 const& lon1,
|
||||
T1 const& lat1,
|
||||
T2 const& lon2,
|
||||
T2 const& lat2,
|
||||
Spheroid const& spheroid)
|
||||
{
|
||||
result_type result;
|
||||
|
||||
// coordinates in radians
|
||||
|
||||
if ( math::equals(lon1, lon2) && math::equals(lat1, lat2) )
|
||||
{
|
||||
return result;
|
||||
}
|
||||
|
||||
CT const c0 = 0;
|
||||
CT const c1 = 1;
|
||||
CT const c2 = 2;
|
||||
CT const c4 = 4;
|
||||
|
||||
CT const pi_half = math::pi<CT>() / c2;
|
||||
CT const f = formula::flattening<CT>(spheroid);
|
||||
CT const one_minus_f = c1 - f;
|
||||
|
||||
// CT const tan_theta1 = one_minus_f * tan(lat1);
|
||||
// CT const tan_theta2 = one_minus_f * tan(lat2);
|
||||
// CT const theta1 = atan(tan_theta1);
|
||||
// CT const theta2 = atan(tan_theta2);
|
||||
|
||||
CT const theta1 = math::equals(lat1, pi_half) ? lat1 :
|
||||
math::equals(lat1, -pi_half) ? lat1 :
|
||||
atan(one_minus_f * tan(lat1));
|
||||
CT const theta2 = math::equals(lat2, pi_half) ? lat2 :
|
||||
math::equals(lat2, -pi_half) ? lat2 :
|
||||
atan(one_minus_f * tan(lat2));
|
||||
|
||||
CT const theta_m = (theta1 + theta2) / c2;
|
||||
CT const d_theta_m = (theta2 - theta1) / c2;
|
||||
CT const d_lambda = lon2 - lon1;
|
||||
CT const d_lambda_m = d_lambda / c2;
|
||||
|
||||
CT const sin_theta_m = sin(theta_m);
|
||||
CT const cos_theta_m = cos(theta_m);
|
||||
CT const sin_d_theta_m = sin(d_theta_m);
|
||||
CT const cos_d_theta_m = cos(d_theta_m);
|
||||
CT const sin2_theta_m = math::sqr(sin_theta_m);
|
||||
CT const cos2_theta_m = math::sqr(cos_theta_m);
|
||||
CT const sin2_d_theta_m = math::sqr(sin_d_theta_m);
|
||||
CT const cos2_d_theta_m = math::sqr(cos_d_theta_m);
|
||||
CT const sin_d_lambda_m = sin(d_lambda_m);
|
||||
CT const sin2_d_lambda_m = math::sqr(sin_d_lambda_m);
|
||||
|
||||
CT const H = cos2_theta_m - sin2_d_theta_m;
|
||||
CT const L = sin2_d_theta_m + H * sin2_d_lambda_m;
|
||||
CT const cos_d = c1 - c2 * L;
|
||||
CT const d = acos(cos_d);
|
||||
CT const sin_d = sin(d);
|
||||
|
||||
CT const one_minus_L = c1 - L;
|
||||
|
||||
if ( math::equals(sin_d, c0)
|
||||
|| math::equals(L, c0)
|
||||
|| math::equals(one_minus_L, c0) )
|
||||
{
|
||||
return result;
|
||||
}
|
||||
|
||||
CT const U = c2 * sin2_theta_m * cos2_d_theta_m / one_minus_L;
|
||||
CT const V = c2 * sin2_d_theta_m * cos2_theta_m / L;
|
||||
CT const X = U + V;
|
||||
CT const Y = U - V;
|
||||
CT const T = d / sin_d;
|
||||
CT const D = c4 * math::sqr(T);
|
||||
CT const E = c2 * cos_d;
|
||||
CT const A = D * E;
|
||||
CT const B = c2 * D;
|
||||
CT const C = T - (A - E) / c2;
|
||||
|
||||
CT const f_sqr = math::sqr(f);
|
||||
CT const f_sqr_per_64 = f_sqr / CT(64);
|
||||
|
||||
if ( BOOST_GEOMETRY_CONDITION(EnableDistance) )
|
||||
{
|
||||
CT const n1 = X * (A + C*X);
|
||||
CT const n2 = Y * (B + E*Y);
|
||||
CT const n3 = D*X*Y;
|
||||
|
||||
CT const delta1d = f * (T*X-Y) / c4;
|
||||
CT const delta2d = f_sqr_per_64 * (n1 - n2 + n3);
|
||||
|
||||
CT const a = get_radius<0>(spheroid);
|
||||
|
||||
//result.distance = a * sin_d * (T - delta1d);
|
||||
result.distance = a * sin_d * (T - delta1d + delta2d);
|
||||
}
|
||||
|
||||
if ( BOOST_GEOMETRY_CONDITION(CalcAzimuths) )
|
||||
{
|
||||
// NOTE: if both cos_latX == 0 then below we'd have 0 * INF
|
||||
// it's a situation when the endpoints are on the poles +-90 deg
|
||||
// in this case the azimuth could either be 0 or +-pi
|
||||
// but above always 0 is returned
|
||||
|
||||
CT const F = c2*Y-E*(c4-X);
|
||||
CT const M = CT(32)*T-(CT(20)*T-A)*X-(B+c4)*Y;
|
||||
CT const G = f*T/c2 + f_sqr_per_64 * M;
|
||||
|
||||
// TODO:
|
||||
// If d_lambda is close to 90 or -90 deg then tan(d_lambda) is big
|
||||
// and F is small. The result is not accurate.
|
||||
// In the edge case the result may be 2 orders of magnitude less
|
||||
// accurate than Andoyer's.
|
||||
CT const tan_d_lambda = tan(d_lambda);
|
||||
CT const Q = -(F*G*tan_d_lambda) / c4;
|
||||
CT const d_lambda_m_p = (d_lambda + Q) / c2;
|
||||
CT const tan_d_lambda_m_p = tan(d_lambda_m_p);
|
||||
|
||||
CT const v = atan2(cos_d_theta_m, sin_theta_m * tan_d_lambda_m_p);
|
||||
CT const u = atan2(-sin_d_theta_m, cos_theta_m * tan_d_lambda_m_p);
|
||||
|
||||
CT const pi = math::pi<CT>();
|
||||
|
||||
if (BOOST_GEOMETRY_CONDITION(EnableAzimuth))
|
||||
{
|
||||
CT alpha1 = v + u;
|
||||
if (alpha1 > pi)
|
||||
{
|
||||
alpha1 -= c2 * pi;
|
||||
}
|
||||
|
||||
result.azimuth = alpha1;
|
||||
}
|
||||
|
||||
if (BOOST_GEOMETRY_CONDITION(EnableReverseAzimuth))
|
||||
{
|
||||
CT alpha2 = pi - (v - u);
|
||||
if (alpha2 > pi)
|
||||
{
|
||||
alpha2 -= c2 * pi;
|
||||
}
|
||||
|
||||
result.reverse_azimuth = alpha2;
|
||||
}
|
||||
}
|
||||
|
||||
if (BOOST_GEOMETRY_CONDITION(CalcQuantities))
|
||||
{
|
||||
typedef differential_quantities<CT, EnableReducedLength, EnableGeodesicScale, 2> quantities;
|
||||
quantities::apply(lon1, lat1, lon2, lat2,
|
||||
result.azimuth, result.reverse_azimuth,
|
||||
get_radius<2>(spheroid), f,
|
||||
result.reduced_length, result.geodesic_scale);
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
};
|
||||
|
||||
}}} // namespace boost::geometry::formula
|
||||
|
||||
|
||||
#endif // BOOST_GEOMETRY_FORMULAS_THOMAS_INVERSE_HPP
|
||||
@@ -0,0 +1,148 @@
|
||||
// Boost.Geometry
|
||||
|
||||
// Copyright (c) 2016-2017 Oracle and/or its affiliates.
|
||||
|
||||
// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
|
||||
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
|
||||
|
||||
// Use, modification and distribution is subject to the Boost Software License,
|
||||
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
|
||||
// http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#ifndef BOOST_GEOMETRY_FORMULAS_MAXIMUM_LATITUDE_HPP
|
||||
#define BOOST_GEOMETRY_FORMULAS_MAXIMUM_LATITUDE_HPP
|
||||
|
||||
#include <boost/geometry/core/srs.hpp>
|
||||
#include <boost/geometry/formulas/flattening.hpp>
|
||||
#include <boost/geometry/formulas/spherical.hpp>
|
||||
#include <boost/mpl/assert.hpp>
|
||||
|
||||
namespace boost { namespace geometry { namespace formula
|
||||
{
|
||||
|
||||
/*!
|
||||
\brief Algorithm to compute the vertex latitude of a geodesic segment. Vertex is
|
||||
a point on the geodesic that maximizes (or minimizes) the latitude.
|
||||
\author See
|
||||
[Wood96] Wood - Vertex Latitudes on Ellipsoid Geodesics, SIAM Rev., 38(4),
|
||||
637–644, 1996
|
||||
*/
|
||||
|
||||
template <typename CT>
|
||||
class vertex_latitude_on_sphere
|
||||
{
|
||||
|
||||
public:
|
||||
template<typename T1, typename T2>
|
||||
static inline CT apply(T1 const& lat1,
|
||||
T2 const& alp1)
|
||||
{
|
||||
return std::acos( math::abs(cos(lat1) * sin(alp1)) );
|
||||
}
|
||||
};
|
||||
|
||||
template <typename CT>
|
||||
class vertex_latitude_on_spheroid
|
||||
{
|
||||
|
||||
public:
|
||||
/*
|
||||
* formula based on paper
|
||||
* [Wood96] Wood - Vertex Latitudes on Ellipsoid Geodesics, SIAM Rev., 38(4),
|
||||
* 637–644, 1996
|
||||
template <typename T1, typename T2, typename Spheroid>
|
||||
static inline CT apply(T1 const& lat1,
|
||||
T2 const& alp1,
|
||||
Spheroid const& spheroid)
|
||||
{
|
||||
CT const f = formula::flattening<CT>(spheroid);
|
||||
|
||||
CT const e2 = f * (CT(2) - f);
|
||||
CT const sin_alp1 = sin(alp1);
|
||||
CT const sin2_lat1 = math::sqr(sin(lat1));
|
||||
CT const cos2_lat1 = CT(1) - sin2_lat1;
|
||||
|
||||
CT const e2_sin2 = CT(1) - e2 * sin2_lat1;
|
||||
CT const cos2_sin2 = cos2_lat1 * math::sqr(sin_alp1);
|
||||
CT const vertex_lat = std::asin( math::sqrt((e2_sin2 - cos2_sin2)
|
||||
/ (e2_sin2 - e2 * cos2_sin2)));
|
||||
return vertex_lat;
|
||||
}
|
||||
*/
|
||||
|
||||
// simpler formula based on Clairaut relation for spheroids
|
||||
template <typename T1, typename T2, typename Spheroid>
|
||||
static inline CT apply(T1 const& lat1,
|
||||
T2 const& alp1,
|
||||
Spheroid const& spheroid)
|
||||
{
|
||||
CT const f = formula::flattening<CT>(spheroid);
|
||||
|
||||
CT const one_minus_f = (CT(1) - f);
|
||||
|
||||
//get the reduced latitude
|
||||
CT const bet1 = atan( one_minus_f * tan(lat1) );
|
||||
|
||||
//apply Clairaut relation
|
||||
CT const betv = vertex_latitude_on_sphere<CT>::apply(bet1, alp1);
|
||||
|
||||
//return the spheroid latitude
|
||||
return atan( tan(betv) / one_minus_f );
|
||||
}
|
||||
|
||||
/*
|
||||
template <typename T>
|
||||
inline static void sign_adjustment(CT lat1, CT lat2, CT vertex_lat, T& vrt_result)
|
||||
{
|
||||
// signbit returns a non-zero value (true) if the sign is negative;
|
||||
// and zero (false) otherwise.
|
||||
bool sign = std::signbit(std::abs(lat1) > std::abs(lat2) ? lat1 : lat2);
|
||||
|
||||
vrt_result.north = sign ? std::max(lat1, lat2) : vertex_lat;
|
||||
vrt_result.south = sign ? vertex_lat * CT(-1) : std::min(lat1, lat2);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
inline static bool vertex_on_segment(CT alp1, CT alp2, CT lat1, CT lat2, T& vrt_result)
|
||||
{
|
||||
CT const half_pi = math::pi<CT>() / CT(2);
|
||||
|
||||
// if the segment does not contain the vertex of the geodesic
|
||||
// then return the endpoint of max (min) latitude
|
||||
if ((alp1 < half_pi && alp2 < half_pi)
|
||||
|| (alp1 > half_pi && alp2 > half_pi))
|
||||
{
|
||||
vrt_result.north = std::max(lat1, lat2);
|
||||
vrt_result.south = std::min(lat1, lat2);
|
||||
return false;
|
||||
}
|
||||
return true;
|
||||
}
|
||||
*/
|
||||
};
|
||||
|
||||
|
||||
template <typename CT, typename CS_Tag>
|
||||
struct vertex_latitude
|
||||
{
|
||||
BOOST_MPL_ASSERT_MSG
|
||||
(
|
||||
false, NOT_IMPLEMENTED_FOR_THIS_COORDINATE_SYSTEM, (types<CS_Tag>)
|
||||
);
|
||||
|
||||
};
|
||||
|
||||
template <typename CT>
|
||||
struct vertex_latitude<CT, spherical_equatorial_tag>
|
||||
: vertex_latitude_on_sphere<CT>
|
||||
{};
|
||||
|
||||
template <typename CT>
|
||||
struct vertex_latitude<CT, geographic_tag>
|
||||
: vertex_latitude_on_spheroid<CT>
|
||||
{};
|
||||
|
||||
|
||||
}}} // namespace boost::geometry::formula
|
||||
|
||||
#endif // BOOST_GEOMETRY_FORMULAS_MAXIMUM_LATITUDE_HPP
|
||||
@@ -0,0 +1,184 @@
|
||||
// Boost.Geometry
|
||||
|
||||
// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
|
||||
|
||||
// This file was modified by Oracle on 2014, 2016.
|
||||
// Modifications copyright (c) 2014-2016 Oracle and/or its affiliates.
|
||||
|
||||
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
|
||||
|
||||
// Use, modification and distribution is subject to the Boost Software License,
|
||||
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
|
||||
// http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#ifndef BOOST_GEOMETRY_FORMULAS_VINCENTY_DIRECT_HPP
|
||||
#define BOOST_GEOMETRY_FORMULAS_VINCENTY_DIRECT_HPP
|
||||
|
||||
|
||||
#include <boost/math/constants/constants.hpp>
|
||||
|
||||
#include <boost/geometry/core/radius.hpp>
|
||||
#include <boost/geometry/core/srs.hpp>
|
||||
|
||||
#include <boost/geometry/util/condition.hpp>
|
||||
#include <boost/geometry/util/math.hpp>
|
||||
|
||||
#include <boost/geometry/formulas/differential_quantities.hpp>
|
||||
#include <boost/geometry/formulas/flattening.hpp>
|
||||
#include <boost/geometry/formulas/result_direct.hpp>
|
||||
|
||||
|
||||
#ifndef BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS
|
||||
#define BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS 1000
|
||||
#endif
|
||||
|
||||
|
||||
namespace boost { namespace geometry { namespace formula
|
||||
{
|
||||
|
||||
/*!
|
||||
\brief The solution of the direct problem of geodesics on latlong coordinates, after Vincenty, 1975
|
||||
\author See
|
||||
- http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
|
||||
- http://www.icsm.gov.au/gda/gdav2.3.pdf
|
||||
\author Adapted from various implementations to get it close to the original document
|
||||
- http://www.movable-type.co.uk/scripts/LatLongVincenty.html
|
||||
- http://exogen.case.edu/projects/geopy/source/geopy.distance.html
|
||||
- http://futureboy.homeip.net/fsp/colorize.fsp?fileName=navigation.frink
|
||||
|
||||
*/
|
||||
template <
|
||||
typename CT,
|
||||
bool EnableCoordinates = true,
|
||||
bool EnableReverseAzimuth = false,
|
||||
bool EnableReducedLength = false,
|
||||
bool EnableGeodesicScale = false
|
||||
>
|
||||
class vincenty_direct
|
||||
{
|
||||
static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale;
|
||||
static const bool CalcCoordinates = EnableCoordinates || CalcQuantities;
|
||||
static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcQuantities;
|
||||
|
||||
public:
|
||||
typedef result_direct<CT> result_type;
|
||||
|
||||
template <typename T, typename Dist, typename Azi, typename Spheroid>
|
||||
static inline result_type apply(T const& lo1,
|
||||
T const& la1,
|
||||
Dist const& distance,
|
||||
Azi const& azimuth12,
|
||||
Spheroid const& spheroid)
|
||||
{
|
||||
result_type result;
|
||||
|
||||
CT const lon1 = lo1;
|
||||
CT const lat1 = la1;
|
||||
|
||||
if ( math::equals(distance, Dist(0)) || distance < Dist(0) )
|
||||
{
|
||||
result.lon2 = lon1;
|
||||
result.lat2 = lat1;
|
||||
return result;
|
||||
}
|
||||
|
||||
CT const radius_a = CT(get_radius<0>(spheroid));
|
||||
CT const radius_b = CT(get_radius<2>(spheroid));
|
||||
CT const flattening = formula::flattening<CT>(spheroid);
|
||||
|
||||
CT const sin_azimuth12 = sin(azimuth12);
|
||||
CT const cos_azimuth12 = cos(azimuth12);
|
||||
|
||||
// U: reduced latitude, defined by tan U = (1-f) tan phi
|
||||
CT const one_min_f = CT(1) - flattening;
|
||||
CT const tan_U1 = one_min_f * tan(lat1);
|
||||
CT const sigma1 = atan2(tan_U1, cos_azimuth12); // (1)
|
||||
|
||||
// may be calculated from tan using 1 sqrt()
|
||||
CT const U1 = atan(tan_U1);
|
||||
CT const sin_U1 = sin(U1);
|
||||
CT const cos_U1 = cos(U1);
|
||||
|
||||
CT const sin_alpha = cos_U1 * sin_azimuth12; // (2)
|
||||
CT const sin_alpha_sqr = math::sqr(sin_alpha);
|
||||
CT const cos_alpha_sqr = CT(1) - sin_alpha_sqr;
|
||||
|
||||
CT const b_sqr = radius_b * radius_b;
|
||||
CT const u_sqr = cos_alpha_sqr * (radius_a * radius_a - b_sqr) / b_sqr;
|
||||
CT const A = CT(1) + (u_sqr/CT(16384)) * (CT(4096) + u_sqr*(CT(-768) + u_sqr*(CT(320) - u_sqr*CT(175)))); // (3)
|
||||
CT const B = (u_sqr/CT(1024))*(CT(256) + u_sqr*(CT(-128) + u_sqr*(CT(74) - u_sqr*CT(47)))); // (4)
|
||||
|
||||
CT s_div_bA = distance / (radius_b * A);
|
||||
CT sigma = s_div_bA; // (7)
|
||||
|
||||
CT previous_sigma;
|
||||
CT sin_sigma;
|
||||
CT cos_sigma;
|
||||
CT cos_2sigma_m;
|
||||
CT cos_2sigma_m_sqr;
|
||||
|
||||
int counter = 0; // robustness
|
||||
|
||||
do
|
||||
{
|
||||
previous_sigma = sigma;
|
||||
|
||||
CT const two_sigma_m = CT(2) * sigma1 + sigma; // (5)
|
||||
|
||||
sin_sigma = sin(sigma);
|
||||
cos_sigma = cos(sigma);
|
||||
CT const sin_sigma_sqr = math::sqr(sin_sigma);
|
||||
cos_2sigma_m = cos(two_sigma_m);
|
||||
cos_2sigma_m_sqr = math::sqr(cos_2sigma_m);
|
||||
|
||||
CT const delta_sigma = B * sin_sigma * (cos_2sigma_m
|
||||
+ (B/CT(4)) * ( cos_sigma * (CT(-1) + CT(2)*cos_2sigma_m_sqr)
|
||||
- (B/CT(6) * cos_2sigma_m * (CT(-3)+CT(4)*sin_sigma_sqr) * (CT(-3)+CT(4)*cos_2sigma_m_sqr)) )); // (6)
|
||||
|
||||
sigma = s_div_bA + delta_sigma; // (7)
|
||||
|
||||
++counter; // robustness
|
||||
|
||||
} while ( geometry::math::abs(previous_sigma - sigma) > CT(1e-12)
|
||||
//&& geometry::math::abs(sigma) < pi
|
||||
&& counter < BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS ); // robustness
|
||||
|
||||
if (BOOST_GEOMETRY_CONDITION(CalcCoordinates))
|
||||
{
|
||||
result.lat2
|
||||
= atan2( sin_U1 * cos_sigma + cos_U1 * sin_sigma * cos_azimuth12,
|
||||
one_min_f * math::sqrt(sin_alpha_sqr + math::sqr(sin_U1 * sin_sigma - cos_U1 * cos_sigma * cos_azimuth12))); // (8)
|
||||
|
||||
CT const lambda = atan2( sin_sigma * sin_azimuth12,
|
||||
cos_U1 * cos_sigma - sin_U1 * sin_sigma * cos_azimuth12); // (9)
|
||||
CT const C = (flattening/CT(16)) * cos_alpha_sqr * ( CT(4) + flattening * ( CT(4) - CT(3) * cos_alpha_sqr ) ); // (10)
|
||||
CT const L = lambda - (CT(1) - C) * flattening * sin_alpha
|
||||
* ( sigma + C * sin_sigma * ( cos_2sigma_m + C * cos_sigma * ( CT(-1) + CT(2) * cos_2sigma_m_sqr ) ) ); // (11)
|
||||
|
||||
result.lon2 = lon1 + L;
|
||||
}
|
||||
|
||||
if (BOOST_GEOMETRY_CONDITION(CalcRevAzimuth))
|
||||
{
|
||||
result.reverse_azimuth
|
||||
= atan2(sin_alpha, -sin_U1 * sin_sigma + cos_U1 * cos_sigma * cos_azimuth12); // (12)
|
||||
}
|
||||
|
||||
if (BOOST_GEOMETRY_CONDITION(CalcQuantities))
|
||||
{
|
||||
typedef differential_quantities<CT, EnableReducedLength, EnableGeodesicScale, 2> quantities;
|
||||
quantities::apply(lon1, lat1, result.lon2, result.lat2,
|
||||
azimuth12, result.reverse_azimuth,
|
||||
radius_b, flattening,
|
||||
result.reduced_length, result.geodesic_scale);
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
}}} // namespace boost::geometry::formula
|
||||
|
||||
|
||||
#endif // BOOST_GEOMETRY_FORMULAS_VINCENTY_DIRECT_HPP
|
||||
@@ -0,0 +1,224 @@
|
||||
// Boost.Geometry
|
||||
|
||||
// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
|
||||
|
||||
// This file was modified by Oracle on 2014, 2016, 2017.
|
||||
// Modifications copyright (c) 2014-2017 Oracle and/or its affiliates.
|
||||
|
||||
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
|
||||
|
||||
// Use, modification and distribution is subject to the Boost Software License,
|
||||
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
|
||||
// http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#ifndef BOOST_GEOMETRY_FORMULAS_VINCENTY_INVERSE_HPP
|
||||
#define BOOST_GEOMETRY_FORMULAS_VINCENTY_INVERSE_HPP
|
||||
|
||||
|
||||
#include <boost/math/constants/constants.hpp>
|
||||
|
||||
#include <boost/geometry/core/radius.hpp>
|
||||
#include <boost/geometry/core/srs.hpp>
|
||||
|
||||
#include <boost/geometry/util/condition.hpp>
|
||||
#include <boost/geometry/util/math.hpp>
|
||||
|
||||
#include <boost/geometry/formulas/differential_quantities.hpp>
|
||||
#include <boost/geometry/formulas/flattening.hpp>
|
||||
#include <boost/geometry/formulas/result_inverse.hpp>
|
||||
|
||||
|
||||
#ifndef BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS
|
||||
#define BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS 1000
|
||||
#endif
|
||||
|
||||
|
||||
namespace boost { namespace geometry { namespace formula
|
||||
{
|
||||
|
||||
/*!
|
||||
\brief The solution of the inverse problem of geodesics on latlong coordinates, after Vincenty, 1975
|
||||
\author See
|
||||
- http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
|
||||
- http://www.icsm.gov.au/gda/gda-v_2.4.pdf
|
||||
\author Adapted from various implementations to get it close to the original document
|
||||
- http://www.movable-type.co.uk/scripts/LatLongVincenty.html
|
||||
- http://exogen.case.edu/projects/geopy/source/geopy.distance.html
|
||||
- http://futureboy.homeip.net/fsp/colorize.fsp?fileName=navigation.frink
|
||||
|
||||
*/
|
||||
template <
|
||||
typename CT,
|
||||
bool EnableDistance,
|
||||
bool EnableAzimuth,
|
||||
bool EnableReverseAzimuth = false,
|
||||
bool EnableReducedLength = false,
|
||||
bool EnableGeodesicScale = false
|
||||
>
|
||||
struct vincenty_inverse
|
||||
{
|
||||
static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale;
|
||||
static const bool CalcAzimuths = EnableAzimuth || EnableReverseAzimuth || CalcQuantities;
|
||||
static const bool CalcFwdAzimuth = EnableAzimuth || CalcQuantities;
|
||||
static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcQuantities;
|
||||
|
||||
public:
|
||||
typedef result_inverse<CT> result_type;
|
||||
|
||||
template <typename T1, typename T2, typename Spheroid>
|
||||
static inline result_type apply(T1 const& lon1,
|
||||
T1 const& lat1,
|
||||
T2 const& lon2,
|
||||
T2 const& lat2,
|
||||
Spheroid const& spheroid)
|
||||
{
|
||||
result_type result;
|
||||
|
||||
if (math::equals(lat1, lat2) && math::equals(lon1, lon2))
|
||||
{
|
||||
return result;
|
||||
}
|
||||
|
||||
CT const c1 = 1;
|
||||
CT const c2 = 2;
|
||||
CT const c3 = 3;
|
||||
CT const c4 = 4;
|
||||
CT const c16 = 16;
|
||||
CT const c_e_12 = CT(1e-12);
|
||||
|
||||
CT const pi = geometry::math::pi<CT>();
|
||||
CT const two_pi = c2 * pi;
|
||||
|
||||
// lambda: difference in longitude on an auxiliary sphere
|
||||
CT L = lon2 - lon1;
|
||||
CT lambda = L;
|
||||
|
||||
if (L < -pi) L += two_pi;
|
||||
if (L > pi) L -= two_pi;
|
||||
|
||||
CT const radius_a = CT(get_radius<0>(spheroid));
|
||||
CT const radius_b = CT(get_radius<2>(spheroid));
|
||||
CT const f = formula::flattening<CT>(spheroid);
|
||||
|
||||
// U: reduced latitude, defined by tan U = (1-f) tan phi
|
||||
CT const one_min_f = c1 - f;
|
||||
CT const tan_U1 = one_min_f * tan(lat1); // above (1)
|
||||
CT const tan_U2 = one_min_f * tan(lat2); // above (1)
|
||||
|
||||
// calculate sin U and cos U using trigonometric identities
|
||||
CT const temp_den_U1 = math::sqrt(c1 + math::sqr(tan_U1));
|
||||
CT const temp_den_U2 = math::sqrt(c1 + math::sqr(tan_U2));
|
||||
// cos = 1 / sqrt(1 + tan^2)
|
||||
CT const cos_U1 = c1 / temp_den_U1;
|
||||
CT const cos_U2 = c1 / temp_den_U2;
|
||||
// sin = tan / sqrt(1 + tan^2)
|
||||
// sin = tan * cos
|
||||
CT const sin_U1 = tan_U1 * cos_U1;
|
||||
CT const sin_U2 = tan_U2 * cos_U2;
|
||||
|
||||
// calculate sin U and cos U directly
|
||||
//CT const U1 = atan(tan_U1);
|
||||
//CT const U2 = atan(tan_U2);
|
||||
//cos_U1 = cos(U1);
|
||||
//cos_U2 = cos(U2);
|
||||
//sin_U1 = tan_U1 * cos_U1; // sin(U1);
|
||||
//sin_U2 = tan_U2 * cos_U2; // sin(U2);
|
||||
|
||||
CT previous_lambda;
|
||||
CT sin_lambda;
|
||||
CT cos_lambda;
|
||||
CT sin_sigma;
|
||||
CT sin_alpha;
|
||||
CT cos2_alpha;
|
||||
CT cos_2sigma_m;
|
||||
CT cos2_2sigma_m;
|
||||
CT sigma;
|
||||
|
||||
int counter = 0; // robustness
|
||||
|
||||
do
|
||||
{
|
||||
previous_lambda = lambda; // (13)
|
||||
sin_lambda = sin(lambda);
|
||||
cos_lambda = cos(lambda);
|
||||
sin_sigma = math::sqrt(math::sqr(cos_U2 * sin_lambda) + math::sqr(cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_lambda)); // (14)
|
||||
CT cos_sigma = sin_U1 * sin_U2 + cos_U1 * cos_U2 * cos_lambda; // (15)
|
||||
sin_alpha = cos_U1 * cos_U2 * sin_lambda / sin_sigma; // (17)
|
||||
cos2_alpha = c1 - math::sqr(sin_alpha);
|
||||
cos_2sigma_m = math::equals(cos2_alpha, 0) ? 0 : cos_sigma - c2 * sin_U1 * sin_U2 / cos2_alpha; // (18)
|
||||
cos2_2sigma_m = math::sqr(cos_2sigma_m);
|
||||
|
||||
CT C = f/c16 * cos2_alpha * (c4 + f * (c4 - c3 * cos2_alpha)); // (10)
|
||||
sigma = atan2(sin_sigma, cos_sigma); // (16)
|
||||
lambda = L + (c1 - C) * f * sin_alpha *
|
||||
(sigma + C * sin_sigma * (cos_2sigma_m + C * cos_sigma * (-c1 + c2 * cos2_2sigma_m))); // (11)
|
||||
|
||||
++counter; // robustness
|
||||
|
||||
} while ( geometry::math::abs(previous_lambda - lambda) > c_e_12
|
||||
&& geometry::math::abs(lambda) < pi
|
||||
&& counter < BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS ); // robustness
|
||||
|
||||
if ( BOOST_GEOMETRY_CONDITION(EnableDistance) )
|
||||
{
|
||||
// Oops getting hard here
|
||||
// (again, problem is that ttmath cannot divide by doubles, which is OK)
|
||||
CT const c1 = 1;
|
||||
CT const c2 = 2;
|
||||
CT const c3 = 3;
|
||||
CT const c4 = 4;
|
||||
CT const c6 = 6;
|
||||
CT const c47 = 47;
|
||||
CT const c74 = 74;
|
||||
CT const c128 = 128;
|
||||
CT const c256 = 256;
|
||||
CT const c175 = 175;
|
||||
CT const c320 = 320;
|
||||
CT const c768 = 768;
|
||||
CT const c1024 = 1024;
|
||||
CT const c4096 = 4096;
|
||||
CT const c16384 = 16384;
|
||||
|
||||
//CT sqr_u = cos2_alpha * (math::sqr(radius_a) - math::sqr(radius_b)) / math::sqr(radius_b); // above (1)
|
||||
CT sqr_u = cos2_alpha * ( math::sqr(radius_a / radius_b) - c1 ); // above (1)
|
||||
|
||||
CT A = c1 + sqr_u/c16384 * (c4096 + sqr_u * (-c768 + sqr_u * (c320 - c175 * sqr_u))); // (3)
|
||||
CT B = sqr_u/c1024 * (c256 + sqr_u * ( -c128 + sqr_u * (c74 - c47 * sqr_u))); // (4)
|
||||
CT const cos_sigma = cos(sigma);
|
||||
CT const sin2_sigma = math::sqr(sin_sigma);
|
||||
CT delta_sigma = B * sin_sigma * (cos_2sigma_m + (B/c4) * (cos_sigma* (-c1 + c2 * cos2_2sigma_m)
|
||||
- (B/c6) * cos_2sigma_m * (-c3 + c4 * sin2_sigma) * (-c3 + c4 * cos2_2sigma_m))); // (6)
|
||||
|
||||
result.distance = radius_b * A * (sigma - delta_sigma); // (19)
|
||||
}
|
||||
|
||||
if ( BOOST_GEOMETRY_CONDITION(CalcAzimuths) )
|
||||
{
|
||||
if (BOOST_GEOMETRY_CONDITION(CalcFwdAzimuth))
|
||||
{
|
||||
result.azimuth = atan2(cos_U2 * sin_lambda, cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_lambda); // (20)
|
||||
}
|
||||
|
||||
if (BOOST_GEOMETRY_CONDITION(CalcRevAzimuth))
|
||||
{
|
||||
result.reverse_azimuth = atan2(cos_U1 * sin_lambda, -sin_U1 * cos_U2 + cos_U1 * sin_U2 * cos_lambda); // (21)
|
||||
}
|
||||
}
|
||||
|
||||
if (BOOST_GEOMETRY_CONDITION(CalcQuantities))
|
||||
{
|
||||
typedef differential_quantities<CT, EnableReducedLength, EnableGeodesicScale, 2> quantities;
|
||||
quantities::apply(lon1, lat1, lon2, lat2,
|
||||
result.azimuth, result.reverse_azimuth,
|
||||
radius_b, f,
|
||||
result.reduced_length, result.geodesic_scale);
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
};
|
||||
|
||||
}}} // namespace boost::geometry::formula
|
||||
|
||||
|
||||
#endif // BOOST_GEOMETRY_FORMULAS_VINCENTY_INVERSE_HPP
|
||||
Reference in New Issue
Block a user