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/util/math.hpp> 44 45 46 namespace boost { namespace geometry { namespace projections 47 { 48 namespace detail 49 { 50 template <typename T> 51 struct mdist 52 { 53 static const int static_size = 20; 54 55 int nb; 56 T es; 57 T E; 58 T b[static_size]; 59 }; 60 61 template <typename T> proj_mdist_ini(T const & es,mdist<T> & b)62 inline bool proj_mdist_ini(T const& es, mdist<T>& b) 63 { 64 T numf, numfi, twon1, denf, denfi, ens, t, twon; 65 T den, El, Es; 66 T E[mdist<T>::static_size]; 67 int i, j; 68 69 /* generate E(e^2) and its terms E[] */ 70 ens = es; 71 numf = twon1 = denfi = 1.; 72 denf = 1.; 73 twon = 4.; 74 Es = El = E[0] = 1.; 75 for (i = 1; i < mdist<T>::static_size ; ++i) 76 { 77 numf *= (twon1 * twon1); 78 den = twon * denf * denf * twon1; 79 t = numf/den; 80 E[i] = t * ens; 81 Es -= E[i]; 82 ens *= es; 83 twon *= 4.; 84 denf *= ++denfi; 85 twon1 += 2.; 86 if (Es == El) /* jump out if no change */ 87 break; 88 El = Es; 89 } 90 b.nb = i - 1; 91 b.es = es; 92 b.E = Es; 93 /* generate b_n coefficients--note: collapse with prefix ratios */ 94 b.b[0] = Es = 1. - Es; 95 numf = denf = 1.; 96 numfi = 2.; 97 denfi = 3.; 98 for (j = 1; j < i; ++j) 99 { 100 Es -= E[j]; 101 numf *= numfi; 102 denf *= denfi; 103 b.b[j] = Es * numf / denf; 104 numfi += 2.; 105 denfi += 2.; 106 } 107 return true; 108 } 109 110 template <typename T> proj_mdist(T const & phi,T const & sphi,T const & cphi,mdist<T> const & b)111 inline T proj_mdist(T const& phi, T const& sphi, T const& cphi, mdist<T> const& b) 112 { 113 T sc, sum, sphi2, D; 114 int i; 115 116 sc = sphi * cphi; 117 sphi2 = sphi * sphi; 118 D = phi * b.E - b.es * sc / sqrt(1. - b.es * sphi2); 119 sum = b.b[i = b.nb]; 120 while (i) sum = b.b[--i] + sphi2 * sum; 121 return(D + sc * sum); 122 } 123 124 template <typename T> proj_inv_mdist(T const & dist,mdist<T> const & b)125 inline T proj_inv_mdist(T const& dist, mdist<T> const& b) 126 { 127 static const T TOL = 1e-14; 128 T s, t, phi, k; 129 int i; 130 131 k = 1./(1.- b.es); 132 i = mdist<T>::static_size; 133 phi = dist; 134 while ( i-- ) { 135 s = sin(phi); 136 t = 1. - b.es * s * s; 137 phi -= t = (proj_mdist(phi, s, cos(phi), b) - dist) * 138 (t * sqrt(t)) * k; 139 if (geometry::math::abs(t) < TOL) /* that is no change */ 140 return phi; 141 } 142 /* convergence failed */ 143 BOOST_THROW_EXCEPTION( projection_exception(error_non_conv_inv_meri_dist) ); 144 } 145 } // namespace detail 146 147 }}} // namespace boost::geometry::projections 148 149 #endif 150