1 // Boost.Geometry - gis-projections (based on PROJ4) 2 3 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands. 4 5 // This file was modified by Oracle on 2017, 2018. 6 // Modifications copyright (c) 2017, 2018, Oracle and/or its affiliates. 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 18 // Last updated version of proj: 5.0.0 19 20 // Original copyright notice: 21 22 // Permission is hereby granted, free of charge, to any person obtaining a 23 // copy of this software and associated documentation files (the "Software"), 24 // to deal in the Software without restriction, including without limitation 25 // the rights to use, copy, modify, merge, publish, distribute, sublicense, 26 // and/or sell copies of the Software, and to permit persons to whom the 27 // Software is furnished to do so, subject to the following conditions: 28 29 // The above copyright notice and this permission notice shall be included 30 // in all copies or substantial portions of the Software. 31 32 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 33 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 34 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 35 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 36 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 37 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 38 // DEALINGS IN THE SOFTWARE. 39 40 #ifndef BOOST_GEOMETRY_PROJECTIONS_VANDG_HPP 41 #define BOOST_GEOMETRY_PROJECTIONS_VANDG_HPP 42 43 #include <boost/geometry/util/math.hpp> 44 45 #include <boost/geometry/srs/projections/impl/base_static.hpp> 46 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp> 47 #include <boost/geometry/srs/projections/impl/projects.hpp> 48 #include <boost/geometry/srs/projections/impl/factory_entry.hpp> 49 50 namespace boost { namespace geometry 51 { 52 53 namespace projections 54 { 55 #ifndef DOXYGEN_NO_DETAIL 56 namespace detail { namespace vandg 57 { 58 59 static const double tolerance = 1.e-10; 60 61 template <typename T> C2_27()62 inline T C2_27() { return .07407407407407407407407407407407; } 63 template <typename T> PI4_3()64 inline T PI4_3() { return boost::math::constants::four_thirds_pi<T>(); } 65 template <typename T> TPISQ()66 inline T TPISQ() { return 19.739208802178717237668981999752; } 67 template <typename T> HPISQ()68 inline T HPISQ() { return 4.9348022005446793094172454999381; } 69 70 // template class, using CRTP to implement forward/inverse 71 template <typename T, typename Parameters> 72 struct base_vandg_spheroid 73 : public base_t_fi<base_vandg_spheroid<T, Parameters>, T, Parameters> 74 { base_vandg_spheroidboost::geometry::projections::detail::vandg::base_vandg_spheroid75 inline base_vandg_spheroid(const Parameters& par) 76 : base_t_fi<base_vandg_spheroid<T, Parameters>, T, Parameters>(*this, par) 77 {} 78 79 // FORWARD(s_forward) spheroid 80 // Project coordinates from geographic (lon, lat) to cartesian (x, y) fwdboost::geometry::projections::detail::vandg::base_vandg_spheroid81 inline void fwd(T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const 82 { 83 static const T half_pi = detail::half_pi<T>(); 84 static const T pi = detail::pi<T>(); 85 86 T al, al2, g, g2, p2; 87 88 p2 = fabs(lp_lat / half_pi); 89 if ((p2 - tolerance) > 1.) { 90 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) ); 91 } 92 if (p2 > 1.) 93 p2 = 1.; 94 if (fabs(lp_lat) <= tolerance) { 95 xy_x = lp_lon; 96 xy_y = 0.; 97 } else if (fabs(lp_lon) <= tolerance || fabs(p2 - 1.) < tolerance) { 98 xy_x = 0.; 99 xy_y = pi * tan(.5 * asin(p2)); 100 if (lp_lat < 0.) xy_y = -xy_y; 101 } else { 102 al = .5 * fabs(pi / lp_lon - lp_lon / pi); 103 al2 = al * al; 104 g = sqrt(1. - p2 * p2); 105 g = g / (p2 + g - 1.); 106 g2 = g * g; 107 p2 = g * (2. / p2 - 1.); 108 p2 = p2 * p2; 109 xy_x = g - p2; g = p2 + al2; 110 xy_x = pi * (al * xy_x + sqrt(al2 * xy_x * xy_x - g * (g2 - p2))) / g; 111 if (lp_lon < 0.) xy_x = -xy_x; 112 xy_y = fabs(xy_x / pi); 113 xy_y = 1. - xy_y * (xy_y + 2. * al); 114 if (xy_y < -tolerance) { 115 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) ); 116 } 117 if (xy_y < 0.) 118 xy_y = 0.; 119 else 120 xy_y = sqrt(xy_y) * (lp_lat < 0. ? -pi : pi); 121 } 122 } 123 124 // INVERSE(s_inverse) spheroid 125 // Project coordinates from cartesian (x, y) to geographic (lon, lat) invboost::geometry::projections::detail::vandg::base_vandg_spheroid126 inline void inv(T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const 127 { 128 static const T half_pi = detail::half_pi<T>(); 129 static const T pi = detail::pi<T>(); 130 static const T pi_sqr = detail::pi_sqr<T>(); 131 static const T third = detail::third<T>(); 132 static const T two_pi = detail::two_pi<T>(); 133 134 static const T C2_27 = vandg::C2_27<T>(); 135 static const T PI4_3 = vandg::PI4_3<T>(); 136 static const T TPISQ = vandg::TPISQ<T>(); 137 static const T HPISQ = vandg::HPISQ<T>(); 138 139 T t, c0, c1, c2, c3, al, r2, r, m, d, ay, x2, y2; 140 141 x2 = xy_x * xy_x; 142 if ((ay = fabs(xy_y)) < tolerance) { 143 lp_lat = 0.; 144 t = x2 * x2 + TPISQ * (x2 + HPISQ); 145 lp_lon = fabs(xy_x) <= tolerance ? 0. : 146 .5 * (x2 - pi_sqr + sqrt(t)) / xy_x; 147 return; 148 } 149 y2 = xy_y * xy_y; 150 r = x2 + y2; r2 = r * r; 151 c1 = - pi * ay * (r + pi_sqr); 152 c3 = r2 + two_pi * (ay * r + pi * (y2 + pi * (ay + half_pi))); 153 c2 = c1 + pi_sqr * (r - 3. * y2); 154 c0 = pi * ay; 155 c2 /= c3; 156 al = c1 / c3 - third * c2 * c2; 157 m = 2. * sqrt(-third * al); 158 d = C2_27 * c2 * c2 * c2 + (c0 * c0 - third * c2 * c1) / c3; 159 if (((t = fabs(d = 3. * d / (al * m))) - tolerance) <= 1.) { 160 d = t > 1. ? (d > 0. ? 0. : pi) : acos(d); 161 lp_lat = pi * (m * cos(d * third + PI4_3) - third * c2); 162 if (xy_y < 0.) lp_lat = -lp_lat; 163 t = r2 + TPISQ * (x2 - y2 + HPISQ); 164 lp_lon = fabs(xy_x) <= tolerance ? 0. : 165 .5 * (r - pi_sqr + (t <= 0. ? 0. : sqrt(t))) / xy_x; 166 } else { 167 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) ); 168 } 169 } 170 get_nameboost::geometry::projections::detail::vandg::base_vandg_spheroid171 static inline std::string get_name() 172 { 173 return "vandg_spheroid"; 174 } 175 176 }; 177 178 // van der Grinten (I) 179 template <typename Parameters> setup_vandg(Parameters & par)180 inline void setup_vandg(Parameters& par) 181 { 182 par.es = 0.; 183 } 184 185 }} // namespace detail::vandg 186 #endif // doxygen 187 188 /*! 189 \brief van der Grinten (I) projection 190 \ingroup projections 191 \tparam Geographic latlong point type 192 \tparam Cartesian xy point type 193 \tparam Parameters parameter type 194 \par Projection characteristics 195 - Miscellaneous 196 - Spheroid 197 \par Example 198 \image html ex_vandg.gif 199 */ 200 template <typename T, typename Parameters> 201 struct vandg_spheroid : public detail::vandg::base_vandg_spheroid<T, Parameters> 202 { 203 template <typename Params> vandg_spheroidboost::geometry::projections::vandg_spheroid204 inline vandg_spheroid(Params const& , Parameters const& par) 205 : detail::vandg::base_vandg_spheroid<T, Parameters>(par) 206 { 207 detail::vandg::setup_vandg(this->m_par); 208 } 209 }; 210 211 #ifndef DOXYGEN_NO_DETAIL 212 namespace detail 213 { 214 215 // Static projection BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::spar::proj_vandg,vandg_spheroid,vandg_spheroid)216 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::spar::proj_vandg, vandg_spheroid, vandg_spheroid) 217 218 // Factory entry(s) 219 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(vandg_entry, vandg_spheroid) 220 221 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(vandg_init) 222 { 223 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(vandg, vandg_entry) 224 } 225 226 } // namespace detail 227 #endif // doxygen 228 229 } // namespace projections 230 231 }} // namespace boost::geometry 232 233 #endif // BOOST_GEOMETRY_PROJECTIONS_VANDG_HPP 234 235