mesomie.hpp 5.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116
  1. #ifndef SRC_MESOMIE_HPP_
  2. #define SRC_MESOMIE_HPP_
  3. //******************************************************************************
  4. // Copyright (C) 2009-2022 Ovidio Pena <ovidio@bytesfall.com>
  5. // Copyright (C) 2013-2022 Konstantin Ladutenko <kostyfisik@gmail.com>
  6. //
  7. // This file is part of scattnlay
  8. //
  9. // This program is free software: you can redistribute it and/or modify
  10. // it under the terms of the GNU General Public License as published by
  11. // the Free Software Foundation, either version 3 of the License, or
  12. // (at your option) any later version.
  13. //
  14. // This program is distributed in the hope that it will be useful,
  15. // but WITHOUT ANY WARRANTY; without even the implied warranty of
  16. // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  17. // GNU General Public License for more details.
  18. //
  19. // The only additional remark is that we expect that all publications
  20. // describing work using this software, or all commercial products
  21. // using it, cite at least one of the following references:
  22. // [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by
  23. // a multilayered sphere," Computer Physics Communications,
  24. // vol. 180, Nov. 2009, pp. 2348-2354.
  25. // [2] K. Ladutenko, U. Pal, A. Rivera, and O. Pena-Rodriguez, "Mie
  26. // calculation of electromagnetic near-field for a multilayered
  27. // sphere," Computer Physics Communications, vol. 214, May 2017,
  28. // pp. 225-230.
  29. //
  30. // You should have received a copy of the GNU General Public License
  31. // along with this program. If not, see <http://www.gnu.org/licenses/>.
  32. //******************************************************************************
  33. //******************************************************************************
  34. // This class implements the algorithm for a multilayered sphere described by:
  35. // [1] W. Yang, "Improved recursive algorithm for light scattering by a
  36. // multilayered sphere,” Applied Optics, vol. 42, Mar. 2003, pp.
  37. // 1710-1720.
  38. //
  39. // You can find the description of all the used equations in:
  40. // [2] O. Pena and U. Pal, "Scattering of electromagnetic radiation by
  41. // a multilayered sphere," Computer Physics Communications,
  42. // vol. 180, Nov. 2009, pp. 2348-2354.
  43. // [3] K. Ladutenko, U. Pal, A. Rivera, and O. Pena-Rodriguez, "Mie
  44. // calculation of electromagnetic near-field for a multilayered
  45. // sphere," Computer Physics Communications, vol. 214, May 2017,
  46. // pp. 225-230.
  47. //
  48. // Hereinafter all equations numbers refer to [2]
  49. //******************************************************************************
  50. #include <iomanip>
  51. #include <iostream>
  52. #include <stdexcept>
  53. #include <vector>
  54. #include "nmie.hpp"
  55. #include "special-functions-impl.hpp"
  56. namespace nmie {
  57. template <typename FloatType>
  58. //******************************************************************************
  59. void MesoMie<FloatType>::calc_ab(int nmax,
  60. FloatType R,
  61. FloatType xd,
  62. FloatType xm,
  63. FloatType eps_d,
  64. FloatType eps_m,
  65. FloatType d_parallel,
  66. FloatType d_perp) {
  67. an_.resize(nmax, static_cast<FloatType>(0.0));
  68. bn_.resize(nmax, static_cast<FloatType>(0.0));
  69. std::vector<std::complex<FloatType>> //
  70. D1_xd(nmax + 1), D3_xd(nmax + 1), //
  71. D1_xm(nmax + 1), D3_xm(nmax + 1), //
  72. Psi_xd(nmax + 1), Zeta_xd(nmax + 1), //
  73. Psi_xm(nmax + 1), Zeta_xm(nmax + 1);
  74. evalPsiZetaD1D3(std::complex<FloatType>(xd), Psi_xd, Zeta_xd, D1_xd, D3_xd);
  75. evalPsiZetaD1D3(std::complex<FloatType>(xm), Psi_xm, Zeta_xm, D1_xm, D3_xm);
  76. for (int n = 1; n <= nmax; n++) {
  77. an_[n] = Psi_xd[n] *
  78. ( //
  79. eps_m * xd * D1_xd[n] - eps_d * xm * D1_xm[n] + //
  80. (eps_m - eps_d) * //
  81. ( //
  82. n * (n + 1) * d_perp + //
  83. xd * D1_xd[n] * xm * D1_xm[n] * d_parallel //
  84. ) / //
  85. R //
  86. ) / //
  87. Zeta_xd[n] *
  88. ( //
  89. eps_m * xd * D3_xd[n] - eps_d * xm * D1_xm[n] + //
  90. (eps_m - eps_d) * //
  91. ( //
  92. n * (n + 1) * d_perp + //
  93. xd * D3_xd[n] * xm * D1_xm[n] * d_parallel //
  94. ) / //
  95. R //
  96. );
  97. bn_[n] = Psi_xd[n] *
  98. ( //
  99. xd * D1_xd[n] - xm * D1_xm[n] + //
  100. (xm * xm - xd * xd) * d_parallel / R //
  101. ) / //
  102. Zeta_xd[n] *
  103. ( //
  104. xd * D3_xd[n] - xm * D1_xm[n] + //
  105. (xm * xm - xd * xd) * d_parallel / R //
  106. ); //
  107. }
  108. }
  109. } // namespace nmie
  110. #endif