proj_mdist.hpp 4.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151
  1. // Boost.Geometry (aka GGL, Generic Geometry Library)
  2. // This file is manually converted from PROJ4
  3. // Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands.
  4. // This file was modified by Oracle on 2018.
  5. // Modifications copyright (c) 2018, Oracle and/or its affiliates.
  6. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
  7. // Use, modification and distribution is subject to the Boost Software License,
  8. // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  9. // http://www.boost.org/LICENSE_1_0.txt)
  10. // This file is converted from PROJ4, http://trac.osgeo.org/proj
  11. // PROJ4 is originally written by Gerald Evenden (then of the USGS)
  12. // PROJ4 is maintained by Frank Warmerdam
  13. // PROJ4 is converted to Geometry Library by Barend Gehrels (Geodan, Amsterdam)
  14. // Original copyright notice:
  15. // Permission is hereby granted, free of charge, to any person obtaining a
  16. // copy of this software and associated documentation files (the "Software"),
  17. // to deal in the Software without restriction, including without limitation
  18. // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  19. // and/or sell copies of the Software, and to permit persons to whom the
  20. // Software is furnished to do so, subject to the following conditions:
  21. // The above copyright notice and this permission notice shall be included
  22. // in all copies or substantial portions of the Software.
  23. // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  24. // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  25. // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  26. // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  27. // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  28. // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  29. // DEALINGS IN THE SOFTWARE.
  30. #ifndef BOOST_GEOMETRY_PROJECTIONS_PROJ_MDIST_HPP
  31. #define BOOST_GEOMETRY_PROJECTIONS_PROJ_MDIST_HPP
  32. #include <boost/geometry/srs/projections/exception.hpp>
  33. #include <boost/geometry/srs/projections/impl/pj_strerrno.hpp>
  34. #include <boost/geometry/util/math.hpp>
  35. namespace boost { namespace geometry { namespace projections
  36. {
  37. namespace detail
  38. {
  39. template <typename T>
  40. struct mdist
  41. {
  42. static const int static_size = 20;
  43. T es;
  44. T E;
  45. T b[static_size];
  46. int nb;
  47. };
  48. template <typename T>
  49. inline bool proj_mdist_ini(T const& es, mdist<T>& b)
  50. {
  51. T numf, numfi, twon1, denf, denfi, ens, t, twon;
  52. T den, El, Es;
  53. T E[mdist<T>::static_size];
  54. int i, j;
  55. /* generate E(e^2) and its terms E[] */
  56. ens = es;
  57. numf = twon1 = denfi = 1.;
  58. denf = 1.;
  59. twon = 4.;
  60. Es = El = E[0] = 1.;
  61. for (i = 1; i < mdist<T>::static_size ; ++i)
  62. {
  63. numf *= (twon1 * twon1);
  64. den = twon * denf * denf * twon1;
  65. t = numf/den;
  66. E[i] = t * ens;
  67. Es -= E[i];
  68. ens *= es;
  69. twon *= 4.;
  70. denf *= ++denfi;
  71. twon1 += 2.;
  72. if (Es == El) /* jump out if no change */
  73. break;
  74. El = Es;
  75. }
  76. b.nb = i - 1;
  77. b.es = es;
  78. b.E = Es;
  79. /* generate b_n coefficients--note: collapse with prefix ratios */
  80. b.b[0] = Es = 1. - Es;
  81. numf = denf = 1.;
  82. numfi = 2.;
  83. denfi = 3.;
  84. for (j = 1; j < i; ++j)
  85. {
  86. Es -= E[j];
  87. numf *= numfi;
  88. denf *= denfi;
  89. b.b[j] = Es * numf / denf;
  90. numfi += 2.;
  91. denfi += 2.;
  92. }
  93. return true;
  94. }
  95. template <typename T>
  96. inline T proj_mdist(T const& phi, T const& sphi, T const& cphi, mdist<T> const& b)
  97. {
  98. T sc, sum, sphi2, D;
  99. int i;
  100. sc = sphi * cphi;
  101. sphi2 = sphi * sphi;
  102. D = phi * b.E - b.es * sc / sqrt(1. - b.es * sphi2);
  103. sum = b.b[i = b.nb];
  104. while (i) sum = b.b[--i] + sphi2 * sum;
  105. return(D + sc * sum);
  106. }
  107. template <typename T>
  108. inline T proj_inv_mdist(T const& dist, mdist<T> const& b)
  109. {
  110. static const T TOL = 1e-14;
  111. T s, t, phi, k;
  112. int i;
  113. k = 1./(1.- b.es);
  114. i = mdist<T>::static_size;
  115. phi = dist;
  116. while ( i-- ) {
  117. s = sin(phi);
  118. t = 1. - b.es * s * s;
  119. phi -= t = (proj_mdist(phi, s, cos(phi), b) - dist) *
  120. (t * sqrt(t)) * k;
  121. if (geometry::math::abs(t) < TOL) /* that is no change */
  122. return phi;
  123. }
  124. /* convergence failed */
  125. BOOST_THROW_EXCEPTION( projection_exception(error_non_conv_inv_meri_dist) );
  126. }
  127. } // namespace detail
  128. }}} // namespace boost::geometry::projections
  129. #endif