]> git.proxmox.com Git - ceph.git/blame - ceph/src/boost/boost/geometry/srs/projections/proj/aitoff.hpp
import new upstream nautilus stable release 14.2.8
[ceph.git] / ceph / src / boost / boost / geometry / srs / projections / proj / aitoff.hpp
CommitLineData
92f5a8d4 1// Boost.Geometry - gis-projections (based on PROJ4)
11fdf7f2
TL
2
3// Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
4
92f5a8d4
TL
5// This file was modified by Oracle on 2017, 2018, 2019.
6// Modifications copyright (c) 2017-2019, Oracle and/or its affiliates.
11fdf7f2
TL
7// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
8
9// Use, modification and distribution is subject to the Boost Software License,
10// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
11// http://www.boost.org/LICENSE_1_0.txt)
12
13// This file is converted from PROJ4, http://trac.osgeo.org/proj
14// PROJ4 is originally written by Gerald Evenden (then of the USGS)
15// PROJ4 is maintained by Frank Warmerdam
16// PROJ4 is converted to Boost.Geometry by Barend Gehrels
17
92f5a8d4 18// Last updated version of proj: 5.0.0
11fdf7f2
TL
19
20// Original copyright notice:
21
22// Purpose: Implementation of the aitoff (Aitoff) and wintri (Winkel Tripel)
92f5a8d4
TL
23// projections.
24// Author: Gerald Evenden (1995)
25// Drazen Tutic, Lovro Gradiser (2015) - add inverse
26// Thomas Knudsen (2016) - revise/add regression tests
11fdf7f2
TL
27// Copyright (c) 1995, Gerald Evenden
28
29// Permission is hereby granted, free of charge, to any person obtaining a
30// copy of this software and associated documentation files (the "Software"),
31// to deal in the Software without restriction, including without limitation
32// the rights to use, copy, modify, merge, publish, distribute, sublicense,
33// and/or sell copies of the Software, and to permit persons to whom the
34// Software is furnished to do so, subject to the following conditions:
35
36// The above copyright notice and this permission notice shall be included
37// in all copies or substantial portions of the Software.
38
39// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
40// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
41// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
42// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
43// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
44// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
45// DEALINGS IN THE SOFTWARE.
46
92f5a8d4
TL
47#ifndef BOOST_GEOMETRY_PROJECTIONS_AITOFF_HPP
48#define BOOST_GEOMETRY_PROJECTIONS_AITOFF_HPP
49
11fdf7f2 50#include <boost/core/ignore_unused.hpp>
11fdf7f2
TL
51
52#include <boost/geometry/srs/projections/impl/base_static.hpp>
53#include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
11fdf7f2 54#include <boost/geometry/srs/projections/impl/factory_entry.hpp>
92f5a8d4
TL
55#include <boost/geometry/srs/projections/impl/pj_param.hpp>
56#include <boost/geometry/srs/projections/impl/projects.hpp>
11fdf7f2 57
92f5a8d4 58#include <boost/geometry/util/math.hpp>
11fdf7f2 59
92f5a8d4 60namespace boost { namespace geometry
11fdf7f2 61{
11fdf7f2
TL
62
63namespace projections
64{
65 #ifndef DOXYGEN_NO_DETAIL
66 namespace detail { namespace aitoff
67 {
92f5a8d4
TL
68 enum mode_type {
69 mode_aitoff = 0,
70 mode_winkel_tripel = 1
71 };
72
11fdf7f2
TL
73 template <typename T>
74 struct par_aitoff
75 {
76 T cosphi1;
92f5a8d4 77 mode_type mode;
11fdf7f2
TL
78 };
79
92f5a8d4
TL
80 template <typename T, typename Parameters>
81 struct base_aitoff_spheroid
11fdf7f2 82 {
92f5a8d4 83 par_aitoff<T> m_proj_parm;
11fdf7f2
TL
84
85 // FORWARD(s_forward) spheroid
86 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
92f5a8d4 87 inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
11fdf7f2 88 {
92f5a8d4 89 T c, d;
11fdf7f2
TL
90
91 if((d = acos(cos(lp_lat) * cos(c = 0.5 * lp_lon)))) {/* basic Aitoff */
92 xy_x = 2. * d * cos(lp_lat) * sin(c) * (xy_y = 1. / sin(d));
93 xy_y *= d * sin(lp_lat);
94 } else
95 xy_x = xy_y = 0.;
92f5a8d4 96 if (this->m_proj_parm.mode == mode_winkel_tripel) { /* Winkel Tripel */
11fdf7f2
TL
97 xy_x = (xy_x + lp_lon * this->m_proj_parm.cosphi1) * 0.5;
98 xy_y = (xy_y + lp_lat) * 0.5;
99 }
100 }
101 /***********************************************************************************
102 *
103 * Inverse functions added by Drazen Tutic and Lovro Gradiser based on paper:
104 *
105 * I.Özbug Biklirici and Cengizhan Ipbüker. A General Algorithm for the Inverse
106 * Transformation of Map Projections Using Jacobian Matrices. In Proceedings of the
107 * Third International Symposium Mathematical & Computational Applications,
108 * pages 175{182, Turkey, September 2002.
109 *
92f5a8d4 110 * Expected accuracy is defined by epsilon = 1e-12. Should be appropriate for
11fdf7f2
TL
111 * most applications of Aitoff and Winkel Tripel projections.
112 *
113 * Longitudes of 180W and 180E can be mixed in solution obtained.
114 *
115 * Inverse for Aitoff projection in poles is undefined, longitude value of 0 is assumed.
116 *
117 * Contact : dtutic@geof.hr
118 * Date: 2015-02-16
119 *
120 ************************************************************************************/
121
122 // INVERSE(s_inverse) sphere
123 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
92f5a8d4 124 inline void inv(Parameters const& , T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
11fdf7f2 125 {
92f5a8d4
TL
126 static const T pi = detail::pi<T>();
127 static const T two_pi = detail::two_pi<T>();
128 static const T epsilon = 1e-12;
11fdf7f2 129
92f5a8d4
TL
130 int iter, max_iter = 10, round = 0, max_round = 20;
131 T D, C, f1, f2, f1p, f1l, f2p, f2l, dp, dl, sl, sp, cp, cl, x, y;
11fdf7f2 132
92f5a8d4
TL
133 if ((fabs(xy_x) < epsilon) && (fabs(xy_y) < epsilon )) {
134 lp_lat = 0.; lp_lon = 0.;
135 return;
136 }
11fdf7f2
TL
137
138 /* intial values for Newton-Raphson method */
139 lp_lat = xy_y; lp_lon = xy_x;
140 do {
141 iter = 0;
142 do {
143 sl = sin(lp_lon * 0.5); cl = cos(lp_lon * 0.5);
144 sp = sin(lp_lat); cp = cos(lp_lat);
145 D = cp * cl;
92f5a8d4
TL
146 C = 1. - D * D;
147 D = acos(D) / math::pow(C, T(1.5));
148 f1 = 2. * D * C * cp * sl;
149 f2 = D * C * sp;
150 f1p = 2.* (sl * cl * sp * cp / C - D * sp * sl);
151 f1l = cp * cp * sl * sl / C + D * cp * cl * sp * sp;
152 f2p = sp * sp * cl / C + D * sl * sl * cp;
153 f2l = 0.5 * (sp * cp * sl / C - D * sp * cp * cp * sl * cl);
154 if (this->m_proj_parm.mode == mode_winkel_tripel) { /* Winkel Tripel */
11fdf7f2
TL
155 f1 = 0.5 * (f1 + lp_lon * this->m_proj_parm.cosphi1);
156 f2 = 0.5 * (f2 + lp_lat);
157 f1p *= 0.5;
158 f1l = 0.5 * (f1l + this->m_proj_parm.cosphi1);
159 f2p = 0.5 * (f2p + 1.);
160 f2l *= 0.5;
161 }
162 f1 -= xy_x; f2 -= xy_y;
163 dl = (f2 * f1p - f1 * f2p) / (dp = f1p * f2l - f2p * f1l);
164 dp = (f1 * f2l - f2 * f1l) / dp;
92f5a8d4 165 dl = fmod(dl, pi); /* set to interval [-M_PI, M_PI] */
11fdf7f2 166 lp_lat -= dp; lp_lon -= dl;
92f5a8d4
TL
167 } while ((fabs(dp) > epsilon || fabs(dl) > epsilon) && (iter++ < max_iter));
168 if (lp_lat > two_pi) lp_lat -= 2.*(lp_lat-two_pi); /* correct if symmetrical solution for Aitoff */
169 if (lp_lat < -two_pi) lp_lat -= 2.*(lp_lat+two_pi); /* correct if symmetrical solution for Aitoff */
170 if ((fabs(fabs(lp_lat) - two_pi) < epsilon) && (!this->m_proj_parm.mode)) lp_lon = 0.; /* if pole in Aitoff, return longitude of 0 */
11fdf7f2
TL
171
172 /* calculate x,y coordinates with solution obtained */
92f5a8d4 173 if((D = acos(cos(lp_lat) * cos(C = 0.5 * lp_lon))) != 0.0) {/* Aitoff */
11fdf7f2
TL
174 x = 2. * D * cos(lp_lat) * sin(C) * (y = 1. / sin(D));
175 y *= D * sin(lp_lat);
176 } else
177 x = y = 0.;
92f5a8d4 178 if (this->m_proj_parm.mode == mode_winkel_tripel) { /* Winkel Tripel */
11fdf7f2
TL
179 x = (x + lp_lon * this->m_proj_parm.cosphi1) * 0.5;
180 y = (y + lp_lat) * 0.5;
181 }
182 /* if too far from given values of x,y, repeat with better approximation of phi,lam */
92f5a8d4 183 } while (((fabs(xy_x-x) > epsilon) || (fabs(xy_y-y) > epsilon)) && (round++ < max_round));
11fdf7f2 184
92f5a8d4
TL
185 if (iter == max_iter && round == max_round)
186 {
187 BOOST_THROW_EXCEPTION( projection_exception(error_non_convergent) );
188 //fprintf(stderr, "Warning: Accuracy of 1e-12 not reached. Last increments: dlat=%e and dlon=%e\n", dp, dl);
189 }
11fdf7f2
TL
190 }
191
192 static inline std::string get_name()
193 {
194 return "aitoff_spheroid";
195 }
196
197 };
198
92f5a8d4
TL
199 template <typename Parameters>
200 inline void setup(Parameters& par)
11fdf7f2 201 {
11fdf7f2
TL
202 par.es = 0.;
203 }
204
205
206 // Aitoff
207 template <typename Parameters, typename T>
208 inline void setup_aitoff(Parameters& par, par_aitoff<T>& proj_parm)
209 {
92f5a8d4
TL
210 proj_parm.mode = mode_aitoff;
211 setup(par);
11fdf7f2
TL
212 }
213
214 // Winkel Tripel
92f5a8d4
TL
215 template <typename Params, typename Parameters, typename T>
216 inline void setup_wintri(Params& params, Parameters& par, par_aitoff<T>& proj_parm)
11fdf7f2 217 {
92f5a8d4
TL
218 static const T two_div_pi = detail::two_div_pi<T>();
219
220 T phi1;
11fdf7f2 221
92f5a8d4
TL
222 proj_parm.mode = mode_winkel_tripel;
223 if (pj_param_r<srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1, phi1)) {
224 if ((proj_parm.cosphi1 = cos(phi1)) == 0.)
225 BOOST_THROW_EXCEPTION( projection_exception(error_lat_larger_than_90) );
11fdf7f2 226 } else /* 50d28' or phi1=acos(2/pi) */
92f5a8d4
TL
227 proj_parm.cosphi1 = two_div_pi;
228 setup(par);
11fdf7f2
TL
229 }
230
231 }} // namespace detail::aitoff
232 #endif // doxygen
233
234 /*!
235 \brief Aitoff projection
236 \ingroup projections
237 \tparam Geographic latlong point type
238 \tparam Cartesian xy point type
239 \tparam Parameters parameter type
240 \par Projection characteristics
241 - Miscellaneous
242 - Spheroid
243 \par Example
244 \image html ex_aitoff.gif
245 */
92f5a8d4
TL
246 template <typename T, typename Parameters>
247 struct aitoff_spheroid : public detail::aitoff::base_aitoff_spheroid<T, Parameters>
11fdf7f2 248 {
92f5a8d4
TL
249 template <typename Params>
250 inline aitoff_spheroid(Params const& , Parameters & par)
11fdf7f2 251 {
92f5a8d4 252 detail::aitoff::setup_aitoff(par, this->m_proj_parm);
11fdf7f2
TL
253 }
254 };
255
256 /*!
257 \brief Winkel Tripel projection
258 \ingroup projections
259 \tparam Geographic latlong point type
260 \tparam Cartesian xy point type
261 \tparam Parameters parameter type
262 \par Projection characteristics
263 - Miscellaneous
264 - Spheroid
265 \par Projection parameters
266 - lat_1: Latitude of first standard parallel (degrees)
267 \par Example
268 \image html ex_wintri.gif
269 */
92f5a8d4
TL
270 template <typename T, typename Parameters>
271 struct wintri_spheroid : public detail::aitoff::base_aitoff_spheroid<T, Parameters>
11fdf7f2 272 {
92f5a8d4
TL
273 template <typename Params>
274 inline wintri_spheroid(Params const& params, Parameters & par)
11fdf7f2 275 {
92f5a8d4 276 detail::aitoff::setup_wintri(params, par, this->m_proj_parm);
11fdf7f2
TL
277 }
278 };
279
280 #ifndef DOXYGEN_NO_DETAIL
281 namespace detail
282 {
283
284 // Static projection
92f5a8d4
TL
285 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_aitoff, aitoff_spheroid)
286 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_wintri, wintri_spheroid)
11fdf7f2
TL
287
288 // Factory entry(s)
92f5a8d4
TL
289 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(aitoff_entry, aitoff_spheroid)
290 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(wintri_entry, wintri_spheroid)
11fdf7f2 291
92f5a8d4 292 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(aitoff_init)
11fdf7f2 293 {
92f5a8d4
TL
294 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(aitoff, aitoff_entry)
295 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(wintri, wintri_entry)
11fdf7f2
TL
296 }
297
298 } // namespace detail
299 #endif // doxygen
300
301} // namespace projections
302
303}} // namespace boost::geometry
304
305#endif // BOOST_GEOMETRY_PROJECTIONS_AITOFF_HPP
306