123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218 |
- #ifndef BOOST_GEOMETRY_PROJECTIONS_EQDC_HPP
- #define BOOST_GEOMETRY_PROJECTIONS_EQDC_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_msfn.hpp>
- #include <boost/geometry/srs/projections/impl/pj_param.hpp>
- #include <boost/geometry/srs/projections/impl/projects.hpp>
- #include <boost/geometry/util/math.hpp>
- #include <boost/math/special_functions/hypot.hpp>
- namespace boost { namespace geometry
- {
- namespace projections
- {
- #ifndef DOXYGEN_NO_DETAIL
- namespace detail { namespace eqdc
- {
- static const double epsilon10 = 1.e-10;
- template <typename T>
- struct par_eqdc
- {
- T phi1;
- T phi2;
- T n;
- T rho0;
- T c;
- detail::en<T> en;
- bool ellips;
- };
- template <typename T, typename Parameters>
- struct base_eqdc_ellipsoid
- {
- par_eqdc<T> m_proj_parm;
-
-
- inline void fwd(Parameters const& , T lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
- {
- T rho = 0.0;
- rho = this->m_proj_parm.c - (this->m_proj_parm.ellips ? pj_mlfn(lp_lat, sin(lp_lat),
- cos(lp_lat), this->m_proj_parm.en) : lp_lat);
- xy_x = rho * sin( lp_lon *= this->m_proj_parm.n );
- xy_y = this->m_proj_parm.rho0 - rho * cos(lp_lon);
- }
-
-
- inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
- {
- static T const 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 = this->m_proj_parm.c - rho;
- if (this->m_proj_parm.ellips)
- lp_lat = pj_inv_mlfn(lp_lat, par.es, this->m_proj_parm.en);
- 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 "eqdc_ellipsoid";
- }
- };
-
- template <typename Params, typename Parameters, typename T>
- inline void setup_eqdc(Params const& params, Parameters& par, par_eqdc<T>& proj_parm)
- {
- T cosphi, sinphi;
- int secant;
- proj_parm.phi1 = pj_get_param_r<T, srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1);
- proj_parm.phi2 = pj_get_param_r<T, srs::spar::lat_2>(params, "lat_2", srs::dpar::lat_2);
- if (fabs(proj_parm.phi1 + proj_parm.phi2) < epsilon10)
- BOOST_THROW_EXCEPTION( projection_exception(error_conic_lat_equal) );
- proj_parm.en = pj_enfn<T>(par.es);
- 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.)) ) {
- double ml1, m1;
- m1 = pj_msfn(sinphi, cosphi, par.es);
- ml1 = pj_mlfn(proj_parm.phi1, sinphi, cosphi, proj_parm.en);
- if (secant) {
- sinphi = sin(proj_parm.phi2);
- cosphi = cos(proj_parm.phi2);
- proj_parm.n = (m1 - pj_msfn(sinphi, cosphi, par.es)) /
- (pj_mlfn(proj_parm.phi2, sinphi, cosphi, proj_parm.en) - ml1);
- }
- proj_parm.c = ml1 + m1 / proj_parm.n;
- proj_parm.rho0 = proj_parm.c - pj_mlfn(par.phi0, sin(par.phi0),
- cos(par.phi0), proj_parm.en);
- } else {
- if (secant)
- proj_parm.n = (cosphi - cos(proj_parm.phi2)) / (proj_parm.phi2 - proj_parm.phi1);
- proj_parm.c = proj_parm.phi1 + cos(proj_parm.phi1) / proj_parm.n;
- proj_parm.rho0 = proj_parm.c - par.phi0;
- }
- }
- }}
- #endif
-
- template <typename T, typename Parameters>
- struct eqdc_ellipsoid : public detail::eqdc::base_eqdc_ellipsoid<T, Parameters>
- {
- template <typename Params>
- inline eqdc_ellipsoid(Params const& params, Parameters const& par)
- {
- detail::eqdc::setup_eqdc(params, par, this->m_proj_parm);
- }
- };
- #ifndef DOXYGEN_NO_DETAIL
- namespace detail
- {
-
- BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_eqdc, eqdc_ellipsoid)
-
- BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(eqdc_entry, eqdc_ellipsoid)
-
- BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(eqdc_init)
- {
- BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(eqdc, eqdc_entry);
- }
- }
- #endif
- }
- }}
- #endif
|