nmie.cc 46 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981
  1. //**********************************************************************************//
  2. // Copyright (C) 2009-2015 Ovidio Pena <ovidio@bytesfall.com> //
  3. // //
  4. // This file is part of scattnlay //
  5. // //
  6. // This program is free software: you can redistribute it and/or modify //
  7. // it under the terms of the GNU General Public License as published by //
  8. // the Free Software Foundation, either version 3 of the License, or //
  9. // (at your option) any later version. //
  10. // //
  11. // This program is distributed in the hope that it will be useful, //
  12. // but WITHOUT ANY WARRANTY; without even the implied warranty of //
  13. // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
  14. // GNU General Public License for more details. //
  15. // //
  16. // The only additional remark is that we expect that all publications //
  17. // describing work using this software, or all commercial products //
  18. // using it, cite the following reference: //
  19. // [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by //
  20. // a multilayered sphere," Computer Physics Communications, //
  21. // vol. 180, Nov. 2009, pp. 2348-2354. //
  22. // //
  23. // You should have received a copy of the GNU General Public License //
  24. // along with this program. If not, see <http://www.gnu.org/licenses/>. //
  25. //**********************************************************************************//
  26. //**********************************************************************************//
  27. // This library implements the algorithm for a multilayered sphere described by: //
  28. // [1] W. Yang, "Improved recursive algorithm for light scattering by a //
  29. // multilayered sphere,” Applied Optics, vol. 42, Mar. 2003, pp. 1710-1720. //
  30. // //
  31. // You can find the description of all the used equations in: //
  32. // [2] O. Pena and U. Pal, "Scattering of electromagnetic radiation by //
  33. // a multilayered sphere," Computer Physics Communications, //
  34. // vol. 180, Nov. 2009, pp. 2348-2354. //
  35. // //
  36. // Hereinafter all equations numbers refer to [2] //
  37. //**********************************************************************************//
  38. #include <math.h>
  39. #include <stdlib.h>
  40. #include <stdio.h>
  41. #include "nmie.h"
  42. #define round(x) ((x) >= 0 ? (int)((x) + 0.5):(int)((x) - 0.5))
  43. // Calculate Nstop - equation (17)
  44. int Nstop(double xL) {
  45. int result;
  46. if (xL <= 8) {
  47. result = round(xL + 4*pow(xL, 1/3) + 1);
  48. } else if (xL <= 4200) {
  49. result = round(xL + 4.05*pow(xL, 1/3) + 2);
  50. } else {
  51. result = round(xL + 4*pow(xL, 1/3) + 2);
  52. }
  53. return result;
  54. }
  55. //**********************************************************************************//
  56. int Nmax(int L, int fl, int pl,
  57. std::vector<double> x,
  58. std::vector<std::complex<double> > m) {
  59. int i, result, ri, riM1;
  60. result = Nstop(x[L - 1]);
  61. for (i = fl; i < L; i++) {
  62. if (i > pl) {
  63. ri = round(std::abs(x[i]*m[i]));
  64. } else {
  65. ri = 0;
  66. }
  67. if (result < ri) {
  68. result = ri;
  69. }
  70. if ((i > fl) && ((i - 1) > pl)) {
  71. riM1 = round(std::abs(x[i - 1]* m[i]));
  72. } else {
  73. riM1 = 0;
  74. }
  75. if (result < riM1) {
  76. result = riM1;
  77. }
  78. }
  79. return result + 15;
  80. }
  81. //**********************************************************************************//
  82. // This function calculates the spherical Bessel (jn) and Hankel (h1n) functions //
  83. // and their derivatives for a given complex value z. See pag. 87 B&H. //
  84. // //
  85. // Input parameters: //
  86. // z: Real argument to evaluate jn and h1n //
  87. // nmax: Maximum number of terms to calculate jn and h1n //
  88. // //
  89. // Output parameters: //
  90. // jn, h1n: Spherical Bessel and Hankel functions //
  91. // jnp, h1np: Derivatives of the spherical Bessel and Hankel functions //
  92. // //
  93. // The implementation follows the algorithm by I.J. Thompson and A.R. Barnett, //
  94. // Comp. Phys. Comm. 47 (1987) 245-257. //
  95. // //
  96. // Complex spherical Bessel functions from n=0..nmax-1 for z in the upper half //
  97. // plane (Im(z) > -3). //
  98. // //
  99. // j[n] = j/n(z) Regular solution: j[0]=sin(z)/z //
  100. // j'[n] = d[j/n(z)]/dz //
  101. // h1[n] = h[0]/n(z) Irregular Hankel function: //
  102. // h1'[n] = d[h[0]/n(z)]/dz h1[0] = j0(z) + i*y0(z) //
  103. // = (sin(z)-i*cos(z))/z //
  104. // = -i*exp(i*z)/z //
  105. // Using complex CF1, and trigonometric forms for n=0 solutions. //
  106. //**********************************************************************************//
  107. int sbesjh(std::complex<double> z, int nmax, std::vector<std::complex<double> >& jn, std::vector<std::complex<double> >& jnp, std::vector<std::complex<double> >& h1n, std::vector<std::complex<double> >& h1np) {
  108. const int limit = 20000;
  109. double const accur = 1.0e-12;
  110. double const tm30 = 1e-30;
  111. int n;
  112. double absc;
  113. std::complex<double> zi, w;
  114. std::complex<double> pl, f, b, d, c, del, jn0, xdb, h1nldb, h1nbdb;
  115. absc = std::abs(std::real(z)) + std::abs(std::imag(z));
  116. if ((absc < accur) || (std::imag(z) < -3.0)) {
  117. return -1;
  118. }
  119. zi = 1.0/z;
  120. w = zi + zi;
  121. pl = double(nmax)*zi;
  122. f = pl + zi;
  123. b = f + f + zi;
  124. d = 0.0;
  125. c = f;
  126. for (n = 0; n < limit; n++) {
  127. d = b - d;
  128. c = b - 1.0/c;
  129. absc = std::abs(std::real(d)) + std::abs(std::imag(d));
  130. if (absc < tm30) {
  131. d = tm30;
  132. }
  133. absc = std::abs(std::real(c)) + std::abs(std::imag(c));
  134. if (absc < tm30) {
  135. c = tm30;
  136. }
  137. d = 1.0/d;
  138. del = d*c;
  139. f = f*del;
  140. b += w;
  141. absc = std::abs(std::real(del - 1.0)) + std::abs(std::imag(del - 1.0));
  142. if (absc < accur) {
  143. // We have obtained the desired accuracy
  144. break;
  145. }
  146. }
  147. if (absc > accur) {
  148. // We were not able to obtain the desired accuracy
  149. return -2;
  150. }
  151. jn[nmax - 1] = tm30;
  152. jnp[nmax - 1] = f*jn[nmax - 1];
  153. // Downward recursion to n=0 (N.B. Coulomb Functions)
  154. for (n = nmax - 2; n >= 0; n--) {
  155. jn[n] = pl*jn[n + 1] + jnp[n + 1];
  156. jnp[n] = pl*jn[n] - jn[n + 1];
  157. pl = pl - zi;
  158. }
  159. // Calculate the n=0 Bessel Functions
  160. jn0 = zi*std::sin(z);
  161. h1n[0] = std::exp(std::complex<double>(0.0, 1.0)*z)*zi*(-std::complex<double>(0.0, 1.0));
  162. h1np[0] = h1n[0]*(std::complex<double>(0.0, 1.0) - zi);
  163. // Rescale j[n], j'[n], converting to spherical Bessel functions.
  164. // Recur h1[n], h1'[n] as spherical Bessel functions.
  165. w = 1.0/jn[0];
  166. pl = zi;
  167. for (n = 0; n < nmax; n++) {
  168. jn[n] = jn0*(w*jn[n]);
  169. jnp[n] = jn0*(w*jnp[n]) - zi*jn[n];
  170. if (n != 0) {
  171. h1n[n] = (pl - zi)*h1n[n - 1] - h1np[n - 1];
  172. // check if hankel is increasing (upward stable)
  173. if (std::abs(h1n[n]) < std::abs(h1n[n - 1])) {
  174. xdb = z;
  175. h1nldb = h1n[n];
  176. h1nbdb = h1n[n - 1];
  177. }
  178. pl += zi;
  179. h1np[n] = -(pl*h1n[n]) + h1n[n - 1];
  180. }
  181. }
  182. // success
  183. return 0;
  184. }
  185. //**********************************************************************************//
  186. // This function calculates the spherical Bessel functions (bj and by) and the //
  187. // logarithmic derivative (bd) for a given complex value z. See pag. 87 B&H. //
  188. // //
  189. // Input parameters: //
  190. // z: Complex argument to evaluate bj, by and bd //
  191. // nmax: Maximum number of terms to calculate bj, by and bd //
  192. // //
  193. // Output parameters: //
  194. // bj, by: Spherical Bessel functions //
  195. // bd: Logarithmic derivative //
  196. //**********************************************************************************//
  197. void sphericalBessel(std::complex<double> z, int nmax, std::vector<double>& bj, std::vector<double>& by, std::vector<double>& bd) {
  198. std::vector<std::complex<double> > jn, jnp, h1n, h1np;
  199. jn.resize(nmax);
  200. jnp.resize(nmax);
  201. h1n.resize(nmax);
  202. h1np.resize(nmax);
  203. int ifail = sbesjh(z, nmax, jn, jnp, h1n, h1np);
  204. for (int n = 0; n < nmax; n++) {
  205. bj[n] = jn[n];
  206. by[n] = (h1n[n] - jn[n])/std::complex<double>(0.0, 1.0);
  207. bd[n] = jnp[n]/jn[n] + 1.0/z;
  208. }
  209. }
  210. // external scattering field = incident + scattered
  211. // BH p.92 (4.37), 94 (4.45), 95 (4.50)
  212. // assume: medium is non-absorbing; refim = 0; Uabs = 0
  213. int fieldExt(int nmax, double Rho, double Phi, double Theta, std::vector<double> Pi, std::vector<double> Tau,
  214. std::vector<std::complex<double> > an, std::vector<std::complex<double> > bn,
  215. std::vector<std::complex<double> >& E, std::vector<std::complex<double> >& H) {
  216. int i, n;
  217. double rn = 0.0;
  218. std::complex<double> zn, xxip, encap;
  219. std::vector<std::complex<double> > vm3o1n, vm3e1n, vn3o1n, vn3e1n;
  220. vm3o1n.resize(3);
  221. vm3e1n.resize(3);
  222. vn3o1n.resize(3);
  223. vn3e1n.resize(3);
  224. std::vector<std::complex<double> > Ei, Hi, Es, Hs;
  225. Ei.resize(3);
  226. Hi.resize(3);
  227. Es.resize(3);
  228. Hs.resize(3);
  229. for (i = 0; i < 3; i++) {
  230. Ei[i] = std::complex<double>(0.0, 0.0);
  231. Hi[i] = std::complex<double>(0.0, 0.0);
  232. Es[i] = std::complex<double>(0.0, 0.0);
  233. Hs[i] = std::complex<double>(0.0, 0.0);
  234. }
  235. std::vector<std::complex<double> > bj, by, bd;
  236. bj.resize(nmax);
  237. by.resize(nmax);
  238. bd.resize(nmax);
  239. // Calculate spherical Bessel and Hankel functions
  240. sphericalBessel(Rho, nmax, bj, by, bd);
  241. for (n = 0; n < nmax; n++) {
  242. rn = double(n + 1);
  243. zn = bj[n] + std::complex<double>(0.0, 1.0)*by[n];
  244. xxip = Rho*(bj[n] + std::complex<double>(0.0, 1.0)*by[n]) - rn*zn;
  245. vm3o1n[0] = std::complex<double>(0.0, 0.0);
  246. vm3o1n[1] = std::cos(Phi)*Pi[n]*zn;
  247. vm3o1n[2] = -(std::sin(Phi)*Tau[n]*zn);
  248. vm3e1n[0] = std::complex<double>(0.0, 0.0);
  249. vm3e1n[1] = -(std::sin(Phi)*Pi[n]*zn);
  250. vm3e1n[2] = -(std::cos(Phi)*Tau[n]*zn);
  251. vn3o1n[0] = std::sin(Phi)*rn*(rn + 1.0)*std::sin(Theta)*Pi[n]*zn/Rho;
  252. vn3o1n[1] = std::sin(Phi)*Tau[n]*xxip/Rho;
  253. vn3o1n[2] = std::cos(Phi)*Pi[n]*xxip/Rho;
  254. vn3e1n[0] = std::cos(Phi)*rn*(rn + 1.0)*std::sin(Theta)*Pi[n]*zn/Rho;
  255. vn3e1n[1] = std::cos(Phi)*Tau[n]*xxip/Rho;
  256. vn3e1n[2] = -(std::sin(Phi)*Pi[n]*xxip/Rho);
  257. // scattered field: BH p.94 (4.45)
  258. encap = std::pow(std::complex<double>(0.0, 1.0), rn)*(2.0*rn + 1.0)/(rn*(rn + 1.0));
  259. for (i = 0; i < 3; i++) {
  260. Es[i] = Es[i] + encap*(std::complex<double>(0.0, 1.0)*an[n]*vn3e1n[i] - bn[n]*vm3o1n[i]);
  261. Hs[i] = Hs[i] + encap*(std::complex<double>(0.0, 1.0)*bn[n]*vn3o1n[i] + an[n]*vm3e1n[i]);
  262. }
  263. }
  264. // incident E field: BH p.89 (4.21); cf. p.92 (4.37), p.93 (4.38)
  265. // basis unit vectors = er, etheta, ephi
  266. std::complex<double> eifac = std::exp(std::complex<double>(0.0, 1.0)*Rho*std::cos(Theta));
  267. Ei[0] = eifac*std::sin(Theta)*std::cos(Phi);
  268. Ei[1] = eifac*std::cos(Theta)*std::cos(Phi);
  269. Ei[2] = -(eifac*std::sin(Phi));
  270. // magnetic field
  271. std::complex<double> hffact = std::complex<double>(refmed, 0.0)/(cc*mu);
  272. for (i = 0; i < 3; i++) {
  273. Hs[i] = hffact*Hs[i];
  274. }
  275. // incident H field: BH p.26 (2.43), p.89 (4.21)
  276. std::complex<double> hffacta = hffact;
  277. std::complex<double> hifac = eifac*hffacta;
  278. Hi[0] = hifac*std::sin(Theta)*std::sin(Phi);
  279. Hi[1] = hifac*std::cos(Theta)*std::sin(Phi);
  280. Hi[2] = hifac*std::cos(Phi);
  281. for (i = 0; i < 3; i++) {
  282. // electric field E [V m-1] = EF*E0
  283. E[i] = Ei[i] + Es[i];
  284. H[i] = Hi[i] + Hs[i];
  285. }
  286. }
  287. // Calculate an - equation (5)
  288. std::complex<double> calc_an(int n, double XL, std::complex<double> Ha, std::complex<double> mL,
  289. std::complex<double> PsiXL, std::complex<double> ZetaXL,
  290. std::complex<double> PsiXLM1, std::complex<double> ZetaXLM1) {
  291. std::complex<double> Num = (Ha/mL + n/XL)*PsiXL - PsiXLM1;
  292. std::complex<double> Denom = (Ha/mL + n/XL)*ZetaXL - ZetaXLM1;
  293. return Num/Denom;
  294. }
  295. // Calculate bn - equation (6)
  296. std::complex<double> calc_bn(int n, double XL, std::complex<double> Hb, std::complex<double> mL,
  297. std::complex<double> PsiXL, std::complex<double> ZetaXL,
  298. std::complex<double> PsiXLM1, std::complex<double> ZetaXLM1) {
  299. std::complex<double> Num = (mL*Hb + n/XL)*PsiXL - PsiXLM1;
  300. std::complex<double> Denom = (mL*Hb + n/XL)*ZetaXL - ZetaXLM1;
  301. return Num/Denom;
  302. }
  303. // Calculates S1 - equation (25a)
  304. std::complex<double> calc_S1(int n, std::complex<double> an, std::complex<double> bn,
  305. double Pi, double Tau) {
  306. return double(n + n + 1)*(Pi*an + Tau*bn)/double(n*n + n);
  307. }
  308. // Calculates S2 - equation (25b) (it's the same as (25a), just switches Pi and Tau)
  309. std::complex<double> calc_S2(int n, std::complex<double> an, std::complex<double> bn,
  310. double Pi, double Tau) {
  311. return calc_S1(n, an, bn, Tau, Pi);
  312. }
  313. //**********************************************************************************//
  314. // This function calculates the Riccati-Bessel functions (Psi and Zeta) for a //
  315. // real argument (x). //
  316. // Equations (20a) - (21b) //
  317. // //
  318. // Input parameters: //
  319. // x: Real argument to evaluate Psi and Zeta //
  320. // nmax: Maximum number of terms to calculate Psi and Zeta //
  321. // //
  322. // Output parameters: //
  323. // Psi, Zeta: Riccati-Bessel functions //
  324. //**********************************************************************************//
  325. void calcPsiZeta(double x, int nmax,
  326. std::vector<std::complex<double> > D1,
  327. std::vector<std::complex<double> > D3,
  328. std::vector<std::complex<double> >& Psi,
  329. std::vector<std::complex<double> >& Zeta) {
  330. int n;
  331. //Upward recurrence for Psi and Zeta - equations (20a) - (21b)
  332. Psi[0] = std::complex<double>(sin(x), 0);
  333. Zeta[0] = std::complex<double>(sin(x), -cos(x));
  334. for (n = 1; n <= nmax; n++) {
  335. Psi[n] = Psi[n - 1]*(n/x - D1[n - 1]);
  336. Zeta[n] = Zeta[n - 1]*(n/x - D3[n - 1]);
  337. }
  338. }
  339. //**********************************************************************************//
  340. // This function calculates the logarithmic derivatives of the Riccati-Bessel //
  341. // functions (D1 and D3) for a complex argument (z). //
  342. // Equations (16a), (16b) and (18a) - (18d) //
  343. // //
  344. // Input parameters: //
  345. // z: Complex argument to evaluate D1 and D3 //
  346. // nmax: Maximum number of terms to calculate D1 and D3 //
  347. // //
  348. // Output parameters: //
  349. // D1, D3: Logarithmic derivatives of the Riccati-Bessel functions //
  350. //**********************************************************************************//
  351. void calcD1D3(std::complex<double> z, int nmax,
  352. std::vector<std::complex<double> >& D1,
  353. std::vector<std::complex<double> >& D3) {
  354. int n;
  355. std::vector<std::complex<double> > PsiZeta;
  356. PsiZeta.resize(nmax + 1);
  357. // Downward recurrence for D1 - equations (16a) and (16b)
  358. D1[nmax] = std::complex<double>(0.0, 0.0);
  359. for (n = nmax; n > 0; n--) {
  360. D1[n - 1] = double(n)/z - 1.0/(D1[n] + double(n)/z);
  361. }
  362. // Upward recurrence for PsiZeta and D3 - equations (18a) - (18d)
  363. PsiZeta[0] = 0.5*(1.0 - std::complex<double>(cos(2.0*z.real()), sin(2.0*z.real()))*exp(-2.0*z.imag()));
  364. D3[0] = std::complex<double>(0.0, 1.0);
  365. for (n = 1; n <= nmax; n++) {
  366. PsiZeta[n] = PsiZeta[n - 1]*(double(n)/z - D1[n - 1])*(double(n)/z- D3[n - 1]);
  367. D3[n] = D1[n] + std::complex<double>(0.0, 1.0)/PsiZeta[n];
  368. }
  369. }
  370. //**********************************************************************************//
  371. // This function calculates Pi and Tau for all values of Theta. //
  372. // Equations (26a) - (26c) //
  373. // //
  374. // Input parameters: //
  375. // nmax: Maximum number of terms to calculate Pi and Tau //
  376. // nTheta: Number of scattering angles //
  377. // Theta: Array containing all the scattering angles where the scattering //
  378. // amplitudes will be calculated //
  379. // //
  380. // Output parameters: //
  381. // Pi, Tau: Angular functions Pi and Tau, as defined in equations (26a) - (26c) //
  382. //**********************************************************************************//
  383. void calcPiTau(int nmax, int nTheta, std::vector<double> Theta,
  384. std::vector< std::vector<double> >& Pi,
  385. std::vector< std::vector<double> >& Tau) {
  386. int n, t;
  387. //****************************************************//
  388. // Equations (26a) - (26c) //
  389. //****************************************************//
  390. for (t = 0; t < nTheta; t++) {
  391. for (n = 0; n < nmax; n++) {
  392. if (n == 0) {
  393. // Initialize Pi and Tau
  394. Pi[t][n] = 1.0;
  395. Tau[t][n] = (n + 1)*cos(Theta[t]);
  396. } else {
  397. // Calculate the actual values
  398. Pi[t][n] = ((n == 1) ? ((n + n + 1)*cos(Theta[t])*Pi[t][n - 1]/n)
  399. : (((n + n + 1)*cos(Theta[t])*Pi[t][n - 1] - (n + 1)*Pi[t][n - 2])/n));
  400. Tau[t][n] = (n + 1)*cos(Theta[t])*Pi[t][n] - (n + 2)*Pi[t][n - 1];
  401. }
  402. }
  403. }
  404. }
  405. //**********************************************************************************//
  406. // This function calculates the scattering coefficients required to calculate //
  407. // both the near- and far-field parameters. //
  408. // //
  409. // Input parameters: //
  410. // L: Number of layers //
  411. // pl: Index of PEC layer. If there is none just send -1 //
  412. // x: Array containing the size parameters of the layers [0..L-1] //
  413. // m: Array containing the relative refractive indexes of the layers [0..L-1] //
  414. // nmax: Maximum number of multipolar expansion terms to be used for the //
  415. // calculations. Only used if you know what you are doing, otherwise set //
  416. // this parameter to -1 and the function will calculate it. //
  417. // //
  418. // Output parameters: //
  419. // an, bn: Complex scattering amplitudes //
  420. // //
  421. // Return value: //
  422. // Number of multipolar expansion terms used for the calculations //
  423. //**********************************************************************************//
  424. int ScattCoeffs(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m, int nmax,
  425. std::vector<std::complex<double> >& an, std::vector<std::complex<double> >& bn) {
  426. //************************************************************************//
  427. // Calculate the index of the first layer. It can be either 0 (default) //
  428. // or the index of the outermost PEC layer. In the latter case all layers //
  429. // below the PEC are discarded. //
  430. //************************************************************************//
  431. int fl = (pl > 0) ? pl : 0;
  432. if (nmax <= 0) {
  433. nmax = Nmax(L, fl, pl, x, m);
  434. }
  435. std::complex<double> z1, z2;
  436. std::complex<double> Num, Denom;
  437. std::complex<double> G1, G2;
  438. std::complex<double> Temp;
  439. int n, l;
  440. //**************************************************************************//
  441. // Note that since Fri, Nov 14, 2014 all arrays start from 0 (zero), which //
  442. // means that index = layer number - 1 or index = n - 1. The only exception //
  443. // are the arrays for representing D1, D3 and Q because they need a value //
  444. // for the index 0 (zero), hence it is important to consider this shift //
  445. // between different arrays. The change was done to optimize memory usage. //
  446. //**************************************************************************//
  447. // Allocate memory to the arrays
  448. std::vector<std::vector<std::complex<double> > > D1_mlxl, D1_mlxlM1;
  449. D1_mlxl.resize(L);
  450. D1_mlxlM1.resize(L);
  451. std::vector<std::vector<std::complex<double> > > D3_mlxl, D3_mlxlM1;
  452. D3_mlxl.resize(L);
  453. D3_mlxlM1.resize(L);
  454. std::vector<std::vector<std::complex<double> > > Q;
  455. Q.resize(L);
  456. std::vector<std::vector<std::complex<double> > > Ha, Hb;
  457. Ha.resize(L);
  458. Hb.resize(L);
  459. for (l = 0; l < L; l++) {
  460. D1_mlxl[l].resize(nmax + 1);
  461. D1_mlxlM1[l].resize(nmax + 1);
  462. D3_mlxl[l].resize(nmax + 1);
  463. D3_mlxlM1[l].resize(nmax + 1);
  464. Q[l].resize(nmax + 1);
  465. Ha[l].resize(nmax);
  466. Hb[l].resize(nmax);
  467. }
  468. an.resize(nmax);
  469. bn.resize(nmax);
  470. std::vector<std::complex<double> > D1XL, D3XL;
  471. D1XL.resize(nmax + 1);
  472. D3XL.resize(nmax + 1);
  473. std::vector<std::complex<double> > PsiXL, ZetaXL;
  474. PsiXL.resize(nmax + 1);
  475. ZetaXL.resize(nmax + 1);
  476. //*************************************************//
  477. // Calculate D1 and D3 for z1 in the first layer //
  478. //*************************************************//
  479. if (fl == pl) { // PEC layer
  480. for (n = 0; n <= nmax; n++) {
  481. D1_mlxl[fl][n] = std::complex<double>(0.0, -1.0);
  482. D3_mlxl[fl][n] = std::complex<double>(0.0, 1.0);
  483. }
  484. } else { // Regular layer
  485. z1 = x[fl]* m[fl];
  486. // Calculate D1 and D3
  487. calcD1D3(z1, nmax, D1_mlxl[fl], D3_mlxl[fl]);
  488. }
  489. //******************************************************************//
  490. // Calculate Ha and Hb in the first layer - equations (7a) and (8a) //
  491. //******************************************************************//
  492. for (n = 0; n < nmax; n++) {
  493. Ha[fl][n] = D1_mlxl[fl][n + 1];
  494. Hb[fl][n] = D1_mlxl[fl][n + 1];
  495. }
  496. //*****************************************************//
  497. // Iteration from the second layer to the last one (L) //
  498. //*****************************************************//
  499. for (l = fl + 1; l < L; l++) {
  500. //************************************************************//
  501. //Calculate D1 and D3 for z1 and z2 in the layers fl+1..L //
  502. //************************************************************//
  503. z1 = x[l]*m[l];
  504. z2 = x[l - 1]*m[l];
  505. //Calculate D1 and D3 for z1
  506. calcD1D3(z1, nmax, D1_mlxl[l], D3_mlxl[l]);
  507. //Calculate D1 and D3 for z2
  508. calcD1D3(z2, nmax, D1_mlxlM1[l], D3_mlxlM1[l]);
  509. //*********************************************//
  510. //Calculate Q, Ha and Hb in the layers fl+1..L //
  511. //*********************************************//
  512. // Upward recurrence for Q - equations (19a) and (19b)
  513. Num = exp(-2.0*(z1.imag() - z2.imag()))*std::complex<double>(cos(-2.0*z2.real()) - exp(-2.0*z2.imag()), sin(-2.0*z2.real()));
  514. Denom = std::complex<double>(cos(-2.0*z1.real()) - exp(-2.0*z1.imag()), sin(-2.0*z1.real()));
  515. Q[l][0] = Num/Denom;
  516. for (n = 1; n <= nmax; n++) {
  517. Num = (z1*D1_mlxl[l][n] + double(n))*(double(n) - z1*D3_mlxl[l][n - 1]);
  518. Denom = (z2*D1_mlxlM1[l][n] + double(n))*(double(n) - z2*D3_mlxlM1[l][n - 1]);
  519. Q[l][n] = (((x[l - 1]*x[l - 1])/(x[l]*x[l])* Q[l][n - 1])*Num)/Denom;
  520. }
  521. // Upward recurrence for Ha and Hb - equations (7b), (8b) and (12) - (15)
  522. for (n = 1; n <= nmax; n++) {
  523. //Ha
  524. if ((l - 1) == pl) { // The layer below the current one is a PEC layer
  525. G1 = -D1_mlxlM1[l][n];
  526. G2 = -D3_mlxlM1[l][n];
  527. } else {
  528. G1 = (m[l]*Ha[l - 1][n - 1]) - (m[l - 1]*D1_mlxlM1[l][n]);
  529. G2 = (m[l]*Ha[l - 1][n - 1]) - (m[l - 1]*D3_mlxlM1[l][n]);
  530. }
  531. Temp = Q[l][n]*G1;
  532. Num = (G2*D1_mlxl[l][n]) - (Temp*D3_mlxl[l][n]);
  533. Denom = G2 - Temp;
  534. Ha[l][n - 1] = Num/Denom;
  535. //Hb
  536. if ((l - 1) == pl) { // The layer below the current one is a PEC layer
  537. G1 = Hb[l - 1][n - 1];
  538. G2 = Hb[l - 1][n - 1];
  539. } else {
  540. G1 = (m[l - 1]*Hb[l - 1][n - 1]) - (m[l]*D1_mlxlM1[l][n]);
  541. G2 = (m[l - 1]*Hb[l - 1][n - 1]) - (m[l]*D3_mlxlM1[l][n]);
  542. }
  543. Temp = Q[l][n]*G1;
  544. Num = (G2*D1_mlxl[l][n]) - (Temp* D3_mlxl[l][n]);
  545. Denom = (G2- Temp);
  546. Hb[l][n - 1] = (Num/ Denom);
  547. }
  548. }
  549. //**************************************//
  550. //Calculate D1, D3, Psi and Zeta for XL //
  551. //**************************************//
  552. // Calculate D1XL and D3XL
  553. calcD1D3(x[L - 1], nmax, D1XL, D3XL);
  554. // Calculate PsiXL and ZetaXL
  555. calcPsiZeta(x[L - 1], nmax, D1XL, D3XL, PsiXL, ZetaXL);
  556. //*********************************************************************//
  557. // Finally, we calculate the scattering coefficients (an and bn) and //
  558. // the angular functions (Pi and Tau). Note that for these arrays the //
  559. // first layer is 0 (zero), in future versions all arrays will follow //
  560. // this convention to save memory. (13 Nov, 2014) //
  561. //*********************************************************************//
  562. for (n = 0; n < nmax; n++) {
  563. //********************************************************************//
  564. //Expressions for calculating an and bn coefficients are not valid if //
  565. //there is only one PEC layer (ie, for a simple PEC sphere). //
  566. //********************************************************************//
  567. if (pl < (L - 1)) {
  568. 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]);
  569. 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]);
  570. } else {
  571. an[n] = calc_an(n + 1, x[L - 1], std::complex<double>(0.0, 0.0), std::complex<double>(1.0, 0.0), PsiXL[n + 1], ZetaXL[n + 1], PsiXL[n], ZetaXL[n]);
  572. bn[n] = PsiXL[n + 1]/ZetaXL[n + 1];
  573. }
  574. }
  575. return nmax;
  576. }
  577. //**********************************************************************************//
  578. // This function calculates the actual scattering parameters and amplitudes //
  579. // //
  580. // Input parameters: //
  581. // L: Number of layers //
  582. // pl: Index of PEC layer. If there is none just send -1 //
  583. // x: Array containing the size parameters of the layers [0..L-1] //
  584. // m: Array containing the relative refractive indexes of the layers [0..L-1] //
  585. // nTheta: Number of scattering angles //
  586. // Theta: Array containing all the scattering angles where the scattering //
  587. // amplitudes will be calculated //
  588. // nmax: Maximum number of multipolar expansion terms to be used for the //
  589. // calculations. Only used if you know what you are doing, otherwise set //
  590. // this parameter to -1 and the function will calculate it //
  591. // //
  592. // Output parameters: //
  593. // Qext: Efficiency factor for extinction //
  594. // Qsca: Efficiency factor for scattering //
  595. // Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca) //
  596. // Qbk: Efficiency factor for backscattering //
  597. // Qpr: Efficiency factor for the radiation pressure //
  598. // g: Asymmetry factor (g = (Qext-Qpr)/Qsca) //
  599. // Albedo: Single scattering albedo (Albedo = Qsca/Qext) //
  600. // S1, S2: Complex scattering amplitudes //
  601. // //
  602. // Return value: //
  603. // Number of multipolar expansion terms used for the calculations //
  604. //**********************************************************************************//
  605. int nMie(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m,
  606. int nTheta, std::vector<double> Theta, int nmax,
  607. double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo,
  608. std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
  609. int i, n, t;
  610. std::vector<std::complex<double> > an, bn;
  611. std::complex<double> Qbktmp;
  612. // Calculate scattering coefficients
  613. nmax = ScattCoeffs(L, pl, x, m, nmax, an, bn);
  614. std::vector< std::vector<double> > Pi, Tau;
  615. Pi.resize(nTheta);
  616. Tau.resize(nTheta);
  617. for (t = 0; t < nTheta; t++) {
  618. Pi[t].resize(nmax);
  619. Tau[t].resize(nmax);
  620. }
  621. calcPiTau(nmax, nTheta, Theta, Pi, Tau);
  622. double x2 = x[L - 1]*x[L - 1];
  623. // Initialize the scattering parameters
  624. *Qext = 0;
  625. *Qsca = 0;
  626. *Qabs = 0;
  627. *Qbk = 0;
  628. Qbktmp = std::complex<double>(0.0, 0.0);
  629. *Qpr = 0;
  630. *g = 0;
  631. *Albedo = 0;
  632. // Initialize the scattering amplitudes
  633. for (t = 0; t < nTheta; t++) {
  634. S1[t] = std::complex<double>(0.0, 0.0);
  635. S2[t] = std::complex<double>(0.0, 0.0);
  636. }
  637. // By using downward recurrence we avoid loss of precision due to float rounding errors
  638. // See: https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html
  639. // http://en.wikipedia.org/wiki/Loss_of_significance
  640. for (i = nmax - 2; i >= 0; i--) {
  641. n = i + 1;
  642. // Equation (27)
  643. *Qext += (n + n + 1)*(an[i].real() + bn[i].real());
  644. // Equation (28)
  645. *Qsca += (n + n + 1)*(an[i].real()*an[i].real() + an[i].imag()*an[i].imag() + bn[i].real()*bn[i].real() + bn[i].imag()*bn[i].imag());
  646. // Equation (29)
  647. // We must check carefully this equation. If we remove the typecast to double then the result changes. Which is the correct one??? Ovidio (2014/12/10)
  648. *Qpr += ((n*(n + 2)/(n + 1))*((an[i]*std::conj(an[n]) + bn[i]*std::conj(bn[n])).real()) + ((double)(n + n + 1)/(n*(n + 1)))*(an[i]*std::conj(bn[i])).real());
  649. // Equation (33)
  650. Qbktmp = Qbktmp + (double)(n + n + 1)*(1 - 2*(n % 2))*(an[i]- bn[i]);
  651. //****************************************************//
  652. // Calculate the scattering amplitudes (S1 and S2) //
  653. // Equations (25a) - (25b) //
  654. //****************************************************//
  655. for (t = 0; t < nTheta; t++) {
  656. S1[t] += calc_S1(n, an[i], bn[i], Pi[t][i], Tau[t][i]);
  657. S2[t] += calc_S2(n, an[i], bn[i], Pi[t][i], Tau[t][i]);
  658. }
  659. }
  660. *Qext = 2*(*Qext)/x2; // Equation (27)
  661. *Qsca = 2*(*Qsca)/x2; // Equation (28)
  662. *Qpr = *Qext - 4*(*Qpr)/x2; // Equation (29)
  663. *Qabs = *Qext - *Qsca; // Equation (30)
  664. *Albedo = *Qsca / *Qext; // Equation (31)
  665. *g = (*Qext - *Qpr) / *Qsca; // Equation (32)
  666. *Qbk = (Qbktmp.real()*Qbktmp.real() + Qbktmp.imag()*Qbktmp.imag())/x2; // Equation (33)
  667. return nmax;
  668. }
  669. //**********************************************************************************//
  670. // This function is just a wrapper to call the full 'nMie' function with fewer //
  671. // parameters, it is here mainly for compatibility with older versions of the //
  672. // program. Also, you can use it if you neither have a PEC layer nor want to define //
  673. // any limit for the maximum number of terms. //
  674. // //
  675. // Input parameters: //
  676. // L: Number of layers //
  677. // x: Array containing the size parameters of the layers [0..L-1] //
  678. // m: Array containing the relative refractive indexes of the layers [0..L-1] //
  679. // nTheta: Number of scattering angles //
  680. // Theta: Array containing all the scattering angles where the scattering //
  681. // amplitudes will be calculated //
  682. // //
  683. // Output parameters: //
  684. // Qext: Efficiency factor for extinction //
  685. // Qsca: Efficiency factor for scattering //
  686. // Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca) //
  687. // Qbk: Efficiency factor for backscattering //
  688. // Qpr: Efficiency factor for the radiation pressure //
  689. // g: Asymmetry factor (g = (Qext-Qpr)/Qsca) //
  690. // Albedo: Single scattering albedo (Albedo = Qsca/Qext) //
  691. // S1, S2: Complex scattering amplitudes //
  692. // //
  693. // Return value: //
  694. // Number of multipolar expansion terms used for the calculations //
  695. //**********************************************************************************//
  696. int nMie(int L, std::vector<double> x, std::vector<std::complex<double> > m,
  697. int nTheta, std::vector<double> Theta,
  698. double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo,
  699. std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
  700. return nMie(L, -1, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
  701. }
  702. //**********************************************************************************//
  703. // This function is just a wrapper to call the full 'nMie' function with fewer //
  704. // parameters, it is useful if you want to include a PEC layer but not a limit //
  705. // for the maximum number of terms. //
  706. // //
  707. // Input parameters: //
  708. // L: Number of layers //
  709. // pl: Index of PEC layer. If there is none just send -1 //
  710. // x: Array containing the size parameters of the layers [0..L-1] //
  711. // m: Array containing the relative refractive indexes of the layers [0..L-1] //
  712. // nTheta: Number of scattering angles //
  713. // Theta: Array containing all the scattering angles where the scattering //
  714. // amplitudes will be calculated //
  715. // //
  716. // Output parameters: //
  717. // Qext: Efficiency factor for extinction //
  718. // Qsca: Efficiency factor for scattering //
  719. // Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca) //
  720. // Qbk: Efficiency factor for backscattering //
  721. // Qpr: Efficiency factor for the radiation pressure //
  722. // g: Asymmetry factor (g = (Qext-Qpr)/Qsca) //
  723. // Albedo: Single scattering albedo (Albedo = Qsca/Qext) //
  724. // S1, S2: Complex scattering amplitudes //
  725. // //
  726. // Return value: //
  727. // Number of multipolar expansion terms used for the calculations //
  728. //**********************************************************************************//
  729. int nMie(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m,
  730. int nTheta, std::vector<double> Theta,
  731. double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo,
  732. std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
  733. return nMie(L, pl, x, m, nTheta, Theta, -1, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
  734. }
  735. //**********************************************************************************//
  736. // This function is just a wrapper to call the full 'nMie' function with fewer //
  737. // parameters, it is useful if you want to include a limit for the maximum number //
  738. // of terms but not a PEC layer. //
  739. // //
  740. // Input parameters: //
  741. // L: Number of layers //
  742. // x: Array containing the size parameters of the layers [0..L-1] //
  743. // m: Array containing the relative refractive indexes of the layers [0..L-1] //
  744. // nTheta: Number of scattering angles //
  745. // Theta: Array containing all the scattering angles where the scattering //
  746. // amplitudes will be calculated //
  747. // nmax: Maximum number of multipolar expansion terms to be used for the //
  748. // calculations. Only used if you know what you are doing, otherwise set //
  749. // this parameter to -1 and the function will calculate it //
  750. // //
  751. // Output parameters: //
  752. // Qext: Efficiency factor for extinction //
  753. // Qsca: Efficiency factor for scattering //
  754. // Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca) //
  755. // Qbk: Efficiency factor for backscattering //
  756. // Qpr: Efficiency factor for the radiation pressure //
  757. // g: Asymmetry factor (g = (Qext-Qpr)/Qsca) //
  758. // Albedo: Single scattering albedo (Albedo = Qsca/Qext) //
  759. // S1, S2: Complex scattering amplitudes //
  760. // //
  761. // Return value: //
  762. // Number of multipolar expansion terms used for the calculations //
  763. //**********************************************************************************//
  764. int nMie(int L, std::vector<double> x, std::vector<std::complex<double> > m,
  765. int nTheta, std::vector<double> Theta, int nmax,
  766. double *Qext, double *Qsca, double *Qabs, double *Qbk, double *Qpr, double *g, double *Albedo,
  767. std::vector<std::complex<double> >& S1, std::vector<std::complex<double> >& S2) {
  768. return nMie(L, -1, x, m, nTheta, Theta, nmax, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2);
  769. }
  770. //**********************************************************************************//
  771. // This function calculates complex electric and magnetic field in the surroundings //
  772. // and inside (TODO) the particle. //
  773. // //
  774. // Input parameters: //
  775. // L: Number of layers //
  776. // pl: Index of PEC layer. If there is none just send 0 (zero) //
  777. // x: Array containing the size parameters of the layers [0..L-1] //
  778. // m: Array containing the relative refractive indexes of the layers [0..L-1] //
  779. // nmax: Maximum number of multipolar expansion terms to be used for the //
  780. // calculations. Only used if you know what you are doing, otherwise set //
  781. // this parameter to 0 (zero) and the function will calculate it. //
  782. // ncoord: Number of coordinate points //
  783. // Coords: Array containing all coordinates where the complex electric and //
  784. // magnetic fields will be calculated //
  785. // //
  786. // Output parameters: //
  787. // E, H: Complex electric and magnetic field at the provided coordinates //
  788. // //
  789. // Return value: //
  790. // Number of multipolar expansion terms used for the calculations //
  791. //**********************************************************************************//
  792. int nField(int L, int pl, std::vector<double> x, std::vector<std::complex<double> > m, int nmax,
  793. int ncoord, std::vector<double> Xp, std::vector<double> Yp, std::vector<double> Zp,
  794. std::vector<std::vector<std::complex<double> > >& E, std::vector<std::vector<std::complex<double> > >& H) {
  795. int i, n, c;
  796. std::vector<std::complex<double> > an, bn;
  797. // Calculate scattering coefficients
  798. nmax = ScattCoeffs(L, pl, x, m, nmax, an, bn);
  799. std::vector<double> Rho, Phi, Theta;
  800. Rho.resize(ncoord);
  801. Phi.resize(ncoord);
  802. Theta.resize(ncoord);
  803. for (c = 0; c < ncoord; c++) {
  804. // Convert to spherical coordinates
  805. Rho = sqrt(Xp[c]*Xp[c] + Yp[c]*Yp[c] + Zp[c]*Zp[c]);
  806. if (Rho < 1e-3) {
  807. Rho = 1e-3;
  808. }
  809. Phi = acos(Xp[c]/sqrt(Xp[c]*Xp[c] + Yp[c]*Yp[c]));
  810. Theta = acos(Xp[c]/Rho[c]);
  811. }
  812. std::vector< std::vector<double> > Pi, Tau;
  813. Pi.resize(ncoord);
  814. Tau.resize(ncoord);
  815. for (c = 0; c < ncoord; c++) {
  816. Pi[c].resize(nmax);
  817. Tau[c].resize(nmax);
  818. }
  819. calcPiTau(nmax, ncoord, Theta, Pi, Tau);
  820. for (c = 0; c < ncoord; c++) {
  821. //*******************************************************//
  822. // external scattering field = incident + scattered //
  823. // BH p.92 (4.37), 94 (4.45), 95 (4.50) //
  824. // assume: medium is non-absorbing; refim = 0; Uabs = 0 //
  825. //*******************************************************//
  826. // Initialize the fields
  827. E[c] = std::complex<double>(0.0, 0.0);
  828. H[c] = std::complex<double>(0.0, 0.0);
  829. // Firstly the easiest case, we want the field outside the particle
  830. if (Rho >= x[L - 1]) {
  831. fieldExt(nmax, Rho, Phi, Theta, Pi[c], Tau[c], an, bn, E[c], H[c]);
  832. }
  833. }
  834. return nmax;
  835. }