123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151 |
- #ifndef BOOST_GEOMETRY_PROJECTIONS_PROJ_MDIST_HPP
- #define BOOST_GEOMETRY_PROJECTIONS_PROJ_MDIST_HPP
- #include <boost/geometry/srs/projections/exception.hpp>
- #include <boost/geometry/srs/projections/impl/pj_strerrno.hpp>
- #include <boost/geometry/util/math.hpp>
- namespace boost { namespace geometry { namespace projections
- {
- namespace detail
- {
- template <typename T>
- struct mdist
- {
- static const int static_size = 20;
- T es;
- T E;
- T b[static_size];
- int nb;
- };
- template <typename T>
- inline bool proj_mdist_ini(T const& es, mdist<T>& b)
- {
- T numf, numfi, twon1, denf, denfi, ens, t, twon;
- T den, El, Es;
- T E[mdist<T>::static_size];
- int i, j;
-
- ens = es;
- numf = twon1 = denfi = 1.;
- denf = 1.;
- twon = 4.;
- Es = El = E[0] = 1.;
- for (i = 1; i < mdist<T>::static_size ; ++i)
- {
- numf *= (twon1 * twon1);
- den = twon * denf * denf * twon1;
- t = numf/den;
- E[i] = t * ens;
- Es -= E[i];
- ens *= es;
- twon *= 4.;
- denf *= ++denfi;
- twon1 += 2.;
- if (Es == El)
- break;
- El = Es;
- }
- b.nb = i - 1;
- b.es = es;
- b.E = Es;
-
- b.b[0] = Es = 1. - Es;
- numf = denf = 1.;
- numfi = 2.;
- denfi = 3.;
- for (j = 1; j < i; ++j)
- {
- Es -= E[j];
- numf *= numfi;
- denf *= denfi;
- b.b[j] = Es * numf / denf;
- numfi += 2.;
- denfi += 2.;
- }
- return true;
- }
- template <typename T>
- inline T proj_mdist(T const& phi, T const& sphi, T const& cphi, mdist<T> const& b)
- {
- T sc, sum, sphi2, D;
- int i;
- sc = sphi * cphi;
- sphi2 = sphi * sphi;
- D = phi * b.E - b.es * sc / sqrt(1. - b.es * sphi2);
- sum = b.b[i = b.nb];
- while (i) sum = b.b[--i] + sphi2 * sum;
- return(D + sc * sum);
- }
- template <typename T>
- inline T proj_inv_mdist(T const& dist, mdist<T> const& b)
- {
- static const T TOL = 1e-14;
- T s, t, phi, k;
- int i;
- k = 1./(1.- b.es);
- i = mdist<T>::static_size;
- phi = dist;
- while ( i-- ) {
- s = sin(phi);
- t = 1. - b.es * s * s;
- phi -= t = (proj_mdist(phi, s, cos(phi), b) - dist) *
- (t * sqrt(t)) * k;
- if (geometry::math::abs(t) < TOL)
- return phi;
- }
-
- BOOST_THROW_EXCEPTION( projection_exception(error_non_conv_inv_meri_dist) );
- }
- }
- }}}
- #endif
|