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> proj_mdist_ini(T const & es,mdist<T> & b)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> proj_mdist(T const & phi,T const & sphi,T const & cphi,mdist<T> const & b)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> proj_inv_mdist(T const & dist,mdist<T> const & b)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 152