boost/geometry/srs/projections/impl/proj_mdist.hpp

   1 // Boost.Geometry (aka GGL, Generic Geometry Library)
   2 // This file is manually converted from PROJ4
   3
   4 // Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands.
   5
   6 // This file was modified by Oracle on 2018.
   7 // Modifications copyright (c) 2018, Oracle and/or its affiliates.
   8 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
   9
  10 // Use, modification and distribution is subject to the Boost Software License,
  11 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  12 // http://www.boost.org/LICENSE_1_0.txt)
  13
  14 // This file is converted from PROJ4, http://trac.osgeo.org/proj
  15 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
  16 // PROJ4 is maintained by Frank Warmerdam
  17 // PROJ4 is converted to Geometry Library by Barend Gehrels (Geodan, Amsterdam)
  18
  19 // Original copyright notice:
  20
  21 // Permission is hereby granted, free of charge, to any person obtaining a
  22 // copy of this software and associated documentation files (the "Software"),
  23 // to deal in the Software without restriction, including without limitation
  24 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  25 // and/or sell copies of the Software, and to permit persons to whom the
  26 // Software is furnished to do so, subject to the following conditions:
  27
  28 // The above copyright notice and this permission notice shall be included
  29 // in all copies or substantial portions of the Software.
  30
  31 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  32 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  33 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  34 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  35 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  36 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  37 // DEALINGS IN THE SOFTWARE.
  38
  39 #ifndef BOOST_GEOMETRY_PROJECTIONS_PROJ_MDIST_HPP
  40 #define BOOST_GEOMETRY_PROJECTIONS_PROJ_MDIST_HPP
  41
  42
  43 #include <boost/geometry/srs/projections/exception.hpp>
  44 #include <boost/geometry/srs/projections/impl/pj_strerrno.hpp>
  45 #include <boost/geometry/util/math.hpp>
  46
  47
  48 namespace boost { namespace geometry { namespace projections
  49 {
  50 namespace detail
  51 {
  52     template <typename T>
  53     struct mdist
  54     {
  55         static const int static_size = 20;
  56
  57         T es;
  58         T E;
  59         T b[static_size];
  60         int nb;
  61     };
  62
  63     template <typename T>
  64     inline bool proj_mdist_ini(T const& es, mdist<T>& b)
  65     {
  66         T numf, numfi, twon1, denf, denfi, ens, t, twon;
  67         T den, El, Es;
  68         T E[mdist<T>::static_size];
  69         int i, j;
  70
  71         /* generate E(e^2) and its terms E[] */
  72         ens = es;
  73         numf = twon1 = denfi = 1.;
  74         denf = 1.;
  75         twon = 4.;
  76         Es = El = E[0] = 1.;
  77         for (i = 1; i < mdist<T>::static_size ; ++i)
  78         {
  79             numf *= (twon1 * twon1);
  80             den = twon * denf * denf * twon1;
  81             t = numf/den;
  82             E[i] = t * ens;
  83             Es -= E[i];
  84             ens *= es;
  85             twon *= 4.;
  86             denf *= ++denfi;
  87             twon1 += 2.;
  88             if (Es == El) /* jump out if no change */
  89                 break;
  90             El = Es;
  91         }
  92         b.nb = i - 1;
  93         b.es = es;
  94         b.E = Es;
  95         /* generate b_n coefficients--note: collapse with prefix ratios */
  96         b.b[0] = Es = 1. - Es;
  97         numf = denf = 1.;
  98         numfi = 2.;
  99         denfi = 3.;
 100         for (j = 1; j < i; ++j)
 101         {
 102             Es -= E[j];
 103             numf *= numfi;
 104             denf *= denfi;
 105             b.b[j] = Es * numf / denf;
 106             numfi += 2.;
 107             denfi += 2.;
 108         }
 109         return true;
 110     }
 111
 112     template <typename T>
 113     inline T proj_mdist(T const& phi, T const& sphi, T const& cphi, mdist<T> const& b)
 114     {
 115         T sc, sum, sphi2, D;
 116         int i;
 117
 118         sc = sphi * cphi;
 119         sphi2 = sphi * sphi;
 120         D = phi * b.E - b.es * sc / sqrt(1. - b.es * sphi2);
 121         sum = b.b[i = b.nb];
 122         while (i) sum = b.b[--i] + sphi2 * sum;
 123         return(D + sc * sum);
 124     }
 125
 126     template <typename T>
 127     inline T proj_inv_mdist(T const& dist, mdist<T> const& b)
 128     {
 129         static const T TOL = 1e-14;
 130         T s, t, phi, k;
 131         int i;
 132
 133         k = 1./(1.- b.es);
 134         i = mdist<T>::static_size;
 135         phi = dist;
 136         while ( i-- ) {
 137             s = sin(phi);
 138             t = 1. - b.es * s * s;
 139             phi -= t = (proj_mdist(phi, s, cos(phi), b) - dist) *
 140                 (t * sqrt(t)) * k;
 141             if (geometry::math::abs(t) < TOL) /* that is no change */
 142                 return phi;
 143         }
 144             /* convergence failed */
 145         BOOST_THROW_EXCEPTION( projection_exception(error_non_conv_inv_meri_dist) );
 146     }
 147 } // namespace detail
 148
 149 }}} // namespace boost::geometry::projections
 150
 151 #endif