projections/proj/bipc.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_BIPC_HPP
0041 #define BOOST_GEOMETRY_PROJECTIONS_BIPC_HPP
0042
0043 #include <boost/geometry/srs/projections/impl/base_static.hpp>
0044 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
0045 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
0046 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
0047 #include <boost/geometry/srs/projections/impl/projects.hpp>
0048
0049 #include <boost/geometry/util/math.hpp>
0050
0051 #include <boost/math/special_functions/hypot.hpp>
0052
0053 namespace boost { namespace geometry
0054 {
0055
0056 namespace projections
0057 {
0058     #ifndef DOXYGEN_NO_DETAIL
0059     namespace detail { namespace bipc
0060     {
0061
0062             static const double epsilon = 1e-10;
0063             static const double epsilon10 = 1e-10;
0064             static const double one_plus_eps = 1.000000001;
0065             static const int n_iter = 10;
0066             static const double lamB = -.34894976726250681539;
0067             static const double n = .63055844881274687180;
0068             static const double F = 1.89724742567461030582;
0069             static const double Azab = .81650043674686363166;
0070             static const double Azba = 1.82261843856185925133;
0071             static const double const_T = 1.27246578267089012270;
0072             static const double rhoc = 1.20709121521568721927;
0073             static const double cAzc = .69691523038678375519;
0074             static const double sAzc = .71715351331143607555;
0075             static const double C45 = .70710678118654752469;
0076             static const double S45 = .70710678118654752410;
0077             static const double C20 = .93969262078590838411;
0078             static const double S20 = -.34202014332566873287;
0079             static const double R110 = 1.91986217719376253360;
0080             static const double R104 = 1.81514242207410275904;
0081
0082             struct par_bipc
0083             {
0084                 bool   noskew;
0085             };
0086
0087             template <typename T, typename Parameters>
0088             struct base_bipc_spheroid
0089             {
0090                 par_bipc m_proj_parm;
0091
0092                 // FORWARD(s_forward)  spheroid
0093                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0094                 inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
0095                 {
0096                     static const T half_pi = detail::half_pi<T>();
0097                     static const T pi = detail::pi<T>();
0098
0099                     T cphi, sphi, tphi, t, al, Az, z, Av, cdlam, sdlam, r;
0100                     int tag;
0101
0102                     cphi = cos(lp_lat);
0103                     sphi = sin(lp_lat);
0104                     cdlam = cos(sdlam = lamB - lp_lon);
0105                     sdlam = sin(sdlam);
0106                     if (fabs(fabs(lp_lat) - half_pi) < epsilon10) {
0107                         Az = lp_lat < 0. ? pi : 0.;
0108                         tphi = HUGE_VAL;
0109                     } else {
0110                         tphi = sphi / cphi;
0111                         Az = atan2(sdlam , C45 * (tphi - cdlam));
0112                     }
0113                     if( (tag = (Az > Azba)) ) {
0114                         cdlam = cos(sdlam = lp_lon + R110);
0115                         sdlam = sin(sdlam);
0116                         z = S20 * sphi + C20 * cphi * cdlam;
0117                         if (fabs(z) > 1.) {
0118                             if (fabs(z) > one_plus_eps)
0119                                 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0120                             else
0121                                 z = z < 0. ? -1. : 1.;
0122                         } else
0123                             z = acos(z);
0124                         if (tphi != HUGE_VAL)
0125                             Az = atan2(sdlam, (C20 * tphi - S20 * cdlam));
0126                         Av = Azab;
0127                         xy_y = rhoc;
0128                     } else {
0129                         z = S45 * (sphi + cphi * cdlam);
0130                         if (fabs(z) > 1.) {
0131                             if (fabs(z) > one_plus_eps)
0132                                 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0133                             else
0134                                 z = z < 0. ? -1. : 1.;
0135                         } else
0136                             z = acos(z);
0137                         Av = Azba;
0138                         xy_y = -rhoc;
0139                     }
0140                     if (z < 0.) {
0141                         BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0142                     }
0143                     r = F * (t = math::pow(tan(T(0.5) * z), n));
0144                     if ((al = .5 * (R104 - z)) < 0.) {
0145                         BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0146                     }
0147                     al = (t + math::pow(al, n)) / const_T;
0148                     if (fabs(al) > 1.) {
0149                         if (fabs(al) > one_plus_eps)
0150                             BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0151                         else
0152                             al = al < 0. ? -1. : 1.;
0153                     } else
0154                         al = acos(al);
0155                     if (fabs(t = n * (Av - Az)) < al)
0156                         r /= cos(al + (tag ? t : -t));
0157                     xy_x = r * sin(t);
0158                     xy_y += (tag ? -r : r) * cos(t);
0159                     if (this->m_proj_parm.noskew) {
0160                         t = xy_x;
0161                         xy_x = -xy_x * cAzc - xy_y * sAzc;
0162                         xy_y = -xy_y * cAzc + t * sAzc;
0163                     }
0164                 }
0165
0166                 // INVERSE(s_inverse)  spheroid
0167                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0168                 inline void inv(Parameters const& , T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
0169                 {
0170                     T t, r, rp, rl, al, z, fAz, Az, s, c, Av;
0171                     int neg, i;
0172
0173                     if (this->m_proj_parm.noskew) {
0174                         t = xy_x;
0175                         xy_x = -xy_x * cAzc + xy_y * sAzc;
0176                         xy_y = -xy_y * cAzc - t * sAzc;
0177                     }
0178                     if( (neg = (xy_x < 0.)) ) {
0179                         xy_y = rhoc - xy_y;
0180                         s = S20;
0181                         c = C20;
0182                         Av = Azab;
0183                     } else {
0184                         xy_y += rhoc;
0185                         s = S45;
0186                         c = C45;
0187                         Av = Azba;
0188                     }
0189                     rl = rp = r = boost::math::hypot(xy_x, xy_y);
0190                     fAz = fabs(Az = atan2(xy_x, xy_y));
0191                     for (i = n_iter; i ; --i) {
0192                         z = 2. * atan(math::pow(r / F,T(1) / n));
0193                         al = acos((math::pow(tan(T(0.5) * z), n) +
0194                            math::pow(tan(T(0.5) * (R104 - z)), n)) / const_T);
0195                         if (fAz < al)
0196                             r = rp * cos(al + (neg ? Az : -Az));
0197                         if (fabs(rl - r) < epsilon)
0198                             break;
0199                         rl = r;
0200                     }
0201                     if (! i)
0202                         BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0203                     Az = Av - Az / n;
0204                     lp_lat = asin(s * cos(z) + c * sin(z) * cos(Az));
0205                     lp_lon = atan2(sin(Az), c / tan(z) - s * cos(Az));
0206                     if (neg)
0207                         lp_lon -= R110;
0208                     else
0209                         lp_lon = lamB - lp_lon;
0210                 }
0211
0212                 static inline std::string get_name()
0213                 {
0214                     return "bipc_spheroid";
0215                 }
0216
0217             };
0218
0219             // Bipolar conic of western hemisphere
0220             template <typename Params, typename Parameters>
0221             inline void setup_bipc(Params const& params, Parameters& par, par_bipc& proj_parm)
0222             {
0223                 proj_parm.noskew = pj_get_param_b<srs::spar::ns>(params, "ns", srs::dpar::ns);
0224                 par.es = 0.;
0225             }
0226
0227     }} // namespace detail::bipc
0228     #endif // doxygen
0229
0230     /*!
0231         \brief Bipolar conic of western hemisphere projection
0232         \ingroup projections
0233         \tparam Geographic latlong point type
0234         \tparam Cartesian xy point type
0235         \tparam Parameters parameter type
0236         \par Projection characteristics
0237          - Conic
0238          - Spheroid
0239         \par Projection parameters
0240          - ns (boolean)
0241         \par Example
0242         \image html ex_bipc.gif
0243     */
0244     template <typename T, typename Parameters>
0245     struct bipc_spheroid : public detail::bipc::base_bipc_spheroid<T, Parameters>
0246     {
0247         template <typename Params>
0248         inline bipc_spheroid(Params const& params, Parameters & par)
0249         {
0250             detail::bipc::setup_bipc(params, par, this->m_proj_parm);
0251         }
0252     };
0253
0254     #ifndef DOXYGEN_NO_DETAIL
0255     namespace detail
0256     {
0257
0258         // Static projection
0259         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_bipc, bipc_spheroid)
0260
0261         // Factory entry(s)
0262         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(bipc_entry, bipc_spheroid)
0263
0264         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(bipc_init)
0265         {
0266             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(bipc, bipc_entry)
0267         }
0268
0269     } // namespace detail
0270     #endif // doxygen
0271
0272 } // namespace projections
0273
0274 }} // namespace boost::geometry
0275
0276 #endif // BOOST_GEOMETRY_PROJECTIONS_BIPC_HPP
0277