proj_mdist.hpp 4.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150
  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/util/math.hpp>
  33. namespace boost { namespace geometry { namespace projections
  34. {
  35. namespace detail
  36. {
  37. template <typename T>
  38. struct mdist
  39. {
  40. static const int static_size = 20;
  41. int nb;
  42. T es;
  43. T E;
  44. T b[static_size];
  45. };
  46. template <typename T>
  47. inline bool proj_mdist_ini(T const& es, mdist<T>& b)
  48. {
  49. T numf, numfi, twon1, denf, denfi, ens, t, twon;
  50. T den, El, Es;
  51. T E[mdist<T>::static_size];
  52. int i, j;
  53. /* generate E(e^2) and its terms E[] */
  54. ens = es;
  55. numf = twon1 = denfi = 1.;
  56. denf = 1.;
  57. twon = 4.;
  58. Es = El = E[0] = 1.;
  59. for (i = 1; i < mdist<T>::static_size ; ++i)
  60. {
  61. numf *= (twon1 * twon1);
  62. den = twon * denf * denf * twon1;
  63. t = numf/den;
  64. E[i] = t * ens;
  65. Es -= E[i];
  66. ens *= es;
  67. twon *= 4.;
  68. denf *= ++denfi;
  69. twon1 += 2.;
  70. if (Es == El) /* jump out if no change */
  71. break;
  72. El = Es;
  73. }
  74. b.nb = i - 1;
  75. b.es = es;
  76. b.E = Es;
  77. /* generate b_n coefficients--note: collapse with prefix ratios */
  78. b.b[0] = Es = 1. - Es;
  79. numf = denf = 1.;
  80. numfi = 2.;
  81. denfi = 3.;
  82. for (j = 1; j < i; ++j)
  83. {
  84. Es -= E[j];
  85. numf *= numfi;
  86. denf *= denfi;
  87. b.b[j] = Es * numf / denf;
  88. numfi += 2.;
  89. denfi += 2.;
  90. }
  91. return true;
  92. }
  93. template <typename T>
  94. inline T proj_mdist(T const& phi, T const& sphi, T const& cphi, mdist<T> const& b)
  95. {
  96. T sc, sum, sphi2, D;
  97. int i;
  98. sc = sphi * cphi;
  99. sphi2 = sphi * sphi;
  100. D = phi * b.E - b.es * sc / sqrt(1. - b.es * sphi2);
  101. sum = b.b[i = b.nb];
  102. while (i) sum = b.b[--i] + sphi2 * sum;
  103. return(D + sc * sum);
  104. }
  105. template <typename T>
  106. inline T proj_inv_mdist(T const& dist, mdist<T> const& b)
  107. {
  108. static const T TOL = 1e-14;
  109. T s, t, phi, k;
  110. int i;
  111. k = 1./(1.- b.es);
  112. i = mdist<T>::static_size;
  113. phi = dist;
  114. while ( i-- ) {
  115. s = sin(phi);
  116. t = 1. - b.es * s * s;
  117. phi -= t = (proj_mdist(phi, s, cos(phi), b) - dist) *
  118. (t * sqrt(t)) * k;
  119. if (geometry::math::abs(t) < TOL) /* that is no change */
  120. return phi;
  121. }
  122. /* convergence failed */
  123. BOOST_THROW_EXCEPTION( projection_exception(error_non_conv_inv_meri_dist) );
  124. }
  125. } // namespace detail
  126. }}} // namespace boost::geometry::projections
  127. #endif