[ceph.git] / ceph / src / boost / boost / geometry / srs / projections / proj / aea.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

 // Author: Gerald Evenden (1995)
 //         Thomas Knudsen (2016) - revise/add regression tests

// Last updated version of proj: 5.0.0

// Original copyright notice:

// Purpose:  Implementation of the aea (Albers Equal Area) projection.
// Author:   Gerald Evenden
// Copyright (c) 1995, Gerald Evenden

// 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_AEA_HPP
#define BOOST_GEOMETRY_PROJECTIONS_AEA_HPP

#include <boost/core/ignore_unused.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_mlfn.hpp>
#include <boost/geometry/srs/projections/impl/pj_msfn.hpp>
#include <boost/geometry/srs/projections/impl/pj_param.hpp>
#include <boost/geometry/srs/projections/impl/pj_qsfn.hpp>


namespace boost { namespace geometry
{

namespace projections
{
    #ifndef DOXYGEN_NO_DETAIL
    namespace detail { namespace aea
    {

            static const double epsilon10 = 1.e-10;
            static const double tolerance7 = 1.e-7;
            static const double epsilon = 1.0e-7;
            static const double tolerance = 1.0e-10;
            static const int n_iter = 15;

            template <typename T>
            struct par_aea
            {
                T    ec;
                T    n;
                T    c;
                T    dd;
                T    n2;
                T    rho0;
                T    phi1;
                T    phi2;
                detail::en<T> en;
                bool ellips;
            };

            /* determine latitude angle phi-1 */
            template <typename T>
            inline T phi1_(T const& qs, T const& Te, T const& Tone_es)
            {
                int i;
                T Phi, sinpi, cospi, con, com, dphi;

                Phi = asin (.5 * qs);
                if (Te < epsilon)
                    return( Phi );
                i = n_iter;
                do {
                    sinpi = sin (Phi);
                    cospi = cos (Phi);
                    con = Te * sinpi;
                    com = 1. - con * con;
                    dphi = .5 * com * com / cospi * (qs / Tone_es -
                       sinpi / com + .5 / Te * log ((1. - con) /
                       (1. + con)));
                    Phi += dphi;
                } while (fabs(dphi) > tolerance && --i);
                return( i ? Phi : HUGE_VAL );
            }

            template <typename T, typename Parameters>
            struct base_aea_ellipsoid
            {
                par_aea<T> m_proj_parm;

                // FORWARD(e_forward)  ellipsoid & spheroid
                // Project coordinates from geographic (lon, lat) to cartesian (x, y)
                inline void fwd(Parameters const& par, T lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
                {
                    T rho = this->m_proj_parm.c - (this->m_proj_parm.ellips
                                                                    ? this->m_proj_parm.n * pj_qsfn(sin(lp_lat), par.e, par.one_es)
                                                                    : this->m_proj_parm.n2 * sin(lp_lat));
                    if (rho < 0.)
                        BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
                    rho = this->m_proj_parm.dd * sqrt(rho);
                    xy_x = rho * sin( lp_lon *= this->m_proj_parm.n );
                    xy_y = this->m_proj_parm.rho0 - rho * cos(lp_lon);
                }

                // INVERSE(e_inverse)  ellipsoid & 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 rho = 0.0;
                    if( (rho = boost::math::hypot(xy_x, xy_y = this->m_proj_parm.rho0 - xy_y)) != 0.0 ) {
                        if (this->m_proj_parm.n < 0.) {
                            rho = -rho;
                            xy_x = -xy_x;
                            xy_y = -xy_y;
                        }
                        lp_lat =  rho / this->m_proj_parm.dd;
                        if (this->m_proj_parm.ellips) {
                            lp_lat = (this->m_proj_parm.c - lp_lat * lp_lat) / this->m_proj_parm.n;
                            if (fabs(this->m_proj_parm.ec - fabs(lp_lat)) > tolerance7) {
                                if ((lp_lat = phi1_(lp_lat, par.e, par.one_es)) == HUGE_VAL)
                                    BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
                            } else
                                lp_lat = lp_lat < 0. ? -half_pi : half_pi;
                        } else if (fabs(lp_lat = (this->m_proj_parm.c - lp_lat * lp_lat) / this->m_proj_parm.n2) <= 1.)
                            lp_lat = asin(lp_lat);
                        else
                            lp_lat = lp_lat < 0. ? -half_pi : half_pi;
                        lp_lon = atan2(xy_x, xy_y) / this->m_proj_parm.n;
                    } else {
                        lp_lon = 0.;
                        lp_lat = this->m_proj_parm.n > 0. ? half_pi : - half_pi;
                    }
                }

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

            };

            template <typename Parameters, typename T>
            inline void setup(Parameters const& par, par_aea<T>& proj_parm) 
            {
                T cosphi, sinphi;
                int secant;

                if (fabs(proj_parm.phi1 + proj_parm.phi2) < epsilon10)
                    BOOST_THROW_EXCEPTION( projection_exception(error_conic_lat_equal) );
                proj_parm.n = sinphi = sin(proj_parm.phi1);
                cosphi = cos(proj_parm.phi1);
                secant = fabs(proj_parm.phi1 - proj_parm.phi2) >= epsilon10;
                if( (proj_parm.ellips = (par.es > 0.))) {
                    T ml1, m1;

                    proj_parm.en = pj_enfn<T>(par.es);
                    m1 = pj_msfn(sinphi, cosphi, par.es);
                    ml1 = pj_qsfn(sinphi, par.e, par.one_es);
                    if (secant) { /* secant cone */
                        T ml2, m2;

                        sinphi = sin(proj_parm.phi2);
                        cosphi = cos(proj_parm.phi2);
                        m2 = pj_msfn(sinphi, cosphi, par.es);
                        ml2 = pj_qsfn(sinphi, par.e, par.one_es);
                        if (ml2 == ml1)
                            BOOST_THROW_EXCEPTION( projection_exception(0) );

                        proj_parm.n = (m1 * m1 - m2 * m2) / (ml2 - ml1);
                    }
                    proj_parm.ec = 1. - .5 * par.one_es * log((1. - par.e) /
                        (1. + par.e)) / par.e;
                    proj_parm.c = m1 * m1 + proj_parm.n * ml1;
                    proj_parm.dd = 1. / proj_parm.n;
                    proj_parm.rho0 = proj_parm.dd * sqrt(proj_parm.c - proj_parm.n * pj_qsfn(sin(par.phi0),
                        par.e, par.one_es));
                } else {
                    if (secant) proj_parm.n = .5 * (proj_parm.n + sin(proj_parm.phi2));
                    proj_parm.n2 = proj_parm.n + proj_parm.n;
                    proj_parm.c = cosphi * cosphi + proj_parm.n2 * sinphi;
                    proj_parm.dd = 1. / proj_parm.n;
                    proj_parm.rho0 = proj_parm.dd * sqrt(proj_parm.c - proj_parm.n2 * sin(par.phi0));
                }
            }


            // Albers Equal Area
            template <typename Params, typename Parameters, typename T>
            inline void setup_aea(Params const& params, Parameters const& par, par_aea<T>& proj_parm)
            {
                proj_parm.phi1 = 0.0;
                proj_parm.phi2 = 0.0;
                bool is_phi1_set = pj_param_r<srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1, proj_parm.phi1);
                bool is_phi2_set = pj_param_r<srs::spar::lat_2>(params, "lat_2", srs::dpar::lat_2, proj_parm.phi2);

                // Boost.Geometry specific, set default parameters manually
                if (! is_phi1_set || ! is_phi2_set) {
                    bool const use_defaults = ! pj_get_param_b<srs::spar::no_defs>(params, "no_defs", srs::dpar::no_defs);
                    if (use_defaults) {
                        if (!is_phi1_set)
                            proj_parm.phi1 = 29.5;
                        if (!is_phi2_set)
                            proj_parm.phi2 = 45.5;
                    }
                }

                setup(par, proj_parm);
            }

            // Lambert Equal Area Conic
            template <typename Params, typename Parameters, typename T>
            inline void setup_leac(Params const& params, Parameters const& par, par_aea<T>& proj_parm)
            {
                static const T half_pi = detail::half_pi<T>();

                proj_parm.phi2 = pj_get_param_r<T, srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1);
                proj_parm.phi1 = pj_get_param_b<srs::spar::south>(params, "south", srs::dpar::south) ? -half_pi : half_pi;
                setup(par, proj_parm);
            }

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

    /*!
        \brief Albers Equal Area projection
        \ingroup projections
        \tparam Geographic latlong point type
        \tparam Cartesian xy point type
        \tparam Parameters parameter type
        \par Projection characteristics
         - Conic
         - Spheroid
         - Ellipsoid
        \par Projection parameters
         - lat_1: Latitude of first standard parallel (degrees)
         - lat_2: Latitude of second standard parallel (degrees)
        \par Example
        \image html ex_aea.gif
    */
    template <typename T, typename Parameters>
    struct aea_ellipsoid : public detail::aea::base_aea_ellipsoid<T, Parameters>
    {
        template <typename Params>
        inline aea_ellipsoid(Params const& params, Parameters const& par)
        {
            detail::aea::setup_aea(params, par, this->m_proj_parm);
        }
    };

    /*!
        \brief Lambert Equal Area Conic projection
        \ingroup projections
        \tparam Geographic latlong point type
        \tparam Cartesian xy point type
        \tparam Parameters parameter type
        \par Projection characteristics
         - Conic
         - Spheroid
         - Ellipsoid
        \par Projection parameters
         - lat_1: Latitude of first standard parallel (degrees)
         - south: Denotes southern hemisphere UTM zone (boolean)
        \par Example
        \image html ex_leac.gif
    */
    template <typename T, typename Parameters>
    struct leac_ellipsoid : public detail::aea::base_aea_ellipsoid<T, Parameters>
    {
        template <typename Params>
        inline leac_ellipsoid(Params const& params, Parameters const& par)
        {
            detail::aea::setup_leac(params, par, this->m_proj_parm);
        }
    };

    #ifndef DOXYGEN_NO_DETAIL
    namespace detail
    {

        // Static projection
        BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_aea, aea_ellipsoid)
        BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_leac, leac_ellipsoid)

        // Factory entry(s)
        BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(aea_entry, aea_ellipsoid)
        BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(leac_entry, leac_ellipsoid)

        BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(aea_init)
        {
            BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(aea, aea_entry)
            BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(leac, leac_entry)
        }

    } // namespace detail
    #endif // doxygen

} // namespace projections

}} // namespace boost::geometry

#endif // BOOST_GEOMETRY_PROJECTIONS_AEA_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 TL	18	// Author: Gerald Evenden (1995)
	19	// Thomas Knudsen (2016) - revise/add regression tests
	20
	21	// Last updated version of proj: 5.0.0
11fdf7f2 TL	22
	23	// Original copyright notice:
	24
	25	// Purpose: Implementation of the aea (Albers Equal Area) projection.
	26	// Author: Gerald Evenden
	27	// Copyright (c) 1995, Gerald Evenden
	28
	29	// Permission is hereby granted, free of charge, to any person obtaining a
	30	// copy of this software and associated documentation files (the "Software"),
	31	// to deal in the Software without restriction, including without limitation
	32	// the rights to use, copy, modify, merge, publish, distribute, sublicense,
	33	// and/or sell copies of the Software, and to permit persons to whom the
	34	// Software is furnished to do so, subject to the following conditions:
	35
	36	// The above copyright notice and this permission notice shall be included
	37	// in all copies or substantial portions of the Software.
	38
	39	// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
	40	// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
	41	// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
	42	// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
	43	// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
	44	// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
	45	// DEALINGS IN THE SOFTWARE.
	46
92f5a8d4 TL	47	#ifndef BOOST_GEOMETRY_PROJECTIONS_AEA_HPP
	48	#define BOOST_GEOMETRY_PROJECTIONS_AEA_HPP
	49
11fdf7f2 TL	50	#include <boost/core/ignore_unused.hpp>
	51	#include <boost/geometry/util/math.hpp>
	52	#include <boost/math/special_functions/hypot.hpp>
	53
	54	#include <boost/geometry/srs/projections/impl/base_static.hpp>
	55	#include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
	56	#include <boost/geometry/srs/projections/impl/projects.hpp>
	57	#include <boost/geometry/srs/projections/impl/factory_entry.hpp>
	58	#include <boost/geometry/srs/projections/impl/pj_mlfn.hpp>
	59	#include <boost/geometry/srs/projections/impl/pj_msfn.hpp>
92f5a8d4	60	#include <boost/geometry/srs/projections/impl/pj_param.hpp>
11fdf7f2 TL	61	#include <boost/geometry/srs/projections/impl/pj_qsfn.hpp>
	62
	63
	64	namespace boost { namespace geometry
	65	{
	66
11fdf7f2 TL	67	namespace projections
	68	{
	69	#ifndef DOXYGEN_NO_DETAIL
	70	namespace detail { namespace aea
	71	{
	72
92f5a8d4 TL	73	static const double epsilon10 = 1.e-10;
	74	static const double tolerance7 = 1.e-7;
	75	static const double epsilon = 1.0e-7;
	76	static const double tolerance = 1.0e-10;
	77	static const int n_iter = 15;
11fdf7f2 TL	78
	79	template <typename T>
	80	struct par_aea
	81	{
	82	T ec;
	83	T n;
	84	T c;
	85	T dd;
	86	T n2;
	87	T rho0;
	88	T phi1;
	89	T phi2;
92f5a8d4 TL	90	detail::en<T> en;
92f5a8d4 TL	91	bool ellips;
11fdf7f2 TL	92	};
	93
	94	/* determine latitude angle phi-1 */
	95	template <typename T>
	96	inline T phi1_(T const& qs, T const& Te, T const& Tone_es)
	97	{
	98	int i;
	99	T Phi, sinpi, cospi, con, com, dphi;
	100
	101	Phi = asin (.5 * qs);
92f5a8d4	102	if (Te < epsilon)
11fdf7f2	103	return( Phi );
92f5a8d4	104	i = n_iter;
11fdf7f2 TL	105	do {
	106	sinpi = sin (Phi);
	107	cospi = cos (Phi);
	108	con = Te * sinpi;
	109	com = 1. - con * con;
	110	dphi = .5 * com * com / cospi * (qs / Tone_es -
	111	sinpi / com + .5 / Te * log ((1. - con) /
	112	(1. + con)));
	113	Phi += dphi;
92f5a8d4	114	} while (fabs(dphi) > tolerance && --i);
11fdf7f2 TL	115	return( i ? Phi : HUGE_VAL );
	116	}
	117
92f5a8d4 TL	118	template <typename T, typename Parameters>
92f5a8d4 TL	119	struct base_aea_ellipsoid
11fdf7f2	120	{
92f5a8d4	121	par_aea<T> m_proj_parm;
11fdf7f2 TL	122
	123	// FORWARD(e_forward) ellipsoid & spheroid
	124	// Project coordinates from geographic (lon, lat) to cartesian (x, y)
92f5a8d4	125	inline void fwd(Parameters const& par, T lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
11fdf7f2	126	{
92f5a8d4 TL	127	T rho = this->m_proj_parm.c - (this->m_proj_parm.ellips
	128	? this->m_proj_parm.n * pj_qsfn(sin(lp_lat), par.e, par.one_es)
	129	: this->m_proj_parm.n2 * sin(lp_lat));
	130	if (rho < 0.)
	131	BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
11fdf7f2 TL	132	rho = this->m_proj_parm.dd * sqrt(rho);
	133	xy_x = rho * sin( lp_lon *= this->m_proj_parm.n );
	134	xy_y = this->m_proj_parm.rho0 - rho * cos(lp_lon);
	135	}
	136
	137	// INVERSE(e_inverse) ellipsoid & spheroid
	138	// Project coordinates from cartesian (x, y) to geographic (lon, lat)
92f5a8d4	139	inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
11fdf7f2	140	{
92f5a8d4	141	static const T half_pi = detail::half_pi<T>();
11fdf7f2	142
92f5a8d4	143	T rho = 0.0;
11fdf7f2 TL	144	if( (rho = boost::math::hypot(xy_x, xy_y = this->m_proj_parm.rho0 - xy_y)) != 0.0 ) {
	145	if (this->m_proj_parm.n < 0.) {
	146	rho = -rho;
	147	xy_x = -xy_x;
	148	xy_y = -xy_y;
	149	}
	150	lp_lat = rho / this->m_proj_parm.dd;
	151	if (this->m_proj_parm.ellips) {
	152	lp_lat = (this->m_proj_parm.c - lp_lat * lp_lat) / this->m_proj_parm.n;
92f5a8d4 TL	153	if (fabs(this->m_proj_parm.ec - fabs(lp_lat)) > tolerance7) {
	154	if ((lp_lat = phi1_(lp_lat, par.e, par.one_es)) == HUGE_VAL)
	155	BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
11fdf7f2	156	} else
92f5a8d4	157	lp_lat = lp_lat < 0. ? -half_pi : half_pi;
11fdf7f2 TL	158	} else if (fabs(lp_lat = (this->m_proj_parm.c - lp_lat * lp_lat) / this->m_proj_parm.n2) <= 1.)
	159	lp_lat = asin(lp_lat);
	160	else
92f5a8d4	161	lp_lat = lp_lat < 0. ? -half_pi : half_pi;
11fdf7f2 TL	162	lp_lon = atan2(xy_x, xy_y) / this->m_proj_parm.n;
	163	} else {
	164	lp_lon = 0.;
92f5a8d4	165	lp_lat = this->m_proj_parm.n > 0. ? half_pi : - half_pi;
11fdf7f2 TL	166	}
	167	}
	168
	169	static inline std::string get_name()
	170	{
	171	return "aea_ellipsoid";
	172	}
	173
	174	};
	175
	176	template <typename Parameters, typename T>
92f5a8d4	177	inline void setup(Parameters const& par, par_aea<T>& proj_parm)
11fdf7f2 TL	178	{
	179	T cosphi, sinphi;
	180	int secant;
	181
92f5a8d4 TL	182	if (fabs(proj_parm.phi1 + proj_parm.phi2) < epsilon10)
92f5a8d4 TL	183	BOOST_THROW_EXCEPTION( projection_exception(error_conic_lat_equal) );
11fdf7f2 TL	184	proj_parm.n = sinphi = sin(proj_parm.phi1);
11fdf7f2 TL	185	cosphi = cos(proj_parm.phi1);
92f5a8d4	186	secant = fabs(proj_parm.phi1 - proj_parm.phi2) >= epsilon10;
11fdf7f2 TL	187	if( (proj_parm.ellips = (par.es > 0.))) {
	188	T ml1, m1;
	189
92f5a8d4	190	proj_parm.en = pj_enfn<T>(par.es);
11fdf7f2 TL	191	m1 = pj_msfn(sinphi, cosphi, par.es);
	192	ml1 = pj_qsfn(sinphi, par.e, par.one_es);
	193	if (secant) { /* secant cone */
	194	T ml2, m2;
	195
	196	sinphi = sin(proj_parm.phi2);
	197	cosphi = cos(proj_parm.phi2);
	198	m2 = pj_msfn(sinphi, cosphi, par.es);
	199	ml2 = pj_qsfn(sinphi, par.e, par.one_es);
92f5a8d4 TL	200	if (ml2 == ml1)
	201	BOOST_THROW_EXCEPTION( projection_exception(0) );
	202
11fdf7f2 TL	203	proj_parm.n = (m1 * m1 - m2 * m2) / (ml2 - ml1);
	204	}
	205	proj_parm.ec = 1. - .5 * par.one_es * log((1. - par.e) /
	206	(1. + par.e)) / par.e;
	207	proj_parm.c = m1 * m1 + proj_parm.n * ml1;
	208	proj_parm.dd = 1. / proj_parm.n;
	209	proj_parm.rho0 = proj_parm.dd * sqrt(proj_parm.c - proj_parm.n * pj_qsfn(sin(par.phi0),
	210	par.e, par.one_es));
	211	} else {
	212	if (secant) proj_parm.n = .5 * (proj_parm.n + sin(proj_parm.phi2));
	213	proj_parm.n2 = proj_parm.n + proj_parm.n;
	214	proj_parm.c = cosphi * cosphi + proj_parm.n2 * sinphi;
	215	proj_parm.dd = 1. / proj_parm.n;
	216	proj_parm.rho0 = proj_parm.dd * sqrt(proj_parm.c - proj_parm.n2 * sin(par.phi0));
	217	}
	218	}
	219
	220
	221	// Albers Equal Area
92f5a8d4 TL	222	template <typename Params, typename Parameters, typename T>
92f5a8d4 TL	223	inline void setup_aea(Params const& params, Parameters const& par, par_aea<T>& proj_parm)
11fdf7f2	224	{
92f5a8d4 TL	225	proj_parm.phi1 = 0.0;
	226	proj_parm.phi2 = 0.0;
	227	bool is_phi1_set = pj_param_r<srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1, proj_parm.phi1);
	228	bool is_phi2_set = pj_param_r<srs::spar::lat_2>(params, "lat_2", srs::dpar::lat_2, proj_parm.phi2);
	229
	230	// Boost.Geometry specific, set default parameters manually
	231	if (! is_phi1_set \|\| ! is_phi2_set) {
	232	bool const use_defaults = ! pj_get_param_b<srs::spar::no_defs>(params, "no_defs", srs::dpar::no_defs);
	233	if (use_defaults) {
	234	if (!is_phi1_set)
	235	proj_parm.phi1 = 29.5;
	236	if (!is_phi2_set)
	237	proj_parm.phi2 = 45.5;
	238	}
	239	}
	240
11fdf7f2 TL	241	setup(par, proj_parm);
	242	}
	243
	244	// Lambert Equal Area Conic
92f5a8d4 TL	245	template <typename Params, typename Parameters, typename T>
92f5a8d4 TL	246	inline void setup_leac(Params const& params, Parameters const& par, par_aea<T>& proj_parm)
11fdf7f2	247	{
92f5a8d4	248	static const T half_pi = detail::half_pi<T>();
11fdf7f2	249
92f5a8d4 TL	250	proj_parm.phi2 = pj_get_param_r<T, srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1);
92f5a8d4 TL	251	proj_parm.phi1 = pj_get_param_b<srs::spar::south>(params, "south", srs::dpar::south) ? -half_pi : half_pi;
11fdf7f2 TL	252	setup(par, proj_parm);
	253	}
	254
	255	}} // namespace detail::aea
	256	#endif // doxygen
	257
	258	/*!
	259	\brief Albers Equal Area projection
	260	\ingroup projections
	261	\tparam Geographic latlong point type
	262	\tparam Cartesian xy point type
	263	\tparam Parameters parameter type
	264	\par Projection characteristics
	265	- Conic
	266	- Spheroid
	267	- Ellipsoid
	268	\par Projection parameters
	269	- lat_1: Latitude of first standard parallel (degrees)
	270	- lat_2: Latitude of second standard parallel (degrees)
	271	\par Example
	272	\image html ex_aea.gif
	273	*/
92f5a8d4 TL	274	template <typename T, typename Parameters>
92f5a8d4 TL	275	struct aea_ellipsoid : public detail::aea::base_aea_ellipsoid<T, Parameters>
11fdf7f2	276	{
92f5a8d4 TL	277	template <typename Params>
92f5a8d4 TL	278	inline aea_ellipsoid(Params const& params, Parameters const& par)
11fdf7f2	279	{
92f5a8d4	280	detail::aea::setup_aea(params, par, this->m_proj_parm);
11fdf7f2 TL	281	}
	282	};
	283
	284	/*!
	285	\brief Lambert Equal Area Conic projection
	286	\ingroup projections
	287	\tparam Geographic latlong point type
	288	\tparam Cartesian xy point type
	289	\tparam Parameters parameter type
	290	\par Projection characteristics
	291	- Conic
	292	- Spheroid
	293	- Ellipsoid
	294	\par Projection parameters
	295	- lat_1: Latitude of first standard parallel (degrees)
	296	- south: Denotes southern hemisphere UTM zone (boolean)
	297	\par Example
	298	\image html ex_leac.gif
	299	*/
92f5a8d4 TL	300	template <typename T, typename Parameters>
92f5a8d4 TL	301	struct leac_ellipsoid : public detail::aea::base_aea_ellipsoid<T, Parameters>
11fdf7f2	302	{
92f5a8d4 TL	303	template <typename Params>
92f5a8d4 TL	304	inline leac_ellipsoid(Params const& params, Parameters const& par)
11fdf7f2	305	{
92f5a8d4	306	detail::aea::setup_leac(params, par, this->m_proj_parm);
11fdf7f2 TL	307	}
	308	};
	309
	310	#ifndef DOXYGEN_NO_DETAIL
	311	namespace detail
	312	{
	313
	314	// Static projection
92f5a8d4 TL	315	BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_aea, aea_ellipsoid)
92f5a8d4 TL	316	BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_leac, leac_ellipsoid)
11fdf7f2 TL	317
11fdf7f2 TL	318	// Factory entry(s)
92f5a8d4 TL	319	BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(aea_entry, aea_ellipsoid)
92f5a8d4 TL	320	BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(leac_entry, leac_ellipsoid)
11fdf7f2	321
92f5a8d4	322	BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(aea_init)
11fdf7f2	323	{
92f5a8d4 TL	324	BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(aea, aea_entry)
92f5a8d4 TL	325	BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(leac, leac_entry)
11fdf7f2 TL	326	}
	327
	328	} // namespace detail
	329	#endif // doxygen
	330
	331	} // namespace projections
	332
	333	}} // namespace boost::geometry
	334
	335	#endif // BOOST_GEOMETRY_PROJECTIONS_AEA_HPP
	336