123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249 |
- #ifndef BOOST_GEOMETRY_PROJECTIONS_LCCA_HPP
- #define BOOST_GEOMETRY_PROJECTIONS_LCCA_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>
- namespace boost { namespace geometry
- {
- namespace projections
- {
- #ifndef DOXYGEN_NO_DETAIL
- namespace detail { namespace lcca
- {
- static const int max_iter = 10;
- static const double del_tol = 1e-12;
- template <typename T>
- struct par_lcca
- {
- detail::en<T> en;
- T r0, l, M0;
- T C;
- };
- template <typename T>
- inline T fS(T const& S, T const& C)
- {
- return(S * ( 1. + S * S * C));
- }
- template <typename T>
- inline T fSp(T const& S, T const& C)
- {
- return(1. + 3.* S * S * C);
- }
- template <typename T, typename Parameters>
- struct base_lcca_ellipsoid
- {
- par_lcca<T> m_proj_parm;
-
-
- inline void fwd(Parameters const& par, T lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
- {
- T S, r, dr;
- S = pj_mlfn(lp_lat, sin(lp_lat), cos(lp_lat), this->m_proj_parm.en) - this->m_proj_parm.M0;
- dr = fS(S, this->m_proj_parm.C);
- r = this->m_proj_parm.r0 - dr;
- xy_x = par.k0 * (r * sin( lp_lon *= this->m_proj_parm.l ) );
- xy_y = par.k0 * (this->m_proj_parm.r0 - r * cos(lp_lon) );
- }
-
-
- inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
- {
- T theta, dr, S, dif;
- int i;
- xy_x /= par.k0;
- xy_y /= par.k0;
- theta = atan2(xy_x , this->m_proj_parm.r0 - xy_y);
- dr = xy_y - xy_x * tan(0.5 * theta);
- lp_lon = theta / this->m_proj_parm.l;
- S = dr;
- for (i = max_iter; i ; --i) {
- S -= (dif = (fS(S, this->m_proj_parm.C) - dr) / fSp(S, this->m_proj_parm.C));
- if (fabs(dif) < del_tol) break;
- }
- if (!i) {
- BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
- }
- lp_lat = pj_inv_mlfn(S + this->m_proj_parm.M0, par.es, this->m_proj_parm.en);
- }
- static inline std::string get_name()
- {
- return "lcca_ellipsoid";
- }
- };
-
- template <typename Parameters, typename T>
- inline void setup_lcca(Parameters const& par, par_lcca<T>& proj_parm)
- {
- T s2p0, N0, R0, tan0;
- proj_parm.en = pj_enfn<T>(par.es);
- if (par.phi0 == 0.) {
- BOOST_THROW_EXCEPTION( projection_exception(error_lat_0_is_zero) );
- }
- proj_parm.l = sin(par.phi0);
- proj_parm.M0 = pj_mlfn(par.phi0, proj_parm.l, cos(par.phi0), proj_parm.en);
- s2p0 = proj_parm.l * proj_parm.l;
- R0 = 1. / (1. - par.es * s2p0);
- N0 = sqrt(R0);
- R0 *= par.one_es * N0;
- tan0 = tan(par.phi0);
- proj_parm.r0 = N0 / tan0;
- proj_parm.C = 1. / (6. * R0 * N0);
- }
- }}
- #endif
-
- template <typename T, typename Parameters>
- struct lcca_ellipsoid : public detail::lcca::base_lcca_ellipsoid<T, Parameters>
- {
- template <typename Params>
- inline lcca_ellipsoid(Params const& , Parameters const& par)
- {
- detail::lcca::setup_lcca(par, this->m_proj_parm);
- }
- };
- #ifndef DOXYGEN_NO_DETAIL
- namespace detail
- {
-
- BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_lcca, lcca_ellipsoid)
-
- BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(lcca_entry, lcca_ellipsoid)
- BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(lcca_init)
- {
- BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(lcca, lcca_entry)
- }
- }
- #endif
- }
- }}
- #endif
|