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