field-Ag-flow.py 8.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245
  1. #!/usr/bin/env python
  2. # -*- coding: UTF-8 -*-
  3. #
  4. # Copyright (C) 2009-2015 Ovidio Peña Rodríguez <ovidio@bytesfall.com>
  5. #
  6. # This file is part of python-scattnlay
  7. #
  8. # This program is free software: you can redistribute it and/or modify
  9. # it under the terms of the GNU General Public License as published by
  10. # the Free Software Foundation, either version 3 of the License, or
  11. # (at your option) any later version.
  12. #
  13. # This program is distributed in the hope that it will be useful,
  14. # but WITHOUT ANY WARRANTY; without even the implied warranty of
  15. # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  16. # GNU General Public License for more details.
  17. #
  18. # The only additional remark is that we expect that all publications
  19. # describing work using this software, or all commercial products
  20. # using it, cite the following reference:
  21. # [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by
  22. # a multilayered sphere," Computer Physics Communications,
  23. # vol. 180, Nov. 2009, pp. 2348-2354.
  24. #
  25. # You should have received a copy of the GNU General Public License
  26. # along with this program. If not, see <http://www.gnu.org/licenses/>.
  27. # This test case calculates the electric field in the
  28. # E-k plane, for an spherical Si-Ag-Si nanoparticle. Core radius is 17.74 nm,
  29. # inner layer 23.31nm, outer layer 22.95nm. Working wavelength is 800nm, we use
  30. # silicon epsilon=13.64+i0.047, silver epsilon= -28.05+i1.525
  31. import scattnlay
  32. from scattnlay import fieldnlay
  33. from scattnlay import scattnlay
  34. import numpy as np
  35. import cmath
  36. def get_index(array,value):
  37. idx = (np.abs(array-value)).argmin()
  38. return idx
  39. #Ec = np.resize(Ec, (npts, npts)).T
  40. def GetFlow(scale_x, scale_z, Ec, Hc, a, b, nmax):
  41. # Initial position
  42. flow_x = [a]
  43. flow_z = [b]
  44. x_pos = flow_x[-1]
  45. z_pos = flow_z[-1]
  46. x_idx = get_index(scale_x, x_pos)
  47. z_idx = get_index(scale_z, z_pos)
  48. S=np.cross(Ec[npts*z_idx+x_idx], Hc[npts*z_idx+x_idx]).real
  49. #if (np.linalg.norm(S)> 1e-4):
  50. Snorm_prev=S/np.linalg.norm(S)
  51. Snorm_prev=Snorm_prev.real
  52. for n in range(0, nmax):
  53. #Get the next position
  54. #1. Find Poynting vector and normalize it
  55. x_pos = flow_x[-1]
  56. z_pos = flow_z[-1]
  57. x_idx = get_index(scale_x, x_pos)
  58. z_idx = get_index(scale_z, z_pos)
  59. #print x_idx, z_idx
  60. S=np.cross(Ec[npts*z_idx+x_idx], Hc[npts*z_idx+x_idx])
  61. #if (np.linalg.norm(S)> 1e-4):
  62. Snorm=S/np.linalg.norm(S)
  63. Snorm=Snorm.real
  64. #2. Evaluate displacement = half of the discrete and new position
  65. dpos = abs(scale_z[0]-scale_z[1])/2.0
  66. dx = dpos*Snorm[0]/abs(Snorm[2])
  67. dz = dpos*Snorm[2]/abs(Snorm[2])
  68. x_pos = x_pos+dx
  69. z_pos = z_pos+dz
  70. #3. Save result
  71. flow_x.append(x_pos)
  72. flow_z.append(z_pos)
  73. return flow_x, flow_z
  74. # # a)
  75. #WL=400 #nm
  76. #core_r = WL/20.0
  77. #epsilon_Ag = -2.0 + 10.0j
  78. # # b)
  79. #WL=400 #nm
  80. #core_r = WL/20.0
  81. #epsilon_Ag = -2.0 + 1.0j
  82. # c)
  83. WL=354 #nm
  84. core_r = WL/20.0
  85. epsilon_Ag = -2.0 + 0.28j
  86. # # d)
  87. # WL=367 #nm
  88. # core_r = WL/20.0
  89. # epsilon_Ag = -2.71 + 0.25j
  90. index_Ag = np.sqrt(epsilon_Ag)
  91. # n1 = 1.53413
  92. # n2 = 0.565838 + 7.23262j
  93. nm = 1.0
  94. x = np.ones((1, 1), dtype = np.float64)
  95. x[0, 0] = 2.0*np.pi*core_r/WL
  96. m = np.ones((1, 1), dtype = np.complex128)
  97. m[0, 0] = index_Ag/nm
  98. print "x =", x
  99. print "m =", m
  100. npts = 281
  101. factor=2
  102. scan = np.linspace(-factor*x[0, 0], factor*x[0, 0], npts)
  103. coordX, coordZ = np.meshgrid(scan, scan)
  104. coordX.resize(npts*npts)
  105. coordZ.resize(npts*npts)
  106. coordY = np.zeros(npts*npts, dtype = np.float64)
  107. coord = np.vstack((coordX, coordY, coordZ)).transpose()
  108. #coord = np.vstack((coordY, coordX, coordZ)).transpose()
  109. terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(x, m)
  110. terms, E, H = fieldnlay(x, m, coord)
  111. P = np.array(map(lambda n: np.linalg.norm(np.cross(E[0][n], H[0][n])).real, range(0, len(E[0]))))
  112. Ec = E[0, :, :]
  113. Hc = H[0, :, :]
  114. try:
  115. import matplotlib.pyplot as plt
  116. from matplotlib import cm
  117. from matplotlib.colors import LogNorm
  118. # min_tick = 0.0
  119. # max_tick = 1.0
  120. Eabs_data = np.resize(P, (npts, npts)).T
  121. #Eabs_data = np.resize(Eabs, (npts, npts)).T
  122. #Eabs_data = np.resize(Eangle, (npts, npts)).T
  123. #Eabs_data = np.resize(Habs, (npts, npts)).T
  124. #Eabs_data = np.resize(Hangle, (npts, npts)).T
  125. fig, axs = plt.subplots(1, 2)#, sharey=True, sharex=True)
  126. idxs = np.where(np.abs(coordX) < 1e-10)
  127. #print H[0, idxs][0, :, 1]
  128. axs[0].errorbar(coordZ[idxs]*WL/2.0/np.pi/nm, P[idxs], fmt = 'r', label = 'Poynting vector')
  129. #axs[0].errorbar(coordZ[idxs]*WL/2.0/np.pi/nm, np.linalg.norm(E[0, idxs][0], axis = 1), fmt = 'g', label = 'E')
  130. axs[0].errorbar(coordZ[idxs]*WL/2.0/np.pi/nm, np.linalg.norm(H[0, idxs][0], axis = 1), fmt = 'b', label = 'H')
  131. axs[0].errorbar(coordZ[idxs]*WL/2.0/np.pi/nm, H[0, idxs][0, :, 1].real, fmt = 'k', label = 'H.r')
  132. axs[0].errorbar(coordZ[idxs]*WL/2.0/np.pi/nm, H[0, idxs][0, :, 1].imag, fmt = 'b', label = 'H.i')
  133. axs[0].legend()
  134. #fig.tight_layout()
  135. # Rescale to better show the axes
  136. scale_x = np.linspace(min(coordX)*WL/2.0/np.pi/nm, max(coordX)*WL/2.0/np.pi/nm, npts)
  137. scale_z = np.linspace(min(coordZ)*WL/2.0/np.pi/nm, max(coordZ)*WL/2.0/np.pi/nm, npts)
  138. # Define scale ticks
  139. min_tick = np.amin(Eabs_data)
  140. max_tick = np.amax(Eabs_data)
  141. # scale_ticks = np.power(10.0, np.linspace(np.log10(min_tick), np.log10(max_tick), 6))
  142. scale_ticks = np.linspace(min_tick, max_tick, 11)
  143. # Interpolation can be 'nearest', 'bilinear' or 'bicubic'
  144. #axs[1].set_title('Eabs')
  145. cax = axs[1].imshow(Eabs_data, interpolation = 'nearest', cmap = cm.jet,
  146. origin = 'lower'
  147. #, vmin = min_tick, vmax = max_tick
  148. , extent = (min(scale_x), max(scale_x), min(scale_z), max(scale_z))
  149. #,norm = LogNorm()
  150. )
  151. axs[1].axis("image")
  152. # # Add colorbar
  153. cbar = fig.colorbar(cax, ticks = [a for a in scale_ticks])
  154. cbar.ax.set_yticklabels(['%5.3g' % (a) for a in scale_ticks]) # vertically oriented colorbar
  155. # pos = list(cbar.axs[1].get_position().bounds)
  156. # fig.text(pos[0] - 0.02, 0.925, '|E|/|E$_0$|', fontsize = 14)
  157. plt.xlabel('Z, nm')
  158. plt.ylabel('X, nm')
  159. # This part draws the nanoshell
  160. from matplotlib import patches
  161. s1 = patches.Arc((0, 0), 2.0*core_r, 2.0*core_r, angle=0.0, zorder=2,
  162. theta1=0.0, theta2=360.0, linewidth=1, color='black')
  163. axs[1].add_patch(s1)
  164. from matplotlib.path import Path
  165. #import matplotlib.patches as patches
  166. flow_total = 41
  167. for flow in range(0,flow_total):
  168. flow_x, flow_z = GetFlow(scale_x, scale_z, Ec, Hc,
  169. min(scale_x)+flow*(scale_x[-1]-scale_x[0])/(flow_total-1),
  170. min(scale_z),
  171. npts*6)
  172. verts = np.vstack((flow_z, flow_x)).transpose().tolist()
  173. #codes = [Path.CURVE4]*len(verts)
  174. codes = [Path.LINETO]*len(verts)
  175. codes[0] = Path.MOVETO
  176. path = Path(verts, codes)
  177. patch = patches.PathPatch(path, facecolor='none', lw=1, edgecolor='white')
  178. axs[1].add_patch(patch)
  179. # # Start powerflow lines in the middle of the image
  180. # flow_total = 131
  181. # for flow in range(0,flow_total):
  182. # flow_x, flow_z = GetFlow(scale_x, scale_z, Ec, Hc,
  183. # min(scale_x)+flow*(scale_x[-1]-scale_x[0])/(flow_total-1),
  184. # 15.0, #min(scale_z),
  185. # npts*6)
  186. # verts = np.vstack((flow_z, flow_x)).transpose().tolist()
  187. # #codes = [Path.CURVE4]*len(verts)
  188. # codes = [Path.LINETO]*len(verts)
  189. # codes[0] = Path.MOVETO
  190. # path = Path(verts, codes)
  191. # patch = patches.PathPatch(path, facecolor='none', lw=1, edgecolor='white')
  192. # axs[1].add_patch(patch)
  193. # plt.savefig("SiAgSi.png")
  194. plt.draw()
  195. plt.show()
  196. plt.clf()
  197. plt.close()
  198. finally:
  199. print("Qabs = "+str(Qabs));
  200. #