1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
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
|
