nmie.cc 69 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328
  1. //**********************************************************************************//
  2. // Copyright (C) 2009-2015 Ovidio Pena <ovidio@bytesfall.com> //
  3. // Copyright (C) 2013-2015 Konstantin Ladutenko <kostyfisik@gmail.com> //
  4. // //
  5. // This file is part of scattnlay //
  6. // //
  7. // This program is free software: you can redistribute it and/or modify //
  8. // it under the terms of the GNU General Public License as published by //
  9. // the Free Software Foundation, either version 3 of the License, or //
  10. // (at your option) any later version. //
  11. // //
  12. // This program is distributed in the hope that it will be useful, //
  13. // but WITHOUT ANY WARRANTY; without even the implied warranty of //
  14. // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
  15. // GNU General Public License for more details. //
  16. // //
  17. // The only additional remark is that we expect that all publications //
  18. // describing work using this software, or all commercial products //
  19. // using it, cite the following reference: //
  20. // [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by //
  21. // a multilayered sphere," Computer Physics Communications, //
  22. // vol. 180, Nov. 2009, pp. 2348-2354. //
  23. // //
  24. // You should have received a copy of the GNU General Public License //
  25. // along with this program. If not, see <http://www.gnu.org/licenses/>. //
  26. //**********************************************************************************//
  27. //**********************************************************************************//
  28. // This class implements the algorithm for a multilayered sphere described by: //
  29. // [1] W. Yang, "Improved recursive algorithm for light scattering by a //
  30. // multilayered sphere,” Applied Optics, vol. 42, Mar. 2003, pp. 1710-1720. //
  31. // //
  32. // You can find the description of all the used equations in: //
  33. // [2] O. Pena and U. Pal, "Scattering of electromagnetic radiation by //
  34. // a multilayered sphere," Computer Physics Communications, //
  35. // vol. 180, Nov. 2009, pp. 2348-2354. //
  36. // //
  37. // Hereinafter all equations numbers refer to [2] //
  38. //**********************************************************************************//
  39. #include "nmie.h"
  40. #include <array>
  41. #include <algorithm>
  42. #include <cstdio>
  43. #include <cstdlib>
  44. #include <stdexcept>
  45. #include <vector>
  46. using namespace std;
  47. namespace nmie {
  48. //helpers
  49. template<class T> inline T pow2(const T value) {return value*value;}
  50. int round(double x) {
  51. return x >= 0 ? (int)(x + 0.5):(int)(x - 0.5);
  52. }
  53. //**********************************************************************************//
  54. // This function emulates a C call to calculate the scattering coefficients //
  55. // required to calculate both the near- and far-field parameters. //
  56. // //
  57. // Input parameters: //
  58. // L: Number of layers //
  59. // pl: Index of PEC layer. If there is none just send -1 //
  60. // x: Array containing the size parameters of the layers [0..L-1] //
  61. // m: Array containing the relative refractive indexes of the layers [0..L-1] //
  62. // nmax: Maximum number of multipolar expansion terms to be used for the //
  63. // calculations. Only use it if you know what you are doing, otherwise //
  64. // set this parameter to -1 and the function will calculate it. //
  65. // //
  66. // Output parameters: //
  67. // an, bn: Complex scattering amplitudes //
  68. // //
  69. // Return value: //
  70. // Number of multipolar expansion terms used for the calculations //
  71. //**********************************************************************************//
  72. int ScattCoeffs(const unsigned int L, const int pl, vector<double>& x, vector<complex<double> >& m, const int nmax, vector<complex<double> >& an, vector<complex<double> >& bn) {
  73. if (x.size() != L || m.size() != L)
  74. throw invalid_argument("Declared number of layers do not fit x and m!");
  75. try {
  76. MultiLayerMie ml_mie;
  77. ml_mie.SetAllLayersSize(x);
  78. ml_mie.SetAllLayersIndex(m);
  79. ml_mie.SetPECLayer(pl);
  80. ml_mie.SetMaxTerms(nmax);
  81. ml_mie.calcScattCoeffs();
  82. an = ml_mie.GetAn();
  83. bn = ml_mie.GetBn();
  84. return ml_mie.GetMaxTerms();
  85. } catch(const invalid_argument& ia) {
  86. // Will catch if ml_mie fails or other errors.
  87. cerr << "Invalid argument: " << ia.what() << endl;
  88. throw invalid_argument(ia);
  89. return -1;
  90. }
  91. return 0;
  92. }
  93. //**********************************************************************************//
  94. // This function emulates a C call to calculate the actual scattering parameters //
  95. // and amplitudes. //
  96. // //
  97. // Input parameters: //
  98. // L: Number of layers //
  99. // pl: Index of PEC layer. If there is none just send -1 //
  100. // x: Array containing the size parameters of the layers [0..L-1] //
  101. // m: Array containing the relative refractive indexes of the layers [0..L-1] //
  102. // nTheta: Number of scattering angles //
  103. // Theta: Array containing all the scattering angles where the scattering //
  104. // amplitudes will be calculated //
  105. // nmax: Maximum number of multipolar expansion terms to be used for the //
  106. // calculations. Only use it if you know what you are doing, otherwise //
  107. // set this parameter to -1 and the function will calculate it //
  108. // //
  109. // Output parameters: //
  110. // Qext: Efficiency factor for extinction //
  111. // Qsca: Efficiency factor for scattering //
  112. // Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca) //
  113. // Qbk: Efficiency factor for backscattering //
  114. // Qpr: Efficiency factor for the radiation pressure //
  115. // g: Asymmetry factor (g = (Qext-Qpr)/Qsca) //
  116. // Albedo: Single scattering albedo (Albedo = Qsca/Qext) //
  117. // S1, S2: Complex scattering amplitudes //
  118. // //
  119. // Return value: //
  120. // Number of multipolar expansion terms used for the calculations //
  121. //**********************************************************************************//
  122. int nMie(const unsigned int L, const int pl, vector<double>& x, vector<complex<double> >& m, const unsigned int nTheta, vector<double>& Theta, const int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, vector<complex<double> >& S1, vector<complex<double> >& S2) {
  123. if (x.size() != L || m.size() != L)
  124. throw invalid_argument("Declared number of layers do not fit x and m!");
  125. if (Theta.size() != nTheta)
  126. throw invalid_argument("Declared number of sample for Theta is not correct!");
  127. try {
  128. MultiLayerMie ml_mie;
  129. ml_mie.SetAllLayersSize(x);
  130. ml_mie.SetAllLayersIndex(m);
  131. ml_mie.SetAngles(Theta);
  132. ml_mie.SetPECLayer(pl);
  133. ml_mie.SetMaxTerms(nmax);
  134. ml_mie.RunMieCalculation();
  135. *Qext = ml_mie.GetQext();
  136. *Qsca = ml_mie.GetQsca();
  137. *Qabs = ml_mie.GetQabs();
  138. *Qbk = ml_mie.GetQbk();
  139. *Qpr = ml_mie.GetQpr();
  140. *g = ml_mie.GetAsymmetryFactor();
  141. *Albedo = ml_mie.GetAlbedo();
  142. S1 = ml_mie.GetS1();
  143. S2 = ml_mie.GetS2();
  144. return ml_mie.GetMaxTerms();
  145. } catch(const invalid_argument& ia) {
  146. // Will catch if ml_mie fails or other errors.
  147. cerr << "Invalid argument: " << ia.what() << endl;
  148. throw invalid_argument(ia);
  149. return -1;
  150. }
  151. return 0;
  152. }
  153. //**********************************************************************************//
  154. // This function is just a wrapper to call the full 'nMie' function with fewer //
  155. // parameters, it is here mainly for compatibility with older versions of the //
  156. // program. Also, you can use it if you neither have a PEC layer nor want to define //
  157. // any limit for the maximum number of terms. //
  158. // //
  159. // Input parameters: //
  160. // L: Number of layers //
  161. // x: Array containing the size parameters of the layers [0..L-1] //
  162. // m: Array containing the relative refractive indexes of the layers [0..L-1] //
  163. // nTheta: Number of scattering angles //
  164. // Theta: Array containing all the scattering angles where the scattering //
  165. // amplitudes will be calculated //
  166. // //
  167. // Output parameters: //
  168. // Qext: Efficiency factor for extinction //
  169. // Qsca: Efficiency factor for scattering //
  170. // Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca) //
  171. // Qbk: Efficiency factor for backscattering //
  172. // Qpr: Efficiency factor for the radiation pressure //
  173. // g: Asymmetry factor (g = (Qext-Qpr)/Qsca) //
  174. // Albedo: Single scattering albedo (Albedo = Qsca/Qext) //
  175. // S1, S2: Complex scattering amplitudes //
  176. // //
  177. // Return value: //
  178. // Number of multipolar expansion terms used for the calculations //
  179. //**********************************************************************************//
  180. int nMie(const unsigned int L, vector<double>& x, vector<complex<double> >& m, const unsigned int nTheta, vector<double>& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, vector<complex<double> >& S1, vector<complex<double> >& S2) {
  181. return nmie::nMie(L, -1, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
  182. }
  183. //**********************************************************************************//
  184. // This function is just a wrapper to call the full 'nMie' function with fewer //
  185. // parameters, it is useful if you want to include a PEC layer but not a limit //
  186. // for the maximum number of terms. //
  187. // //
  188. // Input parameters: //
  189. // L: Number of layers //
  190. // pl: Index of PEC layer. If there is none just send -1 //
  191. // x: Array containing the size parameters of the layers [0..L-1] //
  192. // m: Array containing the relative refractive indexes of the layers [0..L-1] //
  193. // nTheta: Number of scattering angles //
  194. // Theta: Array containing all the scattering angles where the scattering //
  195. // amplitudes will be calculated //
  196. // //
  197. // Output parameters: //
  198. // Qext: Efficiency factor for extinction //
  199. // Qsca: Efficiency factor for scattering //
  200. // Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca) //
  201. // Qbk: Efficiency factor for backscattering //
  202. // Qpr: Efficiency factor for the radiation pressure //
  203. // g: Asymmetry factor (g = (Qext-Qpr)/Qsca) //
  204. // Albedo: Single scattering albedo (Albedo = Qsca/Qext) //
  205. // S1, S2: Complex scattering amplitudes //
  206. // //
  207. // Return value: //
  208. // Number of multipolar expansion terms used for the calculations //
  209. //**********************************************************************************//
  210. int nMie(const unsigned int L, const int pl, vector<double>& x, vector<complex<double> >& m, const unsigned int nTheta, vector<double>& Theta, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, vector<complex<double> >& S1, vector<complex<double> >& S2) {
  211. return nmie::nMie(L, pl, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
  212. }
  213. //**********************************************************************************//
  214. // This function is just a wrapper to call the full 'nMie' function with fewer //
  215. // parameters, it is useful if you want to include a limit for the maximum number //
  216. // of terms but not a PEC layer. //
  217. // //
  218. // Input parameters: //
  219. // L: Number of layers //
  220. // x: Array containing the size parameters of the layers [0..L-1] //
  221. // m: Array containing the relative refractive indexes of the layers [0..L-1] //
  222. // nTheta: Number of scattering angles //
  223. // Theta: Array containing all the scattering angles where the scattering //
  224. // amplitudes will be calculated //
  225. // nmax: Maximum number of multipolar expansion terms to be used for the //
  226. // calculations. Only use it if you know what you are doing, otherwise //
  227. // set this parameter to -1 and the function will calculate it //
  228. // //
  229. // Output parameters: //
  230. // Qext: Efficiency factor for extinction //
  231. // Qsca: Efficiency factor for scattering //
  232. // Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca) //
  233. // Qbk: Efficiency factor for backscattering //
  234. // Qpr: Efficiency factor for the radiation pressure //
  235. // g: Asymmetry factor (g = (Qext-Qpr)/Qsca) //
  236. // Albedo: Single scattering albedo (Albedo = Qsca/Qext) //
  237. // S1, S2: Complex scattering amplitudes //
  238. // //
  239. // Return value: //
  240. // Number of multipolar expansion terms used for the calculations //
  241. //**********************************************************************************//
  242. int nMie(const unsigned int L, vector<double>& x, vector<complex<double> >& m, const unsigned int nTheta, vector<double>& Theta, const int nmax, double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo, vector<complex<double> >& S1, vector<complex<double> >& S2) {
  243. return nmie::nMie(L, -1, x, m, nTheta, Theta, nmax, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
  244. }
  245. //**********************************************************************************//
  246. // This function emulates a C call to calculate complex electric and magnetic field //
  247. // in the surroundings and inside (TODO) the particle. //
  248. // //
  249. // Input parameters: //
  250. // L: Number of layers //
  251. // pl: Index of PEC layer. If there is none just send 0 (zero) //
  252. // x: Array containing the size parameters of the layers [0..L-1] //
  253. // m: Array containing the relative refractive indexes of the layers [0..L-1] //
  254. // nmax: Maximum number of multipolar expansion terms to be used for the //
  255. // calculations. Only use it if you know what you are doing, otherwise //
  256. // set this parameter to 0 (zero) and the function will calculate it. //
  257. // ncoord: Number of coordinate points //
  258. // Coords: Array containing all coordinates where the complex electric and //
  259. // magnetic fields will be calculated //
  260. // //
  261. // Output parameters: //
  262. // E, H: Complex electric and magnetic field at the provided coordinates //
  263. // //
  264. // Return value: //
  265. // Number of multipolar expansion terms used for the calculations //
  266. //**********************************************************************************//
  267. int nField(const unsigned int L, const int pl, const vector<double>& x, const vector<complex<double> >& m, const int nmax, const unsigned int ncoord, const vector<double>& Xp_vec, const vector<double>& Yp_vec, const vector<double>& Zp_vec, vector<vector<complex<double> > >& E, vector<vector<complex<double> > >& H) {
  268. if (x.size() != L || m.size() != L)
  269. throw invalid_argument("Declared number of layers do not fit x and m!");
  270. if (Xp_vec.size() != ncoord || Yp_vec.size() != ncoord || Zp_vec.size() != ncoord
  271. || E.size() != ncoord || H.size() != ncoord)
  272. throw invalid_argument("Declared number of coords do not fit Xp, Yp, Zp, E, or H!");
  273. for (auto f:E)
  274. if (f.size() != 3)
  275. throw invalid_argument("Field E is not 3D!");
  276. for (auto f:H)
  277. if (f.size() != 3)
  278. throw invalid_argument("Field H is not 3D!");
  279. try {
  280. MultiLayerMie ml_mie;
  281. //ml_mie.SetPECLayer(pl); // TODO add PEC layer to field plotting
  282. ml_mie.SetAllLayersSize(x);
  283. ml_mie.SetAllLayersIndex(m);
  284. ml_mie.SetFieldCoords({Xp_vec, Yp_vec, Zp_vec});
  285. ml_mie.RunFieldCalculation();
  286. E = ml_mie.GetFieldE();
  287. H = ml_mie.GetFieldH();
  288. return ml_mie.GetMaxTerms();
  289. } catch(const invalid_argument& ia) {
  290. // Will catch if ml_mie fails or other errors.
  291. cerr << "Invalid argument: " << ia.what() << endl;
  292. throw invalid_argument(ia);
  293. return - 1;
  294. }
  295. return 0;
  296. }
  297. // ********************************************************************** //
  298. // Returns previously calculated Qext //
  299. // ********************************************************************** //
  300. double MultiLayerMie::GetQext() {
  301. if (!isMieCalculated_)
  302. throw invalid_argument("You should run calculations before result request!");
  303. return Qext_;
  304. }
  305. // ********************************************************************** //
  306. // Returns previously calculated Qabs //
  307. // ********************************************************************** //
  308. double MultiLayerMie::GetQabs() {
  309. if (!isMieCalculated_)
  310. throw invalid_argument("You should run calculations before result request!");
  311. return Qabs_;
  312. }
  313. // ********************************************************************** //
  314. // Returns previously calculated Qsca //
  315. // ********************************************************************** //
  316. double MultiLayerMie::GetQsca() {
  317. if (!isMieCalculated_)
  318. throw invalid_argument("You should run calculations before result request!");
  319. return Qsca_;
  320. }
  321. // ********************************************************************** //
  322. // Returns previously calculated Qbk //
  323. // ********************************************************************** //
  324. double MultiLayerMie::GetQbk() {
  325. if (!isMieCalculated_)
  326. throw invalid_argument("You should run calculations before result request!");
  327. return Qbk_;
  328. }
  329. // ********************************************************************** //
  330. // Returns previously calculated Qpr //
  331. // ********************************************************************** //
  332. double MultiLayerMie::GetQpr() {
  333. if (!isMieCalculated_)
  334. throw invalid_argument("You should run calculations before result request!");
  335. return Qpr_;
  336. }
  337. // ********************************************************************** //
  338. // Returns previously calculated assymetry factor //
  339. // ********************************************************************** //
  340. double MultiLayerMie::GetAsymmetryFactor() {
  341. if (!isMieCalculated_)
  342. throw invalid_argument("You should run calculations before result request!");
  343. return asymmetry_factor_;
  344. }
  345. // ********************************************************************** //
  346. // Returns previously calculated Albedo //
  347. // ********************************************************************** //
  348. double MultiLayerMie::GetAlbedo() {
  349. if (!isMieCalculated_)
  350. throw invalid_argument("You should run calculations before result request!");
  351. return albedo_;
  352. }
  353. // ********************************************************************** //
  354. // Returns previously calculated S1 //
  355. // ********************************************************************** //
  356. vector<complex<double> > MultiLayerMie::GetS1() {
  357. if (!isMieCalculated_)
  358. throw invalid_argument("You should run calculations before result request!");
  359. return S1_;
  360. }
  361. // ********************************************************************** //
  362. // Returns previously calculated S2 //
  363. // ********************************************************************** //
  364. vector<complex<double> > MultiLayerMie::GetS2() {
  365. if (!isMieCalculated_)
  366. throw invalid_argument("You should run calculations before result request!");
  367. return S2_;
  368. }
  369. // ********************************************************************** //
  370. // Modify scattering (theta) angles //
  371. // ********************************************************************** //
  372. void MultiLayerMie::SetAngles(const vector<double>& angles) {
  373. isExpCoeffsCalc_ = false;
  374. isScaCoeffsCalc_ = false;
  375. isMieCalculated_ = false;
  376. theta_ = angles;
  377. }
  378. // ********************************************************************** //
  379. // Modify size of all layers //
  380. // ********************************************************************** //
  381. void MultiLayerMie::SetAllLayersSize(const vector<double>& layer_size) {
  382. isExpCoeffsCalc_ = false;
  383. isScaCoeffsCalc_ = false;
  384. isMieCalculated_ = false;
  385. size_param_.clear();
  386. double prev_layer_size = 0.0;
  387. for (auto curr_layer_size : layer_size) {
  388. if (curr_layer_size <= 0.0)
  389. throw invalid_argument("Size parameter should be positive!");
  390. if (prev_layer_size > curr_layer_size)
  391. throw invalid_argument
  392. ("Size parameter for next layer should be larger than the previous one!");
  393. prev_layer_size = curr_layer_size;
  394. size_param_.push_back(curr_layer_size);
  395. }
  396. }
  397. // ********************************************************************** //
  398. // Modify refractive index of all layers //
  399. // ********************************************************************** //
  400. void MultiLayerMie::SetAllLayersIndex(const vector< complex<double> >& index) {
  401. isExpCoeffsCalc_ = false;
  402. isScaCoeffsCalc_ = false;
  403. isMieCalculated_ = false;
  404. refractive_index_ = index;
  405. }
  406. // ********************************************************************** //
  407. // Modify coordinates for field calculation //
  408. // ********************************************************************** //
  409. void MultiLayerMie::SetFieldCoords(const vector< vector<double> >& coords) {
  410. if (coords.size() != 3)
  411. throw invalid_argument("Error! Wrong dimension of field monitor points!");
  412. if (coords[0].size() != coords[1].size() || coords[0].size() != coords[2].size())
  413. throw invalid_argument("Error! Missing coordinates for field monitor points!");
  414. coords_ = coords;
  415. }
  416. // ********************************************************************** //
  417. // ********************************************************************** //
  418. // ********************************************************************** //
  419. void MultiLayerMie::SetPECLayer(int layer_position) {
  420. isExpCoeffsCalc_ = false;
  421. isScaCoeffsCalc_ = false;
  422. isMieCalculated_ = false;
  423. if (layer_position < 0 && layer_position != -1)
  424. throw invalid_argument("Error! Layers are numbered from 0!");
  425. PEC_layer_position_ = layer_position;
  426. }
  427. // ********************************************************************** //
  428. // Set maximun number of terms to be used //
  429. // ********************************************************************** //
  430. void MultiLayerMie::SetMaxTerms(int nmax) {
  431. isExpCoeffsCalc_ = false;
  432. isScaCoeffsCalc_ = false;
  433. isMieCalculated_ = false;
  434. nmax_preset_ = nmax;
  435. }
  436. // ********************************************************************** //
  437. // ********************************************************************** //
  438. // ********************************************************************** //
  439. double MultiLayerMie::GetSizeParameter() {
  440. if (size_param_.size() > 0)
  441. return size_param_.back();
  442. else
  443. return 0;
  444. }
  445. // ********************************************************************** //
  446. // Clear layer information //
  447. // ********************************************************************** //
  448. void MultiLayerMie::ClearLayers() {
  449. isExpCoeffsCalc_ = false;
  450. isScaCoeffsCalc_ = false;
  451. isMieCalculated_ = false;
  452. size_param_.clear();
  453. refractive_index_.clear();
  454. }
  455. // ********************************************************************** //
  456. // ********************************************************************** //
  457. // ********************************************************************** //
  458. // Computational core
  459. // ********************************************************************** //
  460. // ********************************************************************** //
  461. // ********************************************************************** //
  462. // ********************************************************************** //
  463. // Calculate calcNstop - equation (17) //
  464. // ********************************************************************** //
  465. void MultiLayerMie::calcNstop() {
  466. const double& xL = size_param_.back();
  467. if (xL <= 8) {
  468. nmax_ = round(xL + 4.0*pow(xL, 1.0/3.0) + 1);
  469. } else if (xL <= 4200) {
  470. nmax_ = round(xL + 4.05*pow(xL, 1.0/3.0) + 2);
  471. } else {
  472. nmax_ = round(xL + 4.0*pow(xL, 1.0/3.0) + 2);
  473. }
  474. }
  475. // ********************************************************************** //
  476. // Maximum number of terms required for the calculation //
  477. // ********************************************************************** //
  478. void MultiLayerMie::calcNmax(unsigned int first_layer) {
  479. int ri, riM1;
  480. const vector<double>& x = size_param_;
  481. const vector<complex<double> >& m = refractive_index_;
  482. calcNstop(); // Set initial nmax_ value
  483. for (unsigned int i = first_layer; i < x.size(); i++) {
  484. if (static_cast<int>(i) > PEC_layer_position_) // static_cast used to avoid warning
  485. ri = round(abs(x[i]*m[i]));
  486. else
  487. ri = 0;
  488. nmax_ = max(nmax_, ri);
  489. // first layer is pec, if pec is present
  490. if ((i > first_layer) && (static_cast<int>(i - 1) > PEC_layer_position_))
  491. riM1 = round(abs(x[i - 1]* m[i]));
  492. else
  493. riM1 = 0;
  494. nmax_ = max(nmax_, riM1);
  495. }
  496. nmax_ += 15; // Final nmax_ value
  497. }
  498. // ********************************************************************** //
  499. // Calculate an - equation (5) //
  500. // ********************************************************************** //
  501. complex<double> MultiLayerMie::calc_an(int n, double XL, complex<double> Ha, complex<double> mL,
  502. complex<double> PsiXL, complex<double> ZetaXL,
  503. complex<double> PsiXLM1, complex<double> ZetaXLM1) {
  504. complex<double> Num = (Ha/mL + n/XL)*PsiXL - PsiXLM1;
  505. complex<double> Denom = (Ha/mL + n/XL)*ZetaXL - ZetaXLM1;
  506. return Num/Denom;
  507. }
  508. // ********************************************************************** //
  509. // Calculate bn - equation (6) //
  510. // ********************************************************************** //
  511. complex<double> MultiLayerMie::calc_bn(int n, double XL, complex<double> Hb, complex<double> mL,
  512. complex<double> PsiXL, complex<double> ZetaXL,
  513. complex<double> PsiXLM1, complex<double> ZetaXLM1) {
  514. complex<double> Num = (mL*Hb + n/XL)*PsiXL - PsiXLM1;
  515. complex<double> Denom = (mL*Hb + n/XL)*ZetaXL - ZetaXLM1;
  516. return Num/Denom;
  517. }
  518. // ********************************************************************** //
  519. // Calculates S1 - equation (25a) //
  520. // ********************************************************************** //
  521. complex<double> MultiLayerMie::calc_S1(int n, complex<double> an, complex<double> bn,
  522. double Pi, double Tau) {
  523. return double(n + n + 1)*(Pi*an + Tau*bn)/double(n*n + n);
  524. }
  525. // ********************************************************************** //
  526. // Calculates S2 - equation (25b) (it's the same as (25a), just switches //
  527. // Pi and Tau) //
  528. // ********************************************************************** //
  529. complex<double> MultiLayerMie::calc_S2(int n, complex<double> an, complex<double> bn,
  530. double Pi, double Tau) {
  531. return calc_S1(n, an, bn, Tau, Pi);
  532. }
  533. //**********************************************************************************//
  534. // This function calculates the logarithmic derivatives of the Riccati-Bessel //
  535. // functions (D1 and D3) for a complex argument (z). //
  536. // Equations (16a), (16b) and (18a) - (18d) //
  537. // //
  538. // Input parameters: //
  539. // z: Complex argument to evaluate D1 and D3 //
  540. // nmax_: Maximum number of terms to calculate D1 and D3 //
  541. // //
  542. // Output parameters: //
  543. // D1, D3: Logarithmic derivatives of the Riccati-Bessel functions //
  544. //**********************************************************************************//
  545. void MultiLayerMie::calcD1D3(const complex<double> z,
  546. vector<complex<double> >& D1,
  547. vector<complex<double> >& D3) {
  548. // Downward recurrence for D1 - equations (16a) and (16b)
  549. D1[nmax_] = complex<double>(0.0, 0.0);
  550. const complex<double> zinv = complex<double>(1.0, 0.0)/z;
  551. for (int n = nmax_; n > 0; n--) {
  552. D1[n - 1] = static_cast<double>(n)*zinv - 1.0/(D1[n] + static_cast<double>(n)*zinv);
  553. }
  554. <<<<<<< HEAD
  555. <<<<<<< HEAD
  556. <<<<<<< HEAD
  557. <<<<<<< HEAD
  558. if (abs(D1[0]) > 1.0e15) {
  559. throw invalid_argument("Unstable D1! Please, try to change input parameters!\n");
  560. //printf("Warning: Potentially unstable D1! Please, try to change input parameters!\n");
  561. }
  562. =======
  563. if (std::abs(D1[0]) > 1.0e8)
  564. // throw std::invalid_argument("Unstable D1! Please, try to change input parameters!\n");
  565. printf("Warning: Potentially unstable D1! Please, try to change input parameters!\n");
  566. >>>>>>> parent of b4f83e4... Throw if D1 seems to be unstable
  567. =======
  568. if (std::abs(D1[0]) > 1.0e8)
  569. // throw std::invalid_argument("Unstable D1! Please, try to change input parameters!\n");
  570. printf("Warning: Potentially unstable D1! Please, try to change input parameters!\n");
  571. >>>>>>> parent of b4f83e4... Throw if D1 seems to be unstable
  572. =======
  573. if (std::abs(D1[0]) > 1.0e8)
  574. // throw std::invalid_argument("Unstable D1! Please, try to change input parameters!\n");
  575. printf("Warning: Potentially unstable D1! Please, try to change input parameters!\n");
  576. >>>>>>> parent of b4f83e4... Throw if D1 seems to be unstable
  577. =======
  578. if (std::abs(D1[0]) > 1.0e8)
  579. // throw std::invalid_argument("Unstable D1! Please, try to change input parameters!\n");
  580. printf("Warning: Potentially unstable D1! Please, try to change input parameters!\n");
  581. >>>>>>> parent of b4f83e4... Throw if D1 seems to be unstable
  582. // Upward recurrence for PsiZeta and D3 - equations (18a) - (18d)
  583. PsiZeta_[0] = 0.5*(1.0 - complex<double>(cos(2.0*z.real()), sin(2.0*z.real()))
  584. *exp(-2.0*z.imag()));
  585. D3[0] = complex<double>(0.0, 1.0);
  586. for (int n = 1; n <= nmax_; n++) {
  587. PsiZeta_[n] = PsiZeta_[n - 1]*(static_cast<double>(n)*zinv - D1[n - 1])
  588. *(static_cast<double>(n)*zinv - D3[n - 1]);
  589. D3[n] = D1[n] + complex<double>(0.0, 1.0)/PsiZeta_[n];
  590. }
  591. }
  592. //**********************************************************************************//
  593. // This function calculates the Riccati-Bessel functions (Psi and Zeta) for a //
  594. // complex argument (z). //
  595. // Equations (20a) - (21b) //
  596. // //
  597. // Input parameters: //
  598. // z: Complex argument to evaluate Psi and Zeta //
  599. // nmax: Maximum number of terms to calculate Psi and Zeta //
  600. // //
  601. // Output parameters: //
  602. // Psi, Zeta: Riccati-Bessel functions //
  603. //**********************************************************************************//
  604. void MultiLayerMie::calcPsiZeta(complex<double> z,
  605. vector<complex<double> >& Psi,
  606. vector<complex<double> >& Zeta) {
  607. complex<double> c_i(0.0, 1.0);
  608. vector<complex<double> > D1(nmax_ + 1), D3(nmax_ + 1);
  609. // First, calculate the logarithmic derivatives
  610. calcD1D3(z, D1, D3);
  611. // Now, use the upward recurrence to calculate Psi and Zeta - equations (20a) - (21b)
  612. Psi[0] = sin(z);
  613. Zeta[0] = sin(z) - c_i*cos(z);
  614. for (int n = 1; n <= nmax_; n++) {
  615. Psi[n] = Psi[n - 1]*(static_cast<double>(n)/z - D1[n - 1]);
  616. Zeta[n] = Zeta[n - 1]*(static_cast<double>(n)/z - D3[n - 1]);
  617. }
  618. }
  619. //**********************************************************************************//
  620. // This function calculates Pi and Tau for a given value of cos(Theta). //
  621. // Equations (26a) - (26c) //
  622. // //
  623. // Input parameters: //
  624. // nmax_: Maximum number of terms to calculate Pi and Tau //
  625. // nTheta: Number of scattering angles //
  626. // Theta: Array containing all the scattering angles where the scattering //
  627. // amplitudes will be calculated //
  628. // //
  629. // Output parameters: //
  630. // Pi, Tau: Angular functions Pi and Tau, as defined in equations (26a) - (26c) //
  631. //**********************************************************************************//
  632. void MultiLayerMie::calcPiTau(const double& costheta,
  633. vector<double>& Pi, vector<double>& Tau) {
  634. int i;
  635. //****************************************************//
  636. // Equations (26a) - (26c) //
  637. //****************************************************//
  638. // Initialize Pi and Tau
  639. Pi[0] = 1.0; // n=1
  640. Tau[0] = costheta;
  641. // Calculate the actual values
  642. if (nmax_ > 1) {
  643. Pi[1] = 3*costheta*Pi[0]; //n=2
  644. Tau[1] = 2*costheta*Pi[1] - 3*Pi[0];
  645. for (i = 2; i < nmax_; i++) { //n=[3..nmax_]
  646. Pi[i] = ((i + i + 1)*costheta*Pi[i - 1] - (i + 1)*Pi[i - 2])/i;
  647. Tau[i] = (i + 1)*costheta*Pi[i] - (i + 2)*Pi[i - 1];
  648. }
  649. }
  650. } // end of MultiLayerMie::calcPiTau(...)
  651. //**********************************************************************************//
  652. // This function calculates vector spherical harmonics (eq. 4.50, p. 95 BH), //
  653. // required to calculate the near-field parameters. //
  654. // //
  655. // Input parameters: //
  656. // Rho: Radial distance //
  657. // Phi: Azimuthal angle //
  658. // Theta: Polar angle //
  659. // rn: Either the spherical Ricatti-Bessel function of first or third kind //
  660. // Dn: Logarithmic derivative of rn //
  661. // Pi, Tau: Angular functions Pi and Tau //
  662. // n: Order of vector spherical harmonics //
  663. // //
  664. // Output parameters: //
  665. // Mo1n, Me1n, No1n, Ne1n: Complex vector spherical harmonics //
  666. //**********************************************************************************//
  667. void MultiLayerMie::calcSpherHarm(const complex<double> Rho, const double Theta, const double Phi,
  668. const complex<double>& rn, const complex<double>& Dn,
  669. const double& Pi, const double& Tau, const double& n,
  670. vector<complex<double> >& Mo1n, vector<complex<double> >& Me1n,
  671. vector<complex<double> >& No1n, vector<complex<double> >& Ne1n) {
  672. // using eq 4.50 in BH
  673. complex<double> c_zero(0.0, 0.0);
  674. Mo1n[0] = c_zero;
  675. Mo1n[1] = cos(Phi)*Pi*rn/Rho;
  676. Mo1n[2] = -sin(Phi)*Tau*rn/Rho;
  677. Me1n[0] = c_zero;
  678. Me1n[1] = -sin(Phi)*Pi*rn/Rho;
  679. Me1n[2] = -cos(Phi)*Tau*rn/Rho;
  680. No1n[0] = sin(Phi)*(n*n + n)*sin(Theta)*Pi*rn/Rho/Rho;
  681. No1n[1] = sin(Phi)*Tau*Dn*rn/Rho;
  682. No1n[2] = cos(Phi)*Pi*Dn*rn/Rho;
  683. Ne1n[0] = cos(Phi)*(n*n + n)*sin(Theta)*Pi*rn/Rho/Rho;
  684. Ne1n[1] = cos(Phi)*Tau*Dn*rn/Rho;
  685. Ne1n[2] = -sin(Phi)*Pi*Dn*rn/Rho;
  686. } // end of MultiLayerMie::calcSpherHarm(...)
  687. //**********************************************************************************//
  688. // This function calculates the scattering coefficients required to calculate //
  689. // both the near- and far-field parameters. //
  690. // //
  691. // Input parameters: //
  692. // L: Number of layers //
  693. // pl: Index of PEC layer. If there is none just send -1 //
  694. // x: Array containing the size parameters of the layers [0..L-1] //
  695. // m: Array containing the relative refractive indexes of the layers [0..L-1] //
  696. // nmax: Maximum number of multipolar expansion terms to be used for the //
  697. // calculations. Only use it if you know what you are doing, otherwise //
  698. // set this parameter to -1 and the function will calculate it. //
  699. // //
  700. // Output parameters: //
  701. // an, bn: Complex scattering amplitudes //
  702. // //
  703. // Return value: //
  704. // Number of multipolar expansion terms used for the calculations //
  705. //**********************************************************************************//
  706. void MultiLayerMie::calcScattCoeffs() {
  707. isScaCoeffsCalc_ = false;
  708. const vector<double>& x = size_param_;
  709. const vector<complex<double> >& m = refractive_index_;
  710. const int& pl = PEC_layer_position_;
  711. const int L = refractive_index_.size();
  712. //************************************************************************//
  713. // Calculate the index of the first layer. It can be either 0 (default) //
  714. // or the index of the outermost PEC layer. In the latter case all layers //
  715. // below the PEC are discarded. //
  716. // ***********************************************************************//
  717. int fl = (pl > 0) ? pl : 0;
  718. if (nmax_preset_ <= 0) calcNmax(fl);
  719. else nmax_ = nmax_preset_;
  720. complex<double> z1, z2;
  721. //**************************************************************************//
  722. // Note that since Fri, Nov 14, 2014 all arrays start from 0 (zero), which //
  723. // means that index = layer number - 1 or index = n - 1. The only exception //
  724. // are the arrays for representing D1, D3 and Q because they need a value //
  725. // for the index 0 (zero), hence it is important to consider this shift //
  726. // between different arrays. The change was done to optimize memory usage. //
  727. //**************************************************************************//
  728. // Allocate memory to the arrays
  729. vector<complex<double> > D1_mlxl(nmax_ + 1), D1_mlxlM1(nmax_ + 1),
  730. D3_mlxl(nmax_ + 1), D3_mlxlM1(nmax_ + 1);
  731. vector<vector<complex<double> > > Q(L), Ha(L), Hb(L);
  732. for (int l = 0; l < L; l++) {
  733. Q[l].resize(nmax_ + 1);
  734. Ha[l].resize(nmax_);
  735. Hb[l].resize(nmax_);
  736. }
  737. an_.resize(nmax_);
  738. bn_.resize(nmax_);
  739. PsiZeta_.resize(nmax_ + 1);
  740. vector<complex<double> > PsiXL(nmax_ + 1), ZetaXL(nmax_ + 1);
  741. //*************************************************//
  742. // Calculate D1 and D3 for z1 in the first layer //
  743. //*************************************************//
  744. if (fl == pl) { // PEC layer
  745. for (int n = 0; n <= nmax_; n++) {
  746. D1_mlxl[n] = complex<double>(0.0, - 1.0);
  747. D3_mlxl[n] = complex<double>(0.0, 1.0);
  748. }
  749. } else { // Regular layer
  750. z1 = x[fl]* m[fl];
  751. // Calculate D1 and D3
  752. calcD1D3(z1, D1_mlxl, D3_mlxl);
  753. }
  754. //******************************************************************//
  755. // Calculate Ha and Hb in the first layer - equations (7a) and (8a) //
  756. //******************************************************************//
  757. for (int n = 0; n < nmax_; n++) {
  758. Ha[fl][n] = D1_mlxl[n + 1];
  759. Hb[fl][n] = D1_mlxl[n + 1];
  760. }
  761. //*****************************************************//
  762. // Iteration from the second layer to the last one (L) //
  763. //*****************************************************//
  764. complex<double> Temp, Num, Denom;
  765. complex<double> G1, G2;
  766. for (int l = fl + 1; l < L; l++) {
  767. //************************************************************//
  768. //Calculate D1 and D3 for z1 and z2 in the layers fl + 1..L //
  769. //************************************************************//
  770. z1 = x[l]*m[l];
  771. z2 = x[l - 1]*m[l];
  772. //Calculate D1 and D3 for z1
  773. calcD1D3(z1, D1_mlxl, D3_mlxl);
  774. //Calculate D1 and D3 for z2
  775. calcD1D3(z2, D1_mlxlM1, D3_mlxlM1);
  776. //*************************************************//
  777. //Calculate Q, Ha and Hb in the layers fl + 1..L //
  778. //*************************************************//
  779. // Upward recurrence for Q - equations (19a) and (19b)
  780. Num = exp(-2.0*(z1.imag() - z2.imag()))
  781. *complex<double>(cos(-2.0*z2.real()) - exp(-2.0*z2.imag()), sin(-2.0*z2.real()));
  782. Denom = complex<double>(cos(-2.0*z1.real()) - exp(-2.0*z1.imag()), sin(-2.0*z1.real()));
  783. Q[l][0] = Num/Denom;
  784. for (int n = 1; n <= nmax_; n++) {
  785. Num = (z1*D1_mlxl[n] + double(n))*(double(n) - z1*D3_mlxl[n - 1]);
  786. Denom = (z2*D1_mlxlM1[n] + double(n))*(double(n) - z2*D3_mlxlM1[n - 1]);
  787. Q[l][n] = ((pow2(x[l - 1]/x[l])* Q[l][n - 1])*Num)/Denom;
  788. }
  789. // Upward recurrence for Ha and Hb - equations (7b), (8b) and (12) - (15)
  790. for (int n = 1; n <= nmax_; n++) {
  791. //Ha
  792. if ((l - 1) == pl) { // The layer below the current one is a PEC layer
  793. G1 = -D1_mlxlM1[n];
  794. G2 = -D3_mlxlM1[n];
  795. } else {
  796. G1 = (m[l]*Ha[l - 1][n - 1]) - (m[l - 1]*D1_mlxlM1[n]);
  797. G2 = (m[l]*Ha[l - 1][n - 1]) - (m[l - 1]*D3_mlxlM1[n]);
  798. } // end of if PEC
  799. Temp = Q[l][n]*G1;
  800. Num = (G2*D1_mlxl[n]) - (Temp*D3_mlxl[n]);
  801. Denom = G2 - Temp;
  802. Ha[l][n - 1] = Num/Denom;
  803. //Hb
  804. if ((l - 1) == pl) { // The layer below the current one is a PEC layer
  805. G1 = Hb[l - 1][n - 1];
  806. G2 = Hb[l - 1][n - 1];
  807. } else {
  808. G1 = (m[l - 1]*Hb[l - 1][n - 1]) - (m[l]*D1_mlxlM1[n]);
  809. G2 = (m[l - 1]*Hb[l - 1][n - 1]) - (m[l]*D3_mlxlM1[n]);
  810. } // end of if PEC
  811. Temp = Q[l][n]*G1;
  812. Num = (G2*D1_mlxl[n]) - (Temp* D3_mlxl[n]);
  813. Denom = (G2- Temp);
  814. Hb[l][n - 1] = (Num/ Denom);
  815. } // end of for Ha and Hb terms
  816. } // end of for layers iteration
  817. //**************************************//
  818. //Calculate Psi and Zeta for XL //
  819. //**************************************//
  820. // Calculate PsiXL and ZetaXL
  821. calcPsiZeta(x[L - 1], PsiXL, ZetaXL);
  822. //*********************************************************************//
  823. // Finally, we calculate the scattering coefficients (an and bn) and //
  824. // the angular functions (Pi and Tau). Note that for these arrays the //
  825. // first layer is 0 (zero), in future versions all arrays will follow //
  826. // this convention to save memory. (13 Nov, 2014) //
  827. //*********************************************************************//
  828. for (int n = 0; n < nmax_; n++) {
  829. //********************************************************************//
  830. //Expressions for calculating an and bn coefficients are not valid if //
  831. //there is only one PEC layer (ie, for a simple PEC sphere). //
  832. //********************************************************************//
  833. if (pl < (L - 1)) {
  834. an_[n] = calc_an(n + 1, x[L - 1], Ha[L - 1][n], m[L - 1], PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
  835. bn_[n] = calc_bn(n + 1, x[L - 1], Hb[L - 1][n], m[L - 1], PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
  836. } else {
  837. an_[n] = calc_an(n + 1, x[L - 1], complex<double>(0.0, 0.0), complex<double>(1.0, 0.0), PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
  838. bn_[n] = PsiXL[n + 1]/ZetaXL[n + 1];
  839. }
  840. } // end of for an and bn terms
  841. isScaCoeffsCalc_ = true;
  842. } // end of MultiLayerMie::calcScattCoeffs()
  843. //**********************************************************************************//
  844. // This function calculates the actual scattering parameters and amplitudes //
  845. // //
  846. // Input parameters: //
  847. // L: Number of layers //
  848. // pl: Index of PEC layer. If there is none just send -1 //
  849. // x: Array containing the size parameters of the layers [0..L-1] //
  850. // m: Array containing the relative refractive indexes of the layers [0..L-1] //
  851. // nTheta: Number of scattering angles //
  852. // Theta: Array containing all the scattering angles where the scattering //
  853. // amplitudes will be calculated //
  854. // nmax_: Maximum number of multipolar expansion terms to be used for the //
  855. // calculations. Only use it if you know what you are doing, otherwise //
  856. // set this parameter to -1 and the function will calculate it //
  857. // //
  858. // Output parameters: //
  859. // Qext: Efficiency factor for extinction //
  860. // Qsca: Efficiency factor for scattering //
  861. // Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca) //
  862. // Qbk: Efficiency factor for backscattering //
  863. // Qpr: Efficiency factor for the radiation pressure //
  864. // g: Asymmetry factor (g = (Qext-Qpr)/Qsca) //
  865. // Albedo: Single scattering albedo (Albedo = Qsca/Qext) //
  866. // S1, S2: Complex scattering amplitudes //
  867. // //
  868. // Return value: //
  869. // Number of multipolar expansion terms used for the calculations //
  870. //**********************************************************************************//
  871. void MultiLayerMie::RunMieCalculation() {
  872. if (size_param_.size() != refractive_index_.size())
  873. throw invalid_argument("Each size parameter should have only one index!");
  874. if (size_param_.size() == 0)
  875. throw invalid_argument("Initialize model first!");
  876. const vector<double>& x = size_param_;
  877. isExpCoeffsCalc_ = false;
  878. isScaCoeffsCalc_ = false;
  879. isMieCalculated_ = false;
  880. // Calculate scattering coefficients
  881. calcScattCoeffs();
  882. if (!isScaCoeffsCalc_) // TODO seems to be unreachable
  883. throw invalid_argument("Calculation of scattering coefficients failed!");
  884. // Initialize the scattering parameters
  885. Qext_ = 0.0;
  886. Qsca_ = 0.0;
  887. Qabs_ = 0.0;
  888. Qbk_ = 0.0;
  889. Qpr_ = 0.0;
  890. asymmetry_factor_ = 0.0;
  891. albedo_ = 0.0;
  892. // Initialize the scattering amplitudes
  893. vector<complex<double> > tmp1(theta_.size(),complex<double>(0.0, 0.0));
  894. S1_.swap(tmp1);
  895. S2_ = S1_;
  896. vector<double> Pi(nmax_), Tau(nmax_);
  897. complex<double> Qbktmp(0.0, 0.0);
  898. vector< complex<double> > Qbktmp_ch(nmax_ - 1, Qbktmp);
  899. // By using downward recurrence we avoid loss of precision due to float rounding errors
  900. // See: https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html
  901. // http://en.wikipedia.org/wiki/Loss_of_significance
  902. for (int i = nmax_ - 2; i >= 0; i--) {
  903. const int n = i + 1;
  904. // Equation (27)
  905. Qext_ += (n + n + 1.0)*(an_[i].real() + bn_[i].real());
  906. // Equation (28)
  907. Qsca_ += (n + n + 1.0)*(an_[i].real()*an_[i].real() + an_[i].imag()*an_[i].imag()
  908. + bn_[i].real()*bn_[i].real() + bn_[i].imag()*bn_[i].imag());
  909. // Equation (29)
  910. Qpr_ += ((n*(n + 2)/(n + 1))*((an_[i]*conj(an_[n]) + bn_[i]*conj(bn_[n])).real())
  911. + ((double)(n + n + 1)/(n*(n + 1)))*(an_[i]*conj(bn_[i])).real());
  912. // Equation (33)
  913. Qbktmp += (double)(n + n + 1)*(1 - 2*(n % 2))*(an_[i]- bn_[i]);
  914. // Calculate the scattering amplitudes (S1 and S2) //
  915. // Equations (25a) - (25b) //
  916. for (unsigned int t = 0; t < theta_.size(); t++) {
  917. calcPiTau(cos(theta_[t]), Pi, Tau);
  918. S1_[t] += calc_S1(n, an_[i], bn_[i], Pi[i], Tau[i]);
  919. S2_[t] += calc_S2(n, an_[i], bn_[i], Pi[i], Tau[i]);
  920. }
  921. }
  922. double x2 = pow2(x.back());
  923. Qext_ = 2.0*(Qext_)/x2; // Equation (27)
  924. Qsca_ = 2.0*(Qsca_)/x2; // Equation (28)
  925. Qpr_ = Qext_ - 4.0*(Qpr_)/x2; // Equation (29)
  926. Qabs_ = Qext_ - Qsca_; // Equation (30)
  927. albedo_ = Qsca_/Qext_; // Equation (31)
  928. asymmetry_factor_ = (Qext_ - Qpr_)/Qsca_; // Equation (32)
  929. Qbk_ = (Qbktmp.real()*Qbktmp.real() + Qbktmp.imag()*Qbktmp.imag())/x2; // Equation (33)
  930. isMieCalculated_ = true;
  931. }
  932. //**********************************************************************************//
  933. // This function calculates the expansion coefficients inside the particle, //
  934. // required to calculate the near-field parameters. //
  935. // //
  936. // Input parameters: //
  937. // L: Number of layers //
  938. // pl: Index of PEC layer. If there is none just send -1 //
  939. // x: Array containing the size parameters of the layers [0..L-1] //
  940. // m: Array containing the relative refractive indexes of the layers [0..L-1] //
  941. // nmax: Maximum number of multipolar expansion terms to be used for the //
  942. // calculations. Only use it if you know what you are doing, otherwise //
  943. // set this parameter to -1 and the function will calculate it. //
  944. // //
  945. // Output parameters: //
  946. // aln, bln, cln, dln: Complex scattering amplitudes inside the particle //
  947. // //
  948. // Return value: //
  949. // Number of multipolar expansion terms used for the calculations //
  950. //**********************************************************************************//
  951. void MultiLayerMie::calcExpanCoeffs() {
  952. if (!isScaCoeffsCalc_)
  953. throw invalid_argument("(ExpanCoeffs) You should calculate external coefficients first!");
  954. isExpCoeffsCalc_ = false;
  955. complex<double> c_one(1.0, 0.0), c_zero(0.0, 0.0);
  956. const int L = refractive_index_.size();
  957. aln_.resize(L + 1);
  958. bln_.resize(L + 1);
  959. cln_.resize(L + 1);
  960. dln_.resize(L + 1);
  961. for (int l = 0; l <= L; l++) {
  962. aln_[l].resize(nmax_);
  963. bln_[l].resize(nmax_);
  964. cln_[l].resize(nmax_);
  965. dln_[l].resize(nmax_);
  966. }
  967. // Yang, paragraph under eq. A3
  968. // a^(L + 1)_n = a_n, d^(L + 1) = 1 ...
  969. for (int n = 0; n < nmax_; n++) {
  970. aln_[L][n] = an_[n];
  971. bln_[L][n] = bn_[n];
  972. cln_[L][n] = c_one;
  973. dln_[L][n] = c_one;
  974. }
  975. vector<complex<double> > D1z(nmax_ + 1), D1z1(nmax_ + 1), D3z(nmax_ + 1), D3z1(nmax_ + 1);
  976. vector<complex<double> > Psiz(nmax_ + 1), Psiz1(nmax_ + 1), Zetaz(nmax_ + 1), Zetaz1(nmax_ + 1);
  977. complex<double> denomZeta, denomPsi, T1, T2, T3, T4;
  978. auto& m = refractive_index_;
  979. vector< complex<double> > m1(L);
  980. for (int l = 0; l < L - 1; l++) m1[l] = m[l + 1];
  981. m1[L - 1] = complex<double> (1.0, 0.0);
  982. complex<double> z, z1;
  983. for (int l = L - 1; l >= 0; l--) {
  984. z = size_param_[l]*m[l];
  985. z1 = size_param_[l]*m1[l];
  986. calcD1D3(z, D1z, D3z);
  987. calcD1D3(z1, D1z1, D3z1);
  988. calcPsiZeta(z, Psiz, Zetaz);
  989. calcPsiZeta(z1, Psiz1, Zetaz1);
  990. for (int n = 0; n < nmax_; n++) {
  991. int n1 = n + 1;
  992. denomZeta = Zetaz[n1]*(D1z[n1] - D3z[n1]);
  993. denomPsi = Psiz[n1]*(D1z[n1] - D3z[n1]);
  994. T1 = aln_[l + 1][n]*Zetaz1[n1] - dln_[l + 1][n]*Psiz1[n1];
  995. T2 = (bln_[l + 1][n]*Zetaz1[n1] - cln_[l + 1][n]*Psiz1[n1])*m[l]/m1[l];
  996. T3 = (dln_[l + 1][n]*D1z1[n1]*Psiz1[n1] - aln_[l + 1][n]*D3z1[n1]*Zetaz1[n1])*m[l]/m1[l];
  997. T4 = cln_[l + 1][n]*D1z1[n1]*Psiz1[n1] - bln_[l + 1][n]*D3z1[n1]*Zetaz1[n1];
  998. // aln
  999. aln_[l][n] = (D1z[n1]*T1 + T3)/denomZeta;
  1000. // bln
  1001. bln_[l][n] = (D1z[n1]*T2 + T4)/denomZeta;
  1002. // cln
  1003. cln_[l][n] = (D3z[n1]*T2 + T4)/denomPsi;
  1004. // dln
  1005. dln_[l][n] = (D3z[n1]*T1 + T3)/denomPsi;
  1006. } // end of all n
  1007. } // end of all l
  1008. // Check the result and change aln_[0][n] and aln_[0][n] for exact zero
  1009. for (int n = 0; n < nmax_; ++n) {
  1010. if (abs(aln_[0][n]) < 1e-10) aln_[0][n] = 0.0;
  1011. else {
  1012. //throw invalid_argument("Unstable calculation of aln_[0][n]!");
  1013. printf("Warning: Potentially unstable calculation of aln (aln[0][%i] = %g, %gi)\n", n, aln_[0][n].real(), aln_[0][n].imag());
  1014. aln_[0][n] = 0.0;
  1015. }
  1016. if (abs(bln_[0][n]) < 1e-10) bln_[0][n] = 0.0;
  1017. else {
  1018. //throw invalid_argument("Unstable calculation of bln_[0][n]!");
  1019. printf("Warning: Potentially unstable calculation of bln (bln[0][%i] = %g, %gi)\n", n, bln_[0][n].real(), bln_[0][n].imag());
  1020. bln_[0][n] = 0.0;
  1021. }
  1022. }
  1023. isExpCoeffsCalc_ = true;
  1024. } // end of void MultiLayerMie::calcExpanCoeffs()
  1025. //**********************************************************************************//
  1026. // This function calculates the electric (E) and magnetic (H) fields inside and //
  1027. // around the particle. //
  1028. // //
  1029. // Input parameters (coordinates of the point): //
  1030. // Rho: Radial distance //
  1031. // Phi: Azimuthal angle //
  1032. // Theta: Polar angle //
  1033. // //
  1034. // Output parameters: //
  1035. // E, H: Complex electric and magnetic fields //
  1036. //**********************************************************************************//
  1037. void MultiLayerMie::calcField(const double Rho, const double Theta, const double Phi,
  1038. vector<complex<double> >& E, vector<complex<double> >& H) {
  1039. complex<double> c_zero(0.0, 0.0), c_i(0.0, 1.0), c_one(1.0, 0.0);
  1040. vector<complex<double> > ipow = {c_one, c_i, -c_one, -c_i}; // Vector containing precomputed integer powers of i to avoid computation
  1041. vector<complex<double> > M3o1n(3), M3e1n(3), N3o1n(3), N3e1n(3);
  1042. vector<complex<double> > M1o1n(3), M1e1n(3), N1o1n(3), N1e1n(3);
  1043. vector<complex<double> > Psi(nmax_ + 1), D1n(nmax_ + 1), Zeta(nmax_ + 1), D3n(nmax_ + 1);
  1044. vector<double> Pi(nmax_), Tau(nmax_);
  1045. int l = 0; // Layer number
  1046. complex<double> ml;
  1047. // Initialize E and H
  1048. for (int i = 0; i < 3; i++) {
  1049. E[i] = c_zero;
  1050. H[i] = c_zero;
  1051. }
  1052. if (Rho > size_param_.back()) {
  1053. l = size_param_.size();
  1054. ml = c_one;
  1055. } else {
  1056. for (int i = size_param_.size() - 1; i >= 0 ; i--) {
  1057. if (Rho <= size_param_[i]) {
  1058. l = i;
  1059. }
  1060. }
  1061. ml = refractive_index_[l];
  1062. }
  1063. // Calculate logarithmic derivative of the Ricatti-Bessel functions
  1064. calcD1D3(Rho*ml, D1n, D3n);
  1065. // Calculate Ricatti-Bessel functions
  1066. calcPsiZeta(Rho*ml, Psi, Zeta);
  1067. // Calculate angular functions Pi and Tau
  1068. calcPiTau(cos(Theta), Pi, Tau);
  1069. for (int n = nmax_ - 2; n >= 0; n--) {
  1070. int n1 = n + 1;
  1071. double rn = static_cast<double>(n1);
  1072. // using BH 4.12 and 4.50
  1073. calcSpherHarm(Rho*ml, Theta, Phi, Psi[n1], D1n[n1], Pi[n], Tau[n], rn, M1o1n, M1e1n, N1o1n, N1e1n);
  1074. calcSpherHarm(Rho*ml, Theta, Phi, Zeta[n1], D3n[n1], Pi[n], Tau[n], rn, M3o1n, M3e1n, N3o1n, N3e1n);
  1075. // Total field in the lth layer: eqs. (1) and (2) in Yang, Appl. Opt., 42 (2003) 1710-1720
  1076. complex<double> En = ipow[n1 % 4]*(rn + rn + 1.0)/(rn*rn + rn);
  1077. for (int i = 0; i < 3; i++) {
  1078. // electric field E [V m - 1] = EF*E0
  1079. E[i] += En*(cln_[l][n]*M1o1n[i] - c_i*dln_[l][n]*N1e1n[i]
  1080. + c_i*aln_[l][n]*N3e1n[i] - bln_[l][n]*M3o1n[i]);
  1081. H[i] += En*(-dln_[l][n]*M1e1n[i] - c_i*cln_[l][n]*N1o1n[i]
  1082. + c_i*bln_[l][n]*N3o1n[i] + aln_[l][n]*M3e1n[i]);
  1083. }
  1084. } // end of for all n
  1085. // magnetic field
  1086. complex<double> hffact = ml/(cc_*mu_);
  1087. for (int i = 0; i < 3; i++) {
  1088. H[i] = hffact*H[i];
  1089. }
  1090. } // end of MultiLayerMie::calcField(...)
  1091. //**********************************************************************************//
  1092. // This function calculates complex electric and magnetic field in the surroundings //
  1093. // and inside the particle. //
  1094. // //
  1095. // Input parameters: //
  1096. // L: Number of layers //
  1097. // pl: Index of PEC layer. If there is none just send 0 (zero) //
  1098. // x: Array containing the size parameters of the layers [0..L-1] //
  1099. // m: Array containing the relative refractive indexes of the layers [0..L-1] //
  1100. // nmax: Maximum number of multipolar expansion terms to be used for the //
  1101. // calculations. Only use it if you know what you are doing, otherwise //
  1102. // set this parameter to 0 (zero) and the function will calculate it. //
  1103. // ncoord: Number of coordinate points //
  1104. // Coords: Array containing all coordinates where the complex electric and //
  1105. // magnetic fields will be calculated //
  1106. // //
  1107. // Output parameters: //
  1108. // E, H: Complex electric and magnetic field at the provided coordinates //
  1109. // //
  1110. // Return value: //
  1111. // Number of multipolar expansion terms used for the calculations //
  1112. //**********************************************************************************//
  1113. void MultiLayerMie::RunFieldCalculation() {
  1114. double Rho, Theta, Phi;
  1115. // Calculate scattering coefficients an_ and bn_
  1116. calcScattCoeffs();
  1117. // Calculate expansion coefficients aln_, bln_, cln_, and dln_
  1118. calcExpanCoeffs();
  1119. long total_points = coords_[0].size();
  1120. E_.resize(total_points);
  1121. H_.resize(total_points);
  1122. for (auto& f : E_) f.resize(3);
  1123. for (auto& f : H_) f.resize(3);
  1124. for (int point = 0; point < total_points; point++) {
  1125. const double& Xp = coords_[0][point];
  1126. const double& Yp = coords_[1][point];
  1127. const double& Zp = coords_[2][point];
  1128. // Convert to spherical coordinates
  1129. Rho = sqrt(pow2(Xp) + pow2(Yp) + pow2(Zp));
  1130. // If Rho=0 then Theta is undefined. Just set it to zero to avoid problems
  1131. Theta = (Rho > 0.0) ? acos(Zp/Rho) : 0.0;
  1132. // If Xp=Yp=0 then Phi is undefined. Just set it to zero to avoid problems
  1133. if (Xp == 0.0)
  1134. Phi = (Yp != 0.0) ? asin(Yp/sqrt(pow2(Xp) + pow2(Yp))) : 0.0;
  1135. else
  1136. Phi = acos(Xp/sqrt(pow2(Xp) + pow2(Yp)));
  1137. // Avoid convergence problems due to Rho too small
  1138. if (Rho < 1e-5) Rho = 1e-5;
  1139. //*******************************************************//
  1140. // external scattering field = incident + scattered //
  1141. // BH p.92 (4.37), 94 (4.45), 95 (4.50) //
  1142. // assume: medium is non-absorbing; refim = 0; Uabs = 0 //
  1143. //*******************************************************//
  1144. // This array contains the fields in spherical coordinates
  1145. vector<complex<double> > Es(3), Hs(3);
  1146. // Do the actual calculation of electric and magnetic field
  1147. calcField(Rho, Theta, Phi, Es, Hs);
  1148. { //Now, convert the fields back to cartesian coordinates
  1149. E_[point][0] = sin(Theta)*cos(Phi)*Es[0] + cos(Theta)*cos(Phi)*Es[1] - sin(Phi)*Es[2];
  1150. E_[point][1] = sin(Theta)*sin(Phi)*Es[0] + cos(Theta)*sin(Phi)*Es[1] + cos(Phi)*Es[2];
  1151. E_[point][2] = cos(Theta)*Es[0] - sin(Theta)*Es[1];
  1152. H_[point][0] = sin(Theta)*cos(Phi)*Hs[0] + cos(Theta)*cos(Phi)*Hs[1] - sin(Phi)*Hs[2];
  1153. H_[point][1] = sin(Theta)*sin(Phi)*Hs[0] + cos(Theta)*sin(Phi)*Hs[1] + cos(Phi)*Hs[2];
  1154. H_[point][2] = cos(Theta)*Hs[0] - sin(Theta)*Hs[1];
  1155. }
  1156. } // end of for all field coordinates
  1157. } // end of MultiLayerMie::RunFieldCalculation()
  1158. } // end of namespace nmie