main.py 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356
  1. #!/usr/bin/env python
  2. # -*- coding: UTF-8 -*-
  3. #
  4. # Copyright (C) 2009-2021 Ovidio Peña Rodríguez <ovidio@bytesfall.com>
  5. # Copyright (C) 2013-2021 Konstantin Ladutenko <kostyfisik@gmail.com>
  6. #
  7. # This file is part of scattnlay
  8. #
  9. # This program is free software: you can redistribute it and/or modify
  10. # it under the terms of the GNU General Public License as published by
  11. # the Free Software Foundation, either version 3 of the License, or
  12. # (at your option) any later version.
  13. #
  14. # This program is distributed in the hope that it will be useful,
  15. # but WITHOUT ANY WARRANTY; without even the implied warranty of
  16. # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  17. # GNU General Public License for more details.
  18. #
  19. # The only additional remark is that we expect that all publications
  20. # describing work using this software, or all commercial products
  21. # using it, cite at least one of the following references:
  22. # [1] O. Peña and U. Pal, "Scattering of electromagnetic radiation by
  23. # a multilayered sphere," Computer Physics Communications,
  24. # vol. 180, Nov. 2009, pp. 2348-2354.
  25. # [2] K. Ladutenko, U. Pal, A. Rivera, and O. Peña-Rodríguez, "Mie
  26. # calculation of electromagnetic near-field for a multilayered
  27. # sphere," Computer Physics Communications, vol. 214, May 2017,
  28. # pp. 225-230.
  29. #
  30. # You should have received a copy of the GNU General Public License
  31. # along with this program. If not, see <http://www.gnu.org/licenses/>.
  32. import numpy as np
  33. from scattnlay_mp import mie_mp as mie_mp_
  34. from scattnlay_dp import mie_dp
  35. mie = mie_dp()
  36. mie_mp = mie_mp_()
  37. def scattcoeffs_(x, m, nmax=-1, pl=-1, mp=False):
  38. if mp:
  39. from scattnlay_mp import mie_mp as mie_
  40. else:
  41. from scattnlay_dp import mie_dp as mie_
  42. # from scattnlay_mp import mie_mp as mie_
  43. mie = mie_()
  44. mie.SetLayersSize(x)
  45. mie.SetLayersIndex(m)
  46. mie.SetPECLayer(pl)
  47. mie.SetMaxTerms(nmax)
  48. mie.calcScattCoeffs()
  49. terms = mie.GetMaxTerms()
  50. a = mie.GetAn()
  51. b = mie.GetBn()
  52. return terms, a, b
  53. def scattcoeffs(x, m, nmax=-1, pl=-1, mp=False):
  54. """
  55. scattcoeffs(x, m[, nmax, pl, mp])
  56. Calculate the scattering coefficients required to calculate both the
  57. near- and far-field parameters.
  58. x: Size parameters (1D or 2D ndarray)
  59. m: Relative refractive indices (1D or 2D ndarray)
  60. nmax: Maximum number of multipolar expansion terms to be used for the
  61. calculations. Only use it if you know what you are doing, otherwise
  62. set this parameter to -1 and the function will calculate it.
  63. pl: Index of PEC layer. If there is none just send -1.
  64. mp: Use multiple (True) or double (False) precision.
  65. Returns: (terms, an, bn)
  66. with
  67. terms: Number of multipolar expansion terms used for the calculations
  68. an, bn: Complex scattering coefficients
  69. """
  70. if len(m.shape) != 1 and len(m.shape) != 2:
  71. raise ValueError('The relative refractive index (m) should be a 1-D or 2-D NumPy array.')
  72. if len(x.shape) == 1:
  73. if len(m.shape) == 1:
  74. return scattcoeffs_(x, m, nmax=nmax, pl=pl, mp=mp)
  75. else:
  76. raise ValueError('The number of of dimensions for the relative refractive index (m) and for the size parameter (x) must be equal.')
  77. elif len(x.shape) != 2:
  78. raise ValueError('The size parameter (x) should be a 1-D or 2-D NumPy array.')
  79. # Repeat the same m for all wavelengths
  80. if len(m.shape) == 1:
  81. m = np.repeat(m[np.newaxis, :], x.shape[0], axis=0)
  82. if nmax == -1:
  83. nstore = 0
  84. else:
  85. nstore = nmax
  86. terms = np.zeros((x.shape[0]), dtype=int)
  87. an = np.zeros((0, nstore), dtype=complex)
  88. bn = np.zeros((0, nstore), dtype=complex)
  89. for i, xi in enumerate(x):
  90. terms[i], a, b = scattcoeffs_(xi, m[i], nmax=nmax, pl=pl, mp=mp)
  91. if terms[i] > nstore:
  92. nstore = terms[i]
  93. an.resize((an.shape[0], nstore))
  94. bn.resize((bn.shape[0], nstore))
  95. an = np.vstack((an, a))
  96. bn = np.vstack((bn, b))
  97. return terms, an, bn
  98. def expancoeffs_(x, m, nmax=-1, pl=-1, mp=False):
  99. if mp:
  100. from scattnlay_mp import mie_mp as mie_
  101. else:
  102. from scattnlay_dp import mie_dp as mie_
  103. # from scattnlay_mp import mie_mp as mie_
  104. mie = mie_()
  105. mie.SetLayersSize(x)
  106. mie.SetLayersIndex(m)
  107. mie.SetPECLayer(pl)
  108. mie.SetMaxTerms(nmax)
  109. mie.calcScattCoeffs()
  110. mie.calcExpanCoeffs()
  111. terms = mie.GetMaxTerms()
  112. an = mie.GetLayerAn()
  113. bn = mie.GetLayerBn()
  114. cn = mie.GetLayerCn()
  115. dn = mie.GetLayerDn()
  116. return terms, an, bn, cn, dn
  117. # TODO verify that expancoeffs() is really working
  118. def expancoeffs(x, m, nmax=-1, pl=-1, mp=False):
  119. """
  120. expancoeffs(x, m[, nmax, pl, mp])
  121. Calculate the scattering coefficients required to calculate both the
  122. near- and far-field parameters.
  123. x: Size parameters (1D or 2D ndarray)
  124. m: Relative refractive indices (1D or 2D ndarray)
  125. nmax: Maximum number of multipolar expansion terms to be used for the
  126. calculations. Only use it if you know what you are doing, otherwise
  127. set this parameter to -1 and the function will calculate it.
  128. pl: Index of PEC layer. If there is none just send -1.
  129. mp: Use multiple (True) or double (False) precision.
  130. Returns: (terms, an, bn, cn, dn)
  131. with
  132. terms: Number of multipolar expansion terms used for the calculations
  133. an, bn, cn, dn: Complex expansion coefficients of each layer
  134. """
  135. if len(m.shape) != 1 and len(m.shape) != 2:
  136. raise ValueError('The relative refractive index (m) should be a 1-D or 2-D NumPy array.')
  137. if len(x.shape) == 1:
  138. if len(m.shape) == 1:
  139. return expancoeffs_(x, m, nmax=nmax, pl=pl, mp=mp)
  140. else:
  141. raise ValueError('The number of of dimensions for the relative refractive index (m) and for the size parameter (x) must be equal.')
  142. elif len(x.shape) != 2:
  143. raise ValueError('The size parameter (x) should be a 1-D or 2-D NumPy array.')
  144. # Repeat the same m for all wavelengths
  145. if len(m.shape) == 1:
  146. m = np.repeat(m[np.newaxis, :], x.shape[0], axis=0)
  147. if nmax == -1:
  148. nstore = 0
  149. else:
  150. nstore = nmax
  151. terms = np.zeros((x.shape[0]), dtype=int)
  152. an = np.zeros((0, x.shape[1]+1, nstore), dtype=complex)
  153. bn = np.zeros((0, x.shape[1]+1, nstore), dtype=complex)
  154. cn = np.zeros((0, x.shape[1]+1, nstore), dtype=complex)
  155. dn = np.zeros((0, x.shape[1]+1, nstore), dtype=complex)
  156. for i, xi in enumerate(x):
  157. terms[i], a, b, c, d = expancoeffs_(xi, m[i], nmax=nmax, pl=pl, mp=mp)
  158. if terms[i] > nstore:
  159. nstore = terms[i]
  160. an.resize((an.shape[0], an.shape[1], nstore))
  161. bn.resize((bn.shape[0], bn.shape[1], nstore))
  162. cn.resize((cn.shape[0], cn.shape[1], nstore))
  163. dn.resize((dn.shape[0], dn.shape[1], nstore))
  164. an = np.vstack((an, [a]))
  165. bn = np.vstack((bn, [b]))
  166. cn = np.vstack((cn, [c]))
  167. dn = np.vstack((dn, [d]))
  168. return terms, an, bn, cn, dn
  169. def scattnlay_(x, m, theta=np.zeros(0, dtype=float), nmax=-1, pl=-1, mp=False):
  170. if mp:
  171. from scattnlay_mp import mie_mp as mie_
  172. else:
  173. from scattnlay_dp import mie_dp as mie_
  174. mie = mie_()
  175. mie.SetLayersSize(x)
  176. mie.SetLayersIndex(m)
  177. mie.SetAngles(theta)
  178. mie.SetPECLayer(pl)
  179. mie.SetMaxTerms(nmax)
  180. mie.RunMieCalculation()
  181. Qext = mie.GetQext()
  182. Qsca = mie.GetQsca()
  183. Qabs = mie.GetQabs()
  184. Qbk = mie.GetQbk()
  185. Qpr = mie.GetQpr()
  186. g = mie.GetAsymmetryFactor()
  187. Albedo = mie.GetAlbedo()
  188. terms = mie.GetMaxTerms()
  189. S1 = mie.GetS1()
  190. S2 = mie.GetS2()
  191. return terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2
  192. def scattnlay(x, m, theta=np.zeros(0, dtype=float), nmax=-1, pl=-1, mp=False):
  193. """
  194. scattnlay(x, m[, theta, nmax, pl, mp])
  195. Calculate the actual scattering parameters and amplitudes.
  196. x: Size parameters (1D or 2D ndarray)
  197. m: Relative refractive indices (1D or 2D ndarray)
  198. theta: Scattering angles where the scattering amplitudes will be
  199. calculated (optional, 1D ndarray)
  200. nmax: Maximum number of multipolar expansion terms to be used for the
  201. calculations. Only use it if you know what you are doing.
  202. pl: Index of PEC layer. If there is none just send -1.
  203. mp: Use multiple (True) or double (False) precision.
  204. Returns: (terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2)
  205. with
  206. terms: Number of multipolar expansion terms used for the calculations
  207. Qext: Efficiency factor for extinction
  208. Qsca: Efficiency factor for scattering
  209. Qabs: Efficiency factor for absorption (Qabs = Qext - Qsca)
  210. Qbk: Efficiency factor for backscattering
  211. Qpr: Efficiency factor for the radiation pressure
  212. g: Asymmetry factor (g = (Qext-Qpr)/Qsca)
  213. Albedo: Single scattering albedo (Albedo = Qsca/Qext)
  214. S1, S2: Complex scattering amplitudes
  215. """
  216. if len(m.shape) != 1 and len(m.shape) != 2:
  217. raise ValueError('The relative refractive index (m) should be a 1-D or 2-D NumPy array.')
  218. if len(x.shape) == 1:
  219. if len(m.shape) == 1:
  220. return scattnlay_(x, m, theta, nmax=nmax, pl=pl, mp=mp)
  221. else:
  222. raise ValueError('The number of of dimensions for the relative refractive index (m) and for the size parameter (x) must be equal.')
  223. elif len(x.shape) != 2:
  224. raise ValueError('The size parameter (x) should be a 1-D or 2-D NumPy array.')
  225. if len(theta.shape) != 1:
  226. raise ValueError('The scattering angles (theta) should be a 1-D NumPy array.')
  227. # Repeat the same m for all wavelengths
  228. if len(m.shape) == 1:
  229. m = np.repeat(m[np.newaxis, :], x.shape[0], axis=0)
  230. terms = np.zeros((x.shape[0]), dtype=int)
  231. Qext = np.zeros((x.shape[0]), dtype=float)
  232. Qsca = np.zeros((x.shape[0]), dtype=float)
  233. Qabs = np.zeros((x.shape[0]), dtype=float)
  234. Qbk = np.zeros((x.shape[0]), dtype=float)
  235. Qpr = np.zeros((x.shape[0]), dtype=float)
  236. g = np.zeros((x.shape[0]), dtype=float)
  237. Albedo = np.zeros((x.shape[0]), dtype=float)
  238. S1 = np.zeros((x.shape[0], theta.shape[0]), dtype=complex)
  239. S2 = np.zeros((x.shape[0], theta.shape[0]), dtype=complex)
  240. for i, xi in enumerate(x):
  241. terms[i], Qext[i], Qsca[i], Qabs[i], Qbk[i], Qpr[i], g[i], Albedo[i], S1[i], S2[i] = scattnlay_(xi, m[i], theta, nmax=nmax, pl=pl, mp=mp)
  242. return terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2
  243. def fieldnlay_(x, m, xp, yp, zp, nmax=-1, pl=-1, mp=False):
  244. if mp:
  245. from scattnlay_mp import mie_mp as mie_
  246. else:
  247. from scattnlay_dp import mie_dp as mie_
  248. # from scattnlay_mp import mie_mp as mie_
  249. mie = mie_()
  250. mie.SetLayersSize(x)
  251. mie.SetLayersIndex(m)
  252. mie.SetPECLayer(pl)
  253. mie.SetMaxTerms(nmax)
  254. mie.SetFieldCoords(xp, yp, zp)
  255. mie.RunFieldCalculation()
  256. terms = mie.GetMaxTerms()
  257. E = mie.GetFieldE()
  258. H = mie.GetFieldH()
  259. return terms, E, H
  260. def fieldnlay(x, m, xp, yp, zp, nmax=-1, pl=-1, mp=False):
  261. """
  262. fieldnlay(x, m, xp, yp, zp[, nmax, pl, mp])
  263. Calculate the actual scattering parameters and amplitudes.
  264. x: Size parameters (1D or 2D ndarray)
  265. m: Relative refractive indices (1D or 2D ndarray)
  266. xp: Array containing all X coordinates to calculate the complex
  267. electric and magnetic fields (1D* ndarray)
  268. yp: Array containing all Y coordinates to calculate the complex
  269. electric and magnetic fields (1D* ndarray)
  270. zp: Array containing all Z coordinates to calculate the complex
  271. electric and magnetic fields (1D* ndarray)
  272. nmax: Maximum number of multipolar expansion terms to be used for the
  273. calculations. Only use it if you know what you are doing.
  274. pl: Index of PEC layer. If there is none just send -1.
  275. mp: Use multiple (True) or double (False) precision.
  276. Returns: (terms, E, H)
  277. with
  278. terms: Number of multipolar expansion terms used for the calculations
  279. E, H: Complex electric and magnetic field at the provided coordinates
  280. *Note: We assume that the coordinates are referred to the first wavelength
  281. (or structure) and correct it for the following ones
  282. """
  283. if len(m.shape) != 1 and len(m.shape) != 2:
  284. raise ValueError('The relative refractive index (m) should be a 1-D or 2-D NumPy array.')
  285. if len(x.shape) == 1:
  286. if len(m.shape) == 1:
  287. return fieldnlay_(x, m, xp, yp, zp, nmax=nmax, pl=pl, mp=mp)
  288. else:
  289. raise ValueError('The number of of dimensions for the relative refractive index (m) and for the size parameter (x) must be equal.')
  290. elif len(x.shape) != 2:
  291. raise ValueError('The size parameter (x) should be a 1-D or 2-D NumPy array.')
  292. # Repeat the same m for all wavelengths
  293. if len(m.shape) == 1:
  294. m = np.repeat(m[np.newaxis, :], x.shape[0], axis=0)
  295. terms = np.zeros((x.shape[0]), dtype=int)
  296. E = np.zeros((x.shape[0], xp.shape[0], 3), dtype=complex)
  297. H = np.zeros((x.shape[0], xp.shape[0], 3), dtype=complex)
  298. for i, xi in enumerate(x):
  299. # (2020/05/12) We assume that the coordinates are referred to the first wavelength
  300. # (or structure) and correct it for the following ones
  301. terms[i], E[i], H[i] = fieldnlay_(xi, m[i], xp*xi[-1]/x[0, -1], yp*xi[-1]/x[0, -1], zp*xi[-1]/x[0, -1], nmax=nmax, pl=pl, mp=mp)
  302. return terms, E, H