[ceph.git] / ceph / src / boost / boost / geometry / srs / projections / proj / laea.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_LAEA_HPP
#define BOOST_GEOMETRY_PROJECTIONS_LAEA_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/projects.hpp>
#include <boost/geometry/srs/projections/impl/factory_entry.hpp>
#include <boost/geometry/srs/projections/impl/pj_auth.hpp>
#include <boost/geometry/srs/projections/impl/pj_qsfn.hpp>

namespace boost { namespace geometry
{

namespace projections
{
    #ifndef DOXYGEN_NO_DETAIL
    namespace detail { namespace laea
    {
            static const double epsilon10 = 1.e-10;

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

            template <typename T>
            struct par_laea
            {
                T   sinb1;
                T   cosb1;
                T   xmf;
                T   ymf;
                T   mmf;
                T   qp;
                T   dd;
                T   rq;
                detail::apa<T> apa;
                mode_type mode;
            };

            template <typename T, typename Parameters>
            struct base_laea_ellipsoid
            {
                par_laea<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 const& lp_lat, T& xy_x, T& xy_y) const
                {
                    static const T half_pi = detail::half_pi<T>();

                    T coslam, sinlam, sinphi, q, sinb=0.0, cosb=0.0, b=0.0;

                    coslam = cos(lp_lon);
                    sinlam = sin(lp_lon);
                    sinphi = sin(lp_lat);
                    q = pj_qsfn(sinphi, par.e, par.one_es);

                    if (this->m_proj_parm.mode == obliq || this->m_proj_parm.mode == equit) {
                        sinb = q / this->m_proj_parm.qp;
                        cosb = sqrt(1. - sinb * sinb);
                    }

                    switch (this->m_proj_parm.mode) {
                    case obliq:
                        b = 1. + this->m_proj_parm.sinb1 * sinb + this->m_proj_parm.cosb1 * cosb * coslam;
                        break;
                    case equit:
                        b = 1. + cosb * coslam;
                        break;
                    case n_pole:
                        b = half_pi + lp_lat;
                        q = this->m_proj_parm.qp - q;
                        break;
                    case s_pole:
                        b = lp_lat - half_pi;
                        q = this->m_proj_parm.qp + q;
                        break;
                    }
                    if (fabs(b) < epsilon10) {
                        BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
                    }

                    switch (this->m_proj_parm.mode) {
                    case obliq:
                        b = sqrt(2. / b);
                        xy_y = this->m_proj_parm.ymf * b * (this->m_proj_parm.cosb1 * sinb - this->m_proj_parm.sinb1 * cosb * coslam);
                        goto eqcon;
                        break;
                    case equit:
                        b = sqrt(2. / (1. + cosb * coslam));
                        xy_y = b * sinb * this->m_proj_parm.ymf;
                eqcon:
                        xy_x = this->m_proj_parm.xmf * b * cosb * sinlam;
                        break;
                    case n_pole:
                    case s_pole:
                        if (q >= 0.) {
                            b = sqrt(q);
                            xy_x = b * sinlam;
                            xy_y = coslam * (this->m_proj_parm.mode == s_pole ? b : -b);
                        } else
                            xy_x = xy_y = 0.;
                        break;
                    }
                }

                // 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
                {
                    T cCe, sCe, q, rho, ab=0.0;

                    switch (this->m_proj_parm.mode) {
                    case equit:
                    case obliq:
                        xy_x /= this->m_proj_parm.dd;
                        xy_y *=  this->m_proj_parm.dd;
                        rho = boost::math::hypot(xy_x, xy_y);
                        if (rho < epsilon10) {
                            lp_lon = 0.;
                            lp_lat = par.phi0;
                            return;
                        }
                        sCe = 2. * asin(.5 * rho / this->m_proj_parm.rq);
                        cCe = cos(sCe);
                        sCe = sin(sCe);
                        xy_x *= sCe;
                        if (this->m_proj_parm.mode == obliq) {
                            ab = cCe * this->m_proj_parm.sinb1 + xy_y * sCe * this->m_proj_parm.cosb1 / rho;
                            xy_y = rho * this->m_proj_parm.cosb1 * cCe - xy_y * this->m_proj_parm.sinb1 * sCe;
                        } else {
                            ab = xy_y * sCe / rho;
                            xy_y = rho * cCe;
                        }
                        break;
                    case n_pole:
                        xy_y = -xy_y;
                        BOOST_FALLTHROUGH;
                    case s_pole:
                        q = (xy_x * xy_x + xy_y * xy_y);
                        if (q == 0.0) {
                            lp_lon = 0.;
                            lp_lat = par.phi0;
                            return;
                        }
                        ab = 1. - q / this->m_proj_parm.qp;
                        if (this->m_proj_parm.mode == s_pole)
                            ab = - ab;
                        break;
                    }
                    lp_lon = atan2(xy_x, xy_y);
                    lp_lat = pj_authlat(asin(ab), this->m_proj_parm.apa);
                }

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

            };

            template <typename T, typename Parameters>
            struct base_laea_spheroid
            {
                par_laea<T> m_proj_parm;

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

                    T  coslam, cosphi, sinphi;

                    sinphi = sin(lp_lat);
                    cosphi = cos(lp_lat);
                    coslam = cos(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.sinb1 * sinphi + this->m_proj_parm.cosb1 * cosphi * coslam;
                oblcon:
                        if (xy_y <= epsilon10) {
                            BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
                        }
                        xy_y = sqrt(2. / xy_y);
                        xy_x = xy_y * cosphi * sin(lp_lon);
                        xy_y *= this->m_proj_parm.mode == equit ? sinphi :
                           this->m_proj_parm.cosb1 * sinphi - this->m_proj_parm.sinb1 * cosphi * coslam;
                        break;
                    case n_pole:
                        coslam = -coslam;
                        BOOST_FALLTHROUGH;
                    case s_pole:
                        if (fabs(lp_lat + par.phi0) < epsilon10) {
                            BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
                        }
                        xy_y = fourth_pi - lp_lat * .5;
                        xy_y = 2. * (this->m_proj_parm.mode == s_pole ? cos(xy_y) : sin(xy_y));
                        xy_x = xy_y * sin(lp_lon);
                        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 xy_x, T xy_y, T& lp_lon, T& lp_lat) const
                {
                    static const T half_pi = detail::half_pi<T>();

                    T  cosz=0.0, rh, sinz=0.0;

                    rh = boost::math::hypot(xy_x, xy_y);
                    if ((lp_lat = rh * .5 ) > 1.) {
                        BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
                    }
                    lp_lat = 2. * asin(lp_lat);
                    if (this->m_proj_parm.mode == obliq || this->m_proj_parm.mode == equit) {
                        sinz = sin(lp_lat);
                        cosz = cos(lp_lat);
                    }
                    switch (this->m_proj_parm.mode) {
                    case equit:
                        lp_lat = fabs(rh) <= epsilon10 ? 0. : asin(xy_y * sinz / rh);
                        xy_x *= sinz;
                        xy_y = cosz * rh;
                        break;
                    case obliq:
                        lp_lat = fabs(rh) <= epsilon10 ? par.phi0 :
                           asin(cosz * this->m_proj_parm.sinb1 + xy_y * sinz * this->m_proj_parm.cosb1 / rh);
                        xy_x *= sinz * this->m_proj_parm.cosb1;
                        xy_y = (cosz - sin(lp_lat) * this->m_proj_parm.sinb1) * rh;
                        break;
                    case n_pole:
                        xy_y = -xy_y;
                        lp_lat = half_pi - lp_lat;
                        break;
                    case s_pole:
                        lp_lat -= half_pi;
                        break;
                    }
                    lp_lon = (xy_y == 0. && (this->m_proj_parm.mode == equit || this->m_proj_parm.mode == obliq)) ?
                        0. : atan2(xy_x, xy_y);
                }

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

            };

            // Lambert Azimuthal Equal Area
            template <typename Parameters, typename T>
            inline void setup_laea(Parameters& par, par_laea<T>& proj_parm)
            {
                static const T half_pi = detail::half_pi<T>();

                T t;

                t = fabs(par.phi0);
                if (fabs(t - half_pi) < epsilon10)
                    proj_parm.mode = par.phi0 < 0. ? s_pole : n_pole;
                else if (fabs(t) < epsilon10)
                    proj_parm.mode = equit;
                else
                    proj_parm.mode = obliq;
                if (par.es != 0.0) {
                    double sinphi;

                    par.e = sqrt(par.es); // TODO : Isn't it already set?
                    proj_parm.qp = pj_qsfn(1., par.e, par.one_es);
                    proj_parm.mmf = .5 / (1. - par.es);
                    proj_parm.apa = pj_authset<T>(par.es);
                    switch (proj_parm.mode) {
                    case n_pole:
                    case s_pole:
                        proj_parm.dd = 1.;
                        break;
                    case equit:
                        proj_parm.dd = 1. / (proj_parm.rq = sqrt(.5 * proj_parm.qp));
                        proj_parm.xmf = 1.;
                        proj_parm.ymf = .5 * proj_parm.qp;
                        break;
                    case obliq:
                        proj_parm.rq = sqrt(.5 * proj_parm.qp);
                        sinphi = sin(par.phi0);
                        proj_parm.sinb1 = pj_qsfn(sinphi, par.e, par.one_es) / proj_parm.qp;
                        proj_parm.cosb1 = sqrt(1. - proj_parm.sinb1 * proj_parm.sinb1);
                        proj_parm.dd = cos(par.phi0) / (sqrt(1. - par.es * sinphi * sinphi) *
                           proj_parm.rq * proj_parm.cosb1);
                        proj_parm.ymf = (proj_parm.xmf = proj_parm.rq) / proj_parm.dd;
                        proj_parm.xmf *= proj_parm.dd;
                        break;
                    }
                } else {
                    if (proj_parm.mode == obliq) {
                        proj_parm.sinb1 = sin(par.phi0);
                        proj_parm.cosb1 = cos(par.phi0);
                    }
                }
            }

    }} // namespace laea
    #endif // doxygen

    /*!
        \brief Lambert Azimuthal Equal Area projection
        \ingroup projections
        \tparam Geographic latlong point type
        \tparam Cartesian xy point type
        \tparam Parameters parameter type
        \par Projection characteristics
         - Azimuthal
         - Spheroid
         - Ellipsoid
        \par Example
        \image html ex_laea.gif
    */
    template <typename T, typename Parameters>
    struct laea_ellipsoid : public detail::laea::base_laea_ellipsoid<T, Parameters>
    {
        template <typename Params>
        inline laea_ellipsoid(Params const& , Parameters & par)
        {
            detail::laea::setup_laea(par, this->m_proj_parm);
        }
    };

    /*!
        \brief Lambert Azimuthal Equal Area projection
        \ingroup projections
        \tparam Geographic latlong point type
        \tparam Cartesian xy point type
        \tparam Parameters parameter type
        \par Projection characteristics
         - Azimuthal
         - Spheroid
         - Ellipsoid
        \par Example
        \image html ex_laea.gif
    */
    template <typename T, typename Parameters>
    struct laea_spheroid : public detail::laea::base_laea_spheroid<T, Parameters>
    {
        template <typename Params>
        inline laea_spheroid(Params const& , Parameters & par)
        {
            detail::laea::setup_laea(par, this->m_proj_parm);
        }
    };

    #ifndef DOXYGEN_NO_DETAIL
    namespace detail
    {

        // Static projection
        BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_laea, laea_spheroid, laea_ellipsoid)

        // Factory entry(s)
        BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(laea_entry, laea_spheroid, laea_ellipsoid)
        
        BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(laea_init)
        {
            BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(laea, laea_entry)
        }

    } // namespace detail
    #endif // doxygen

} // namespace projections

}} // namespace boost::geometry

#endif // BOOST_GEOMETRY_PROJECTIONS_LAEA_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_LAEA_HPP
	41	#define BOOST_GEOMETRY_PROJECTIONS_LAEA_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>
	49	#include <boost/geometry/srs/projections/impl/projects.hpp>
	50	#include <boost/geometry/srs/projections/impl/factory_entry.hpp>
	51	#include <boost/geometry/srs/projections/impl/pj_auth.hpp>
	52	#include <boost/geometry/srs/projections/impl/pj_qsfn.hpp>
	53
	54	namespace boost { namespace geometry
	55	{
	56
11fdf7f2 TL	57	namespace projections
	58	{
	59	#ifndef DOXYGEN_NO_DETAIL
	60	namespace detail { namespace laea
	61	{
92f5a8d4 TL	62	static const double epsilon10 = 1.e-10;
	63
	64	enum mode_type {
	65	n_pole = 0,
	66	s_pole = 1,
	67	equit = 2,
	68	obliq = 3
	69	};
11fdf7f2 TL	70
	71	template <typename T>
	72	struct par_laea
	73	{
	74	T sinb1;
	75	T cosb1;
	76	T xmf;
	77	T ymf;
	78	T mmf;
	79	T qp;
	80	T dd;
	81	T rq;
92f5a8d4 TL	82	detail::apa<T> apa;
92f5a8d4 TL	83	mode_type mode;
11fdf7f2 TL	84	};
11fdf7f2 TL	85
92f5a8d4 TL	86	template <typename T, typename Parameters>
92f5a8d4 TL	87	struct base_laea_ellipsoid
11fdf7f2	88	{
92f5a8d4	89	par_laea<T> m_proj_parm;
11fdf7f2 TL	90
	91	// FORWARD(e_forward) ellipsoid
	92	// Project coordinates from geographic (lon, lat) to cartesian (x, y)
92f5a8d4	93	inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
11fdf7f2	94	{
92f5a8d4	95	static const T half_pi = detail::half_pi<T>();
11fdf7f2	96
92f5a8d4	97	T coslam, sinlam, sinphi, q, sinb=0.0, cosb=0.0, b=0.0;
11fdf7f2 TL	98
	99	coslam = cos(lp_lon);
	100	sinlam = sin(lp_lon);
	101	sinphi = sin(lp_lat);
92f5a8d4 TL	102	q = pj_qsfn(sinphi, par.e, par.one_es);
	103
	104	if (this->m_proj_parm.mode == obliq \|\| this->m_proj_parm.mode == equit) {
11fdf7f2 TL	105	sinb = q / this->m_proj_parm.qp;
	106	cosb = sqrt(1. - sinb * sinb);
	107	}
92f5a8d4	108
11fdf7f2	109	switch (this->m_proj_parm.mode) {
92f5a8d4	110	case obliq:
11fdf7f2 TL	111	b = 1. + this->m_proj_parm.sinb1 * sinb + this->m_proj_parm.cosb1 * cosb * coslam;
11fdf7f2 TL	112	break;
92f5a8d4	113	case equit:
11fdf7f2 TL	114	b = 1. + cosb * coslam;
11fdf7f2 TL	115	break;
92f5a8d4 TL	116	case n_pole:
92f5a8d4 TL	117	b = half_pi + lp_lat;
11fdf7f2 TL	118	q = this->m_proj_parm.qp - q;
11fdf7f2 TL	119	break;
92f5a8d4 TL	120	case s_pole:
92f5a8d4 TL	121	b = lp_lat - half_pi;
11fdf7f2 TL	122	q = this->m_proj_parm.qp + q;
	123	break;
	124	}
92f5a8d4 TL	125	if (fabs(b) < epsilon10) {
	126	BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
	127	}
	128
11fdf7f2	129	switch (this->m_proj_parm.mode) {
92f5a8d4 TL	130	case obliq:
	131	b = sqrt(2. / b);
	132	xy_y = this->m_proj_parm.ymf * b * (this->m_proj_parm.cosb1 * sinb - this->m_proj_parm.sinb1 * cosb * coslam);
11fdf7f2 TL	133	goto eqcon;
11fdf7f2 TL	134	break;
92f5a8d4 TL	135	case equit:
	136	b = sqrt(2. / (1. + cosb * coslam));
	137	xy_y = b * sinb * this->m_proj_parm.ymf;
11fdf7f2 TL	138	eqcon:
	139	xy_x = this->m_proj_parm.xmf * b * cosb * sinlam;
	140	break;
92f5a8d4 TL	141	case n_pole:
92f5a8d4 TL	142	case s_pole:
11fdf7f2	143	if (q >= 0.) {
92f5a8d4 TL	144	b = sqrt(q);
	145	xy_x = b * sinlam;
	146	xy_y = coslam * (this->m_proj_parm.mode == s_pole ? b : -b);
11fdf7f2 TL	147	} else
	148	xy_x = xy_y = 0.;
	149	break;
	150	}
	151	}
	152
	153	// INVERSE(e_inverse) ellipsoid
	154	// Project coordinates from cartesian (x, y) to geographic (lon, lat)
92f5a8d4	155	inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
11fdf7f2	156	{
92f5a8d4	157	T cCe, sCe, q, rho, ab=0.0;
11fdf7f2 TL	158
11fdf7f2 TL	159	switch (this->m_proj_parm.mode) {
92f5a8d4 TL	160	case equit:
	161	case obliq:
	162	xy_x /= this->m_proj_parm.dd;
	163	xy_y *= this->m_proj_parm.dd;
	164	rho = boost::math::hypot(xy_x, xy_y);
	165	if (rho < epsilon10) {
11fdf7f2	166	lp_lon = 0.;
92f5a8d4	167	lp_lat = par.phi0;
11fdf7f2 TL	168	return;
11fdf7f2 TL	169	}
92f5a8d4 TL	170	sCe = 2. * asin(.5 * rho / this->m_proj_parm.rq);
	171	cCe = cos(sCe);
	172	sCe = sin(sCe);
	173	xy_x *= sCe;
	174	if (this->m_proj_parm.mode == obliq) {
	175	ab = cCe * this->m_proj_parm.sinb1 + xy_y * sCe * this->m_proj_parm.cosb1 / rho;
11fdf7f2 TL	176	xy_y = rho * this->m_proj_parm.cosb1 * cCe - xy_y * this->m_proj_parm.sinb1 * sCe;
11fdf7f2 TL	177	} else {
92f5a8d4	178	ab = xy_y * sCe / rho;
11fdf7f2 TL	179	xy_y = rho * cCe;
	180	}
	181	break;
92f5a8d4	182	case n_pole:
11fdf7f2 TL	183	xy_y = -xy_y;
11fdf7f2 TL	184	BOOST_FALLTHROUGH;
92f5a8d4 TL	185	case s_pole:
	186	q = (xy_x * xy_x + xy_y * xy_y);
	187	if (q == 0.0) {
11fdf7f2	188	lp_lon = 0.;
92f5a8d4	189	lp_lat = par.phi0;
11fdf7f2 TL	190	return;
11fdf7f2 TL	191	}
11fdf7f2	192	ab = 1. - q / this->m_proj_parm.qp;
92f5a8d4	193	if (this->m_proj_parm.mode == s_pole)
11fdf7f2 TL	194	ab = - ab;
	195	break;
	196	}
	197	lp_lon = atan2(xy_x, xy_y);
	198	lp_lat = pj_authlat(asin(ab), this->m_proj_parm.apa);
	199	}
	200
	201	static inline std::string get_name()
	202	{
	203	return "laea_ellipsoid";
	204	}
	205
	206	};
	207
92f5a8d4 TL	208	template <typename T, typename Parameters>
92f5a8d4 TL	209	struct base_laea_spheroid
11fdf7f2	210	{
92f5a8d4	211	par_laea<T> m_proj_parm;
11fdf7f2 TL	212
	213	// FORWARD(s_forward) spheroid
	214	// Project coordinates from geographic (lon, lat) to cartesian (x, y)
92f5a8d4	215	inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
11fdf7f2	216	{
92f5a8d4	217	static const T fourth_pi = detail::fourth_pi<T>();
11fdf7f2	218
92f5a8d4	219	T coslam, cosphi, sinphi;
11fdf7f2 TL	220
	221	sinphi = sin(lp_lat);
	222	cosphi = cos(lp_lat);
	223	coslam = cos(lp_lon);
	224	switch (this->m_proj_parm.mode) {
92f5a8d4	225	case equit:
11fdf7f2 TL	226	xy_y = 1. + cosphi * coslam;
11fdf7f2 TL	227	goto oblcon;
92f5a8d4	228	case obliq:
11fdf7f2 TL	229	xy_y = 1. + this->m_proj_parm.sinb1 * sinphi + this->m_proj_parm.cosb1 * cosphi * coslam;
11fdf7f2 TL	230	oblcon:
92f5a8d4 TL	231	if (xy_y <= epsilon10) {
	232	BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
	233	}
	234	xy_y = sqrt(2. / xy_y);
	235	xy_x = xy_y * cosphi * sin(lp_lon);
	236	xy_y *= this->m_proj_parm.mode == equit ? sinphi :
11fdf7f2 TL	237	this->m_proj_parm.cosb1 * sinphi - this->m_proj_parm.sinb1 * cosphi * coslam;
11fdf7f2 TL	238	break;
92f5a8d4	239	case n_pole:
11fdf7f2 TL	240	coslam = -coslam;
11fdf7f2 TL	241	BOOST_FALLTHROUGH;
92f5a8d4 TL	242	case s_pole:
	243	if (fabs(lp_lat + par.phi0) < epsilon10) {
	244	BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
	245	}
	246	xy_y = fourth_pi - lp_lat * .5;
	247	xy_y = 2. * (this->m_proj_parm.mode == s_pole ? cos(xy_y) : sin(xy_y));
11fdf7f2 TL	248	xy_x = xy_y * sin(lp_lon);
	249	xy_y *= coslam;
	250	break;
	251	}
	252	}
	253
	254	// INVERSE(s_inverse) spheroid
	255	// Project coordinates from cartesian (x, y) to geographic (lon, lat)
92f5a8d4	256	inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
11fdf7f2	257	{
92f5a8d4	258	static const T half_pi = detail::half_pi<T>();
11fdf7f2	259
92f5a8d4	260	T cosz=0.0, rh, sinz=0.0;
11fdf7f2 TL	261
11fdf7f2 TL	262	rh = boost::math::hypot(xy_x, xy_y);
92f5a8d4 TL	263	if ((lp_lat = rh * .5 ) > 1.) {
	264	BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
	265	}
11fdf7f2	266	lp_lat = 2. * asin(lp_lat);
92f5a8d4	267	if (this->m_proj_parm.mode == obliq \|\| this->m_proj_parm.mode == equit) {
11fdf7f2 TL	268	sinz = sin(lp_lat);
	269	cosz = cos(lp_lat);
	270	}
	271	switch (this->m_proj_parm.mode) {
92f5a8d4 TL	272	case equit:
92f5a8d4 TL	273	lp_lat = fabs(rh) <= epsilon10 ? 0. : asin(xy_y * sinz / rh);
11fdf7f2 TL	274	xy_x *= sinz;
	275	xy_y = cosz * rh;
	276	break;
92f5a8d4 TL	277	case obliq:
92f5a8d4 TL	278	lp_lat = fabs(rh) <= epsilon10 ? par.phi0 :
11fdf7f2 TL	279	asin(cosz * this->m_proj_parm.sinb1 + xy_y * sinz * this->m_proj_parm.cosb1 / rh);
	280	xy_x = sinz this->m_proj_parm.cosb1;
	281	xy_y = (cosz - sin(lp_lat) * this->m_proj_parm.sinb1) * rh;
	282	break;
92f5a8d4	283	case n_pole:
11fdf7f2	284	xy_y = -xy_y;
92f5a8d4	285	lp_lat = half_pi - lp_lat;
11fdf7f2	286	break;
92f5a8d4 TL	287	case s_pole:
92f5a8d4 TL	288	lp_lat -= half_pi;
11fdf7f2 TL	289	break;
11fdf7f2 TL	290	}
92f5a8d4	291	lp_lon = (xy_y == 0. && (this->m_proj_parm.mode == equit \|\| this->m_proj_parm.mode == obliq)) ?
11fdf7f2 TL	292	0. : atan2(xy_x, xy_y);
	293	}
	294
	295	static inline std::string get_name()
	296	{
	297	return "laea_spheroid";
	298	}
	299
	300	};
	301
	302	// Lambert Azimuthal Equal Area
	303	template <typename Parameters, typename T>
	304	inline void setup_laea(Parameters& par, par_laea<T>& proj_parm)
	305	{
92f5a8d4	306	static const T half_pi = detail::half_pi<T>();
11fdf7f2 TL	307
	308	T t;
	309
92f5a8d4 TL	310	t = fabs(par.phi0);
	311	if (fabs(t - half_pi) < epsilon10)
	312	proj_parm.mode = par.phi0 < 0. ? s_pole : n_pole;
	313	else if (fabs(t) < epsilon10)
	314	proj_parm.mode = equit;
11fdf7f2	315	else
92f5a8d4 TL	316	proj_parm.mode = obliq;
92f5a8d4 TL	317	if (par.es != 0.0) {
11fdf7f2 TL	318	double sinphi;
11fdf7f2 TL	319
92f5a8d4	320	par.e = sqrt(par.es); // TODO : Isn't it already set?
11fdf7f2 TL	321	proj_parm.qp = pj_qsfn(1., par.e, par.one_es);
11fdf7f2 TL	322	proj_parm.mmf = .5 / (1. - par.es);
92f5a8d4	323	proj_parm.apa = pj_authset<T>(par.es);
11fdf7f2	324	switch (proj_parm.mode) {
92f5a8d4 TL	325	case n_pole:
92f5a8d4 TL	326	case s_pole:
11fdf7f2 TL	327	proj_parm.dd = 1.;
11fdf7f2 TL	328	break;
92f5a8d4	329	case equit:
11fdf7f2 TL	330	proj_parm.dd = 1. / (proj_parm.rq = sqrt(.5 * proj_parm.qp));
	331	proj_parm.xmf = 1.;
	332	proj_parm.ymf = .5 * proj_parm.qp;
	333	break;
92f5a8d4	334	case obliq:
11fdf7f2 TL	335	proj_parm.rq = sqrt(.5 * proj_parm.qp);
	336	sinphi = sin(par.phi0);
	337	proj_parm.sinb1 = pj_qsfn(sinphi, par.e, par.one_es) / proj_parm.qp;
	338	proj_parm.cosb1 = sqrt(1. - proj_parm.sinb1 * proj_parm.sinb1);
	339	proj_parm.dd = cos(par.phi0) / (sqrt(1. - par.es * sinphi * sinphi) *
	340	proj_parm.rq * proj_parm.cosb1);
	341	proj_parm.ymf = (proj_parm.xmf = proj_parm.rq) / proj_parm.dd;
	342	proj_parm.xmf *= proj_parm.dd;
	343	break;
	344	}
	345	} else {
92f5a8d4	346	if (proj_parm.mode == obliq) {
11fdf7f2 TL	347	proj_parm.sinb1 = sin(par.phi0);
	348	proj_parm.cosb1 = cos(par.phi0);
	349	}
	350	}
	351	}
	352
	353	}} // namespace laea
	354	#endif // doxygen
	355
	356	/*!
	357	\brief Lambert Azimuthal Equal Area projection
	358	\ingroup projections
	359	\tparam Geographic latlong point type
	360	\tparam Cartesian xy point type
	361	\tparam Parameters parameter type
	362	\par Projection characteristics
	363	- Azimuthal
	364	- Spheroid
	365	- Ellipsoid
	366	\par Example
	367	\image html ex_laea.gif
	368	*/
92f5a8d4 TL	369	template <typename T, typename Parameters>
92f5a8d4 TL	370	struct laea_ellipsoid : public detail::laea::base_laea_ellipsoid<T, Parameters>
11fdf7f2	371	{
92f5a8d4 TL	372	template <typename Params>
92f5a8d4 TL	373	inline laea_ellipsoid(Params const& , Parameters & par)
11fdf7f2	374	{
92f5a8d4	375	detail::laea::setup_laea(par, this->m_proj_parm);
11fdf7f2 TL	376	}
	377	};
	378
	379	/*!
	380	\brief Lambert Azimuthal Equal Area projection
	381	\ingroup projections
	382	\tparam Geographic latlong point type
	383	\tparam Cartesian xy point type
	384	\tparam Parameters parameter type
	385	\par Projection characteristics
	386	- Azimuthal
	387	- Spheroid
	388	- Ellipsoid
	389	\par Example
	390	\image html ex_laea.gif
	391	*/
92f5a8d4 TL	392	template <typename T, typename Parameters>
92f5a8d4 TL	393	struct laea_spheroid : public detail::laea::base_laea_spheroid<T, Parameters>
11fdf7f2	394	{
92f5a8d4 TL	395	template <typename Params>
92f5a8d4 TL	396	inline laea_spheroid(Params const& , Parameters & par)
11fdf7f2	397	{
92f5a8d4	398	detail::laea::setup_laea(par, this->m_proj_parm);
11fdf7f2 TL	399	}
	400	};
	401
	402	#ifndef DOXYGEN_NO_DETAIL
	403	namespace detail
	404	{
	405
	406	// Static projection
92f5a8d4	407	BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_laea, laea_spheroid, laea_ellipsoid)
11fdf7f2 TL	408
11fdf7f2 TL	409	// Factory entry(s)
92f5a8d4 TL	410	BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(laea_entry, laea_spheroid, laea_ellipsoid)
	411
	412	BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(laea_init)
11fdf7f2	413	{
92f5a8d4	414	BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(laea, laea_entry)
11fdf7f2 TL	415	}
	416
	417	} // namespace detail
	418	#endif // doxygen
	419
	420	} // namespace projections
	421
	422	}} // namespace boost::geometry
	423
	424	#endif // BOOST_GEOMETRY_PROJECTIONS_LAEA_HPP
	425