nmie-applied-impl.hpp 22 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430
  1. #ifndef SRC_NMIE_APPLIED_IMPL_HPP_
  2. #define SRC_NMIE_APPLIED_IMPL_HPP_
  3. //**********************************************************************************//
  4. // Copyright (C) 2009-2018 Ovidio Pena <ovidio@bytesfall.com> //
  5. // Copyright (C) 2013-2018 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. // @brief Wrapper class around nMie function for ease of use //
  34. // //
  35. //**********************************************************************************//
  36. #include "nmie-applied.hpp"
  37. #include "nmie-precision.hpp"
  38. #include <array>
  39. #include <algorithm>
  40. #include <cstdio>
  41. #include <cstdlib>
  42. #include <stdexcept>
  43. #include <vector>
  44. namespace nmie {
  45. // ********************************************************************** //
  46. // ********************************************************************** //
  47. // ********************************************************************** //
  48. template <typename FloatType>
  49. void MultiLayerMieApplied<FloatType>::GetFailed() {
  50. FloatType faild_x = 9.42477796076938;
  51. //FloatType faild_x = 9.42477796076937;
  52. std::complex<FloatType> z(faild_x, 0.0);
  53. std::vector<int> nmax_local_array = {20, 100, 500, 2500};
  54. for (auto nmax_local : nmax_local_array) {
  55. std::vector<std::complex<FloatType> > D1_failed(nmax_local + 1);
  56. // Downward recurrence for D1 - equations (16a) and (16b)
  57. D1_failed[nmax_local] = std::complex<FloatType>(0.0, 0.0);
  58. const std::complex<FloatType> zinv = std::complex<FloatType>(1.0, 0.0)/z;
  59. for (int n = nmax_local; n > 0; n--) {
  60. D1_failed[n - 1] = FloatType(n)*zinv - 1.0/(D1_failed[n] + FloatType(n)*zinv);
  61. }
  62. printf("Faild D1[0] from reccurence (z = %16.14f, nmax = %d): %g\n",
  63. faild_x, nmax_local, D1_failed[0].real());
  64. }
  65. printf("Faild D1[0] from continued fraction (z = %16.14f): %g\n", faild_x,
  66. calcD1confra(0,z).real());
  67. //D1[nmax_] = calcD1confra(nmax_, z);
  68. }
  69. // ********************************************************************** //
  70. // ********************************************************************** //
  71. // ********************************************************************** //
  72. template <typename FloatType>
  73. void MultiLayerMieApplied<FloatType>::AddTargetLayer(FloatType width, std::complex<FloatType> layer_index) {
  74. this->MarkUncalculated();
  75. if (width <= 0)
  76. throw std::invalid_argument("Layer width should be positive!");
  77. target_width_.push_back(width);
  78. target_index_.push_back(layer_index);
  79. } // end of void MultiLayerMieApplied<FloatType>::AddTargetLayer(...)
  80. // ********************************************************************** //
  81. // ********************************************************************** //
  82. // ********************************************************************** //
  83. template <typename FloatType>
  84. void MultiLayerMieApplied<FloatType>::AddTargetLayerReIm(FloatType width,
  85. FloatType re_layer_index, FloatType im_layer_index) {
  86. this->MarkUncalculated();
  87. if (width <= 0)
  88. throw std::invalid_argument("Layer width should be positive!");
  89. target_width_.push_back(width);
  90. target_index_.push_back(std::complex<FloatType>(re_layer_index,im_layer_index));
  91. } // end of void MultiLayerMieApplied<FloatType>::AddTargetLayer(...)
  92. // ********************************************************************** //
  93. // ********************************************************************** //
  94. // ********************************************************************** //
  95. template <typename FloatType>
  96. void MultiLayerMieApplied<FloatType>::SetTargetPEC(FloatType radius) {
  97. this->MarkUncalculated();
  98. if (target_width_.size() != 0 || target_index_.size() != 0)
  99. throw std::invalid_argument("Error! Define PEC target radius before any other layers!");
  100. // Add layer of any index...
  101. AddTargetLayer(radius, std::complex<FloatType>(0.0, 0.0));
  102. // ... and mark it as PEC
  103. this->SetPECLayer(0);
  104. }
  105. // ********************************************************************** //
  106. // ********************************************************************** //
  107. // ********************************************************************** //
  108. template <typename FloatType>
  109. void MultiLayerMieApplied<FloatType>::SetCoatingIndex(std::vector<std::complex<FloatType> > index) {
  110. this->MarkUncalculated();
  111. coating_index_.clear();
  112. for (auto value : index) coating_index_.push_back(value);
  113. } // end of void MultiLayerMieApplied<FloatType>::SetCoatingIndex(std::vector<complex> index);
  114. // ********************************************************************** //
  115. // ********************************************************************** //
  116. // ********************************************************************** //
  117. template <typename FloatType>
  118. void MultiLayerMieApplied<FloatType>::SetCoatingWidth(std::vector<FloatType> width) {
  119. this->MarkUncalculated();
  120. coating_width_.clear();
  121. for (auto w : width)
  122. if (w <= 0)
  123. throw std::invalid_argument("Coating width should be positive!");
  124. else coating_width_.push_back(w);
  125. }
  126. // end of void MultiLayerMieApplied<FloatType>::SetCoatingWidth(...);
  127. // ********************************************************************** //
  128. // ********************************************************************** //
  129. // ********************************************************************** //
  130. template <typename FloatType>
  131. void MultiLayerMieApplied<FloatType>::SetWidthSP(const std::vector<FloatType>& size_parameter) {
  132. this->MarkUncalculated();
  133. this->size_param_.clear();
  134. FloatType prev_size_parameter = 0.0;
  135. for (auto layer_size_parameter : size_parameter) {
  136. if (layer_size_parameter <= 0.0)
  137. throw std::invalid_argument("Size parameter should be positive!");
  138. if (prev_size_parameter > layer_size_parameter)
  139. throw std::invalid_argument
  140. ("Size parameter for next layer should be larger than the previous one!");
  141. prev_size_parameter = layer_size_parameter;
  142. this->size_param_.push_back(layer_size_parameter);
  143. }
  144. }
  145. // end of void MultiLayerMieApplied<FloatType>::SetWidthSP(...);
  146. // ********************************************************************** //
  147. // ********************************************************************** //
  148. // ********************************************************************** //
  149. template <typename FloatType>
  150. void MultiLayerMieApplied<FloatType>::SetIndexSP(const std::vector< std::complex<FloatType> >& index) {
  151. this->MarkUncalculated();
  152. //refractive_index_.clear();
  153. this->refractive_index_ = index;
  154. // for (auto value : index) refractive_index_.push_back(value);
  155. } // end of void MultiLayerMieApplied<FloatType>::SetIndexSP(...);
  156. // ********************************************************************** //
  157. // ********************************************************************** //
  158. // ********************************************************************** //
  159. template <typename FloatType>
  160. void MultiLayerMieApplied<FloatType>::SetFieldPointsSP(const std::vector< std::vector<FloatType> >& coords_sp) {
  161. if (coords_sp.size() != 3)
  162. throw std::invalid_argument("Error! Wrong dimension of field monitor points!");
  163. if (coords_sp[0].size() != coords_sp[1].size() || coords_sp[0].size() != coords_sp[2].size())
  164. throw std::invalid_argument("Error! Missing coordinates for field monitor points!");
  165. this->coords_ = coords_sp;
  166. // for (int i = 0; i < coords_sp_[0].size(); ++i) {
  167. // printf("%g, %g, %g\n", coords_sp_[0][i], coords_sp_[1][i], coords_sp_[2][i]);
  168. // }
  169. } // end of void MultiLayerMieApplied<FloatType>::SetFieldPointsSP(...)
  170. // ********************************************************************** //
  171. // ********************************************************************** //
  172. // ********************************************************************** //
  173. template <typename FloatType>
  174. void MultiLayerMieApplied<FloatType>::GenerateSizeParameter() {
  175. this->MarkUncalculated();
  176. this->size_param_.clear();
  177. FloatType radius = 0.0;
  178. for (auto width : target_width_) {
  179. radius += width;
  180. this->size_param_.push_back(2*this->PI_*radius/wavelength_);
  181. }
  182. for (auto width : coating_width_) {
  183. radius += width;
  184. this->size_param_.push_back(2*this->PI_*radius/wavelength_);
  185. }
  186. this->total_radius_ = radius;
  187. } // end of void MultiLayerMieApplied<FloatType>::GenerateSizeParameter();
  188. // ********************************************************************** //
  189. // ********************************************************************** //
  190. // ********************************************************************** //
  191. template <typename FloatType>
  192. void MultiLayerMieApplied<FloatType>::GenerateIndex() {
  193. this->MarkUncalculated();
  194. this->refractive_index_.clear();
  195. for (auto index : this->target_index_)
  196. this->refractive_index_.push_back(index);
  197. for (auto index : this->coating_index_)
  198. this->refractive_index_.push_back(index);
  199. } // end of void MultiLayerMieApplied<FloatType>::GenerateIndex();
  200. // ********************************************************************** //
  201. // ********************************************************************** //
  202. // ********************************************************************** //
  203. template <typename FloatType>
  204. FloatType MultiLayerMieApplied<FloatType>::GetTotalRadius() {
  205. if (!this->isMieCalculated()) GenerateSizeParameter();
  206. return this->total_radius_;
  207. } // end of FloatType MultiLayerMieApplied<FloatType>::GetTotalRadius();
  208. // ********************************************************************** //
  209. // ********************************************************************** //
  210. // ********************************************************************** //
  211. template <typename FloatType>
  212. std::vector< std::vector<FloatType> >
  213. MultiLayerMieApplied<FloatType>::GetSpectra(FloatType from_WL, FloatType to_WL, int samples) {
  214. if (!this->isMieCalculated())
  215. throw std::invalid_argument("You should run calculations before result request!");
  216. std::vector< std::vector<FloatType> > spectra;
  217. FloatType step_WL = (to_WL - from_WL)/static_cast<FloatType>(samples);
  218. FloatType wavelength_backup = wavelength_;
  219. long fails = 0;
  220. for (FloatType WL = from_WL; WL < to_WL; WL += step_WL) {
  221. wavelength_ = WL;
  222. try {
  223. RunMieCalculation();
  224. } catch(const std::invalid_argument& ia) {
  225. fails++;
  226. continue;
  227. }
  228. //printf("%3.1f ",WL);
  229. spectra.push_back(std::vector<FloatType>({wavelength_, this->GetQext(),
  230. this->GetQsca(), this->GetQabs(), this->GetQbk()}));
  231. } // end of for each WL in spectra
  232. printf("Spectrum has %li fails\n",fails);
  233. wavelength_ = wavelength_backup;
  234. return spectra;
  235. }
  236. // ********************************************************************** //
  237. // ********************************************************************** //
  238. // ********************************************************************** //
  239. template <typename FloatType>
  240. void MultiLayerMieApplied<FloatType>::ClearTarget() {
  241. this->MarkUncalculated();
  242. this->target_width_.clear();
  243. this->target_index_.clear();
  244. }
  245. // ********************************************************************** //
  246. // ********************************************************************** //
  247. // ********************************************************************** //
  248. template <typename FloatType>
  249. void MultiLayerMieApplied<FloatType>::ClearCoating() {
  250. this->MarkUncalculated();
  251. this->coating_width_.clear();
  252. this->coating_index_.clear();
  253. }
  254. // ********************************************************************** //
  255. // ********************************************************************** //
  256. // ********************************************************************** //
  257. template <typename FloatType>
  258. void MultiLayerMieApplied<FloatType>::ClearLayers() {
  259. this->MarkUncalculated();
  260. this->ClearTarget();
  261. this->ClearCoating();
  262. }
  263. // ********************************************************************** //
  264. // ********************************************************************** //
  265. // ********************************************************************** //
  266. template <typename FloatType>
  267. void MultiLayerMieApplied<FloatType>::ClearAllDesign() {
  268. this->MarkUncalculated();
  269. this->ClearLayers();
  270. this->size_param_.clear();
  271. this->refractive_index_.clear();
  272. }
  273. // ********************************************************************** //
  274. // ********************************************************************** //
  275. // ********************************************************************** //
  276. // Computational core
  277. // ********************************************************************** //
  278. // ********************************************************************** //
  279. // ********************************************************************** //
  280. //**********************************************************************************//
  281. // Function CONFRA ported from MIEV0.f (Wiscombe,1979)
  282. // Ref. to NCAR Technical Notes, Wiscombe, 1979
  283. /*
  284. c Compute Bessel function ratio A-sub-N from its
  285. c continued fraction using Lentz method
  286. c ZINV = Reciprocal of argument of A
  287. c I N T E R N A L V A R I A B L E S
  288. c ------------------------------------
  289. c CAK Term in continued fraction expansion of A (Eq. R25)
  290. c a_k
  291. c CAPT Factor used in Lentz iteration for A (Eq. R27)
  292. c T_k
  293. c CNUMER Numerator in capT (Eq. R28A)
  294. c N_k
  295. c CDENOM Denominator in capT (Eq. R28B)
  296. c D_k
  297. c CDTD Product of two successive denominators of capT factors
  298. c (Eq. R34C)
  299. c xi_1
  300. c CNTN Product of two successive numerators of capT factors
  301. c (Eq. R34B)
  302. c xi_2
  303. c EPS1 Ill-conditioning criterion
  304. c EPS2 Convergence criterion
  305. c KK Subscript k of cAk (Eq. R25B)
  306. c k
  307. c KOUNT Iteration counter (used to prevent infinite looping)
  308. c MAXIT Max. allowed no. of iterations
  309. c MM + 1 and - 1, alternately
  310. */
  311. template <typename FloatType>
  312. std::complex<FloatType> MultiLayerMieApplied<FloatType>::calcD1confra(const int N, const std::complex<FloatType> z) {
  313. // NTMR -> nmax_ - 1 \\TODO nmax_ ?
  314. //int N = nmax_ - 1;
  315. int KK, KOUNT, MAXIT = 10000, MM;
  316. // FloatType EPS1=1.0e-2;
  317. FloatType EPS2=1.0e-8;
  318. std::complex<FloatType> CAK, CAPT, CDENOM, CDTD, CNTN, CNUMER;
  319. std::complex<FloatType> one = std::complex<FloatType>(1.0,0.0);
  320. std::complex<FloatType> ZINV = one/z;
  321. // c ** Eq. R25a
  322. std::complex<FloatType> CONFRA = static_cast<std::complex<FloatType> >(N + 1)*ZINV; //debug ZINV
  323. MM = - 1;
  324. KK = 2*N +3; //debug 3
  325. // c ** Eq. R25b, k=2
  326. CAK = static_cast<std::complex<FloatType> >(MM*KK)*ZINV; //debug -3 ZINV
  327. CDENOM = CAK;
  328. CNUMER = CDENOM + one/CONFRA; //-3zinv+z
  329. KOUNT = 1;
  330. //10 CONTINUE
  331. do { ++KOUNT;
  332. if (KOUNT > MAXIT) {
  333. printf("re(%g):im(%g)\t\n", CONFRA.real(), CONFRA.imag());
  334. throw std::invalid_argument("ConFra--Iteration failed to converge!\n");
  335. }
  336. MM *= - 1; KK += 2; //debug mm=1 kk=5
  337. CAK = static_cast<std::complex<FloatType> >(MM*KK)*ZINV; // ** Eq. R25b //debug 5zinv
  338. // //c ** Eq. R32 Ill-conditioned case -- stride two terms instead of one
  339. // if (std::abs(CNUMER/CAK) >= EPS1 || std::abs(CDENOM/CAK) >= EPS1) {
  340. // //c ** Eq. R34
  341. // CNTN = CAK*CNUMER + 1.0;
  342. // CDTD = CAK*CDENOM + 1.0;
  343. // CONFRA = (CNTN/CDTD)*CONFRA; // ** Eq. R33
  344. // MM *= - 1; KK += 2;
  345. // CAK = static_cast<std::complex<FloatType> >(MM*KK)*ZINV; // ** Eq. R25b
  346. // //c ** Eq. R35
  347. // CNUMER = CAK + CNUMER/CNTN;
  348. // CDENOM = CAK + CDENOM/CDTD;
  349. // ++KOUNT;
  350. // //GO TO 10
  351. // continue;
  352. // } else { //c *** Well-conditioned case
  353. {
  354. CAPT = CNUMER/CDENOM; // ** Eq. R27 //debug (-3zinv + z)/(-3zinv)
  355. // printf("re(%g):im(%g)**\t", CAPT.real(), CAPT.imag());
  356. CONFRA = CAPT*CONFRA; // ** Eq. R26
  357. //if (N == 0) {output=true;printf(" re:");prn(CONFRA.real());printf(" im:"); prn(CONFRA.imag());output=false;};
  358. //c ** Check for convergence; Eq. R31
  359. if (std::abs(CAPT.real() - 1.0) >= EPS2 || std::abs(CAPT.imag()) >= EPS2) {
  360. //c ** Eq. R30
  361. CNUMER = CAK + one/CNUMER;
  362. CDENOM = CAK + one/CDENOM;
  363. continue;
  364. //GO TO 10
  365. } // end of if < eps2
  366. }
  367. break;
  368. } while(1);
  369. //if (N == 0) printf(" return confra for z=(%g,%g)\n", ZINV.real(), ZINV.imag());
  370. return CONFRA;
  371. }
  372. // ********************************************************************** //
  373. // ********************************************************************** //
  374. // ********************************************************************** //
  375. template <typename FloatType>
  376. void MultiLayerMieApplied<FloatType>::ConvertToSP() {
  377. this->MarkUncalculated();
  378. if (target_width_.size() + coating_width_.size() == 0)
  379. return; // Nothing to convert, we suppose that SP was set directly
  380. GenerateSizeParameter();
  381. GenerateIndex();
  382. if (this->size_param_.size() != this->refractive_index_.size())
  383. throw std::invalid_argument("Ivalid conversion of width to size parameter units!/n");
  384. }
  385. // ********************************************************************** //
  386. // ********************************************************************** //
  387. // ********************************************************************** //
  388. template <typename FloatType>
  389. void MultiLayerMieApplied<FloatType>::RunMieCalculation() {
  390. ConvertToSP();
  391. this->MultiLayerMie<FloatType>::RunMieCalculation();
  392. }
  393. // ********************************************************************** //
  394. // ********************************************************************** //
  395. // ********************************************************************** //
  396. template <typename FloatType>
  397. void MultiLayerMieApplied<FloatType>::GetExpanCoeffs( std::vector< std::vector<std::complex<FloatType> > >& aln, std::vector< std::vector<std::complex<FloatType> > >& bln, std::vector< std::vector<std::complex<FloatType> > >& cln, std::vector< std::vector<std::complex<FloatType> > >& dln) {
  398. ConvertToSP(); // Case of call before running full Mie calculation.
  399. // Calculate scattering coefficients an_ and bn_
  400. this->calcScattCoeffs();
  401. // Calculate expansion coefficients aln_, bln_, cln_, and dln_
  402. this->calcExpanCoeffs();
  403. aln = this->aln_;
  404. bln = this->bln_;
  405. cln = this->cln_;
  406. dln = this->dln_;
  407. } // end of void MultiLayerMieApplied<FloatType>::GetExpanCoeffs( ...)
  408. // ********************************************************************** //
  409. // ********************************************************************** //
  410. // ********************************************************************** //
  411. } // end of namespace nmie
  412. #endif // SRC_NMIE_APPLIED_IMPL_HPP_