nmie.cc 46 KB

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