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

   1 #ifndef BOOST_GEOMETRY_PROJECTIONS_LSAT_HPP
   2 #define BOOST_GEOMETRY_PROJECTIONS_LSAT_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, 2018.
  10 // Modifications copyright (c) 2017-2018, 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
  46 #include <boost/geometry/srs/projections/impl/base_static.hpp>
  47 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
  48 #include <boost/geometry/srs/projections/impl/projects.hpp>
  49 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
  50 #include <boost/geometry/srs/projections/impl/aasincos.hpp>
  51
  52 namespace boost { namespace geometry
  53 {
  54
  55 namespace srs { namespace par4
  56 {
  57     struct lsat {};
  58
  59 }} //namespace srs::par4
  60
  61 namespace projections
  62 {
  63     #ifndef DOXYGEN_NO_DETAIL
  64     namespace detail { namespace lsat
  65     {
  66             static const double TOL = 1e-7;
  67             //static const double PI_HALFPI = 4.71238898038468985766;
  68             //static const double TWOPI_HALFPI = 7.85398163397448309610;
  69
  70             template <typename T>
  71             struct par_lsat
  72             {
  73                 T a2, a4, b, c1, c3;
  74                 T q, t, u, w, p22, sa, ca, xj, rlm, rlm2;
  75             };
  76
  77             /* based upon Snyder and Linck, USGS-NMD */
  78             template <typename T>
  79             inline void
  80             seraz0(T lam, T const& mult, par_lsat<T>& proj_parm)
  81             {
  82                 T sdsq, h, s, fc, sd, sq, d__1;
  83
  84                 lam *= geometry::math::d2r<T>();
  85                 sd = sin(lam);
  86                 sdsq = sd * sd;
  87                 s = proj_parm.p22 * proj_parm.sa * cos(lam) * sqrt((1. + proj_parm.t * sdsq) / ((
  88                     1. + proj_parm.w * sdsq) * (1. + proj_parm.q * sdsq)));
  89                 d__1 = 1. + proj_parm.q * sdsq;
  90                 h = sqrt((1. + proj_parm.q * sdsq) / (1. + proj_parm.w * sdsq)) * ((1. +
  91                     proj_parm.w * sdsq) / (d__1 * d__1) - proj_parm.p22 * proj_parm.ca);
  92                 sq = sqrt(proj_parm.xj * proj_parm.xj + s * s);
  93                 proj_parm.b += fc = mult * (h * proj_parm.xj - s * s) / sq;
  94                 proj_parm.a2 += fc * cos(lam + lam);
  95                 proj_parm.a4 += fc * cos(lam * 4.);
  96                 fc = mult * s * (h + proj_parm.xj) / sq;
  97                 proj_parm.c1 += fc * cos(lam);
  98                 proj_parm.c3 += fc * cos(lam * 3.);
  99             }
 100
 101             // template class, using CRTP to implement forward/inverse
 102             template <typename CalculationType, typename Parameters>
 103             struct base_lsat_ellipsoid : public base_t_fi<base_lsat_ellipsoid<CalculationType, Parameters>,
 104                      CalculationType, Parameters>
 105             {
 106
 107                 typedef CalculationType geographic_type;
 108                 typedef CalculationType cartesian_type;
 109
 110                 par_lsat<CalculationType> m_proj_parm;
 111
 112                 inline base_lsat_ellipsoid(const Parameters& par)
 113                     : base_t_fi<base_lsat_ellipsoid<CalculationType, Parameters>,
 114                      CalculationType, Parameters>(*this, par) {}
 115
 116                 // FORWARD(e_forward)  ellipsoid
 117                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
 118                 inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
 119                 {
 120                     static const CalculationType FORTPI = detail::FORTPI<CalculationType>();
 121                     static const CalculationType HALFPI = detail::HALFPI<CalculationType>();
 122                     static const CalculationType PI_HALFPI = detail::PI_HALFPI<CalculationType>();
 123                     static const CalculationType TWOPI_HALFPI = detail::TWOPI_HALFPI<CalculationType>();
 124
 125                     int l, nn;
 126                     CalculationType lamt, xlam, sdsq, c, d, s, lamdp, phidp, lampp, tanph,
 127                                     lamtp, cl, sd, sp, fac, sav, tanphi;
 128
 129                     if (lp_lat > HALFPI)
 130                         lp_lat = HALFPI;
 131                     else if (lp_lat < -HALFPI)
 132                         lp_lat = -HALFPI;
 133                     lampp = lp_lat >= 0. ? HALFPI : PI_HALFPI;
 134                     tanphi = tan(lp_lat);
 135                     for (nn = 0;;) {
 136                         sav = lampp;
 137                         lamtp = lp_lon + this->m_proj_parm.p22 * lampp;
 138                         cl = cos(lamtp);
 139                         if (fabs(cl) < TOL)
 140                             lamtp -= TOL;
 141                         fac = lampp - sin(lampp) * (cl < 0. ? -HALFPI : HALFPI);
 142                         for (l = 50; l; --l) {
 143                             lamt = lp_lon + this->m_proj_parm.p22 * sav;
 144                             if (fabs(c = cos(lamt)) < TOL)
 145                                 lamt -= TOL;
 146                             xlam = (this->m_par.one_es * tanphi * this->m_proj_parm.sa + sin(lamt) * this->m_proj_parm.ca) / c;
 147                             lamdp = atan(xlam) + fac;
 148                             if (fabs(fabs(sav) - fabs(lamdp)) < TOL)
 149                                 break;
 150                             sav = lamdp;
 151                         }
 152                         if (!l || ++nn >= 3 || (lamdp > this->m_proj_parm.rlm && lamdp < this->m_proj_parm.rlm2))
 153                             break;
 154                         if (lamdp <= this->m_proj_parm.rlm)
 155                             lampp = TWOPI_HALFPI;
 156                         else if (lamdp >= this->m_proj_parm.rlm2)
 157                             lampp = HALFPI;
 158                     }
 159                     if (l) {
 160                         sp = sin(lp_lat);
 161                         phidp = aasin((this->m_par.one_es * this->m_proj_parm.ca * sp - this->m_proj_parm.sa * cos(lp_lat) *
 162                             sin(lamt)) / sqrt(1. - this->m_par.es * sp * sp));
 163                         tanph = log(tan(FORTPI + .5 * phidp));
 164                         sd = sin(lamdp);
 165                         sdsq = sd * sd;
 166                         s = this->m_proj_parm.p22 * this->m_proj_parm.sa * cos(lamdp) * sqrt((1. + this->m_proj_parm.t * sdsq)
 167                              / ((1. + this->m_proj_parm.w * sdsq) * (1. + this->m_proj_parm.q * sdsq)));
 168                         d = sqrt(this->m_proj_parm.xj * this->m_proj_parm.xj + s * s);
 169                         xy_x = this->m_proj_parm.b * lamdp + this->m_proj_parm.a2 * sin(2. * lamdp) + this->m_proj_parm.a4 *
 170                             sin(lamdp * 4.) - tanph * s / d;
 171                         xy_y = this->m_proj_parm.c1 * sd + this->m_proj_parm.c3 * sin(lamdp * 3.) + tanph * this->m_proj_parm.xj / d;
 172                     } else
 173                         xy_x = xy_y = HUGE_VAL;
 174                 }
 175
 176                 // INVERSE(e_inverse)  ellipsoid
 177                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 178                 inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
 179                 {
 180                     static const CalculationType FORTPI = detail::FORTPI<CalculationType>();
 181                     static const CalculationType HALFPI = detail::HALFPI<CalculationType>();
 182
 183                     int nn;
 184                     CalculationType lamt, sdsq, s, lamdp, phidp, sppsq, dd, sd, sl, fac, scl, sav, spp;
 185
 186                     lamdp = xy_x / this->m_proj_parm.b;
 187                     nn = 50;
 188                     do {
 189                         sav = lamdp;
 190                         sd = sin(lamdp);
 191                         sdsq = sd * sd;
 192                         s = this->m_proj_parm.p22 * this->m_proj_parm.sa * cos(lamdp) * sqrt((1. + this->m_proj_parm.t * sdsq)
 193                              / ((1. + this->m_proj_parm.w * sdsq) * (1. + this->m_proj_parm.q * sdsq)));
 194                         lamdp = xy_x + xy_y * s / this->m_proj_parm.xj - this->m_proj_parm.a2 * sin(
 195                             2. * lamdp) - this->m_proj_parm.a4 * sin(lamdp * 4.) - s / this->m_proj_parm.xj * (
 196                             this->m_proj_parm.c1 * sin(lamdp) + this->m_proj_parm.c3 * sin(lamdp * 3.));
 197                         lamdp /= this->m_proj_parm.b;
 198                     } while (fabs(lamdp - sav) >= TOL && --nn);
 199                     sl = sin(lamdp);
 200                     fac = exp(sqrt(1. + s * s / this->m_proj_parm.xj / this->m_proj_parm.xj) * (xy_y -
 201                         this->m_proj_parm.c1 * sl - this->m_proj_parm.c3 * sin(lamdp * 3.)));
 202                     phidp = 2. * (atan(fac) - FORTPI);
 203                     dd = sl * sl;
 204                     if (fabs(cos(lamdp)) < TOL)
 205                         lamdp -= TOL;
 206                     spp = sin(phidp);
 207                     sppsq = spp * spp;
 208                     lamt = atan(((1. - sppsq * this->m_par.rone_es) * tan(lamdp) *
 209                         this->m_proj_parm.ca - spp * this->m_proj_parm.sa * sqrt((1. + this->m_proj_parm.q * dd) * (
 210                         1. - sppsq) - sppsq * this->m_proj_parm.u) / cos(lamdp)) / (1. - sppsq
 211                         * (1. + this->m_proj_parm.u)));
 212                     sl = lamt >= 0. ? 1. : -1.;
 213                     scl = cos(lamdp) >= 0. ? 1. : -1;
 214                     lamt -= HALFPI * (1. - scl) * sl;
 215                     lp_lon = lamt - this->m_proj_parm.p22 * lamdp;
 216                     if (fabs(this->m_proj_parm.sa) < TOL)
 217                         lp_lat = aasin(spp / sqrt(this->m_par.one_es * this->m_par.one_es + this->m_par.es * sppsq));
 218                     else
 219                         lp_lat = atan((tan(lamdp) * cos(lamt) - this->m_proj_parm.ca * sin(lamt)) /
 220                             (this->m_par.one_es * this->m_proj_parm.sa));
 221                 }
 222
 223                 static inline std::string get_name()
 224                 {
 225                     return "lsat_ellipsoid";
 226                 }
 227
 228             };
 229
 230             // Space oblique for LANDSAT
 231             template <typename Parameters, typename T>
 232             inline void setup_lsat(Parameters& par, par_lsat<T>& proj_parm)
 233             {
 234                 int land, path;
 235                 T lam, alf, esc, ess;
 236
 237                 land = pj_param(par.params, "ilsat").i;
 238                 if (land <= 0 || land > 5)
 239                     BOOST_THROW_EXCEPTION( projection_exception(-28) );
 240                 path = pj_param(par.params, "ipath").i;
 241                 if (path <= 0 || path > (land <= 3 ? 251 : 233))
 242                     BOOST_THROW_EXCEPTION( projection_exception(-29) );
 243                 if (land <= 3) {
 244                     par.lam0 = geometry::math::d2r<T>() * 128.87 - geometry::math::two_pi<T>() / 251. * path;
 245                     proj_parm.p22 = 103.2669323;
 246                     alf = geometry::math::d2r<T>() * 99.092;
 247                 } else {
 248                     par.lam0 = geometry::math::d2r<T>() * 129.3 - geometry::math::two_pi<T>() / 233. * path;
 249                     proj_parm.p22 = 98.8841202;
 250                     alf = geometry::math::d2r<T>() * 98.2;
 251                 }
 252                 proj_parm.p22 /= 1440.;
 253                 proj_parm.sa = sin(alf);
 254                 proj_parm.ca = cos(alf);
 255                 if (fabs(proj_parm.ca) < 1e-9)
 256                     proj_parm.ca = 1e-9;
 257                 esc = par.es * proj_parm.ca * proj_parm.ca;
 258                 ess = par.es * proj_parm.sa * proj_parm.sa;
 259                 proj_parm.w = (1. - esc) * par.rone_es;
 260                 proj_parm.w = proj_parm.w * proj_parm.w - 1.;
 261                 proj_parm.q = ess * par.rone_es;
 262                 proj_parm.t = ess * (2. - par.es) * par.rone_es * par.rone_es;
 263                 proj_parm.u = esc * par.rone_es;
 264                 proj_parm.xj = par.one_es * par.one_es * par.one_es;
 265                 proj_parm.rlm = geometry::math::pi<T>() * (1. / 248. + .5161290322580645);
 266                 proj_parm.rlm2 = proj_parm.rlm + geometry::math::two_pi<T>();
 267                 proj_parm.a2 = proj_parm.a4 = proj_parm.b = proj_parm.c1 = proj_parm.c3 = 0.;
 268                 seraz0(0., 1., proj_parm);
 269                 for (lam = 9.; lam <= 81.0001; lam += 18.)
 270                     seraz0(lam, 4., proj_parm);
 271                 for (lam = 18; lam <= 72.0001; lam += 18.)
 272                     seraz0(lam, 2., proj_parm);
 273                 seraz0(90., 1., proj_parm);
 274                 proj_parm.a2 /= 30.;
 275                 proj_parm.a4 /= 60.;
 276                 proj_parm.b /= 30.;
 277                 proj_parm.c1 /= 15.;
 278                 proj_parm.c3 /= 45.;
 279             }
 280
 281     }} // namespace detail::lsat
 282     #endif // doxygen
 283
 284     /*!
 285         \brief Space oblique for LANDSAT projection
 286         \ingroup projections
 287         \tparam Geographic latlong point type
 288         \tparam Cartesian xy point type
 289         \tparam Parameters parameter type
 290         \par Projection characteristics
 291          - Cylindrical
 292          - Spheroid
 293          - Ellipsoid
 294         \par Projection parameters
 295          - lsat (integer)
 296          - path (integer)
 297         \par Example
 298         \image html ex_lsat.gif
 299     */
 300     template <typename CalculationType, typename Parameters>
 301     struct lsat_ellipsoid : public detail::lsat::base_lsat_ellipsoid<CalculationType, Parameters>
 302     {
 303         inline lsat_ellipsoid(const Parameters& par) : detail::lsat::base_lsat_ellipsoid<CalculationType, Parameters>(par)
 304         {
 305             detail::lsat::setup_lsat(this->m_par, this->m_proj_parm);
 306         }
 307     };
 308
 309     #ifndef DOXYGEN_NO_DETAIL
 310     namespace detail
 311     {
 312
 313         // Static projection
 314         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::lsat, lsat_ellipsoid, lsat_ellipsoid)
 315
 316         // Factory entry(s)
 317         template <typename CalculationType, typename Parameters>
 318         class lsat_entry : public detail::factory_entry<CalculationType, Parameters>
 319         {
 320             public :
 321                 virtual base_v<CalculationType, Parameters>* create_new(const Parameters& par) const
 322                 {
 323                     return new base_v_fi<lsat_ellipsoid<CalculationType, Parameters>, CalculationType, Parameters>(par);
 324                 }
 325         };
 326
 327         template <typename CalculationType, typename Parameters>
 328         inline void lsat_init(detail::base_factory<CalculationType, Parameters>& factory)
 329         {
 330             factory.add_to_factory("lsat", new lsat_entry<CalculationType, Parameters>);
 331         }
 332
 333     } // namespace detail
 334     #endif // doxygen
 335
 336 } // namespace projections
 337
 338 }} // namespace boost::geometry
 339
 340 #endif // BOOST_GEOMETRY_PROJECTIONS_LSAT_HPP
 341