nmie-applied-impl.hpp 20 KB

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