projections/proj/vandg.hpp

0001 // Boost.Geometry - gis-projections (based on PROJ4)
0002
0003 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
0004
0005 // This file was modified by Oracle on 2017, 2018, 2019.
0006 // Modifications copyright (c) 2017-2019, Oracle and/or its affiliates.
0007 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
0008
0009 // Use, modification and distribution is subject to the Boost Software License,
0010 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0011 // http://www.boost.org/LICENSE_1_0.txt)
0012
0013 // This file is converted from PROJ4, http://trac.osgeo.org/proj
0014 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
0015 // PROJ4 is maintained by Frank Warmerdam
0016 // PROJ4 is converted to Boost.Geometry by Barend Gehrels
0017
0018 // Last updated version of proj: 5.0.0
0019
0020 // Original copyright notice:
0021
0022 // Permission is hereby granted, free of charge, to any person obtaining a
0023 // copy of this software and associated documentation files (the "Software"),
0024 // to deal in the Software without restriction, including without limitation
0025 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
0026 // and/or sell copies of the Software, and to permit persons to whom the
0027 // Software is furnished to do so, subject to the following conditions:
0028
0029 // The above copyright notice and this permission notice shall be included
0030 // in all copies or substantial portions of the Software.
0031
0032 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
0033 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0034 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
0035 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0036 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
0037 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
0038 // DEALINGS IN THE SOFTWARE.
0039
0040 #ifndef BOOST_GEOMETRY_PROJECTIONS_VANDG_HPP
0041 #define BOOST_GEOMETRY_PROJECTIONS_VANDG_HPP
0042
0043 #include <boost/geometry/util/math.hpp>
0044
0045 #include <boost/geometry/srs/projections/impl/base_static.hpp>
0046 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
0047 #include <boost/geometry/srs/projections/impl/projects.hpp>
0048 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
0049
0050 namespace boost { namespace geometry
0051 {
0052
0053 namespace projections
0054 {
0055     #ifndef DOXYGEN_NO_DETAIL
0056     namespace detail { namespace vandg
0057     {
0058
0059             static const double tolerance = 1.e-10;
0060
0061             template <typename T>
0062             inline T C2_27() { return .07407407407407407407407407407407; }
0063             template <typename T>
0064             inline T PI4_3() { return boost::math::constants::four_thirds_pi<T>(); }
0065             template <typename T>
0066             inline T TPISQ() { return 19.739208802178717237668981999752; }
0067             template <typename T>
0068             inline T HPISQ() { return 4.9348022005446793094172454999381; }
0069
0070             template <typename T, typename Parameters>
0071             struct base_vandg_spheroid
0072             {
0073                 // FORWARD(s_forward)  spheroid
0074                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0075                 inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
0076                 {
0077                     static const T half_pi = detail::half_pi<T>();
0078                     static const T pi = detail::pi<T>();
0079
0080                     T  al, al2, g, g2, p2;
0081
0082                     p2 = fabs(lp_lat / half_pi);
0083                     if ((p2 - tolerance) > 1.) {
0084                         BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0085                     }
0086                     if (p2 > 1.)
0087                         p2 = 1.;
0088                     if (fabs(lp_lat) <= tolerance) {
0089                         xy_x = lp_lon;
0090                         xy_y = 0.;
0091                     } else if (fabs(lp_lon) <= tolerance || fabs(p2 - 1.) < tolerance) {
0092                         xy_x = 0.;
0093                         xy_y = pi * tan(.5 * asin(p2));
0094                         if (lp_lat < 0.) xy_y = -xy_y;
0095                     } else {
0096                         al = .5 * fabs(pi / lp_lon - lp_lon / pi);
0097                         al2 = al * al;
0098                         g = sqrt(1. - p2 * p2);
0099                         g = g / (p2 + g - 1.);
0100                         g2 = g * g;
0101                         p2 = g * (2. / p2 - 1.);
0102                         p2 = p2 * p2;
0103                         xy_x = g - p2; g = p2 + al2;
0104                         xy_x = pi * (al * xy_x + sqrt(al2 * xy_x * xy_x - g * (g2 - p2))) / g;
0105                         if (lp_lon < 0.) xy_x = -xy_x;
0106                         xy_y = fabs(xy_x / pi);
0107                         xy_y = 1. - xy_y * (xy_y + 2. * al);
0108                         if (xy_y < -tolerance) {
0109                             BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0110                         }
0111                         if (xy_y < 0.)
0112                             xy_y = 0.;
0113                         else
0114                             xy_y = sqrt(xy_y) * (lp_lat < 0. ? -pi : pi);
0115                     }
0116                 }
0117
0118                 // INVERSE(s_inverse)  spheroid
0119                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0120                 inline void inv(Parameters const& , T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
0121                 {
0122                     static const T half_pi = detail::half_pi<T>();
0123                     static const T pi = detail::pi<T>();
0124                     static const T pi_sqr = detail::pi_sqr<T>();
0125                     static const T third = detail::third<T>();
0126                     static const T two_pi = detail::two_pi<T>();
0127
0128                     static const T C2_27 = vandg::C2_27<T>();
0129                     static const T PI4_3 = vandg::PI4_3<T>();
0130                     static const T TPISQ = vandg::TPISQ<T>();
0131                     static const T HPISQ = vandg::HPISQ<T>();
0132
0133                     T t, c0, c1, c2, c3, al, r2, r, m, d, ay, x2, y2;
0134
0135                     x2 = xy_x * xy_x;
0136                     if ((ay = fabs(xy_y)) < tolerance) {
0137                         lp_lat = 0.;
0138                         t = x2 * x2 + TPISQ * (x2 + HPISQ);
0139                         lp_lon = fabs(xy_x) <= tolerance ? 0. :
0140                            .5 * (x2 - pi_sqr + sqrt(t)) / xy_x;
0141                             return;
0142                     }
0143                     y2 = xy_y * xy_y;
0144                     r = x2 + y2;    r2 = r * r;
0145                     c1 = - pi * ay * (r + pi_sqr);
0146                     c3 = r2 + two_pi * (ay * r + pi * (y2 + pi * (ay + half_pi)));
0147                     c2 = c1 + pi_sqr * (r - 3. *  y2);
0148                     c0 = pi * ay;
0149                     c2 /= c3;
0150                     al = c1 / c3 - third * c2 * c2;
0151                     m = 2. * sqrt(-third * al);
0152                     d = C2_27 * c2 * c2 * c2 + (c0 * c0 - third * c2 * c1) / c3;
0153                     if (((t = fabs(d = 3. * d / (al * m))) - tolerance) <= 1.) {
0154                         d = t > 1. ? (d > 0. ? 0. : pi) : acos(d);
0155                         lp_lat = pi * (m * cos(d * third + PI4_3) - third * c2);
0156                         if (xy_y < 0.) lp_lat = -lp_lat;
0157                         t = r2 + TPISQ * (x2 - y2 + HPISQ);
0158                         lp_lon = fabs(xy_x) <= tolerance ? 0. :
0159                            .5 * (r - pi_sqr + (t <= 0. ? 0. : sqrt(t))) / xy_x;
0160                     } else {
0161                         BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0162                     }
0163                 }
0164
0165                 static inline std::string get_name()
0166                 {
0167                     return "vandg_spheroid";
0168                 }
0169
0170             };
0171
0172             // van der Grinten (I)
0173             template <typename Parameters>
0174             inline void setup_vandg(Parameters& par)
0175             {
0176                 par.es = 0.;
0177             }
0178
0179     }} // namespace detail::vandg
0180     #endif // doxygen
0181
0182     /*!
0183         \brief van der Grinten (I) projection
0184         \ingroup projections
0185         \tparam Geographic latlong point type
0186         \tparam Cartesian xy point type
0187         \tparam Parameters parameter type
0188         \par Projection characteristics
0189          - Miscellaneous
0190          - Spheroid
0191         \par Example
0192         \image html ex_vandg.gif
0193     */
0194     template <typename T, typename Parameters>
0195     struct vandg_spheroid : public detail::vandg::base_vandg_spheroid<T, Parameters>
0196     {
0197         template <typename Params>
0198         inline vandg_spheroid(Params const& , Parameters & par)
0199         {
0200             detail::vandg::setup_vandg(par);
0201         }
0202     };
0203
0204     #ifndef DOXYGEN_NO_DETAIL
0205     namespace detail
0206     {
0207
0208         // Static projection
0209         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_vandg, vandg_spheroid)
0210
0211         // Factory entry(s)
0212         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(vandg_entry, vandg_spheroid)
0213
0214         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(vandg_init)
0215         {
0216             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(vandg, vandg_entry)
0217         }
0218
0219     } // namespace detail
0220     #endif // doxygen
0221
0222 } // namespace projections
0223
0224 }} // namespace boost::geometry
0225
0226 #endif // BOOST_GEOMETRY_PROJECTIONS_VANDG_HPP
0227