[ceph.git] / ceph / src / boost / boost / geometry / srs / projections / proj / stere.hpp

// Boost.Geometry - gis-projections (based on PROJ4)

// Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.

// 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,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)

// This file is converted from PROJ4, http://trac.osgeo.org/proj
// PROJ4 is originally written by Gerald Evenden (then of the USGS)
// PROJ4 is maintained by Frank Warmerdam
// PROJ4 is converted to Boost.Geometry by Barend Gehrels

// Last updated version of proj: 5.0.0

// Original copyright notice:

// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:

// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.

// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.

#ifndef BOOST_GEOMETRY_PROJECTIONS_STERE_HPP
#define BOOST_GEOMETRY_PROJECTIONS_STERE_HPP

#include <boost/config.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/math/special_functions/hypot.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/factory_entry.hpp>
#include <boost/geometry/srs/projections/impl/pj_param.hpp>
#include <boost/geometry/srs/projections/impl/pj_tsfn.hpp>
#include <boost/geometry/srs/projections/impl/projects.hpp>

namespace boost { namespace geometry
{

namespace projections
{
    #ifndef DOXYGEN_NO_DETAIL
    namespace detail { namespace stere
    {
            static const double epsilon10 = 1.e-10;
            static const double tolerance = 1.e-8;
            static const int n_iter = 8;
            static const double conv_tolerance = 1.e-10;

            enum mode_type {
                s_pole = 0,
                n_pole = 1,
                obliq  = 2,
                equit  = 3
            };

            template <typename T>
            struct par_stere
            {
                T   phits;
                T   sinX1;
                T   cosX1;
                T   akm1;
                mode_type mode;
            };

            template <typename T>
            inline T ssfn_(T const& phit, T sinphi, T const& eccen)
            {
                static const T half_pi = detail::half_pi<T>();

                sinphi *= eccen;
                return (tan (.5 * (half_pi + phit)) *
                   math::pow((T(1) - sinphi) / (T(1) + sinphi), T(0.5) * eccen));
            }

            template <typename T, typename Parameters>
            struct base_stere_ellipsoid
            {
                par_stere<T> m_proj_parm;

                // FORWARD(e_forward)  ellipsoid
                // Project coordinates from geographic (lon, lat) to cartesian (x, y)
                inline void fwd(Parameters const& par, T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const
                {
                    static const T half_pi = detail::half_pi<T>();

                    T coslam, sinlam, sinX=0.0, cosX=0.0, X, A = 0.0, sinphi;

                    coslam = cos(lp_lon);
                    sinlam = sin(lp_lon);
                    sinphi = sin(lp_lat);
                    if (this->m_proj_parm.mode == obliq || this->m_proj_parm.mode == equit) {
                        sinX = sin(X = 2. * atan(ssfn_(lp_lat, sinphi, par.e)) - half_pi);
                        cosX = cos(X);
                    }
                    switch (this->m_proj_parm.mode) {
                    case obliq:
                        A = this->m_proj_parm.akm1 / (this->m_proj_parm.cosX1 * (1. + this->m_proj_parm.sinX1 * sinX +
                           this->m_proj_parm.cosX1 * cosX * coslam));
                        xy_y = A * (this->m_proj_parm.cosX1 * sinX - this->m_proj_parm.sinX1 * cosX * coslam);
                        goto xmul; /* but why not just  xy.x = A * cosX; break;  ? */

                    case equit:
                        // TODO: calculate denominator once
                        /* avoid zero division */
                        if (1. + cosX * coslam == 0.0) {
                            xy_y = HUGE_VAL;
                        } else {
                            A = this->m_proj_parm.akm1 / (1. + cosX * coslam);
                            xy_y = A * sinX;
                        }
                xmul:
                        xy_x = A * cosX;
                        break;

                    case s_pole:
                        lp_lat = -lp_lat;
                        coslam = - coslam;
                        sinphi = -sinphi;
                        BOOST_FALLTHROUGH;
                    case n_pole:
                        xy_x = this->m_proj_parm.akm1 * pj_tsfn(lp_lat, sinphi, par.e);
                        xy_y = - xy_x * coslam;
                        break;
                    }

                    xy_x = xy_x * sinlam;
                }

                // INVERSE(e_inverse)  ellipsoid
                // Project coordinates from cartesian (x, y) to geographic (lon, lat)
                inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
                {
                    static const T half_pi = detail::half_pi<T>();

                    T cosphi, sinphi, tp=0.0, phi_l=0.0, rho, halfe=0.0, halfpi=0.0;
                    int i;

                    rho = boost::math::hypot(xy_x, xy_y);
                    switch (this->m_proj_parm.mode) {
                    case obliq:
                    case equit:
                        cosphi = cos( tp = 2. * atan2(rho * this->m_proj_parm.cosX1 , this->m_proj_parm.akm1) );
                        sinphi = sin(tp);
                        if( rho == 0.0 )
                            phi_l = asin(cosphi * this->m_proj_parm.sinX1);
                        else
                            phi_l = asin(cosphi * this->m_proj_parm.sinX1 + (xy_y * sinphi * this->m_proj_parm.cosX1 / rho));

                        tp = tan(.5 * (half_pi + phi_l));
                        xy_x *= sinphi;
                        xy_y = rho * this->m_proj_parm.cosX1 * cosphi - xy_y * this->m_proj_parm.sinX1* sinphi;
                        halfpi = half_pi;
                        halfe = .5 * par.e;
                        break;
                    case n_pole:
                        xy_y = -xy_y;
                        BOOST_FALLTHROUGH;
                    case s_pole:
                        phi_l = half_pi - 2. * atan(tp = - rho / this->m_proj_parm.akm1);
                        halfpi = -half_pi;
                        halfe = -.5 * par.e;
                        break;
                    }
                    for (i = n_iter; i--; phi_l = lp_lat) {
                        sinphi = par.e * sin(phi_l);
                        lp_lat = T(2) * atan(tp * math::pow((T(1)+sinphi)/(T(1)-sinphi), halfe)) - halfpi;
                        if (fabs(phi_l - lp_lat) < conv_tolerance) {
                            if (this->m_proj_parm.mode == s_pole)
                                lp_lat = -lp_lat;
                            lp_lon = (xy_x == 0. && xy_y == 0.) ? 0. : atan2(xy_x, xy_y);
                            return;
                        }
                    }
                    BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
                }

                static inline std::string get_name()
                {
                    return "stere_ellipsoid";
                }

            };

            template <typename T, typename Parameters>
            struct base_stere_spheroid
            {
                par_stere<T> m_proj_parm;

                // FORWARD(s_forward)  spheroid
                // Project coordinates from geographic (lon, lat) to cartesian (x, y)
                inline void fwd(Parameters const& , T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const
                {
                    static const T fourth_pi = detail::fourth_pi<T>();
                    static const T half_pi = detail::half_pi<T>();

                    T  sinphi, cosphi, coslam, sinlam;

                    sinphi = sin(lp_lat);
                    cosphi = cos(lp_lat);
                    coslam = cos(lp_lon);
                    sinlam = sin(lp_lon);
                    switch (this->m_proj_parm.mode) {
                    case equit:
                        xy_y = 1. + cosphi * coslam;
                        goto oblcon;
                    case obliq:
                        xy_y = 1. + this->m_proj_parm.sinX1 * sinphi + this->m_proj_parm.cosX1 * cosphi * coslam;
                oblcon:
                        if (xy_y <= epsilon10) {
                            BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
                        }
                        xy_x = (xy_y = this->m_proj_parm.akm1 / xy_y) * cosphi * sinlam;
                        xy_y *= (this->m_proj_parm.mode == equit) ? sinphi :
                           this->m_proj_parm.cosX1 * sinphi - this->m_proj_parm.sinX1 * cosphi * coslam;
                        break;
                    case n_pole:
                        coslam = - coslam;
                        lp_lat = - lp_lat;
                        BOOST_FALLTHROUGH;
                    case s_pole:
                        if (fabs(lp_lat - half_pi) < tolerance) {
                            BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
                        }
                        xy_x = sinlam * ( xy_y = this->m_proj_parm.akm1 * tan(fourth_pi + .5 * lp_lat) );
                        xy_y *= coslam;
                        break;
                    }
                }

                // INVERSE(s_inverse)  spheroid
                // Project coordinates from cartesian (x, y) to geographic (lon, lat)
                inline void inv(Parameters const& par, T const& xy_x, T xy_y, T& lp_lon, T& lp_lat) const
                {
                    T  c, rh, sinc, cosc;

                    sinc = sin(c = 2. * atan((rh = boost::math::hypot(xy_x, xy_y)) / this->m_proj_parm.akm1));
                    cosc = cos(c);
                    lp_lon = 0.;

                    switch (this->m_proj_parm.mode) {
                    case equit:
                        if (fabs(rh) <= epsilon10)
                            lp_lat = 0.;
                        else
                            lp_lat = asin(xy_y * sinc / rh);
                        if (cosc != 0. || xy_x != 0.)
                            lp_lon = atan2(xy_x * sinc, cosc * rh);
                        break;
                    case obliq:
                        if (fabs(rh) <= epsilon10)
                            lp_lat = par.phi0;
                        else
                            lp_lat = asin(cosc * this->m_proj_parm.sinX1 + xy_y * sinc * this->m_proj_parm.cosX1 / rh);
                        if ((c = cosc - this->m_proj_parm.sinX1 * sin(lp_lat)) != 0. || xy_x != 0.)
                            lp_lon = atan2(xy_x * sinc * this->m_proj_parm.cosX1, c * rh);
                        break;
                    case n_pole:
                        xy_y = -xy_y;
                        BOOST_FALLTHROUGH;
                    case s_pole:
                        if (fabs(rh) <= epsilon10)
                            lp_lat = par.phi0;
                        else
                            lp_lat = asin(this->m_proj_parm.mode == s_pole ? - cosc : cosc);
                        lp_lon = (xy_x == 0. && xy_y == 0.) ? 0. : atan2(xy_x, xy_y);
                        break;
                    }
                }

                static inline std::string get_name()
                {
                    return "stere_spheroid";
                }

            };

            template <typename Parameters, typename T>
            inline void setup(Parameters const& par, par_stere<T>& proj_parm)  /* general initialization */
            {
                static const T fourth_pi = detail::fourth_pi<T>();
                static const T half_pi = detail::half_pi<T>();

                T t;

                if (fabs((t = fabs(par.phi0)) - half_pi) < epsilon10)
                    proj_parm.mode = par.phi0 < 0. ? s_pole : n_pole;
                else
                    proj_parm.mode = t > epsilon10 ? obliq : equit;
                proj_parm.phits = fabs(proj_parm.phits);

                if (par.es != 0.0) {
                    T X;

                    switch (proj_parm.mode) {
                    case n_pole:
                    case s_pole:
                        if (fabs(proj_parm.phits - half_pi) < epsilon10)
                            proj_parm.akm1 = 2. * par.k0 /
                               sqrt(math::pow(T(1)+par.e,T(1)+par.e)*math::pow(T(1)-par.e,T(1)-par.e));
                        else {
                            proj_parm.akm1 = cos(proj_parm.phits) /
                               pj_tsfn(proj_parm.phits, t = sin(proj_parm.phits), par.e);
                            t *= par.e;
                            proj_parm.akm1 /= sqrt(1. - t * t);
                        }
                        break;
                    case equit:
                    case obliq:
                        t = sin(par.phi0);
                        X = 2. * atan(ssfn_(par.phi0, t, par.e)) - half_pi;
                        t *= par.e;
                        proj_parm.akm1 = 2. * par.k0 * cos(par.phi0) / sqrt(1. - t * t);
                        proj_parm.sinX1 = sin(X);
                        proj_parm.cosX1 = cos(X);
                        break;
                    }
                } else {
                    switch (proj_parm.mode) {
                    case obliq:
                        proj_parm.sinX1 = sin(par.phi0);
                        proj_parm.cosX1 = cos(par.phi0);
                        BOOST_FALLTHROUGH;
                    case equit:
                        proj_parm.akm1 = 2. * par.k0;
                        break;
                    case s_pole:
                    case n_pole:
                        proj_parm.akm1 = fabs(proj_parm.phits - half_pi) >= epsilon10 ?
                           cos(proj_parm.phits) / tan(fourth_pi - .5 * proj_parm.phits) :
                           2. * par.k0 ;
                        break;
                    }
                }
            }


            // Stereographic
            template <typename Params, typename Parameters, typename T>
            inline void setup_stere(Params const& params, Parameters const& par, par_stere<T>& proj_parm)
            {
                static const T half_pi = detail::half_pi<T>();

                if (! pj_param_r<srs::spar::lat_ts>(params, "lat_ts", srs::dpar::lat_ts, proj_parm.phits))
                    proj_parm.phits = half_pi;

                setup(par, proj_parm);
            }

            // Universal Polar Stereographic
            template <typename Params, typename Parameters, typename T>
            inline void setup_ups(Params const& params, Parameters& par, par_stere<T>& proj_parm)
            {
                static const T half_pi = detail::half_pi<T>();

                /* International Ellipsoid */
                par.phi0 = pj_get_param_b<srs::spar::south>(params, "south", srs::dpar::south) ? -half_pi: half_pi;
                if (par.es == 0.0) {
                    BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) );
                }
                par.k0 = .994;
                par.x0 = 2000000.;
                par.y0 = 2000000.;
                proj_parm.phits = half_pi;
                par.lam0 = 0.;

                setup(par, proj_parm);
            }

    }} // namespace detail::stere
    #endif // doxygen

    /*!
        \brief Stereographic projection
        \ingroup projections
        \tparam Geographic latlong point type
        \tparam Cartesian xy point type
        \tparam Parameters parameter type
        \par Projection characteristics
         - Azimuthal
         - Spheroid
         - Ellipsoid
        \par Projection parameters
         - lat_ts: Latitude of true scale (degrees)
        \par Example
        \image html ex_stere.gif
    */
    template <typename T, typename Parameters>
    struct stere_ellipsoid : public detail::stere::base_stere_ellipsoid<T, Parameters>
    {
        template <typename Params>
        inline stere_ellipsoid(Params const& params, Parameters const& par)
        {
            detail::stere::setup_stere(params, par, this->m_proj_parm);
        }
    };

    /*!
        \brief Stereographic projection
        \ingroup projections
        \tparam Geographic latlong point type
        \tparam Cartesian xy point type
        \tparam Parameters parameter type
        \par Projection characteristics
         - Azimuthal
         - Spheroid
         - Ellipsoid
        \par Projection parameters
         - lat_ts: Latitude of true scale (degrees)
        \par Example
        \image html ex_stere.gif
    */
    template <typename T, typename Parameters>
    struct stere_spheroid : public detail::stere::base_stere_spheroid<T, Parameters>
    {
        template <typename Params>
        inline stere_spheroid(Params const& params, Parameters const& par)
        {
            detail::stere::setup_stere(params, par, this->m_proj_parm);
        }
    };

    /*!
        \brief Universal Polar Stereographic projection
        \ingroup projections
        \tparam Geographic latlong point type
        \tparam Cartesian xy point type
        \tparam Parameters parameter type
        \par Projection characteristics
         - Azimuthal
         - Spheroid
         - Ellipsoid
        \par Projection parameters
         - south: Denotes southern hemisphere UTM zone (boolean)
        \par Example
        \image html ex_ups.gif
    */
    template <typename T, typename Parameters>
    struct ups_ellipsoid : public detail::stere::base_stere_ellipsoid<T, Parameters>
    {
        template <typename Params>
        inline ups_ellipsoid(Params const& params, Parameters & par)
        {
            detail::stere::setup_ups(params, par, this->m_proj_parm);
        }
    };

    /*!
        \brief Universal Polar Stereographic projection
        \ingroup projections
        \tparam Geographic latlong point type
        \tparam Cartesian xy point type
        \tparam Parameters parameter type
        \par Projection characteristics
         - Azimuthal
         - Spheroid
         - Ellipsoid
        \par Projection parameters
         - south: Denotes southern hemisphere UTM zone (boolean)
        \par Example
        \image html ex_ups.gif
    */
    template <typename T, typename Parameters>
    struct ups_spheroid : public detail::stere::base_stere_spheroid<T, Parameters>
    {
        template <typename Params>
        inline ups_spheroid(Params const& params, Parameters & par)
        {
            detail::stere::setup_ups(params, par, this->m_proj_parm);
        }
    };

    #ifndef DOXYGEN_NO_DETAIL
    namespace detail
    {

        // Static projection
        BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_stere, stere_spheroid, stere_ellipsoid)
        BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_ups, ups_spheroid, ups_ellipsoid)

        // Factory entry(s)
        BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(stere_entry, stere_spheroid, stere_ellipsoid)
        BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(ups_entry, ups_spheroid, ups_ellipsoid)
        
        BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(stere_init)
        {
            BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(stere, stere_entry)
            BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(ups, ups_entry)
        }

    } // namespace detail
    #endif // doxygen

} // namespace projections

}} // namespace boost::geometry

#endif // BOOST_GEOMETRY_PROJECTIONS_STERE_HPP
Commit	Line	Data
92f5a8d4	1	// Boost.Geometry - gis-projections (based on PROJ4)
11fdf7f2 TL	2
	3	// Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
	4
92f5a8d4 TL	5	// This file was modified by Oracle on 2017, 2018, 2019.
92f5a8d4 TL	6	// Modifications copyright (c) 2017-2019, Oracle and/or its affiliates.
11fdf7f2 TL	7	// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
	8
	9	// Use, modification and distribution is subject to the Boost Software License,
	10	// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
	11	// http://www.boost.org/LICENSE_1_0.txt)
	12
	13	// This file is converted from PROJ4, http://trac.osgeo.org/proj
	14	// PROJ4 is originally written by Gerald Evenden (then of the USGS)
	15	// PROJ4 is maintained by Frank Warmerdam
	16	// PROJ4 is converted to Boost.Geometry by Barend Gehrels
	17
92f5a8d4	18	// Last updated version of proj: 5.0.0
11fdf7f2 TL	19
	20	// Original copyright notice:
	21
	22	// Permission is hereby granted, free of charge, to any person obtaining a
	23	// copy of this software and associated documentation files (the "Software"),
	24	// to deal in the Software without restriction, including without limitation
	25	// the rights to use, copy, modify, merge, publish, distribute, sublicense,
	26	// and/or sell copies of the Software, and to permit persons to whom the
	27	// Software is furnished to do so, subject to the following conditions:
	28
	29	// The above copyright notice and this permission notice shall be included
	30	// in all copies or substantial portions of the Software.
	31
	32	// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
	33	// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
	34	// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
	35	// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
	36	// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
	37	// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
	38	// DEALINGS IN THE SOFTWARE.
	39
92f5a8d4 TL	40	#ifndef BOOST_GEOMETRY_PROJECTIONS_STERE_HPP
	41	#define BOOST_GEOMETRY_PROJECTIONS_STERE_HPP
	42
11fdf7f2 TL	43	#include <boost/config.hpp>
	44	#include <boost/geometry/util/math.hpp>
	45	#include <boost/math/special_functions/hypot.hpp>
	46
	47	#include <boost/geometry/srs/projections/impl/base_static.hpp>
	48	#include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
11fdf7f2	49	#include <boost/geometry/srs/projections/impl/factory_entry.hpp>
92f5a8d4	50	#include <boost/geometry/srs/projections/impl/pj_param.hpp>
11fdf7f2	51	#include <boost/geometry/srs/projections/impl/pj_tsfn.hpp>
92f5a8d4	52	#include <boost/geometry/srs/projections/impl/projects.hpp>
11fdf7f2 TL	53
	54	namespace boost { namespace geometry
	55	{
	56
11fdf7f2 TL	57	namespace projections
	58	{
	59	#ifndef DOXYGEN_NO_DETAIL
	60	namespace detail { namespace stere
	61	{
92f5a8d4 TL	62	static const double epsilon10 = 1.e-10;
	63	static const double tolerance = 1.e-8;
	64	static const int n_iter = 8;
	65	static const double conv_tolerance = 1.e-10;
	66
	67	enum mode_type {
	68	s_pole = 0,
	69	n_pole = 1,
	70	obliq = 2,
	71	equit = 3
	72	};
11fdf7f2 TL	73
	74	template <typename T>
	75	struct par_stere
	76	{
	77	T phits;
	78	T sinX1;
	79	T cosX1;
	80	T akm1;
92f5a8d4	81	mode_type mode;
11fdf7f2 TL	82	};
	83
	84	template <typename T>
	85	inline T ssfn_(T const& phit, T sinphi, T const& eccen)
	86	{
92f5a8d4 TL	87	static const T half_pi = detail::half_pi<T>();
92f5a8d4 TL	88
11fdf7f2	89	sinphi *= eccen;
92f5a8d4 TL	90	return (tan (.5 * (half_pi + phit)) *
92f5a8d4 TL	91	math::pow((T(1) - sinphi) / (T(1) + sinphi), T(0.5) * eccen));
11fdf7f2 TL	92	}
11fdf7f2 TL	93
92f5a8d4 TL	94	template <typename T, typename Parameters>
92f5a8d4 TL	95	struct base_stere_ellipsoid
11fdf7f2	96	{
92f5a8d4	97	par_stere<T> m_proj_parm;
11fdf7f2 TL	98
	99	// FORWARD(e_forward) ellipsoid
	100	// Project coordinates from geographic (lon, lat) to cartesian (x, y)
92f5a8d4	101	inline void fwd(Parameters const& par, T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const
11fdf7f2	102	{
92f5a8d4	103	static const T half_pi = detail::half_pi<T>();
11fdf7f2	104
92f5a8d4	105	T coslam, sinlam, sinX=0.0, cosX=0.0, X, A = 0.0, sinphi;
11fdf7f2 TL	106
	107	coslam = cos(lp_lon);
	108	sinlam = sin(lp_lon);
	109	sinphi = sin(lp_lat);
92f5a8d4 TL	110	if (this->m_proj_parm.mode == obliq \|\| this->m_proj_parm.mode == equit) {
92f5a8d4 TL	111	sinX = sin(X = 2. * atan(ssfn_(lp_lat, sinphi, par.e)) - half_pi);
11fdf7f2 TL	112	cosX = cos(X);
	113	}
	114	switch (this->m_proj_parm.mode) {
92f5a8d4	115	case obliq:
11fdf7f2 TL	116	A = this->m_proj_parm.akm1 / (this->m_proj_parm.cosX1 * (1. + this->m_proj_parm.sinX1 * sinX +
	117	this->m_proj_parm.cosX1 * cosX * coslam));
	118	xy_y = A * (this->m_proj_parm.cosX1 * sinX - this->m_proj_parm.sinX1 * cosX * coslam);
92f5a8d4 TL	119	goto xmul; /* but why not just xy.x = A * cosX; break; ? */
	120
	121	case equit:
	122	// TODO: calculate denominator once
	123	/* avoid zero division */
	124	if (1. + cosX * coslam == 0.0) {
	125	xy_y = HUGE_VAL;
	126	} else {
	127	A = this->m_proj_parm.akm1 / (1. + cosX * coslam);
	128	xy_y = A * sinX;
	129	}
11fdf7f2 TL	130	xmul:
	131	xy_x = A * cosX;
	132	break;
92f5a8d4 TL	133
92f5a8d4 TL	134	case s_pole:
11fdf7f2 TL	135	lp_lat = -lp_lat;
	136	coslam = - coslam;
	137	sinphi = -sinphi;
	138	BOOST_FALLTHROUGH;
92f5a8d4 TL	139	case n_pole:
92f5a8d4 TL	140	xy_x = this->m_proj_parm.akm1 * pj_tsfn(lp_lat, sinphi, par.e);
11fdf7f2 TL	141	xy_y = - xy_x * coslam;
	142	break;
	143	}
92f5a8d4	144
11fdf7f2 TL	145	xy_x = xy_x * sinlam;
	146	}
	147
	148	// INVERSE(e_inverse) ellipsoid
	149	// Project coordinates from cartesian (x, y) to geographic (lon, lat)
92f5a8d4	150	inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
11fdf7f2	151	{
92f5a8d4	152	static const T half_pi = detail::half_pi<T>();
11fdf7f2	153
92f5a8d4	154	T cosphi, sinphi, tp=0.0, phi_l=0.0, rho, halfe=0.0, halfpi=0.0;
11fdf7f2 TL	155	int i;
	156
	157	rho = boost::math::hypot(xy_x, xy_y);
	158	switch (this->m_proj_parm.mode) {
92f5a8d4 TL	159	case obliq:
92f5a8d4 TL	160	case equit:
11fdf7f2 TL	161	cosphi = cos( tp = 2. * atan2(rho * this->m_proj_parm.cosX1 , this->m_proj_parm.akm1) );
11fdf7f2 TL	162	sinphi = sin(tp);
92f5a8d4	163	if( rho == 0.0 )
11fdf7f2	164	phi_l = asin(cosphi * this->m_proj_parm.sinX1);
92f5a8d4	165	else
11fdf7f2 TL	166	phi_l = asin(cosphi * this->m_proj_parm.sinX1 + (xy_y * sinphi * this->m_proj_parm.cosX1 / rho));
11fdf7f2 TL	167
92f5a8d4	168	tp = tan(.5 * (half_pi + phi_l));
11fdf7f2 TL	169	xy_x *= sinphi;
11fdf7f2 TL	170	xy_y = rho * this->m_proj_parm.cosX1 * cosphi - xy_y * this->m_proj_parm.sinX1* sinphi;
92f5a8d4 TL	171	halfpi = half_pi;
92f5a8d4 TL	172	halfe = .5 * par.e;
11fdf7f2	173	break;
92f5a8d4	174	case n_pole:
11fdf7f2 TL	175	xy_y = -xy_y;
11fdf7f2 TL	176	BOOST_FALLTHROUGH;
92f5a8d4 TL	177	case s_pole:
	178	phi_l = half_pi - 2. * atan(tp = - rho / this->m_proj_parm.akm1);
	179	halfpi = -half_pi;
	180	halfe = -.5 * par.e;
11fdf7f2 TL	181	break;
11fdf7f2 TL	182	}
92f5a8d4 TL	183	for (i = n_iter; i--; phi_l = lp_lat) {
	184	sinphi = par.e * sin(phi_l);
	185	lp_lat = T(2) * atan(tp * math::pow((T(1)+sinphi)/(T(1)-sinphi), halfe)) - halfpi;
	186	if (fabs(phi_l - lp_lat) < conv_tolerance) {
	187	if (this->m_proj_parm.mode == s_pole)
11fdf7f2 TL	188	lp_lat = -lp_lat;
	189	lp_lon = (xy_x == 0. && xy_y == 0.) ? 0. : atan2(xy_x, xy_y);
	190	return;
	191	}
	192	}
92f5a8d4	193	BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
11fdf7f2 TL	194	}
	195
	196	static inline std::string get_name()
	197	{
	198	return "stere_ellipsoid";
	199	}
	200
	201	};
	202
92f5a8d4 TL	203	template <typename T, typename Parameters>
92f5a8d4 TL	204	struct base_stere_spheroid
11fdf7f2	205	{
92f5a8d4	206	par_stere<T> m_proj_parm;
11fdf7f2 TL	207
	208	// FORWARD(s_forward) spheroid
	209	// Project coordinates from geographic (lon, lat) to cartesian (x, y)
92f5a8d4	210	inline void fwd(Parameters const& , T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const
11fdf7f2	211	{
92f5a8d4 TL	212	static const T fourth_pi = detail::fourth_pi<T>();
92f5a8d4 TL	213	static const T half_pi = detail::half_pi<T>();
11fdf7f2	214
92f5a8d4	215	T sinphi, cosphi, coslam, sinlam;
11fdf7f2 TL	216
	217	sinphi = sin(lp_lat);
	218	cosphi = cos(lp_lat);
	219	coslam = cos(lp_lon);
	220	sinlam = sin(lp_lon);
	221	switch (this->m_proj_parm.mode) {
92f5a8d4	222	case equit:
11fdf7f2 TL	223	xy_y = 1. + cosphi * coslam;
11fdf7f2 TL	224	goto oblcon;
92f5a8d4	225	case obliq:
11fdf7f2 TL	226	xy_y = 1. + this->m_proj_parm.sinX1 * sinphi + this->m_proj_parm.cosX1 * cosphi * coslam;
11fdf7f2 TL	227	oblcon:
92f5a8d4 TL	228	if (xy_y <= epsilon10) {
	229	BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
	230	}
11fdf7f2	231	xy_x = (xy_y = this->m_proj_parm.akm1 / xy_y) * cosphi * sinlam;
92f5a8d4	232	xy_y *= (this->m_proj_parm.mode == equit) ? sinphi :
11fdf7f2 TL	233	this->m_proj_parm.cosX1 * sinphi - this->m_proj_parm.sinX1 * cosphi * coslam;
11fdf7f2 TL	234	break;
92f5a8d4	235	case n_pole:
11fdf7f2 TL	236	coslam = - coslam;
	237	lp_lat = - lp_lat;
	238	BOOST_FALLTHROUGH;
92f5a8d4 TL	239	case s_pole:
	240	if (fabs(lp_lat - half_pi) < tolerance) {
	241	BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
	242	}
	243	xy_x = sinlam * ( xy_y = this->m_proj_parm.akm1 * tan(fourth_pi + .5 * lp_lat) );
11fdf7f2 TL	244	xy_y *= coslam;
	245	break;
	246	}
	247	}
	248
	249	// INVERSE(s_inverse) spheroid
	250	// Project coordinates from cartesian (x, y) to geographic (lon, lat)
92f5a8d4	251	inline void inv(Parameters const& par, T const& xy_x, T xy_y, T& lp_lon, T& lp_lat) const
11fdf7f2	252	{
92f5a8d4	253	T c, rh, sinc, cosc;
11fdf7f2 TL	254
	255	sinc = sin(c = 2. * atan((rh = boost::math::hypot(xy_x, xy_y)) / this->m_proj_parm.akm1));
	256	cosc = cos(c);
	257	lp_lon = 0.;
92f5a8d4	258
11fdf7f2	259	switch (this->m_proj_parm.mode) {
92f5a8d4 TL	260	case equit:
92f5a8d4 TL	261	if (fabs(rh) <= epsilon10)
11fdf7f2 TL	262	lp_lat = 0.;
	263	else
	264	lp_lat = asin(xy_y * sinc / rh);
	265	if (cosc != 0. \|\| xy_x != 0.)
	266	lp_lon = atan2(xy_x * sinc, cosc * rh);
	267	break;
92f5a8d4 TL	268	case obliq:
	269	if (fabs(rh) <= epsilon10)
	270	lp_lat = par.phi0;
11fdf7f2 TL	271	else
	272	lp_lat = asin(cosc * this->m_proj_parm.sinX1 + xy_y * sinc * this->m_proj_parm.cosX1 / rh);
	273	if ((c = cosc - this->m_proj_parm.sinX1 * sin(lp_lat)) != 0. \|\| xy_x != 0.)
	274	lp_lon = atan2(xy_x * sinc * this->m_proj_parm.cosX1, c * rh);
	275	break;
92f5a8d4	276	case n_pole:
11fdf7f2 TL	277	xy_y = -xy_y;
11fdf7f2 TL	278	BOOST_FALLTHROUGH;
92f5a8d4 TL	279	case s_pole:
	280	if (fabs(rh) <= epsilon10)
	281	lp_lat = par.phi0;
11fdf7f2	282	else
92f5a8d4	283	lp_lat = asin(this->m_proj_parm.mode == s_pole ? - cosc : cosc);
11fdf7f2 TL	284	lp_lon = (xy_x == 0. && xy_y == 0.) ? 0. : atan2(xy_x, xy_y);
	285	break;
	286	}
	287	}
	288
	289	static inline std::string get_name()
	290	{
	291	return "stere_spheroid";
	292	}
	293
	294	};
	295
	296	template <typename Parameters, typename T>
92f5a8d4	297	inline void setup(Parameters const& par, par_stere<T>& proj_parm) /* general initialization */
11fdf7f2	298	{
92f5a8d4 TL	299	static const T fourth_pi = detail::fourth_pi<T>();
92f5a8d4 TL	300	static const T half_pi = detail::half_pi<T>();
11fdf7f2 TL	301
	302	T t;
	303
92f5a8d4 TL	304	if (fabs((t = fabs(par.phi0)) - half_pi) < epsilon10)
92f5a8d4 TL	305	proj_parm.mode = par.phi0 < 0. ? s_pole : n_pole;
11fdf7f2	306	else
92f5a8d4	307	proj_parm.mode = t > epsilon10 ? obliq : equit;
11fdf7f2	308	proj_parm.phits = fabs(proj_parm.phits);
92f5a8d4	309
11fdf7f2 TL	310	if (par.es != 0.0) {
	311	T X;
	312
	313	switch (proj_parm.mode) {
92f5a8d4 TL	314	case n_pole:
	315	case s_pole:
	316	if (fabs(proj_parm.phits - half_pi) < epsilon10)
11fdf7f2	317	proj_parm.akm1 = 2. * par.k0 /
92f5a8d4	318	sqrt(math::pow(T(1)+par.e,T(1)+par.e)*math::pow(T(1)-par.e,T(1)-par.e));
11fdf7f2 TL	319	else {
	320	proj_parm.akm1 = cos(proj_parm.phits) /
	321	pj_tsfn(proj_parm.phits, t = sin(proj_parm.phits), par.e);
	322	t *= par.e;
	323	proj_parm.akm1 /= sqrt(1. - t * t);
	324	}
	325	break;
92f5a8d4 TL	326	case equit:
92f5a8d4 TL	327	case obliq:
11fdf7f2	328	t = sin(par.phi0);
92f5a8d4	329	X = 2. * atan(ssfn_(par.phi0, t, par.e)) - half_pi;
11fdf7f2 TL	330	t *= par.e;
	331	proj_parm.akm1 = 2. * par.k0 * cos(par.phi0) / sqrt(1. - t * t);
	332	proj_parm.sinX1 = sin(X);
	333	proj_parm.cosX1 = cos(X);
	334	break;
	335	}
	336	} else {
	337	switch (proj_parm.mode) {
92f5a8d4	338	case obliq:
11fdf7f2 TL	339	proj_parm.sinX1 = sin(par.phi0);
	340	proj_parm.cosX1 = cos(par.phi0);
	341	BOOST_FALLTHROUGH;
92f5a8d4	342	case equit:
11fdf7f2 TL	343	proj_parm.akm1 = 2. * par.k0;
11fdf7f2 TL	344	break;
92f5a8d4 TL	345	case s_pole:
	346	case n_pole:
	347	proj_parm.akm1 = fabs(proj_parm.phits - half_pi) >= epsilon10 ?
	348	cos(proj_parm.phits) / tan(fourth_pi - .5 * proj_parm.phits) :
11fdf7f2 TL	349	2. * par.k0 ;
	350	break;
	351	}
	352	}
	353	}
	354
	355
	356	// Stereographic
92f5a8d4 TL	357	template <typename Params, typename Parameters, typename T>
92f5a8d4 TL	358	inline void setup_stere(Params const& params, Parameters const& par, par_stere<T>& proj_parm)
11fdf7f2	359	{
92f5a8d4 TL	360	static const T half_pi = detail::half_pi<T>();
	361
	362	if (! pj_param_r<srs::spar::lat_ts>(params, "lat_ts", srs::dpar::lat_ts, proj_parm.phits))
	363	proj_parm.phits = half_pi;
11fdf7f2	364
11fdf7f2 TL	365	setup(par, proj_parm);
	366	}
	367
	368	// Universal Polar Stereographic
92f5a8d4 TL	369	template <typename Params, typename Parameters, typename T>
92f5a8d4 TL	370	inline void setup_ups(Params const& params, Parameters& par, par_stere<T>& proj_parm)
11fdf7f2	371	{
92f5a8d4	372	static const T half_pi = detail::half_pi<T>();
11fdf7f2 TL	373
11fdf7f2 TL	374	/* International Ellipsoid */
92f5a8d4 TL	375	par.phi0 = pj_get_param_b<srs::spar::south>(params, "south", srs::dpar::south) ? -half_pi: half_pi;
	376	if (par.es == 0.0) {
	377	BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) );
	378	}
11fdf7f2 TL	379	par.k0 = .994;
	380	par.x0 = 2000000.;
	381	par.y0 = 2000000.;
92f5a8d4	382	proj_parm.phits = half_pi;
11fdf7f2	383	par.lam0 = 0.;
92f5a8d4	384
11fdf7f2 TL	385	setup(par, proj_parm);
	386	}
	387
	388	}} // namespace detail::stere
	389	#endif // doxygen
	390
	391	/*!
	392	\brief Stereographic projection
	393	\ingroup projections
	394	\tparam Geographic latlong point type
	395	\tparam Cartesian xy point type
	396	\tparam Parameters parameter type
	397	\par Projection characteristics
	398	- Azimuthal
	399	- Spheroid
	400	- Ellipsoid
	401	\par Projection parameters
	402	- lat_ts: Latitude of true scale (degrees)
	403	\par Example
	404	\image html ex_stere.gif
	405	*/
92f5a8d4 TL	406	template <typename T, typename Parameters>
92f5a8d4 TL	407	struct stere_ellipsoid : public detail::stere::base_stere_ellipsoid<T, Parameters>
11fdf7f2	408	{
92f5a8d4 TL	409	template <typename Params>
92f5a8d4 TL	410	inline stere_ellipsoid(Params const& params, Parameters const& par)
11fdf7f2	411	{
92f5a8d4	412	detail::stere::setup_stere(params, par, this->m_proj_parm);
11fdf7f2 TL	413	}
	414	};
	415
	416	/*!
	417	\brief Stereographic projection
	418	\ingroup projections
	419	\tparam Geographic latlong point type
	420	\tparam Cartesian xy point type
	421	\tparam Parameters parameter type
	422	\par Projection characteristics
	423	- Azimuthal
	424	- Spheroid
	425	- Ellipsoid
	426	\par Projection parameters
	427	- lat_ts: Latitude of true scale (degrees)
	428	\par Example
	429	\image html ex_stere.gif
	430	*/
92f5a8d4 TL	431	template <typename T, typename Parameters>
92f5a8d4 TL	432	struct stere_spheroid : public detail::stere::base_stere_spheroid<T, Parameters>
11fdf7f2	433	{
92f5a8d4 TL	434	template <typename Params>
92f5a8d4 TL	435	inline stere_spheroid(Params const& params, Parameters const& par)
11fdf7f2	436	{
92f5a8d4	437	detail::stere::setup_stere(params, par, this->m_proj_parm);
11fdf7f2 TL	438	}
	439	};
	440
	441	/*!
	442	\brief Universal Polar Stereographic projection
	443	\ingroup projections
	444	\tparam Geographic latlong point type
	445	\tparam Cartesian xy point type
	446	\tparam Parameters parameter type
	447	\par Projection characteristics
	448	- Azimuthal
	449	- Spheroid
	450	- Ellipsoid
	451	\par Projection parameters
	452	- south: Denotes southern hemisphere UTM zone (boolean)
	453	\par Example
	454	\image html ex_ups.gif
	455	*/
92f5a8d4 TL	456	template <typename T, typename Parameters>
92f5a8d4 TL	457	struct ups_ellipsoid : public detail::stere::base_stere_ellipsoid<T, Parameters>
11fdf7f2	458	{
92f5a8d4 TL	459	template <typename Params>
92f5a8d4 TL	460	inline ups_ellipsoid(Params const& params, Parameters & par)
11fdf7f2	461	{
92f5a8d4	462	detail::stere::setup_ups(params, par, this->m_proj_parm);
11fdf7f2 TL	463	}
	464	};
	465
	466	/*!
	467	\brief Universal Polar Stereographic projection
	468	\ingroup projections
	469	\tparam Geographic latlong point type
	470	\tparam Cartesian xy point type
	471	\tparam Parameters parameter type
	472	\par Projection characteristics
	473	- Azimuthal
	474	- Spheroid
	475	- Ellipsoid
	476	\par Projection parameters
	477	- south: Denotes southern hemisphere UTM zone (boolean)
	478	\par Example
	479	\image html ex_ups.gif
	480	*/
92f5a8d4 TL	481	template <typename T, typename Parameters>
92f5a8d4 TL	482	struct ups_spheroid : public detail::stere::base_stere_spheroid<T, Parameters>
11fdf7f2	483	{
92f5a8d4 TL	484	template <typename Params>
92f5a8d4 TL	485	inline ups_spheroid(Params const& params, Parameters & par)
11fdf7f2	486	{
92f5a8d4	487	detail::stere::setup_ups(params, par, this->m_proj_parm);
11fdf7f2 TL	488	}
	489	};
	490
	491	#ifndef DOXYGEN_NO_DETAIL
	492	namespace detail
	493	{
	494
	495	// Static projection
92f5a8d4 TL	496	BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_stere, stere_spheroid, stere_ellipsoid)
92f5a8d4 TL	497	BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_ups, ups_spheroid, ups_ellipsoid)
11fdf7f2 TL	498
11fdf7f2 TL	499	// Factory entry(s)
92f5a8d4 TL	500	BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(stere_entry, stere_spheroid, stere_ellipsoid)
	501	BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(ups_entry, ups_spheroid, ups_ellipsoid)
	502
	503	BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(stere_init)
11fdf7f2	504	{
92f5a8d4 TL	505	BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(stere, stere_entry)
92f5a8d4 TL	506	BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(ups, ups_entry)
11fdf7f2 TL	507	}
	508
	509	} // namespace detail
	510	#endif // doxygen
	511
	512	} // namespace projections
	513
	514	}} // namespace boost::geometry
	515
	516	#endif // BOOST_GEOMETRY_PROJECTIONS_STERE_HPP
	517