ceph/src/boost/boost/geometry/srs/projections/proj/mod_ster.hpp

   1 #ifndef BOOST_GEOMETRY_PROJECTIONS_MOD_STER_HPP
   2 #define BOOST_GEOMETRY_PROJECTIONS_MOD_STER_HPP
   3
   4 // Boost.Geometry - extensions-gis-projections (based on PROJ4)
   5 // This file is automatically generated. DO NOT EDIT.
   6
   7 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
   8
   9 // This file was modified by Oracle on 2017.
  10 // Modifications copyright (c) 2017, Oracle and/or its affiliates.
  11 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
  12
  13 // Use, modification and distribution is subject to the Boost Software License,
  14 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  15 // http://www.boost.org/LICENSE_1_0.txt)
  16
  17 // This file is converted from PROJ4, http://trac.osgeo.org/proj
  18 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
  19 // PROJ4 is maintained by Frank Warmerdam
  20 // PROJ4 is converted to Boost.Geometry by Barend Gehrels
  21
  22 // Last updated version of proj: 4.9.1
  23
  24 // Original copyright notice:
  25
  26 // Permission is hereby granted, free of charge, to any person obtaining a
  27 // copy of this software and associated documentation files (the "Software"),
  28 // to deal in the Software without restriction, including without limitation
  29 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  30 // and/or sell copies of the Software, and to permit persons to whom the
  31 // Software is furnished to do so, subject to the following conditions:
  32
  33 // The above copyright notice and this permission notice shall be included
  34 // in all copies or substantial portions of the Software.
  35
  36 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  37 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  38 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  39 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  40 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  41 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  42 // DEALINGS IN THE SOFTWARE.
  43
  44 #include <boost/geometry/util/math.hpp>
  45 #include <boost/math/special_functions/hypot.hpp>
  46
  47 #include <boost/geometry/srs/projections/impl/base_static.hpp>
  48 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
  49 #include <boost/geometry/srs/projections/impl/projects.hpp>
  50 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
  51 #include <boost/geometry/srs/projections/impl/aasincos.hpp>
  52 #include <boost/geometry/srs/projections/impl/pj_zpoly1.hpp>
  53
  54 namespace boost { namespace geometry
  55 {
  56
  57 namespace srs { namespace par4
  58 {
  59     struct mil_os {};
  60     struct lee_os {};
  61     struct gs48 {};
  62     struct alsk {};
  63     struct gs50 {};
  64
  65 }} //namespace srs::par4
  66
  67 namespace projections
  68 {
  69     #ifndef DOXYGEN_NO_DETAIL
  70     namespace detail { namespace mod_ster
  71     {
  72
  73             static const double EPSLN = 1e-10;
  74
  75             template <typename T>
  76             struct par_mod_ster
  77             {
  78                 COMPLEX<T> *zcoeff;
  79                 T          cchio, schio;
  80                 int        n;
  81             };
  82
  83             /* based upon Snyder and Linck, USGS-NMD */
  84
  85             // template class, using CRTP to implement forward/inverse
  86             template <typename CalculationType, typename Parameters>
  87             struct base_mod_ster_ellipsoid : public base_t_fi<base_mod_ster_ellipsoid<CalculationType, Parameters>,
  88                      CalculationType, Parameters>
  89             {
  90
  91                 typedef CalculationType geographic_type;
  92                 typedef CalculationType cartesian_type;
  93
  94                 par_mod_ster<CalculationType> m_proj_parm;
  95
  96                 inline base_mod_ster_ellipsoid(const Parameters& par)
  97                     : base_t_fi<base_mod_ster_ellipsoid<CalculationType, Parameters>,
  98                      CalculationType, Parameters>(*this, par) {}
  99
 100                 // FORWARD(e_forward)  ellipsoid
 101                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
 102                 inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
 103                 {
 104                     static const CalculationType HALFPI = detail::HALFPI<CalculationType>();
 105
 106                     CalculationType sinlon, coslon, esphi, chi, schi, cchi, s;
 107                     COMPLEX<CalculationType> p;
 108
 109                     sinlon = sin(lp_lon);
 110                     coslon = cos(lp_lon);
 111                     esphi = this->m_par.e * sin(lp_lat);
 112                     chi = 2. * atan(tan((HALFPI + lp_lat) * .5) *
 113                         pow((1. - esphi) / (1. + esphi), this->m_par.e * .5)) - HALFPI;
 114                     schi = sin(chi);
 115                     cchi = cos(chi);
 116                     s = 2. / (1. + this->m_proj_parm.schio * schi + this->m_proj_parm.cchio * cchi * coslon);
 117                     p.r = s * cchi * sinlon;
 118                     p.i = s * (this->m_proj_parm.cchio * schi - this->m_proj_parm.schio * cchi * coslon);
 119                     p = pj_zpoly1(p, this->m_proj_parm.zcoeff, this->m_proj_parm.n);
 120                     xy_x = p.r;
 121                     xy_y = p.i;
 122                 }
 123
 124                 // INVERSE(e_inverse)  ellipsoid
 125                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 126                 inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
 127                 {
 128                     static const CalculationType HALFPI = detail::HALFPI<CalculationType>();
 129
 130                     int nn;
 131                     COMPLEX<CalculationType> p, fxy, fpxy, dp;
 132                     CalculationType den, rh = 0, z, sinz = 0, cosz = 0, chi, phi = 0, dphi, esphi;
 133
 134                     p.r = xy_x;
 135                     p.i = xy_y;
 136                     for (nn = 20; nn ;--nn) {
 137                         fxy = pj_zpolyd1(p, this->m_proj_parm.zcoeff, this->m_proj_parm.n, &fpxy);
 138                         fxy.r -= xy_x;
 139                         fxy.i -= xy_y;
 140                         den = fpxy.r * fpxy.r + fpxy.i * fpxy.i;
 141                         dp.r = -(fxy.r * fpxy.r + fxy.i * fpxy.i) / den;
 142                         dp.i = -(fxy.i * fpxy.r - fxy.r * fpxy.i) / den;
 143                         p.r += dp.r;
 144                         p.i += dp.i;
 145                         if ((fabs(dp.r) + fabs(dp.i)) <= EPSLN)
 146                             break;
 147                     }
 148                     if (nn) {
 149                         rh = boost::math::hypot(p.r, p.i);
 150                         z = 2. * atan(.5 * rh);
 151                         sinz = sin(z);
 152                         cosz = cos(z);
 153                         lp_lon = this->m_par.lam0;
 154                         if (fabs(rh) <= EPSLN) {
 155                             lp_lat = this->m_par.phi0;
 156                             return;
 157                         }
 158                         chi = aasin(cosz * this->m_proj_parm.schio + p.i * sinz * this->m_proj_parm.cchio / rh);
 159                         phi = chi;
 160                         for (nn = 20; nn ;--nn) {
 161                             esphi = this->m_par.e * sin(phi);
 162                             dphi = 2. * atan(tan((HALFPI + chi) * .5) *
 163                                 pow((1. + esphi) / (1. - esphi), this->m_par.e * .5)) - HALFPI - phi;
 164                             phi += dphi;
 165                             if (fabs(dphi) <= EPSLN)
 166                                 break;
 167                         }
 168                     }
 169                     if (nn) {
 170                         lp_lat = phi;
 171                         lp_lon = atan2(p.r * sinz, rh * this->m_proj_parm.cchio * cosz - p.i *
 172                             this->m_proj_parm.schio * sinz);
 173                     } else
 174                         lp_lon = lp_lat = HUGE_VAL;
 175                 }
 176
 177                 static inline std::string get_name()
 178                 {
 179                     return "mod_ster_ellipsoid";
 180                 }
 181
 182             };
 183
 184             template <typename Parameters, typename T>
 185             inline void setup(Parameters& par, par_mod_ster<T>& proj_parm)  /* general initialization */
 186             {
 187                 T esphi, chio;
 188
 189                 if (par.es) {
 190                     esphi = par.e * sin(par.phi0);
 191                     chio = 2. * atan(tan((geometry::math::half_pi<T>() + par.phi0) * .5) *
 192                         pow((1. - esphi) / (1. + esphi), par.e * .5)) - geometry::math::half_pi<T>();
 193                 } else
 194                     chio = par.phi0;
 195                 proj_parm.schio = sin(chio);
 196                 proj_parm.cchio = cos(chio);
 197             }
 198
 199
 200             // Miller Oblated Stereographic
 201             template <typename Parameters, typename T>
 202             inline void setup_mil_os(Parameters& par, par_mod_ster<T>& proj_parm)
 203             {
 204                 static COMPLEX<T> /* Miller Oblated Stereographic */
 205             AB[] = {
 206                 {0.924500,    0.},
 207                 {0.,            0.},
 208                 {0.019430,    0.}
 209             };
 210
 211                 proj_parm.n = 2;
 212                 par.lam0 = geometry::math::d2r<T>() * 20.;
 213                 par.phi0 = geometry::math::d2r<T>() * 18.;
 214                 proj_parm.zcoeff = AB;
 215                 par.es = 0.;
 216                 setup(par, proj_parm);
 217             }
 218
 219             // Lee Oblated Stereographic
 220             template <typename Parameters, typename T>
 221             inline void setup_lee_os(Parameters& par, par_mod_ster<T>& proj_parm)
 222             {
 223                 static COMPLEX<T> /* Lee Oblated Stereographic */
 224             AB[] = {
 225                 {0.721316,    0.},
 226                 {0.,            0.},
 227                     {-0.0088162,     -0.00617325}
 228             };
 229
 230                 proj_parm.n = 2;
 231                 par.lam0 = geometry::math::d2r<T>() * -165.;
 232                 par.phi0 = geometry::math::d2r<T>() * -10.;
 233                 proj_parm.zcoeff = AB;
 234                 par.es = 0.;
 235                 setup(par, proj_parm);
 236             }
 237
 238             // Mod. Stererographics of 48 U.S.
 239             template <typename Parameters, typename T>
 240             inline void setup_gs48(Parameters& par, par_mod_ster<T>& proj_parm)
 241             {
 242                 static COMPLEX<T> /* 48 United States */
 243             AB[] = {
 244                 {0.98879,    0.},
 245                 {0.,        0.},
 246                 {-0.050909,    0.},
 247                 {0.,        0.},
 248                     {0.075528,    0.}
 249             };
 250
 251                 proj_parm.n = 4;
 252                 par.lam0 = geometry::math::d2r<T>() * -96.;
 253                 par.phi0 = geometry::math::d2r<T>() * -39.;
 254                 proj_parm.zcoeff = AB;
 255                 par.es = 0.;
 256                 par.a = 6370997.;
 257                 setup(par, proj_parm);
 258             }
 259
 260             // Mod. Stererographics of Alaska
 261             template <typename Parameters, typename T>
 262             inline void setup_alsk(Parameters& par, par_mod_ster<T>& proj_parm)
 263             {
 264                 static COMPLEX<T>
 265             ABe[] = { /* Alaska ellipsoid */
 266                 {.9945303,    0.},
 267                 {.0052083,    -.0027404},
 268                 {.0072721,    .0048181},
 269                 {-.0151089,    -.1932526},
 270                 {.0642675,    -.1381226},
 271                 {.3582802,    -.2884586}},
 272             ABs[] = { /* Alaska sphere */
 273                 {.9972523,    0.},
 274                 {.0052513,    -.0041175},
 275                 {.0074606,    .0048125},
 276                 {-.0153783,    -.1968253},
 277                 {.0636871,    -.1408027},
 278                     {.3660976,    -.2937382}
 279             };
 280
 281                 proj_parm.n = 5;
 282                 par.lam0 = geometry::math::d2r<T>() * -152.;
 283                 par.phi0 = geometry::math::d2r<T>() * 64.;
 284                 if (par.es) { /* fixed ellipsoid/sphere */
 285                     proj_parm.zcoeff = ABe;
 286                     par.a = 6378206.4;
 287                     par.e = sqrt(par.es = 0.00676866);
 288                 } else {
 289                     proj_parm.zcoeff = ABs;
 290                     par.a = 6370997.;
 291                 }
 292                 setup(par, proj_parm);
 293             }
 294
 295             // Mod. Stererographics of 50 U.S.
 296             template <typename Parameters, typename T>
 297             inline void setup_gs50(Parameters& par, par_mod_ster<T>& proj_parm)
 298             {
 299                 static COMPLEX<T>
 300             ABe[] = { /* GS50 ellipsoid */
 301                 {.9827497,    0.},
 302                 {.0210669,    .0053804},
 303                 {-.1031415,    -.0571664},
 304                 {-.0323337,    -.0322847},
 305                 {.0502303,    .1211983},
 306                 {.0251805,    .0895678},
 307                 {-.0012315,    -.1416121},
 308                 {.0072202,    -.1317091},
 309                 {-.0194029,    .0759677},
 310                     {-.0210072,    .0834037}
 311             },
 312             ABs[] = { /* GS50 sphere */
 313                 {.9842990,    0.},
 314                 {.0211642,    .0037608},
 315                 {-.1036018,    -.0575102},
 316                 {-.0329095,    -.0320119},
 317                 {.0499471,    .1223335},
 318                 {.0260460,    .0899805},
 319                 {.0007388,    -.1435792},
 320                 {.0075848,    -.1334108},
 321                 {-.0216473,    .0776645},
 322                     {-.0225161,    .0853673}
 323             };
 324
 325                 proj_parm.n = 9;
 326                 par.lam0 = geometry::math::d2r<T>() * -120.;
 327                 par.phi0 = geometry::math::d2r<T>() * 45.;
 328                 if (par.es) { /* fixed ellipsoid/sphere */
 329                     proj_parm.zcoeff = ABe;
 330                     par.a = 6378206.4;
 331                     par.e = sqrt(par.es = 0.00676866);
 332                 } else {
 333                     proj_parm.zcoeff = ABs;
 334                     par.a = 6370997.;
 335                 }
 336                 setup(par, proj_parm);
 337             }
 338
 339     }} // namespace detail::mod_ster
 340     #endif // doxygen
 341
 342     /*!
 343         \brief Miller Oblated Stereographic projection
 344         \ingroup projections
 345         \tparam Geographic latlong point type
 346         \tparam Cartesian xy point type
 347         \tparam Parameters parameter type
 348         \par Projection characteristics
 349          - Azimuthal (mod)
 350         \par Example
 351         \image html ex_mil_os.gif
 352     */
 353     template <typename CalculationType, typename Parameters>
 354     struct mil_os_ellipsoid : public detail::mod_ster::base_mod_ster_ellipsoid<CalculationType, Parameters>
 355     {
 356         inline mil_os_ellipsoid(const Parameters& par) : detail::mod_ster::base_mod_ster_ellipsoid<CalculationType, Parameters>(par)
 357         {
 358             detail::mod_ster::setup_mil_os(this->m_par, this->m_proj_parm);
 359         }
 360     };
 361
 362     /*!
 363         \brief Lee Oblated Stereographic projection
 364         \ingroup projections
 365         \tparam Geographic latlong point type
 366         \tparam Cartesian xy point type
 367         \tparam Parameters parameter type
 368         \par Projection characteristics
 369          - Azimuthal (mod)
 370         \par Example
 371         \image html ex_lee_os.gif
 372     */
 373     template <typename CalculationType, typename Parameters>
 374     struct lee_os_ellipsoid : public detail::mod_ster::base_mod_ster_ellipsoid<CalculationType, Parameters>
 375     {
 376         inline lee_os_ellipsoid(const Parameters& par) : detail::mod_ster::base_mod_ster_ellipsoid<CalculationType, Parameters>(par)
 377         {
 378             detail::mod_ster::setup_lee_os(this->m_par, this->m_proj_parm);
 379         }
 380     };
 381
 382     /*!
 383         \brief Mod. Stererographics of 48 U.S. projection
 384         \ingroup projections
 385         \tparam Geographic latlong point type
 386         \tparam Cartesian xy point type
 387         \tparam Parameters parameter type
 388         \par Projection characteristics
 389          - Azimuthal (mod)
 390         \par Example
 391         \image html ex_gs48.gif
 392     */
 393     template <typename CalculationType, typename Parameters>
 394     struct gs48_ellipsoid : public detail::mod_ster::base_mod_ster_ellipsoid<CalculationType, Parameters>
 395     {
 396         inline gs48_ellipsoid(const Parameters& par) : detail::mod_ster::base_mod_ster_ellipsoid<CalculationType, Parameters>(par)
 397         {
 398             detail::mod_ster::setup_gs48(this->m_par, this->m_proj_parm);
 399         }
 400     };
 401
 402     /*!
 403         \brief Mod. Stererographics of Alaska projection
 404         \ingroup projections
 405         \tparam Geographic latlong point type
 406         \tparam Cartesian xy point type
 407         \tparam Parameters parameter type
 408         \par Projection characteristics
 409          - Azimuthal (mod)
 410         \par Example
 411         \image html ex_alsk.gif
 412     */
 413     template <typename CalculationType, typename Parameters>
 414     struct alsk_ellipsoid : public detail::mod_ster::base_mod_ster_ellipsoid<CalculationType, Parameters>
 415     {
 416         inline alsk_ellipsoid(const Parameters& par) : detail::mod_ster::base_mod_ster_ellipsoid<CalculationType, Parameters>(par)
 417         {
 418             detail::mod_ster::setup_alsk(this->m_par, this->m_proj_parm);
 419         }
 420     };
 421
 422     /*!
 423         \brief Mod. Stererographics of 50 U.S. projection
 424         \ingroup projections
 425         \tparam Geographic latlong point type
 426         \tparam Cartesian xy point type
 427         \tparam Parameters parameter type
 428         \par Projection characteristics
 429          - Azimuthal (mod)
 430         \par Example
 431         \image html ex_gs50.gif
 432     */
 433     template <typename CalculationType, typename Parameters>
 434     struct gs50_ellipsoid : public detail::mod_ster::base_mod_ster_ellipsoid<CalculationType, Parameters>
 435     {
 436         inline gs50_ellipsoid(const Parameters& par) : detail::mod_ster::base_mod_ster_ellipsoid<CalculationType, Parameters>(par)
 437         {
 438             detail::mod_ster::setup_gs50(this->m_par, this->m_proj_parm);
 439         }
 440     };
 441
 442     #ifndef DOXYGEN_NO_DETAIL
 443     namespace detail
 444     {
 445
 446         // Static projection
 447         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::mil_os, mil_os_ellipsoid, mil_os_ellipsoid)
 448         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::lee_os, lee_os_ellipsoid, lee_os_ellipsoid)
 449         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::gs48, gs48_ellipsoid, gs48_ellipsoid)
 450         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::alsk, alsk_ellipsoid, alsk_ellipsoid)
 451         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::gs50, gs50_ellipsoid, gs50_ellipsoid)
 452
 453         // Factory entry(s)
 454         template <typename CalculationType, typename Parameters>
 455         class mil_os_entry : public detail::factory_entry<CalculationType, Parameters>
 456         {
 457             public :
 458                 virtual base_v<CalculationType, Parameters>* create_new(const Parameters& par) const
 459                 {
 460                     return new base_v_fi<mil_os_ellipsoid<CalculationType, Parameters>, CalculationType, Parameters>(par);
 461                 }
 462         };
 463
 464         template <typename CalculationType, typename Parameters>
 465         class lee_os_entry : public detail::factory_entry<CalculationType, Parameters>
 466         {
 467             public :
 468                 virtual base_v<CalculationType, Parameters>* create_new(const Parameters& par) const
 469                 {
 470                     return new base_v_fi<lee_os_ellipsoid<CalculationType, Parameters>, CalculationType, Parameters>(par);
 471                 }
 472         };
 473
 474         template <typename CalculationType, typename Parameters>
 475         class gs48_entry : public detail::factory_entry<CalculationType, Parameters>
 476         {
 477             public :
 478                 virtual base_v<CalculationType, Parameters>* create_new(const Parameters& par) const
 479                 {
 480                     return new base_v_fi<gs48_ellipsoid<CalculationType, Parameters>, CalculationType, Parameters>(par);
 481                 }
 482         };
 483
 484         template <typename CalculationType, typename Parameters>
 485         class alsk_entry : public detail::factory_entry<CalculationType, Parameters>
 486         {
 487             public :
 488                 virtual base_v<CalculationType, Parameters>* create_new(const Parameters& par) const
 489                 {
 490                     return new base_v_fi<alsk_ellipsoid<CalculationType, Parameters>, CalculationType, Parameters>(par);
 491                 }
 492         };
 493
 494         template <typename CalculationType, typename Parameters>
 495         class gs50_entry : public detail::factory_entry<CalculationType, Parameters>
 496         {
 497             public :
 498                 virtual base_v<CalculationType, Parameters>* create_new(const Parameters& par) const
 499                 {
 500                     return new base_v_fi<gs50_ellipsoid<CalculationType, Parameters>, CalculationType, Parameters>(par);
 501                 }
 502         };
 503
 504         template <typename CalculationType, typename Parameters>
 505         inline void mod_ster_init(detail::base_factory<CalculationType, Parameters>& factory)
 506         {
 507             factory.add_to_factory("mil_os", new mil_os_entry<CalculationType, Parameters>);
 508             factory.add_to_factory("lee_os", new lee_os_entry<CalculationType, Parameters>);
 509             factory.add_to_factory("gs48", new gs48_entry<CalculationType, Parameters>);
 510             factory.add_to_factory("alsk", new alsk_entry<CalculationType, Parameters>);
 511             factory.add_to_factory("gs50", new gs50_entry<CalculationType, Parameters>);
 512         }
 513
 514     } // namespace detail
 515     #endif // doxygen
 516
 517 } // namespace projections
 518
 519 }} // namespace boost::geometry
 520
 521 #endif // BOOST_GEOMETRY_PROJECTIONS_MOD_STER_HPP
 522