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

   1 #ifndef BOOST_GEOMETRY_PROJECTIONS_KROVAK_HPP
   2 #define BOOST_GEOMETRY_PROJECTIONS_KROVAK_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 // Purpose:  Implementation of the krovak (Krovak) projection.
  27 // Definition: http://www.ihsenergy.com/epsg/guid7.html#1.4.3
  28 // Author:   Thomas Flemming, tf@ttqv.com
  29 // Copyright (c) 2001, Thomas Flemming, tf@ttqv.com
  30
  31 // Permission is hereby granted, free of charge, to any person obtaining a
  32 // copy of this software and associated documentation files (the "Software"),
  33 // to deal in the Software without restriction, including without limitation
  34 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  35 // and/or sell copies of the Software, and to permit persons to whom the
  36 // Software is furnished to do so, subject to the following conditions:
  37
  38 // The above copyright notice and this permission notice shall be included
  39 // in all copies or substantial portions of the Software.
  40
  41 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  42 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  43 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  44 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  45 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  46 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  47 // DEALINGS IN THE SOFTWARE.
  48
  49 #include <boost/geometry/srs/projections/impl/base_static.hpp>
  50 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
  51 #include <boost/geometry/srs/projections/impl/projects.hpp>
  52 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
  53
  54 namespace boost { namespace geometry
  55 {
  56
  57 namespace srs { namespace par4
  58 {
  59     struct krovak {};
  60
  61 }} //namespace srs::par4
  62
  63 namespace projections
  64 {
  65     #ifndef DOXYGEN_NO_DETAIL
  66     namespace detail { namespace krovak
  67     {
  68             template <typename T>
  69             struct par_krovak
  70             {
  71                 T    C_x;
  72             };
  73
  74             /**
  75                NOTES: According to EPSG the full Krovak projection method should have
  76                       the following parameters.  Within PROJ.4 the azimuth, and pseudo
  77                       standard parallel are hardcoded in the algorithm and can't be
  78                       altered from outside.  The others all have defaults to match the
  79                       common usage with Krovak projection.
  80
  81               lat_0 = latitude of centre of the projection
  82
  83               lon_0 = longitude of centre of the projection
  84
  85               ** = azimuth (true) of the centre line passing through the centre of the projection
  86
  87               ** = latitude of pseudo standard parallel
  88
  89               k  = scale factor on the pseudo standard parallel
  90
  91               x_0 = False Easting of the centre of the projection at the apex of the cone
  92
  93               y_0 = False Northing of the centre of the projection at the apex of the cone
  94
  95              **/
  96
  97             // template class, using CRTP to implement forward/inverse
  98             template <typename CalculationType, typename Parameters>
  99             struct base_krovak_ellipsoid : public base_t_fi<base_krovak_ellipsoid<CalculationType, Parameters>,
 100                      CalculationType, Parameters>
 101             {
 102
 103                 typedef CalculationType geographic_type;
 104                 typedef CalculationType cartesian_type;
 105
 106                 par_krovak<CalculationType> m_proj_parm;
 107
 108                 inline base_krovak_ellipsoid(const Parameters& par)
 109                     : base_t_fi<base_krovak_ellipsoid<CalculationType, Parameters>,
 110                      CalculationType, Parameters>(*this, par) {}
 111
 112                 // FORWARD(e_forward)  ellipsoid
 113                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
 114                 inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
 115                 {
 116                 /* calculate xy from lat/lon */
 117
 118                 /* Constants, identical to inverse transform function */
 119                     CalculationType s45, s90, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a, s0, n;
 120                     CalculationType gfi, u, fi0, deltav, s, d, eps, ro;
 121
 122
 123                     s45 = 0.785398163397448;    /* 45 DEG */
 124                     s90 = 2 * s45;
 125                     fi0 = this->m_par.phi0;    /* Latitude of projection centre 49 DEG 30' */
 126
 127                    /* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must
 128                       be set to 1 here.
 129                       Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.1528128,
 130                       e2=0.006674372230614;
 131                    */
 132                     a =  1; /* 6377397.155; */
 133                     /* e2 = this->m_par.es;*/       /* 0.006674372230614; */
 134                     e2 = 0.006674372230614;
 135                     e = sqrt(e2);
 136
 137                     alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2));
 138
 139                     uq = 1.04216856380474;      /* DU(2, 59, 42, 42.69689) */
 140                     u0 = asin(sin(fi0) / alfa);
 141                     g = pow(   (1. + e * sin(fi0)) / (1. - e * sin(fi0)) , alfa * e / 2.  );
 142
 143                     k = tan( u0 / 2. + s45) / pow  (tan(fi0 / 2. + s45) , alfa) * g;
 144
 145                     k1 = this->m_par.k0;
 146                     n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2));
 147                     s0 = 1.37008346281555;       /* Latitude of pseudo standard parallel 78 DEG 30'00" N */
 148                     n = sin(s0);
 149                     ro0 = k1 * n0 / tan(s0);
 150                     ad = s90 - uq;
 151
 152                 /* Transformation */
 153
 154                     gfi =pow ( ((1. + e * sin(lp_lat)) /
 155                                (1. - e * sin(lp_lat))) , (alfa * e / 2.));
 156
 157                     u= 2. * (atan(k * pow( tan(lp_lat / 2. + s45), alfa) / gfi)-s45);
 158
 159                     deltav = - lp_lon * alfa;
 160
 161                     s = asin(cos(ad) * sin(u) + sin(ad) * cos(u) * cos(deltav));
 162                     d = asin(cos(u) * sin(deltav) / cos(s));
 163                     eps = n * d;
 164                     ro = ro0 * pow(tan(s0 / 2. + s45) , n) / pow(tan(s / 2. + s45) , n)   ;
 165
 166                    /* x and y are reverted! */
 167                     xy_y = ro * cos(eps) / a;
 168                     xy_x = ro * sin(eps) / a;
 169
 170                         if( !pj_param(this->m_par.params, "tczech").i )
 171                       {
 172                         xy_y *= -1.0;
 173                         xy_x *= -1.0;
 174                       }
 175                 }
 176
 177                 // INVERSE(e_inverse)  ellipsoid
 178                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 179                 inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
 180                 {
 181                     /* calculate lat/lon from xy */
 182
 183                 /* Constants, identisch wie in der Umkehrfunktion */
 184                     CalculationType s45, s90, fi0, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a, s0, n;
 185                     CalculationType u, deltav, s, d, eps, ro, fi1, xy0;
 186                     int ok;
 187
 188                     s45 = 0.785398163397448;    /* 45 DEG */
 189                     s90 = 2 * s45;
 190                     fi0 = this->m_par.phi0;    /* Latitude of projection centre 49 DEG 30' */
 191
 192
 193                    /* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must
 194                       be set to 1 here.
 195                       Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.1528128,
 196                       e2=0.006674372230614;
 197                    */
 198                     a = 1; /* 6377397.155; */
 199                     /* e2 = this->m_par.es; */      /* 0.006674372230614; */
 200                     e2 = 0.006674372230614;
 201                     e = sqrt(e2);
 202
 203                     alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2));
 204                     uq = 1.04216856380474;      /* DU(2, 59, 42, 42.69689) */
 205                     u0 = asin(sin(fi0) / alfa);
 206                     g = pow(   (1. + e * sin(fi0)) / (1. - e * sin(fi0)) , alfa * e / 2.  );
 207
 208                     k = tan( u0 / 2. + s45) / pow  (tan(fi0 / 2. + s45) , alfa) * g;
 209
 210                     k1 = this->m_par.k0;
 211                     n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2));
 212                     s0 = 1.37008346281555;       /* Latitude of pseudo standard parallel 78 DEG 30'00" N */
 213                     n = sin(s0);
 214                     ro0 = k1 * n0 / tan(s0);
 215                     ad = s90 - uq;
 216
 217
 218                 /* Transformation */
 219                    /* revert y, x*/
 220                     xy0=xy_x;
 221                     xy_x=xy_y;
 222                     xy_y=xy0;
 223
 224                         if( !pj_param(this->m_par.params, "tczech").i )
 225                       {
 226                         xy_x *= -1.0;
 227                         xy_y *= -1.0;
 228                       }
 229
 230                     ro = sqrt(xy_x * xy_x + xy_y * xy_y);
 231                     eps = atan2(xy_y, xy_x);
 232                     d = eps / sin(s0);
 233                     s = 2. * (atan(  pow(ro0 / ro, 1. / n) * tan(s0 / 2. + s45)) - s45);
 234
 235                     u = asin(cos(ad) * sin(s) - sin(ad) * cos(s) * cos(d));
 236                     deltav = asin(cos(s) * sin(d) / cos(u));
 237
 238                     lp_lon = this->m_par.lam0 - deltav / alfa;
 239
 240                 /* ITERATION FOR lp_lat */
 241                    fi1 = u;
 242
 243                    ok = 0;
 244                    do
 245                    {
 246                        lp_lat = 2. * ( atan( pow( k, -1. / alfa)  *
 247                                             pow( tan(u / 2. + s45) , 1. / alfa)  *
 248                                             pow( (1. + e * sin(fi1)) / (1. - e * sin(fi1)) , e / 2.)
 249                                            )  - s45);
 250
 251                       if (fabs(fi1 - lp_lat) < 0.000000000000001) ok=1;
 252                       fi1 = lp_lat;
 253
 254                    }
 255                    while (ok==0);
 256
 257                    lp_lon -= this->m_par.lam0;
 258                 }
 259
 260                 static inline std::string get_name()
 261                 {
 262                     return "krovak_ellipsoid";
 263                 }
 264
 265             };
 266
 267             // Krovak
 268             template <typename Parameters, typename T>
 269             inline void setup_krovak(Parameters& par, par_krovak<T>& proj_parm)
 270             {
 271                 T ts;
 272                 /* read some Parameters,
 273                  * here Latitude Truescale */
 274
 275                 ts = pj_param(par.params, "rlat_ts").f;
 276                 proj_parm.C_x = ts;
 277
 278                 /* we want Bessel as fixed ellipsoid */
 279                 par.a = 6377397.155;
 280                 par.e = sqrt(par.es = 0.006674372230614);
 281
 282                     /* if latitude of projection center is not set, use 49d30'N */
 283                 if (!pj_param(par.params, "tlat_0").i)
 284                         par.phi0 = 0.863937979737193;
 285
 286                     /* if center long is not set use 42d30'E of Ferro - 17d40' for Ferro */
 287                     /* that will correspond to using longitudes relative to greenwich    */
 288                     /* as input and output, instead of lat/long relative to Ferro */
 289                 if (!pj_param(par.params, "tlon_0").i)
 290                         par.lam0 = 0.7417649320975901 - 0.308341501185665;
 291
 292                     /* if scale not set default to 0.9999 */
 293                 if (!pj_param(par.params, "tk").i)
 294                         par.k0 = 0.9999;
 295
 296                 /* always the same */
 297             }
 298
 299     }} // namespace detail::krovak
 300     #endif // doxygen
 301
 302     /*!
 303         \brief Krovak projection
 304         \ingroup projections
 305         \tparam Geographic latlong point type
 306         \tparam Cartesian xy point type
 307         \tparam Parameters parameter type
 308         \par Projection characteristics
 309          - Pseudocylindrical
 310          - Ellipsoid
 311         \par Projection parameters
 312          - lat_ts: Latitude of true scale (degrees)
 313          - lat_0: Latitude of origin
 314          - lon_0: Central meridian
 315          - k: Scale factor on the pseudo standard parallel
 316         \par Example
 317         \image html ex_krovak.gif
 318     */
 319     template <typename CalculationType, typename Parameters>
 320     struct krovak_ellipsoid : public detail::krovak::base_krovak_ellipsoid<CalculationType, Parameters>
 321     {
 322         inline krovak_ellipsoid(const Parameters& par) : detail::krovak::base_krovak_ellipsoid<CalculationType, Parameters>(par)
 323         {
 324             detail::krovak::setup_krovak(this->m_par, this->m_proj_parm);
 325         }
 326     };
 327
 328     #ifndef DOXYGEN_NO_DETAIL
 329     namespace detail
 330     {
 331
 332         // Static projection
 333         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::krovak, krovak_ellipsoid, krovak_ellipsoid)
 334
 335         // Factory entry(s)
 336         template <typename CalculationType, typename Parameters>
 337         class krovak_entry : public detail::factory_entry<CalculationType, Parameters>
 338         {
 339             public :
 340                 virtual base_v<CalculationType, Parameters>* create_new(const Parameters& par) const
 341                 {
 342                     return new base_v_fi<krovak_ellipsoid<CalculationType, Parameters>, CalculationType, Parameters>(par);
 343                 }
 344         };
 345
 346         template <typename CalculationType, typename Parameters>
 347         inline void krovak_init(detail::base_factory<CalculationType, Parameters>& factory)
 348         {
 349             factory.add_to_factory("krovak", new krovak_entry<CalculationType, Parameters>);
 350         }
 351
 352     } // namespace detail
 353     #endif // doxygen
 354
 355 } // namespace projections
 356
 357 }} // namespace boost::geometry
 358
 359 #endif // BOOST_GEOMETRY_PROJECTIONS_KROVAK_HPP
 360