]> git.proxmox.com Git - ceph.git/blobdiff - ceph/src/boost/boost/geometry/srs/projections/proj/omerc.hpp
import new upstream nautilus stable release 14.2.8
[ceph.git] / ceph / src / boost / boost / geometry / srs / projections / proj / omerc.hpp
index 4da6871d13b5277ce5dee3efae4a079c4b3da8b2..9b9f5170fd9e0292a09ece9df6a8aa3b66558bc4 100644 (file)
@@ -1,13 +1,9 @@
-#ifndef BOOST_GEOMETRY_PROJECTIONS_OMERC_HPP
-#define BOOST_GEOMETRY_PROJECTIONS_OMERC_HPP
-
-// Boost.Geometry - extensions-gis-projections (based on PROJ4)
-// This file is automatically generated. DO NOT EDIT.
+// Boost.Geometry - gis-projections (based on PROJ4)
 
 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
 
-// This file was modified by Oracle on 2017.
-// Modifications copyright (c) 2017, Oracle and/or its affiliates.
+// This file was modified by Oracle on 2017, 2018, 2019.
+// Modifications copyright (c) 2017-2019, Oracle and/or its affiliates.
 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
 
 // Use, modification and distribution is subject to the Boost Software License,
@@ -19,7 +15,7 @@
 // PROJ4 is maintained by Frank Warmerdam
 // PROJ4 is converted to Boost.Geometry by Barend Gehrels
 
-// Last updated version of proj: 4.9.1
+// Last updated version of proj: 5.0.0
 
 // Original copyright notice:
 
 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 // DEALINGS IN THE SOFTWARE.
 
+#ifndef BOOST_GEOMETRY_PROJECTIONS_OMERC_HPP
+#define BOOST_GEOMETRY_PROJECTIONS_OMERC_HPP
+
 #include <boost/geometry/util/math.hpp>
 
 #include <boost/geometry/srs/projections/impl/base_static.hpp>
 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
-#include <boost/geometry/srs/projections/impl/projects.hpp>
 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
+#include <boost/geometry/srs/projections/impl/pj_param.hpp>
 #include <boost/geometry/srs/projections/impl/pj_phi2.hpp>
 #include <boost/geometry/srs/projections/impl/pj_tsfn.hpp>
+#include <boost/geometry/srs/projections/impl/projects.hpp>
 
 namespace boost { namespace geometry
 {
 
-namespace srs { namespace par4
-{
-    struct omerc {};
-
-}} //namespace srs::par4
-
 namespace projections
 {
     #ifndef DOXYGEN_NO_DETAIL
     namespace detail { namespace omerc
     {
-            static const double TOL = 1.e-7;
-            static const double EPS = 1.e-10;
-
             template <typename T>
             struct par_omerc
             {
-                T   A, B, E, AB, ArB, BrA, rB, singam, cosgam, sinrot, cosrot;
-                T   v_pole_n, v_pole_s, u_0;
-                int no_rot;
+                T    A, B, E, AB, ArB, BrA, rB, singam, cosgam, sinrot, cosrot;
+                T    v_pole_n, v_pole_s, u_0;
+                bool no_rot;
             };
 
-            // template class, using CRTP to implement forward/inverse
-            template <typename CalculationType, typename Parameters>
-            struct base_omerc_ellipsoid : public base_t_fi<base_omerc_ellipsoid<CalculationType, Parameters>,
-                     CalculationType, Parameters>
-            {
-
-                typedef CalculationType geographic_type;
-                typedef CalculationType cartesian_type;
-
-                par_omerc<CalculationType> m_proj_parm;
+            static const double tolerance = 1.e-7;
+            static const double epsilon = 1.e-10;
 
-                inline base_omerc_ellipsoid(const Parameters& par)
-                    : base_t_fi<base_omerc_ellipsoid<CalculationType, Parameters>,
-                     CalculationType, Parameters>(*this, par) {}
+            template <typename T, typename Parameters>
+            struct base_omerc_ellipsoid
+            {
+                par_omerc<T> m_proj_parm;
 
                 // FORWARD(e_forward)  ellipsoid
                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
-                inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
+                inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
                 {
-                    static const CalculationType HALFPI = detail::HALFPI<CalculationType>();
+                    static const T half_pi = detail::half_pi<T>();
 
-                    CalculationType  Q, S, T, U, V, temp, u, v;
+                    T  s, t, U, V, W, temp, u, v;
 
-                    if (fabs(fabs(lp_lat) - HALFPI) > EPS) {
-                        Q = this->m_proj_parm.E / pow(pj_tsfn(lp_lat, sin(lp_lat), this->m_par.e), this->m_proj_parm.B);
-                        temp = 1. / Q;
-                        S = .5 * (Q - temp);
-                        T = .5 * (Q + temp);
+                    if (fabs(fabs(lp_lat) - half_pi) > epsilon) {
+                        W = this->m_proj_parm.E / math::pow(pj_tsfn(lp_lat, sin(lp_lat), par.e), this->m_proj_parm.B);
+                        temp = 1. / W;
+                        s = .5 * (W - temp);
+                        t = .5 * (W + temp);
                         V = sin(this->m_proj_parm.B * lp_lon);
-                        U = (S * this->m_proj_parm.singam - V * this->m_proj_parm.cosgam) / T;
-                        if (fabs(fabs(U) - 1.0) < EPS)
-                            BOOST_THROW_EXCEPTION( projection_exception(-20) );
+                        U = (s * this->m_proj_parm.singam - V * this->m_proj_parm.cosgam) / t;
+                        if (fabs(fabs(U) - 1.0) < epsilon) {
+                            BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
+                        }
                         v = 0.5 * this->m_proj_parm.ArB * log((1. - U)/(1. + U));
                         temp = cos(this->m_proj_parm.B * lp_lon);
-                                if(fabs(temp) < TOL) {
-                                    u = this->m_proj_parm.A * lp_lon;
-                                } else {
-                                    u = this->m_proj_parm.ArB * atan2((S * this->m_proj_parm.cosgam + V * this->m_proj_parm.singam), temp);
-                                }
+                        if(fabs(temp) < tolerance) {
+                            u = this->m_proj_parm.A * lp_lon;
+                        } else {
+                            u = this->m_proj_parm.ArB * atan2((s * this->m_proj_parm.cosgam + V * this->m_proj_parm.singam), temp);
+                        }
                     } else {
                         v = lp_lat > 0 ? this->m_proj_parm.v_pole_n : this->m_proj_parm.v_pole_s;
                         u = this->m_proj_parm.ArB * lp_lat;
@@ -132,11 +117,11 @@ namespace projections
 
                 // INVERSE(e_inverse)  ellipsoid
                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
-                inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
+                inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
                 {
-                    static const CalculationType HALFPI = detail::HALFPI<CalculationType>();
+                    static const T half_pi = detail::half_pi<T>();
 
-                    CalculationType  u, v, Qp, Sp, Tp, Vp, Up;
+                    T  u, v, Qp, Sp, Tp, Vp, Up;
 
                     if (this->m_proj_parm.no_rot) {
                         v = xy_y;
@@ -150,13 +135,14 @@ namespace projections
                     Tp = .5 * (Qp + 1. / Qp);
                     Vp = sin(this->m_proj_parm.BrA * u);
                     Up = (Vp * this->m_proj_parm.cosgam + Sp * this->m_proj_parm.singam) / Tp;
-                    if (fabs(fabs(Up) - 1.) < EPS) {
+                    if (fabs(fabs(Up) - 1.) < epsilon) {
                         lp_lon = 0.;
-                        lp_lat = Up < 0. ? -HALFPI : HALFPI;
+                        lp_lat = Up < 0. ? -half_pi : half_pi;
                     } else {
                         lp_lat = this->m_proj_parm.E / sqrt((1. + Up) / (1. - Up));
-                        if ((lp_lat = pj_phi2(pow(lp_lat, 1. / this->m_proj_parm.B), this->m_par.e)) == HUGE_VAL)
-                            BOOST_THROW_EXCEPTION( projection_exception(-20) );
+                        if ((lp_lat = pj_phi2(math::pow(lp_lat, T(1) / this->m_proj_parm.B), par.e)) == HUGE_VAL) {
+                            BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
+                        }
                         lp_lon = - this->m_proj_parm.rB * atan2((Sp * this->m_proj_parm.cosgam -
                             Vp * this->m_proj_parm.singam), cos(this->m_proj_parm.BrA * u));
                     }
@@ -170,50 +156,49 @@ namespace projections
             };
 
             // Oblique Mercator
-            template <typename Parameters, typename T>
-            inline void setup_omerc(Parameters& par, par_omerc<T>& proj_parm)
+            template <typename Params, typename Parameters, typename T>
+            inline void setup_omerc(Params const& params, Parameters & par, par_omerc<T>& proj_parm)
             {
-                static const T FORTPI = detail::FORTPI<T>();
-                static const T HALFPI = detail::HALFPI<T>();
-                static const T ONEPI = detail::ONEPI<T>();
-                static const T TWOPI = detail::TWOPI<T>();
+                static const T fourth_pi = detail::fourth_pi<T>();
+                static const T half_pi = detail::half_pi<T>();
+                static const T pi = detail::pi<T>();
+                static const T two_pi = detail::two_pi<T>();
 
                 T con, com, cosph0, D, F, H, L, sinph0, p, J, gamma=0,
-                  gamma0, lamc=0, lam1=0, lam2=0, phi1=0, phi2=0, alpha_c=0.0;
+                  gamma0, lamc=0, lam1=0, lam2=0, phi1=0, phi2=0, alpha_c=0;
                 int alp, gam, no_off = 0;
 
-                proj_parm.no_rot = pj_param(par.params, "tno_rot").i;
-                    if ((alp = pj_param(par.params, "talpha").i) != 0)
-                    alpha_c = pj_param(par.params, "ralpha").f;
-                    if ((gam = pj_param(par.params, "tgamma").i) != 0)
-                    gamma = pj_param(par.params, "rgamma").f;
+                proj_parm.no_rot = pj_get_param_b<srs::spar::no_rot>(params, "no_rot", srs::dpar::no_rot);
+                alp = pj_param_r<srs::spar::alpha>(params, "alpha", srs::dpar::alpha, alpha_c);
+                gam = pj_param_r<srs::spar::gamma>(params, "gamma", srs::dpar::gamma, gamma);
                 if (alp || gam) {
-                    lamc    = pj_param(par.params, "rlonc").f;
-                    no_off =
-                                /* For libproj4 compatability */
-                                pj_param(par.params, "tno_off").i
-                                /* for backward compatibility */
-                                || pj_param(par.params, "tno_uoff").i;
-                    if( no_off )
-                    {
-                        /* Mark the parameter as used, so that the pj_get_def() return them */
-                        pj_param(par.params, "sno_uoff");
-                        pj_param(par.params, "sno_off");
-                    }
+                    lamc = pj_get_param_r<T, srs::spar::lonc>(params, "lonc", srs::dpar::lonc);
+                    // NOTE: This is not needed in Boost.Geometry
+                    //no_off =
+                    //            /* For libproj4 compatability */
+                    //            pj_param_exists(par.params, "no_off")
+                    //            /* for backward compatibility */
+                    //            || pj_param_exists(par.params, "no_uoff");
+                    //if( no_off )
+                    //{
+                    //    /* Mark the parameter as used, so that the pj_get_def() return them */
+                    //    pj_get_param_s(par.params, "no_uoff");
+                    //    pj_get_param_s(par.params, "no_off");
+                    //}
                 } else {
-                    lam1 = pj_param(par.params, "rlon_1").f;
-                    phi1 = pj_param(par.params, "rlat_1").f;
-                    lam2 = pj_param(par.params, "rlon_2").f;
-                    phi2 = pj_param(par.params, "rlat_2").f;
-                    if (fabs(phi1 - phi2) <= TOL ||
-                        (con = fabs(phi1)) <= TOL ||
-                        fabs(con - HALFPI) <= TOL ||
-                        fabs(fabs(par.phi0) - HALFPI) <= TOL ||
-                        fabs(fabs(phi2) - HALFPI) <= TOL)
-                        BOOST_THROW_EXCEPTION( projection_exception(-33) );
+                    lam1 = pj_get_param_r<T, srs::spar::lon_1>(params, "lon_1", srs::dpar::lon_1);
+                    phi1 = pj_get_param_r<T, srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1);
+                    lam2 = pj_get_param_r<T, srs::spar::lon_2>(params, "lon_2", srs::dpar::lon_2);
+                    phi2 = pj_get_param_r<T, srs::spar::lat_2>(params, "lat_2", srs::dpar::lat_2);
+                    if (fabs(phi1 - phi2) <= tolerance ||
+                        (con = fabs(phi1)) <= tolerance ||
+                        fabs(con - half_pi) <= tolerance ||
+                        fabs(fabs(par.phi0) - half_pi) <= tolerance ||
+                        fabs(fabs(phi2) - half_pi) <= tolerance)
+                        BOOST_THROW_EXCEPTION( projection_exception(error_lat_0_or_alpha_eq_90) );
                 }
                 com = sqrt(par.one_es);
-                if (fabs(par.phi0) > EPS) {
+                if (fabs(par.phi0) > epsilon) {
                     sinph0 = sin(par.phi0);
                     cosph0 = cos(par.phi0);
                     con = 1. - par.es * sinph0 * sinph0;
@@ -229,7 +214,7 @@ namespace projections
                             F = -F;
                     }
                     proj_parm.E = F += D;
-                    proj_parm.E *= pow(pj_tsfn(par.phi0, sinph0, par.e), proj_parm.B);
+                    proj_parm.E *= math::pow(pj_tsfn(par.phi0, sinph0, par.e), proj_parm.B);
                 } else {
                     proj_parm.B = 1. / com;
                     proj_parm.A = par.k0;
@@ -237,33 +222,29 @@ namespace projections
                 }
                 if (alp || gam) {
                     if (alp) {
-                        gamma0 = asin(sin(alpha_c) / D);
+                        gamma0 = aasin(sin(alpha_c) / D);
                         if (!gam)
                             gamma = alpha_c;
                     } else
-                        alpha_c = asin(D*sin(gamma0 = gamma));
-                    if ((con = fabs(alpha_c)) <= TOL ||
-                        fabs(con - ONEPI) <= TOL ||
-                        fabs(fabs(par.phi0) - HALFPI) <= TOL)
-                        BOOST_THROW_EXCEPTION( projection_exception(-32) );
-                    par.lam0 = lamc - asin(.5 * (F - 1. / F) *
+                        alpha_c = aasin(D*sin(gamma0 = gamma));
+                    par.lam0 = lamc - aasin(.5 * (F - 1. / F) *
                        tan(gamma0)) / proj_parm.B;
                 } else {
-                    H = pow(pj_tsfn(phi1, sin(phi1), par.e), proj_parm.B);
-                    L = pow(pj_tsfn(phi2, sin(phi2), par.e), proj_parm.B);
+                    H = math::pow(pj_tsfn(phi1, sin(phi1), par.e), proj_parm.B);
+                    L = math::pow(pj_tsfn(phi2, sin(phi2), par.e), proj_parm.B);
                     F = proj_parm.E / H;
                     p = (L - H) / (L + H);
                     J = proj_parm.E * proj_parm.E;
                     J = (J - L * H) / (J + L * H);
-                    if ((con = lam1 - lam2) < -ONEPI)
-                        lam2 -= TWOPI;
-                    else if (con > ONEPI)
-                        lam2 += TWOPI;
+                    if ((con = lam1 - lam2) < -pi)
+                        lam2 -= two_pi;
+                    else if (con > pi)
+                        lam2 += two_pi;
                     par.lam0 = adjlon(.5 * (lam1 + lam2) - atan(
                        J * tan(.5 * proj_parm.B * (lam1 - lam2)) / p) / proj_parm.B);
                     gamma0 = atan(2. * sin(proj_parm.B * adjlon(lam1 - par.lam0)) /
                        (F - 1. / F));
-                    gamma = alpha_c = asin(D * sin(gamma0));
+                    gamma = alpha_c = aasin(D * sin(gamma0));
                 }
                 proj_parm.singam = sin(gamma0);
                 proj_parm.cosgam = cos(gamma0);
@@ -274,13 +255,13 @@ namespace projections
                 if (no_off)
                     proj_parm.u_0 = 0;
                 else {
-                    proj_parm.u_0 = fabs(proj_parm.ArB * atan2(sqrt(D * D - 1.), cos(alpha_c)));
+                    proj_parm.u_0 = fabs(proj_parm.ArB * atan(sqrt(D * D - 1.) / cos(alpha_c)));
                     if (par.phi0 < 0.)
                         proj_parm.u_0 = - proj_parm.u_0;
                 }
                 F = 0.5 * gamma0;
-                proj_parm.v_pole_n = proj_parm.ArB * log(tan(FORTPI - F));
-                proj_parm.v_pole_s = proj_parm.ArB * log(tan(FORTPI + F));
+                proj_parm.v_pole_n = proj_parm.ArB * log(tan(fourth_pi - F));
+                proj_parm.v_pole_s = proj_parm.ArB * log(tan(fourth_pi + F));
             }
 
     }} // namespace detail::omerc
@@ -310,12 +291,13 @@ namespace projections
         \par Example
         \image html ex_omerc.gif
     */
-    template <typename CalculationType, typename Parameters>
-    struct omerc_ellipsoid : public detail::omerc::base_omerc_ellipsoid<CalculationType, Parameters>
+    template <typename T, typename Parameters>
+    struct omerc_ellipsoid : public detail::omerc::base_omerc_ellipsoid<T, Parameters>
     {
-        inline omerc_ellipsoid(const Parameters& par) : detail::omerc::base_omerc_ellipsoid<CalculationType, Parameters>(par)
+        template <typename Params>
+        inline omerc_ellipsoid(Params const& params, Parameters & par)
         {
-            detail::omerc::setup_omerc(this->m_par, this->m_proj_parm);
+            detail::omerc::setup_omerc(params, par, this->m_proj_parm);
         }
     };
 
@@ -324,23 +306,14 @@ namespace projections
     {
 
         // Static projection
-        BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::omerc, omerc_ellipsoid, omerc_ellipsoid)
+        BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_omerc, omerc_ellipsoid)
 
         // Factory entry(s)
-        template <typename CalculationType, typename Parameters>
-        class omerc_entry : public detail::factory_entry<CalculationType, Parameters>
-        {
-            public :
-                virtual base_v<CalculationType, Parameters>* create_new(const Parameters& par) const
-                {
-                    return new base_v_fi<omerc_ellipsoid<CalculationType, Parameters>, CalculationType, Parameters>(par);
-                }
-        };
+        BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(omerc_entry, omerc_ellipsoid)
 
-        template <typename CalculationType, typename Parameters>
-        inline void omerc_init(detail::base_factory<CalculationType, Parameters>& factory)
+        BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(omerc_init)
         {
-            factory.add_to_factory("omerc", new omerc_entry<CalculationType, Parameters>);
+            BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(omerc, omerc_entry)
         }
 
     } // namespace detail