[ceph.git] / ceph / src / boost / boost / geometry / formulas / vertex_longitude.hpp

// Boost.Geometry

// Copyright (c) 2016-2018 Oracle and/or its affiliates.

// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle

// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_GEOMETRY_FORMULAS_MAXIMUM_LONGITUDE_HPP
#define BOOST_GEOMETRY_FORMULAS_MAXIMUM_LONGITUDE_HPP

#include <boost/geometry/formulas/spherical.hpp>
#include <boost/geometry/formulas/flattening.hpp>

#include <boost/mpl/assert.hpp>

#include <boost/math/special_functions/hypot.hpp>

namespace boost { namespace geometry { namespace formula
{

/*!
\brief Algorithm to compute the vertex longitude of a geodesic segment. Vertex is
a point on the geodesic that maximizes (or minimizes) the latitude. The algorithm
is given the vertex latitude.
*/

//Classes for spesific CS

template <typename CT>
class vertex_longitude_on_sphere
{

public:

    template <typename T>
    static inline CT apply(T const& lat1, //segment point 1
                           T const& lat2, //segment point 2
                           T const& lat3, //vertex latitude
                           T const& sin_l12,
                           T const& cos_l12) //lon1 -lon2
    {
        //https://en.wikipedia.org/wiki/Great-circle_navigation#Finding_way-points
        CT const A = sin(lat1) * cos(lat2) * cos(lat3) * sin_l12;
        CT const B = sin(lat1) * cos(lat2) * cos(lat3) * cos_l12
                - cos(lat1) * sin(lat2) * cos(lat3);
        CT lon = atan2(B, A);
        return lon + math::pi<CT>();
    }
};

template <typename CT>
class vertex_longitude_on_spheroid
{
    template<typename T>
    static inline void normalize(T& x, T& y)
    {
        T h = boost::math::hypot(x, y);
        x /= h;
        y /= h;
    }

public:

    template <typename T, typename Spheroid>
    static inline CT apply(T const& lat1, //segment point 1
                           T const& lat2, //segment point 2
                           T const& lat3, //vertex latitude
                           T& alp1,
                           Spheroid const& spheroid)
    {
        // We assume that segment points lay on different side w.r.t.
        // the vertex

        // Constants
        CT const c0 = 0;
        CT const c2 = 2;
        CT const half_pi = math::pi<CT>() / c2;
        if (math::equals(lat1, half_pi)
                || math::equals(lat2, half_pi)
                || math::equals(lat1, -half_pi)
                || math::equals(lat2, -half_pi))
        {
            // one segment point is the pole
            return c0;
        }

        // More constants
        CT const f = flattening<CT>(spheroid);
        CT const pi = math::pi<CT>();
        CT const c1 = 1;
        CT const cminus1 = -1;

        // First, compute longitude on auxiliary sphere

        CT const one_minus_f = c1 - f;
        CT const bet1 = atan(one_minus_f * tan(lat1));
        CT const bet2 = atan(one_minus_f * tan(lat2));
        CT const bet3 = atan(one_minus_f * tan(lat3));

        CT cos_bet1 = cos(bet1);
        CT cos_bet2 = cos(bet2);
        CT const sin_bet1 = sin(bet1);
        CT const sin_bet2 = sin(bet2);
        CT const sin_bet3 = sin(bet3);

        CT omg12 = 0;

        if (bet1 < c0)
        {
            cos_bet1 *= cminus1;
            omg12 += pi;
        }
        if (bet2 < c0)
        {
            cos_bet2 *= cminus1;
            omg12 += pi;
        }

        CT const sin_alp1 = sin(alp1);
        CT const cos_alp1 = math::sqrt(c1 - math::sqr(sin_alp1));

        CT const norm = math::sqrt(math::sqr(cos_alp1) + math::sqr(sin_alp1 * sin_bet1));
        CT const sin_alp0 = sin(atan2(sin_alp1 * cos_bet1, norm));

        BOOST_ASSERT(cos_bet2 != c0);
        CT const sin_alp2 = sin_alp1 * cos_bet1 / cos_bet2;

        CT const cos_alp0 = math::sqrt(c1 - math::sqr(sin_alp0));
        CT const cos_alp2 = math::sqrt(c1 - math::sqr(sin_alp2));

        CT const sig1 = atan2(sin_bet1, cos_alp1 * cos_bet1);
        CT const sig2 = atan2(sin_bet2, -cos_alp2 * cos_bet2); //lat3 is a vertex

        CT const cos_sig1 = cos(sig1);
        CT const sin_sig1 = math::sqrt(c1 - math::sqr(cos_sig1));

        CT const cos_sig2 = cos(sig2);
        CT const sin_sig2 = math::sqrt(c1 - math::sqr(cos_sig2));

        CT const omg1 = atan2(sin_alp0 * sin_sig1, cos_sig1);
        CT const omg2 = atan2(sin_alp0 * sin_sig2, cos_sig2);

        omg12 += omg1 - omg2;

        CT const sin_omg12 = sin(omg12);
        CT const cos_omg12 = cos(omg12);

        CT omg13 = geometry::formula::vertex_longitude_on_sphere<CT>
                ::apply(bet1, bet2, bet3, sin_omg12, cos_omg12);

        if (lat1 * lat2 < c0)//different hemispheres
        {
            if ((lat2 - lat1) * lat3  > c0)// ascending segment
            {
                omg13 = pi - omg13;
            }
        }

        // Second, compute the ellipsoidal longitude

        CT const e2 = f * (c2 - f);
        CT const ep = math::sqrt(e2 / (c1 - e2));
        CT const k2 = math::sqr(ep * cos_alp0);
        CT const sqrt_k2_plus_one = math::sqrt(c1 + k2);
        CT const eps = (sqrt_k2_plus_one - c1) / (sqrt_k2_plus_one + c1);
        CT const eps2 = eps * eps;
        CT const n = f / (c2 - f);

        // sig3 is the length from equator to the vertex
        CT sig3;
        if(sin_bet3 > c0)
        {
            sig3 = half_pi;
        } else {
            sig3 = -half_pi;
        }
        CT const cos_sig3 = 0;
        CT const sin_sig3 = 1;

        CT sig13 = sig3 - sig1;
        if (sig13 > pi)
        {
            sig13 -= 2 * pi;
        }

        // Order 2 approximation
        CT const c1over2 = 0.5;
        CT const c1over4 = 0.25;
        CT const c1over8 = 0.125;
        CT const c1over16 = 0.0625;
        CT const c4 = 4;
        CT const c8 = 8;

        CT const A3 = 1 - (c1over2 - c1over2 * n) * eps - c1over4 * eps2;
        CT const C31 = (c1over4 - c1over4 * n) * eps + c1over8 * eps2;
        CT const C32 = c1over16 * eps2;

        CT const sin2_sig3 = c2 * cos_sig3 * sin_sig3;
        CT const sin4_sig3 = sin_sig3 * (-c4 * cos_sig3
                                         + c8 * cos_sig3 * cos_sig3 * cos_sig3);
        CT const sin2_sig1 = c2 * cos_sig1 * sin_sig1;
        CT const sin4_sig1 = sin_sig1 * (-c4 * cos_sig1
                                         + c8 * cos_sig1 * cos_sig1 * cos_sig1);
        CT const I3 = A3 * (sig13
                            + C31 * (sin2_sig3 - sin2_sig1)
                            + C32 * (sin4_sig3 - sin4_sig1));

        CT const sign = bet3 >= c0
                      ? c1
                      : cminus1;
        
        CT const dlon_max = omg13 - sign * f * sin_alp0 * I3;

        return dlon_max;
    }
};

//CS_tag dispatching

template <typename CT, typename CS_Tag>
struct compute_vertex_lon
{
    BOOST_MPL_ASSERT_MSG
    (
        false, NOT_IMPLEMENTED_FOR_THIS_COORDINATE_SYSTEM, (types<CS_Tag>)
    );

};

template <typename CT>
struct compute_vertex_lon<CT, spherical_equatorial_tag>
{
    template <typename Strategy>
    static inline CT apply(CT const& lat1,
                           CT const& lat2,
                           CT const& vertex_lat,
                           CT const& sin_l12,
                           CT const& cos_l12,
                           CT,
                           Strategy)
    {
        return vertex_longitude_on_sphere<CT>
                ::apply(lat1,
                        lat2,
                        vertex_lat,
                        sin_l12,
                        cos_l12);
    }
};

template <typename CT>
struct compute_vertex_lon<CT, geographic_tag>
{
    template <typename Strategy>
    static inline CT apply(CT const& lat1,
                           CT const& lat2,
                           CT const& vertex_lat,
                           CT,
                           CT,
                           CT& alp1,
                           Strategy const& azimuth_strategy)
    {
        return vertex_longitude_on_spheroid<CT>
                ::apply(lat1,
                        lat2,
                        vertex_lat,
                        alp1,
                        azimuth_strategy.model());
    }
};

// Vertex longitude interface
// Assume that lon1 < lon2 and vertex_lat is the latitude of the vertex

template <typename CT, typename CS_Tag>
class vertex_longitude
{
public :
    template <typename Strategy>
    static inline CT apply(CT& lon1,
                           CT& lat1,
                           CT& lon2,
                           CT& lat2,
                           CT const& vertex_lat,
                           CT& alp1,
                           Strategy const& azimuth_strategy)
    {
        CT const c0 = 0;
        CT pi = math::pi<CT>();

        //Vertex is a segment's point
        if (math::equals(vertex_lat, lat1))
        {
            return lon1;
        }
        if (math::equals(vertex_lat, lat2))
        {
            return lon2;
        }

        //Segment lay on meridian
        if (math::equals(lon1, lon2))
        {
            return (std::max)(lat1, lat2);
        }
        BOOST_ASSERT(lon1 < lon2);

        CT dlon = compute_vertex_lon<CT, CS_Tag>::apply(lat1, lat2,
                                                        vertex_lat,
                                                        sin(lon1 - lon2),
                                                        cos(lon1 - lon2),
                                                        alp1,
                                                        azimuth_strategy);

        CT vertex_lon = std::fmod(lon1 + dlon, 2 * pi);

        if (vertex_lat < c0)
        {
            vertex_lon -= pi;
        }

        if (std::abs(lon1 - lon2) > pi)
        {
            vertex_lon -= pi;
        }

        return vertex_lon;
    }
};

template <typename CT>
class vertex_longitude<CT, cartesian_tag>
{
public :
    template <typename Strategy>
    static inline CT apply(CT& /*lon1*/,
                           CT& /*lat1*/,
                           CT& lon2,
                           CT& /*lat2*/,
                           CT const& /*vertex_lat*/,
                           CT& /*alp1*/,
                           Strategy const& /*azimuth_strategy*/)
    {
        return lon2;
    }
};

}}} // namespace boost::geometry::formula
#endif // BOOST_GEOMETRY_FORMULAS_MAXIMUM_LONGITUDE_HPP
Commit	Line	Data
b32b8144 FG	1	// Boost.Geometry
b32b8144 FG	2
92f5a8d4	3	// Copyright (c) 2016-2018 Oracle and/or its affiliates.
b32b8144 FG	4
	5	// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
	6	// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
	7
	8	// Use, modification and distribution is subject to the Boost Software License,
	9	// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
	10	// http://www.boost.org/LICENSE_1_0.txt)
	11
	12	#ifndef BOOST_GEOMETRY_FORMULAS_MAXIMUM_LONGITUDE_HPP
	13	#define BOOST_GEOMETRY_FORMULAS_MAXIMUM_LONGITUDE_HPP
	14
	15	#include <boost/geometry/formulas/spherical.hpp>
	16	#include <boost/geometry/formulas/flattening.hpp>
11fdf7f2	17
b32b8144 FG	18	#include <boost/mpl/assert.hpp>
	19
	20	#include <boost/math/special_functions/hypot.hpp>
	21
	22	namespace boost { namespace geometry { namespace formula
	23	{
	24
	25	/*!
	26	\brief Algorithm to compute the vertex longitude of a geodesic segment. Vertex is
	27	a point on the geodesic that maximizes (or minimizes) the latitude. The algorithm
	28	is given the vertex latitude.
	29	*/
	30
	31	//Classes for spesific CS
	32
	33	template <typename CT>
	34	class vertex_longitude_on_sphere
	35	{
	36
	37	public:
	38
	39	template <typename T>
	40	static inline CT apply(T const& lat1, //segment point 1
	41	T const& lat2, //segment point 2
	42	T const& lat3, //vertex latitude
	43	T const& sin_l12,
	44	T const& cos_l12) //lon1 -lon2
	45	{
	46	//https://en.wikipedia.org/wiki/Great-circle_navigation#Finding_way-points
	47	CT const A = sin(lat1) * cos(lat2) * cos(lat3) * sin_l12;
	48	CT const B = sin(lat1) * cos(lat2) * cos(lat3) * cos_l12
	49	- cos(lat1) * sin(lat2) * cos(lat3);
	50	CT lon = atan2(B, A);
	51	return lon + math::pi<CT>();
	52	}
	53	};
	54
	55	template <typename CT>
	56	class vertex_longitude_on_spheroid
	57	{
	58	template<typename T>
	59	static inline void normalize(T& x, T& y)
	60	{
	61	T h = boost::math::hypot(x, y);
	62	x /= h;
	63	y /= h;
	64	}
	65
	66	public:
	67
	68	template <typename T, typename Spheroid>
	69	static inline CT apply(T const& lat1, //segment point 1
	70	T const& lat2, //segment point 2
	71	T const& lat3, //vertex latitude
	72	T& alp1,
	73	Spheroid const& spheroid)
	74	{
	75	// We assume that segment points lay on different side w.r.t.
	76	// the vertex
	77
	78	// Constants
	79	CT const c0 = 0;
	80	CT const c2 = 2;
	81	CT const half_pi = math::pi<CT>() / c2;
82	if (math::equals(lat1, half_pi)
83	\|\| math::equals(lat2, half_pi)
84	\|\| math::equals(lat1, -half_pi)
85	\|\| math::equals(lat2, -half_pi))
86	{
87	// one segment point is the pole
88	return c0;
89	}
90
91	// More constants
92	CT const f = flattening<CT>(spheroid);
93	CT const pi = math::pi<CT>();
94	CT const c1 = 1;
95	CT const cminus1 = -1;
96
97	// First, compute longitude on auxiliary sphere
98
99	CT const one_minus_f = c1 - f;
100	CT const bet1 = atan(one_minus_f * tan(lat1));
101	CT const bet2 = atan(one_minus_f * tan(lat2));
102	CT const bet3 = atan(one_minus_f * tan(lat3));
103
104	CT cos_bet1 = cos(bet1);
105	CT cos_bet2 = cos(bet2);
106	CT const sin_bet1 = sin(bet1);
107	CT const sin_bet2 = sin(bet2);
108	CT const sin_bet3 = sin(bet3);
109
110	CT omg12 = 0;
111
112	if (bet1 < c0)
113	{
114	cos_bet1 *= cminus1;
115	omg12 += pi;
116	}
117	if (bet2 < c0)
118	{
119	cos_bet2 *= cminus1;
120	omg12 += pi;
121	}
122
123	CT const sin_alp1 = sin(alp1);
124	CT const cos_alp1 = math::sqrt(c1 - math::sqr(sin_alp1));
125
126	CT const norm = math::sqrt(math::sqr(cos_alp1) + math::sqr(sin_alp1 * sin_bet1));
127	CT const sin_alp0 = sin(atan2(sin_alp1 * cos_bet1, norm));
128
129	BOOST_ASSERT(cos_bet2 != c0);
130	CT const sin_alp2 = sin_alp1 * cos_bet1 / cos_bet2;
131
132	CT const cos_alp0 = math::sqrt(c1 - math::sqr(sin_alp0));
133	CT const cos_alp2 = math::sqrt(c1 - math::sqr(sin_alp2));
134
135	CT const sig1 = atan2(sin_bet1, cos_alp1 * cos_bet1);
136	CT const sig2 = atan2(sin_bet2, -cos_alp2 * cos_bet2); //lat3 is a vertex
137
138	CT const cos_sig1 = cos(sig1);
139	CT const sin_sig1 = math::sqrt(c1 - math::sqr(cos_sig1));
140
141	CT const cos_sig2 = cos(sig2);
142	CT const sin_sig2 = math::sqrt(c1 - math::sqr(cos_sig2));
143
144	CT const omg1 = atan2(sin_alp0 * sin_sig1, cos_sig1);
145	CT const omg2 = atan2(sin_alp0 * sin_sig2, cos_sig2);
146
147	omg12 += omg1 - omg2;
148
149	CT const sin_omg12 = sin(omg12);
150	CT const cos_omg12 = cos(omg12);
151
152	CT omg13 = geometry::formula::vertex_longitude_on_sphere<CT>
153	::apply(bet1, bet2, bet3, sin_omg12, cos_omg12);
154
155	if (lat1 * lat2 < c0)//different hemispheres
156	{
157	if ((lat2 - lat1) * lat3 > c0)// ascending segment
158	{
159	omg13 = pi - omg13;
160	}
161	}
162
163	// Second, compute the ellipsoidal longitude
164
165	CT const e2 = f * (c2 - f);
166	CT const ep = math::sqrt(e2 / (c1 - e2));
167	CT const k2 = math::sqr(ep * cos_alp0);
168	CT const sqrt_k2_plus_one = math::sqrt(c1 + k2);
169	CT const eps = (sqrt_k2_plus_one - c1) / (sqrt_k2_plus_one + c1);
170	CT const eps2 = eps * eps;
171	CT const n = f / (c2 - f);
172
173	// sig3 is the length from equator to the vertex
174	CT sig3;
175	if(sin_bet3 > c0)
176	{
177	sig3 = half_pi;
178	} else {
179	sig3 = -half_pi;
180	}
181	CT const cos_sig3 = 0;
182	CT const sin_sig3 = 1;
183
184	CT sig13 = sig3 - sig1;
185	if (sig13 > pi)
186	{
187	sig13 -= 2 * pi;
188	}
189
190	// Order 2 approximation
191	CT const c1over2 = 0.5;
192	CT const c1over4 = 0.25;
193	CT const c1over8 = 0.125;
194	CT const c1over16 = 0.0625;
195	CT const c4 = 4;
196	CT const c8 = 8;
197
198	CT const A3 = 1 - (c1over2 - c1over2 * n) * eps - c1over4 * eps2;
199	CT const C31 = (c1over4 - c1over4 * n) * eps + c1over8 * eps2;
200	CT const C32 = c1over16 * eps2;
201
202	CT const sin2_sig3 = c2 * cos_sig3 * sin_sig3;
203	CT const sin4_sig3 = sin_sig3 * (-c4 * cos_sig3
204	+ c8 * cos_sig3 * cos_sig3 * cos_sig3);
205	CT const sin2_sig1 = c2 * cos_sig1 * sin_sig1;
206	CT const sin4_sig1 = sin_sig1 * (-c4 * cos_sig1
207	+ c8 * cos_sig1 * cos_sig1 * cos_sig1);
208	CT const I3 = A3 * (sig13
209	+ C31 * (sin2_sig3 - sin2_sig1)
210	+ C32 * (sin4_sig3 - sin4_sig1));
211
11fdf7f2 TL	212	CT const sign = bet3 >= c0
	213	? c1
	214	: cminus1;
	215
b32b8144 FG	216	CT const dlon_max = omg13 - sign * f * sin_alp0 * I3;
	217
	218	return dlon_max;
	219	}
	220	};
	221
	222	//CS_tag dispatching
	223
	224	template <typename CT, typename CS_Tag>
	225	struct compute_vertex_lon
	226	{
	227	BOOST_MPL_ASSERT_MSG
	228	(
	229	false, NOT_IMPLEMENTED_FOR_THIS_COORDINATE_SYSTEM, (types<CS_Tag>)
	230	);
	231
	232	};
	233
	234	template <typename CT>
	235	struct compute_vertex_lon<CT, spherical_equatorial_tag>
	236	{
	237	template <typename Strategy>
	238	static inline CT apply(CT const& lat1,
	239	CT const& lat2,
	240	CT const& vertex_lat,
	241	CT const& sin_l12,
	242	CT const& cos_l12,
	243	CT,
	244	Strategy)
	245	{
	246	return vertex_longitude_on_sphere<CT>
	247	::apply(lat1,
	248	lat2,
	249	vertex_lat,
	250	sin_l12,
	251	cos_l12);
	252	}
	253	};
	254
	255	template <typename CT>
	256	struct compute_vertex_lon<CT, geographic_tag>
	257	{
	258	template <typename Strategy>
	259	static inline CT apply(CT const& lat1,
	260	CT const& lat2,
	261	CT const& vertex_lat,
	262	CT,
	263	CT,
	264	CT& alp1,
	265	Strategy const& azimuth_strategy)
	266	{
	267	return vertex_longitude_on_spheroid<CT>
	268	::apply(lat1,
	269	lat2,
	270	vertex_lat,
	271	alp1,
	272	azimuth_strategy.model());
	273	}
	274	};
	275
	276	// Vertex longitude interface
	277	// Assume that lon1 < lon2 and vertex_lat is the latitude of the vertex
	278
	279	template <typename CT, typename CS_Tag>
280	class vertex_longitude
281	{
282	public :
283	template <typename Strategy>
284	static inline CT apply(CT& lon1,
285	CT& lat1,
286	CT& lon2,
287	CT& lat2,
288	CT const& vertex_lat,
289	CT& alp1,
290	Strategy const& azimuth_strategy)
291	{
292	CT const c0 = 0;
293	CT pi = math::pi<CT>();
294
295	//Vertex is a segment's point
296	if (math::equals(vertex_lat, lat1))
297	{
298	return lon1;
299	}
300	if (math::equals(vertex_lat, lat2))
301	{
302	return lon2;
303	}
304
305	//Segment lay on meridian
306	if (math::equals(lon1, lon2))
307	{
308	return (std::max)(lat1, lat2);
309	}
310	BOOST_ASSERT(lon1 < lon2);
311
312	CT dlon = compute_vertex_lon<CT, CS_Tag>::apply(lat1, lat2,
313	vertex_lat,
314	sin(lon1 - lon2),
315	cos(lon1 - lon2),
316	alp1,
317	azimuth_strategy);
318
319	CT vertex_lon = std::fmod(lon1 + dlon, 2 * pi);
320
321	if (vertex_lat < c0)
322	{
323	vertex_lon -= pi;
324	}
325
326	if (std::abs(lon1 - lon2) > pi)
327	{
328	vertex_lon -= pi;
329	}
330
331	return vertex_lon;
332	}
333	};
334
92f5a8d4 TL	335	template <typename CT>
	336	class vertex_longitude<CT, cartesian_tag>
	337	{
	338	public :
	339	template <typename Strategy>
	340	static inline CT apply(CT& /lon1/,
	341	CT& /lat1/,
	342	CT& lon2,
	343	CT& /lat2/,
	344	CT const& /vertex_lat/,
	345	CT& /alp1/,
	346	Strategy const& /azimuth_strategy/)
	347	{
	348	return lon2;
	349	}
	350	};
	351
b32b8144 FG	352	}}} // namespace boost::geometry::formula
	353	#endif // BOOST_GEOMETRY_FORMULAS_MAXIMUM_LONGITUDE_HPP
	354