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