boost/geometry/srs/projections/proj/nzmg.hpp

   1 // Boost.Geometry - gis-projections (based on PROJ4)
   2
   3 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
   4
   5 // This file was modified by Oracle on 2017, 2018, 2019.
   6 // Modifications copyright (c) 2017-2019, Oracle and/or its affiliates.
   7 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
   8
   9 // Use, modification and distribution is subject to the Boost Software License,
  10 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  11 // http://www.boost.org/LICENSE_1_0.txt)
  12
  13 // This file is converted from PROJ4, http://trac.osgeo.org/proj
  14 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
  15 // PROJ4 is maintained by Frank Warmerdam
  16 // PROJ4 is converted to Boost.Geometry by Barend Gehrels
  17
  18 // Last updated version of proj: 5.0.0
  19
  20 // Original copyright notice:
  21
  22 // Purpose:  Implementation of the nzmg (New Zealand Map Grid) projection.
  23 //           Very loosely based upon DMA code by Bradford W. Drew
  24 // Author:   Gerald Evenden
  25 // Copyright (c) 1995, Gerald Evenden
  26
  27 // Permission is hereby granted, free of charge, to any person obtaining a
  28 // copy of this software and associated documentation files (the "Software"),
  29 // to deal in the Software without restriction, including without limitation
  30 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  31 // and/or sell copies of the Software, and to permit persons to whom the
  32 // Software is furnished to do so, subject to the following conditions:
  33
  34 // The above copyright notice and this permission notice shall be included
  35 // in all copies or substantial portions of the Software.
  36
  37 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  38 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  39 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  40 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  41 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  42 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  43 // DEALINGS IN THE SOFTWARE.
  44
  45 #ifndef BOOST_GEOMETRY_PROJECTIONS_NZMG_HPP
  46 #define BOOST_GEOMETRY_PROJECTIONS_NZMG_HPP
  47
  48 #include <boost/geometry/util/math.hpp>
  49
  50 #include <boost/geometry/srs/projections/impl/base_static.hpp>
  51 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
  52 #include <boost/geometry/srs/projections/impl/projects.hpp>
  53 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
  54 #include <boost/geometry/srs/projections/impl/pj_zpoly1.hpp>
  55
  56 namespace boost { namespace geometry
  57 {
  58
  59 namespace projections
  60 {
  61     #ifndef DOXYGEN_NO_DETAIL
  62     namespace detail { namespace nzmg
  63     {
  64
  65             static const double epsilon = 1e-10;
  66             static const int Nbf = 5;
  67             static const int Ntpsi = 9;
  68             static const int Ntphi = 8;
  69
  70             template <typename T>
  71             inline T sec5_to_rad() { return 0.4848136811095359935899141023; }
  72             template <typename T>
  73             inline T rad_to_sec5() { return 2.062648062470963551564733573; }
  74
  75             template <typename T>
  76             inline const pj_complex<T> * bf()
  77             {
  78                 static const pj_complex<T> result[] = {
  79                     {.7557853228,    0.0},
  80                     {.249204646,    .003371507},
  81                     {-.001541739,    .041058560},
  82                     {-.10162907,    .01727609},
  83                     {-.26623489,    -.36249218},
  84                     {-.6870983,    -1.1651967}
  85                 };
  86                 return result;
  87             }
  88
  89             template <typename T>
  90             inline const T * tphi()
  91             {
  92                 static const T result[] = { 1.5627014243, .5185406398, -.03333098,
  93                                             -.1052906,   -.0368594,     .007317,
  94                                              .01220,      .00394,      -.0013 };
  95                 return result;
  96             }
  97             template <typename T>
  98             inline const T * tpsi()
  99             {
 100                 static const T result[] = { .6399175073, -.1358797613, .063294409, -.02526853, .0117879,
 101                                            -.0055161,     .0026906,   -.001333,     .00067,   -.00034 };
 102                 return result;
 103             }
 104
 105             template <typename T, typename Parameters>
 106             struct base_nzmg_ellipsoid
 107             {
 108                 // FORWARD(e_forward)  ellipsoid
 109                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
 110                 inline void fwd(Parameters const& par, T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const
 111                 {
 112                     static const T rad_to_sec5 = nzmg::rad_to_sec5<T>();
 113
 114                     pj_complex<T> p;
 115                     const T * C;
 116                     int i;
 117
 118                     lp_lat = (lp_lat - par.phi0) * rad_to_sec5;
 119                     for (p.r = *(C = tpsi<T>() + (i = Ntpsi)); i ; --i)
 120                         p.r = *--C + lp_lat * p.r;
 121                     p.r *= lp_lat;
 122                     p.i = lp_lon;
 123                     p = pj_zpoly1(p, bf<T>(), Nbf);
 124                     xy_x = p.i;
 125                     xy_y = p.r;
 126                 }
 127
 128                 // INVERSE(e_inverse)  ellipsoid
 129                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 130                 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
 131                 {
 132                     static const T sec5_to_rad = nzmg::sec5_to_rad<T>();
 133
 134                     int nn, i;
 135                     pj_complex<T> p, f, fp, dp;
 136                     T den;
 137                     const T* C;
 138
 139                     p.r = xy_y;
 140                     p.i = xy_x;
 141                     for (nn = 20; nn ;--nn) {
 142                         f = pj_zpolyd1(p, bf<T>(), Nbf, &fp);
 143                         f.r -= xy_y;
 144                         f.i -= xy_x;
 145                         den = fp.r * fp.r + fp.i * fp.i;
 146                         p.r += dp.r = -(f.r * fp.r + f.i * fp.i) / den;
 147                         p.i += dp.i = -(f.i * fp.r - f.r * fp.i) / den;
 148                         if ((fabs(dp.r) + fabs(dp.i)) <= epsilon)
 149                             break;
 150                     }
 151                     if (nn) {
 152                         lp_lon = p.i;
 153                         for (lp_lat = *(C = tphi<T>() + (i = Ntphi)); i ; --i)
 154                             lp_lat = *--C + p.r * lp_lat;
 155                         lp_lat = par.phi0 + p.r * lp_lat * sec5_to_rad;
 156                     } else
 157                         lp_lon = lp_lat = HUGE_VAL;
 158                 }
 159
 160                 static inline std::string get_name()
 161                 {
 162                     return "nzmg_ellipsoid";
 163                 }
 164
 165             };
 166
 167             // New Zealand Map Grid
 168             template <typename Parameters>
 169             inline void setup_nzmg(Parameters& par)
 170             {
 171                 typedef typename Parameters::type calc_t;
 172                 static const calc_t d2r = geometry::math::d2r<calc_t>();
 173
 174                 /* force to International major axis */
 175                 par.a = 6378388.0;
 176                 par.ra = 1. / par.a;
 177                 par.lam0 = 173. * d2r;
 178                 par.phi0 = -41. * d2r;
 179                 par.x0 = 2510000.;
 180                 par.y0 = 6023150.;
 181             }
 182
 183     }} // namespace detail::nzmg
 184     #endif // doxygen
 185
 186     /*!
 187         \brief New Zealand Map Grid projection
 188         \ingroup projections
 189         \tparam Geographic latlong point type
 190         \tparam Cartesian xy point type
 191         \tparam Parameters parameter type
 192         \par Projection characteristics
 193          - Fixed Earth
 194         \par Example
 195         \image html ex_nzmg.gif
 196     */
 197     template <typename T, typename Parameters>
 198     struct nzmg_ellipsoid : public detail::nzmg::base_nzmg_ellipsoid<T, Parameters>
 199     {
 200         template <typename Params>
 201         inline nzmg_ellipsoid(Params const& , Parameters & par)
 202         {
 203             detail::nzmg::setup_nzmg(par);
 204         }
 205     };
 206
 207     #ifndef DOXYGEN_NO_DETAIL
 208     namespace detail
 209     {
 210
 211         // Static projection
 212         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_nzmg, nzmg_ellipsoid)
 213
 214         // Factory entry(s)
 215         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(nzmg_entry, nzmg_ellipsoid)
 216
 217         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(nzmg_init)
 218         {
 219             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(nzmg, nzmg_entry)
 220         }
 221
 222     } // namespace detail
 223     #endif // doxygen
 224
 225 } // namespace projections
 226
 227 }} // namespace boost::geometry
 228
 229 #endif // BOOST_GEOMETRY_PROJECTIONS_NZMG_HPP
 230