[ceph.git] / ceph / src / boost / boost / geometry / srs / projections / proj / imw_p.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_IMW_P_HPP
#define BOOST_GEOMETRY_PROJECTIONS_IMW_P_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_mlfn.hpp>
#include <boost/geometry/srs/projections/impl/pj_param.hpp>
#include <boost/geometry/srs/projections/impl/projects.hpp>

#include <boost/geometry/util/math.hpp>

namespace boost { namespace geometry
{

namespace projections
{
    #ifndef DOXYGEN_NO_DETAIL
    namespace detail { namespace imw_p
    {

            static const double tolerance = 1e-10;
            static const double epsilon = 1e-10;

            template <typename T>
            struct point_xy { T x, y; }; // specific for IMW_P

            enum mode_type {
                none_is_zero  =  0, /* phi_1 and phi_2 != 0 */
                phi_1_is_zero =  1, /* phi_1 = 0 */
                phi_2_is_zero = -1  /* phi_2 = 0 */
            };

            template <typename T>
            struct par_imw_p
            {
                T    P, Pp, Q, Qp, R_1, R_2, sphi_1, sphi_2, C2;
                T    phi_1, phi_2, lam_1;
                detail::en<T> en;
                mode_type mode;
            };

            template <typename Params, typename T>
            inline int phi12(Params const& params,
                             par_imw_p<T> & proj_parm, T *del, T *sig)
            {
                int err = 0;

                if (!pj_param_r<srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1, proj_parm.phi_1) ||
                    !pj_param_r<srs::spar::lat_2>(params, "lat_2", srs::dpar::lat_2, proj_parm.phi_2)) {
                    err = -41;
                } else {
                    //proj_parm.phi_1 = pj_get_param_r(par.params, "lat_1"); // set above
                    //proj_parm.phi_2 = pj_get_param_r(par.params, "lat_2"); // set above
                    *del = 0.5 * (proj_parm.phi_2 - proj_parm.phi_1);
                    *sig = 0.5 * (proj_parm.phi_2 + proj_parm.phi_1);
                    err = (fabs(*del) < epsilon || fabs(*sig) < epsilon) ? -42 : 0;
                }
                return err;
            }
            template <typename Parameters, typename T>
            inline point_xy<T> loc_for(T const& lp_lam, T const& lp_phi,
                                       Parameters const& par,
                                       par_imw_p<T> const& proj_parm,
                                       T *yc)
            {
                point_xy<T> xy;

                if (lp_phi == 0.0) {
                    xy.x = lp_lam;
                    xy.y = 0.;
                } else {
                    T xa, ya, xb, yb, xc, D, B, m, sp, t, R, C;

                    sp = sin(lp_phi);
                    m = pj_mlfn(lp_phi, sp, cos(lp_phi), proj_parm.en);
                    xa = proj_parm.Pp + proj_parm.Qp * m;
                    ya = proj_parm.P + proj_parm.Q * m;
                    R = 1. / (tan(lp_phi) * sqrt(1. - par.es * sp * sp));
                    C = sqrt(R * R - xa * xa);
                    if (lp_phi < 0.) C = - C;
                    C += ya - R;
                    if (proj_parm.mode == phi_2_is_zero) {
                        xb = lp_lam;
                        yb = proj_parm.C2;
                    } else {
                        t = lp_lam * proj_parm.sphi_2;
                        xb = proj_parm.R_2 * sin(t);
                        yb = proj_parm.C2 + proj_parm.R_2 * (1. - cos(t));
                    }
                    if (proj_parm.mode == phi_1_is_zero) {
                        xc = lp_lam;
                        *yc = 0.;
                    } else {
                        t = lp_lam * proj_parm.sphi_1;
                        xc = proj_parm.R_1 * sin(t);
                        *yc = proj_parm.R_1 * (1. - cos(t));
                    }
                    D = (xb - xc)/(yb - *yc);
                    B = xc + D * (C + R - *yc);
                    xy.x = D * sqrt(R * R * (1 + D * D) - B * B);
                    if (lp_phi > 0)
                        xy.x = - xy.x;
                    xy.x = (B + xy.x) / (1. + D * D);
                    xy.y = sqrt(R * R - xy.x * xy.x);
                    if (lp_phi > 0)
                        xy.y = - xy.y;
                    xy.y += C + R;
                }
                return (xy);
            }
            template <typename Parameters, typename T>
            inline void xy(Parameters const& par, par_imw_p<T> const& proj_parm,
                           T const& phi,
                           T *x, T *y, T *sp, T *R)
            {
                T F;

                *sp = sin(phi);
                *R = 1./(tan(phi) * sqrt(1. - par.es * *sp * *sp ));
                F = proj_parm.lam_1 * *sp;
                *y = *R * (1 - cos(F));
                *x = *R * sin(F);
            }

            template <typename T, typename Parameters>
            struct base_imw_p_ellipsoid
            {
                par_imw_p<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
                {
                    T yc = 0;
                    point_xy<T> xy = loc_for(lp_lon, lp_lat, par, m_proj_parm, &yc);
                    xy_x = xy.x; xy_y = xy.y;
                }

                // INVERSE(e_inverse)  ellipsoid
                // Project coordinates from cartesian (x, y) to geographic (lon, lat)
                inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
                {
                    point_xy<T> t;
                    T yc = 0.0;
                    int i = 0;
                    const int n_max_iter = 1000; /* Arbitrarily choosen number... */

                    lp_lat = this->m_proj_parm.phi_2;
                    lp_lon = xy_x / cos(lp_lat);
                    do {
                        t = loc_for(lp_lon, lp_lat, par, m_proj_parm, &yc);
                        lp_lat = ((lp_lat - this->m_proj_parm.phi_1) * (xy_y - yc) / (t.y - yc)) + this->m_proj_parm.phi_1;
                        lp_lon = lp_lon * xy_x / t.x;
                        i++;
                    } while (i < n_max_iter &&
                             (fabs(t.x - xy_x) > tolerance || fabs(t.y - xy_y) > tolerance));

                    if( i == n_max_iter )
                    {
                        lp_lon = lp_lat = HUGE_VAL;
                    }
                }

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

            };

            // International Map of the World Polyconic
            template <typename Params, typename Parameters, typename T>
            inline void setup_imw_p(Params const& params, Parameters const& par, par_imw_p<T>& proj_parm)
            {
                T del, sig, s, t, x1, x2, T2, y1, m1, m2, y2;
                int err;

                proj_parm.en = pj_enfn<T>(par.es);
                if( (err = phi12(params, proj_parm, &del, &sig)) != 0)
                    BOOST_THROW_EXCEPTION( projection_exception(err) );
                if (proj_parm.phi_2 < proj_parm.phi_1) { /* make sure proj_parm.phi_1 most southerly */
                    del = proj_parm.phi_1;
                    proj_parm.phi_1 = proj_parm.phi_2;
                    proj_parm.phi_2 = del;
                }
                if (pj_param_r<srs::spar::lon_1>(params, "lon_1", srs::dpar::lon_1, proj_parm.lam_1)) {
                    /* empty */
                } else { /* use predefined based upon latitude */
                    sig = fabs(sig * geometry::math::r2d<T>());
                    if (sig <= 60)      sig = 2.;
                    else if (sig <= 76) sig = 4.;
                    else                sig = 8.;
                    proj_parm.lam_1 = sig * geometry::math::d2r<T>();
                }
                proj_parm.mode = none_is_zero;
                if (proj_parm.phi_1 != 0.0)
                    xy(par, proj_parm, proj_parm.phi_1, &x1, &y1, &proj_parm.sphi_1, &proj_parm.R_1);
                else {
                    proj_parm.mode = phi_1_is_zero;
                    y1 = 0.;
                    x1 = proj_parm.lam_1;
                }
                if (proj_parm.phi_2 != 0.0)
                    xy(par, proj_parm, proj_parm.phi_2, &x2, &T2, &proj_parm.sphi_2, &proj_parm.R_2);
                else {
                    proj_parm.mode = phi_2_is_zero;
                    T2 = 0.;
                    x2 = proj_parm.lam_1;
                }
                m1 = pj_mlfn(proj_parm.phi_1, proj_parm.sphi_1, cos(proj_parm.phi_1), proj_parm.en);
                m2 = pj_mlfn(proj_parm.phi_2, proj_parm.sphi_2, cos(proj_parm.phi_2), proj_parm.en);
                t = m2 - m1;
                s = x2 - x1;
                y2 = sqrt(t * t - s * s) + y1;
                proj_parm.C2 = y2 - T2;
                t = 1. / t;
                proj_parm.P = (m2 * y1 - m1 * y2) * t;
                proj_parm.Q = (y2 - y1) * t;
                proj_parm.Pp = (m2 * x1 - m1 * x2) * t;
                proj_parm.Qp = (x2 - x1) * t;
            }

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

    /*!
        \brief International Map of the World Polyconic projection
        \ingroup projections
        \tparam Geographic latlong point type
        \tparam Cartesian xy point type
        \tparam Parameters parameter type
        \par Projection characteristics
         - Mod. Polyconic
         - Ellipsoid
        \par Projection parameters
         - lat_1: Latitude of first standard parallel
         - lat_2: Latitude of second standard parallel
         - lon_1 (degrees)
        \par Example
        \image html ex_imw_p.gif
    */
    template <typename T, typename Parameters>
    struct imw_p_ellipsoid : public detail::imw_p::base_imw_p_ellipsoid<T, Parameters>
    {
        template <typename Params>
        inline imw_p_ellipsoid(Params const& params, Parameters const& par)
        {
            detail::imw_p::setup_imw_p(params, par, this->m_proj_parm);
        }
    };

    #ifndef DOXYGEN_NO_DETAIL
    namespace detail
    {

        // Static projection
        BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_imw_p, imw_p_ellipsoid)

        // Factory entry(s)
        BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(imw_p_entry, imw_p_ellipsoid)
        
        BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(imw_p_init)
        {
            BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(imw_p, imw_p_entry)
        }

    } // namespace detail
    #endif // doxygen

} // namespace projections

}} // namespace boost::geometry

#endif // BOOST_GEOMETRY_PROJECTIONS_IMW_P_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_IMW_P_HPP
92f5a8d4 TL	41	#define BOOST_GEOMETRY_PROJECTIONS_IMW_P_HPP
11fdf7f2 TL	42
	43	#include <boost/geometry/srs/projections/impl/base_static.hpp>
	44	#include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
11fdf7f2 TL	45	#include <boost/geometry/srs/projections/impl/factory_entry.hpp>
11fdf7f2 TL	46	#include <boost/geometry/srs/projections/impl/pj_mlfn.hpp>
92f5a8d4 TL	47	#include <boost/geometry/srs/projections/impl/pj_param.hpp>
92f5a8d4 TL	48	#include <boost/geometry/srs/projections/impl/projects.hpp>
11fdf7f2	49
92f5a8d4	50	#include <boost/geometry/util/math.hpp>
11fdf7f2	51
92f5a8d4	52	namespace boost { namespace geometry
11fdf7f2	53	{
11fdf7f2 TL	54
	55	namespace projections
	56	{
	57	#ifndef DOXYGEN_NO_DETAIL
	58	namespace detail { namespace imw_p
	59	{
	60
92f5a8d4 TL	61	static const double tolerance = 1e-10;
92f5a8d4 TL	62	static const double epsilon = 1e-10;
11fdf7f2 TL	63
11fdf7f2 TL	64	template <typename T>
92f5a8d4 TL	65	struct point_xy { T x, y; }; // specific for IMW_P
	66
	67	enum mode_type {
	68	none_is_zero = 0, /* phi_1 and phi_2 != 0 */
	69	phi_1_is_zero = 1, /* phi_1 = 0 */
	70	phi_2_is_zero = -1 /* phi_2 = 0 */
	71	};
11fdf7f2 TL	72
	73	template <typename T>
	74	struct par_imw_p
	75	{
	76	T P, Pp, Q, Qp, R_1, R_2, sphi_1, sphi_2, C2;
	77	T phi_1, phi_2, lam_1;
92f5a8d4 TL	78	detail::en<T> en;
92f5a8d4 TL	79	mode_type mode;
11fdf7f2 TL	80	};
11fdf7f2 TL	81
92f5a8d4 TL	82	template <typename Params, typename T>
	83	inline int phi12(Params const& params,
	84	par_imw_p<T> & proj_parm, T del, T sig)
11fdf7f2 TL	85	{
	86	int err = 0;
	87
92f5a8d4 TL	88	if (!pj_param_r<srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1, proj_parm.phi_1) \|\|
92f5a8d4 TL	89	!pj_param_r<srs::spar::lat_2>(params, "lat_2", srs::dpar::lat_2, proj_parm.phi_2)) {
11fdf7f2 TL	90	err = -41;
11fdf7f2 TL	91	} else {
92f5a8d4 TL	92	//proj_parm.phi_1 = pj_get_param_r(par.params, "lat_1"); // set above
92f5a8d4 TL	93	//proj_parm.phi_2 = pj_get_param_r(par.params, "lat_2"); // set above
11fdf7f2 TL	94	del = 0.5 (proj_parm.phi_2 - proj_parm.phi_1);
11fdf7f2 TL	95	sig = 0.5 (proj_parm.phi_2 + proj_parm.phi_1);
92f5a8d4	96	err = (fabs(del) < epsilon \|\| fabs(sig) < epsilon) ? -42 : 0;
11fdf7f2 TL	97	}
	98	return err;
	99	}
	100	template <typename Parameters, typename T>
92f5a8d4 TL	101	inline point_xy<T> loc_for(T const& lp_lam, T const& lp_phi,
	102	Parameters const& par,
	103	par_imw_p<T> const& proj_parm,
	104	T *yc)
11fdf7f2	105	{
92f5a8d4	106	point_xy<T> xy;
11fdf7f2	107
92f5a8d4	108	if (lp_phi == 0.0) {
11fdf7f2 TL	109	xy.x = lp_lam;
	110	xy.y = 0.;
	111	} else {
	112	T xa, ya, xb, yb, xc, D, B, m, sp, t, R, C;
	113
	114	sp = sin(lp_phi);
	115	m = pj_mlfn(lp_phi, sp, cos(lp_phi), proj_parm.en);
	116	xa = proj_parm.Pp + proj_parm.Qp * m;
	117	ya = proj_parm.P + proj_parm.Q * m;
	118	R = 1. / (tan(lp_phi) * sqrt(1. - par.es * sp * sp));
	119	C = sqrt(R * R - xa * xa);
	120	if (lp_phi < 0.) C = - C;
	121	C += ya - R;
92f5a8d4	122	if (proj_parm.mode == phi_2_is_zero) {
11fdf7f2 TL	123	xb = lp_lam;
	124	yb = proj_parm.C2;
	125	} else {
	126	t = lp_lam * proj_parm.sphi_2;
	127	xb = proj_parm.R_2 * sin(t);
	128	yb = proj_parm.C2 + proj_parm.R_2 * (1. - cos(t));
	129	}
92f5a8d4	130	if (proj_parm.mode == phi_1_is_zero) {
11fdf7f2 TL	131	xc = lp_lam;
	132	*yc = 0.;
	133	} else {
	134	t = lp_lam * proj_parm.sphi_1;
	135	xc = proj_parm.R_1 * sin(t);
	136	yc = proj_parm.R_1 (1. - cos(t));
	137	}
	138	D = (xb - xc)/(yb - *yc);
	139	B = xc + D * (C + R - *yc);
	140	xy.x = D * sqrt(R * R * (1 + D * D) - B * B);
	141	if (lp_phi > 0)
	142	xy.x = - xy.x;
	143	xy.x = (B + xy.x) / (1. + D * D);
	144	xy.y = sqrt(R * R - xy.x * xy.x);
	145	if (lp_phi > 0)
	146	xy.y = - xy.y;
	147	xy.y += C + R;
	148	}
	149	return (xy);
	150	}
	151	template <typename Parameters, typename T>
92f5a8d4 TL	152	inline void xy(Parameters const& par, par_imw_p<T> const& proj_parm,
	153	T const& phi,
	154	T x, T y, T sp, T R)
11fdf7f2 TL	155	{
	156	T F;
	157
	158	*sp = sin(phi);
	159	R = 1./(tan(phi) sqrt(1. - par.es * sp *sp ));
	160	F = proj_parm.lam_1 * *sp;
	161	y = R * (1 - cos(F));
	162	x = R * sin(F);
	163	}
	164
92f5a8d4 TL	165	template <typename T, typename Parameters>
92f5a8d4 TL	166	struct base_imw_p_ellipsoid
11fdf7f2	167	{
92f5a8d4	168	par_imw_p<T> m_proj_parm;
11fdf7f2 TL	169
	170	// FORWARD(e_forward) ellipsoid
	171	// Project coordinates from geographic (lon, lat) to cartesian (x, y)
92f5a8d4	172	inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
11fdf7f2	173	{
92f5a8d4 TL	174	T yc = 0;
92f5a8d4 TL	175	point_xy<T> xy = loc_for(lp_lon, lp_lat, par, m_proj_parm, &yc);
11fdf7f2 TL	176	xy_x = xy.x; xy_y = xy.y;
	177	}
	178
	179	// INVERSE(e_inverse) ellipsoid
	180	// Project coordinates from cartesian (x, y) to geographic (lon, lat)
92f5a8d4	181	inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
11fdf7f2	182	{
92f5a8d4 TL	183	point_xy<T> t;
	184	T yc = 0.0;
	185	int i = 0;
	186	const int n_max_iter = 1000; /* Arbitrarily choosen number... */
11fdf7f2 TL	187
	188	lp_lat = this->m_proj_parm.phi_2;
	189	lp_lon = xy_x / cos(lp_lat);
	190	do {
92f5a8d4	191	t = loc_for(lp_lon, lp_lat, par, m_proj_parm, &yc);
11fdf7f2 TL	192	lp_lat = ((lp_lat - this->m_proj_parm.phi_1) * (xy_y - yc) / (t.y - yc)) + this->m_proj_parm.phi_1;
11fdf7f2 TL	193	lp_lon = lp_lon * xy_x / t.x;
92f5a8d4 TL	194	i++;
	195	} while (i < n_max_iter &&
	196	(fabs(t.x - xy_x) > tolerance \|\| fabs(t.y - xy_y) > tolerance));
	197
	198	if( i == n_max_iter )
	199	{
	200	lp_lon = lp_lat = HUGE_VAL;
	201	}
11fdf7f2 TL	202	}
	203
	204	static inline std::string get_name()
	205	{
	206	return "imw_p_ellipsoid";
	207	}
	208
	209	};
	210
	211	// International Map of the World Polyconic
92f5a8d4 TL	212	template <typename Params, typename Parameters, typename T>
92f5a8d4 TL	213	inline void setup_imw_p(Params const& params, Parameters const& par, par_imw_p<T>& proj_parm)
11fdf7f2 TL	214	{
11fdf7f2 TL	215	T del, sig, s, t, x1, x2, T2, y1, m1, m2, y2;
92f5a8d4	216	int err;
11fdf7f2	217
92f5a8d4 TL	218	proj_parm.en = pj_enfn<T>(par.es);
	219	if( (err = phi12(params, proj_parm, &del, &sig)) != 0)
	220	BOOST_THROW_EXCEPTION( projection_exception(err) );
11fdf7f2 TL	221	if (proj_parm.phi_2 < proj_parm.phi_1) { /* make sure proj_parm.phi_1 most southerly */
	222	del = proj_parm.phi_1;
	223	proj_parm.phi_1 = proj_parm.phi_2;
	224	proj_parm.phi_2 = del;
	225	}
92f5a8d4 TL	226	if (pj_param_r<srs::spar::lon_1>(params, "lon_1", srs::dpar::lon_1, proj_parm.lam_1)) {
	227	/* empty */
	228	} else { /* use predefined based upon latitude */
11fdf7f2	229	sig = fabs(sig * geometry::math::r2d<T>());
92f5a8d4	230	if (sig <= 60) sig = 2.;
11fdf7f2 TL	231	else if (sig <= 76) sig = 4.;
	232	else sig = 8.;
	233	proj_parm.lam_1 = sig * geometry::math::d2r<T>();
	234	}
92f5a8d4 TL	235	proj_parm.mode = none_is_zero;
	236	if (proj_parm.phi_1 != 0.0)
	237	xy(par, proj_parm, proj_parm.phi_1, &x1, &y1, &proj_parm.sphi_1, &proj_parm.R_1);
11fdf7f2	238	else {
92f5a8d4	239	proj_parm.mode = phi_1_is_zero;
11fdf7f2 TL	240	y1 = 0.;
	241	x1 = proj_parm.lam_1;
	242	}
92f5a8d4 TL	243	if (proj_parm.phi_2 != 0.0)
92f5a8d4 TL	244	xy(par, proj_parm, proj_parm.phi_2, &x2, &T2, &proj_parm.sphi_2, &proj_parm.R_2);
11fdf7f2	245	else {
92f5a8d4	246	proj_parm.mode = phi_2_is_zero;
11fdf7f2 TL	247	T2 = 0.;
	248	x2 = proj_parm.lam_1;
	249	}
	250	m1 = pj_mlfn(proj_parm.phi_1, proj_parm.sphi_1, cos(proj_parm.phi_1), proj_parm.en);
	251	m2 = pj_mlfn(proj_parm.phi_2, proj_parm.sphi_2, cos(proj_parm.phi_2), proj_parm.en);
	252	t = m2 - m1;
	253	s = x2 - x1;
	254	y2 = sqrt(t * t - s * s) + y1;
	255	proj_parm.C2 = y2 - T2;
	256	t = 1. / t;
	257	proj_parm.P = (m2 * y1 - m1 * y2) * t;
	258	proj_parm.Q = (y2 - y1) * t;
	259	proj_parm.Pp = (m2 * x1 - m1 * x2) * t;
	260	proj_parm.Qp = (x2 - x1) * t;
	261	}
	262
	263	}} // namespace detail::imw_p
	264	#endif // doxygen
	265
	266	/*!
	267	\brief International Map of the World Polyconic projection
	268	\ingroup projections
	269	\tparam Geographic latlong point type
	270	\tparam Cartesian xy point type
	271	\tparam Parameters parameter type
	272	\par Projection characteristics
	273	- Mod. Polyconic
	274	- Ellipsoid
	275	\par Projection parameters
	276	- lat_1: Latitude of first standard parallel
	277	- lat_2: Latitude of second standard parallel
	278	- lon_1 (degrees)
	279	\par Example
	280	\image html ex_imw_p.gif
	281	*/
92f5a8d4 TL	282	template <typename T, typename Parameters>
92f5a8d4 TL	283	struct imw_p_ellipsoid : public detail::imw_p::base_imw_p_ellipsoid<T, Parameters>
11fdf7f2	284	{
92f5a8d4 TL	285	template <typename Params>
92f5a8d4 TL	286	inline imw_p_ellipsoid(Params const& params, Parameters const& par)
11fdf7f2	287	{
92f5a8d4	288	detail::imw_p::setup_imw_p(params, par, this->m_proj_parm);
11fdf7f2 TL	289	}
	290	};
	291
	292	#ifndef DOXYGEN_NO_DETAIL
	293	namespace detail
	294	{
	295
	296	// Static projection
92f5a8d4	297	BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_imw_p, imw_p_ellipsoid)
11fdf7f2 TL	298
11fdf7f2 TL	299	// Factory entry(s)
92f5a8d4 TL	300	BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(imw_p_entry, imw_p_ellipsoid)
	301
	302	BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(imw_p_init)
11fdf7f2	303	{
92f5a8d4	304	BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(imw_p, imw_p_entry)
11fdf7f2 TL	305	}
	306
	307	} // namespace detail
	308	#endif // doxygen
	309
	310	} // namespace projections
	311
	312	}} // namespace boost::geometry
	313
	314	#endif // BOOST_GEOMETRY_PROJECTIONS_IMW_P_HPP
	315