fieldplot.py 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310
  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. # Several functions to plot field and streamlines (power flow lines).
  28. import scattnlay
  29. from scattnlay import fieldnlay
  30. from scattnlay import scattnlay
  31. import numpy as np
  32. import cmath
  33. def unit_vector(vector):
  34. """ Returns the unit vector of the vector. """
  35. return vector / np.linalg.norm(vector)
  36. def angle_between(v1, v2):
  37. """ Returns the angle in radians between vectors 'v1' and 'v2'::
  38. >>> angle_between((1, 0, 0), (0, 1, 0))
  39. 1.5707963267948966
  40. >>> angle_between((1, 0, 0), (1, 0, 0))
  41. 0.0
  42. >>> angle_between((1, 0, 0), (-1, 0, 0))
  43. 3.141592653589793
  44. """
  45. v1_u = unit_vector(v1)
  46. v2_u = unit_vector(v2)
  47. angle = np.arccos(np.dot(v1_u, v2_u))
  48. if np.isnan(angle):
  49. if (v1_u == v2_u).all():
  50. return 0.0
  51. else:
  52. return np.pi
  53. return angle
  54. ###############################################################################
  55. def GetFlow3D(x0, y0, z0, max_length, max_angle, x, m):
  56. # Initial position
  57. flow_x = [x0]
  58. flow_y = [y0]
  59. flow_z = [z0]
  60. max_step = x[-1]/3
  61. min_step = x[0]/2000
  62. # max_step = min_step
  63. step = min_step*2.0
  64. if max_step < min_step:
  65. max_step = min_step
  66. coord = np.vstack(([flow_x[-1]], [flow_y[-1]], [flow_z[-1]])).transpose()
  67. terms, E, H = fieldnlay(np.array([x]), np.array([m]), coord)
  68. Ec, Hc = E[0, 0, :], H[0, 0, :]
  69. S = np.cross(Ec, Hc.conjugate()).real
  70. Snorm_prev = S/np.linalg.norm(S)
  71. Sprev = S
  72. length = 0
  73. dpos = step
  74. count = 0
  75. while length < max_length:
  76. count = count + 1
  77. if (count>3000): # Limit length of the absorbed power streamlines
  78. break
  79. if step<max_step:
  80. step = step*2.0
  81. r = np.sqrt(flow_x[-1]**2 + flow_y[-1]**2 + flow_z[-1]**2)
  82. while step > min_step:
  83. #Evaluate displacement from previous poynting vector
  84. dpos = step
  85. dx = dpos*Snorm_prev[0];
  86. dy = dpos*Snorm_prev[1];
  87. dz = dpos*Snorm_prev[2];
  88. #Test the next position not to turn\chang size for more than max_angle
  89. coord = np.vstack(([flow_x[-1]+dx], [flow_y[-1]+dy], [flow_z[-1]+dz])).transpose()
  90. terms, E, H = fieldnlay(np.array([x]), np.array([m]), coord)
  91. Ec, Hc = E[0, 0, :], H[0, 0, :]
  92. Eth = max(np.absolute(Ec))/1e10
  93. Hth = max(np.absolute(Hc))/1e10
  94. for i in xrange(0,len(Ec)):
  95. if abs(Ec[i]) < Eth:
  96. Ec[i] = 0+0j
  97. if abs(Hc[i]) < Hth:
  98. Hc[i] = 0+0j
  99. S = np.cross(Ec, Hc.conjugate()).real
  100. Snorm = S/np.linalg.norm(S)
  101. diff = (S-Sprev)/max(np.linalg.norm(S), np.linalg.norm(Sprev))
  102. if np.linalg.norm(diff)<max_angle:
  103. # angle = angle_between(Snorm, Snorm_prev)
  104. # if abs(angle) < max_angle:
  105. break
  106. step = step/2.0
  107. #3. Save result
  108. Sprev = S
  109. Snorm_prev = Snorm
  110. dx = dpos*Snorm_prev[0];
  111. dy = dpos*Snorm_prev[1];
  112. dz = dpos*Snorm_prev[2];
  113. length = length + step
  114. flow_x.append(flow_x[-1] + dx)
  115. flow_y.append(flow_y[-1] + dy)
  116. flow_z.append(flow_z[-1] + dz)
  117. return np.array(flow_x), np.array(flow_y), np.array(flow_z)
  118. ###############################################################################
  119. def GetField(crossplane, npts, factor, x, m):
  120. """
  121. crossplane: XZ, YZ, XY
  122. npts: number of point in each direction
  123. factor: ratio of plotting size to outer size of the particle
  124. x: size parameters for particle layers
  125. m: relative index values for particle layers
  126. """
  127. scan = np.linspace(-factor*x[-1], factor*x[-1], npts)
  128. zero = np.zeros(npts*npts, dtype = np.float64)
  129. if crossplane=='XZ':
  130. coordX, coordZ = np.meshgrid(scan, scan)
  131. coordX.resize(npts*npts)
  132. coordZ.resize(npts*npts)
  133. coordY = zero
  134. coordPlot1 = coordX
  135. coordPlot2 = coordZ
  136. elif crossplane=='YZ':
  137. coordY, coordZ = np.meshgrid(scan, scan)
  138. coordY.resize(npts*npts)
  139. coordZ.resize(npts*npts)
  140. coordX = zero
  141. coordPlot1 = coordY
  142. coordPlot2 = coordZ
  143. elif crossplane=='XY':
  144. coordX, coordY = np.meshgrid(scan, scan)
  145. coordX.resize(npts*npts)
  146. coordY.resize(npts*npts)
  147. coordZ = zero
  148. coordPlot1 = coordY
  149. coordPlot2 = coordX
  150. coord = np.vstack((coordX, coordY, coordZ)).transpose()
  151. terms, E, H = fieldnlay(np.array([x]), np.array([m]), coord)
  152. Ec = E[0, :, :]
  153. Hc = H[0, :, :]
  154. P=[]
  155. P = np.array(map(lambda n: np.linalg.norm(np.cross(Ec[n], Hc[n])).real, range(0, len(E[0]))))
  156. # for n in range(0, len(E[0])):
  157. # P.append(np.linalg.norm( np.cross(Ec[n], np.conjugate(Hc[n]) ).real/2 ))
  158. return Ec, Hc, P, coordPlot1, coordPlot2
  159. ###############################################################################
  160. def fieldplot(x,m, WL, comment='', WL_units=' ', crossplane='XZ', field_to_plot='Pabs',npts=101, factor=2.1, flow_total=11):
  161. Ec, Hc, P, coordX, coordZ = GetField(crossplane, npts, factor, x, m)
  162. Er = np.absolute(Ec)
  163. Hr = np.absolute(Hc)
  164. try:
  165. import matplotlib.pyplot as plt
  166. from matplotlib import cm
  167. from matplotlib.colors import LogNorm
  168. if field_to_plot == 'Pabs':
  169. Eabs_data = np.resize(P, (npts, npts)).T
  170. label = r'$\operatorname{Re}(E \times H)$'
  171. elif field_to_plot == 'Eabs':
  172. Eabs = np.sqrt(Er[ :, 0]**2 + Er[ :, 1]**2 + Er[ :, 2]**2)
  173. Eabs_data = np.resize(Eabs, (npts, npts)).T
  174. label = r'$|E|$'
  175. elif field_to_plot == 'Habs':
  176. Habs= np.sqrt(Hr[ :, 0]**2 + Hr[ :, 1]**2 + Hr[ :, 2]**2)
  177. Eabs_data = np.resize(Habs, (npts, npts)).T
  178. label = r'$|H|$'
  179. elif field_to_plot == 'angleEx':
  180. Eangle = np.angle(Ec[ :, 0])/np.pi*180
  181. Eabs_data = np.resize(Eangle, (npts, npts)).T
  182. label = r'$arg(E_x)$'
  183. elif field_to_plot == 'angleHy':
  184. Hangle = np.angle(Hc[ :, 1])/np.pi*180
  185. Eabs_data = np.resize(Hangle, (npts, npts)).T
  186. label = r'$arg(H_y)$'
  187. fig, ax = plt.subplots(1,1)
  188. # Rescale to better show the axes
  189. scale_x = np.linspace(min(coordX)*WL/2.0/np.pi, max(coordX)*WL/2.0/np.pi, npts)
  190. scale_z = np.linspace(min(coordZ)*WL/2.0/np.pi, max(coordZ)*WL/2.0/np.pi, npts)
  191. # Define scale ticks
  192. min_tick = np.amin(Eabs_data[~np.isnan(Eabs_data)])
  193. max_tick = np.amax(Eabs_data[~np.isnan(Eabs_data)])
  194. scale_ticks = np.linspace(min_tick, max_tick, 6)
  195. # Interpolation can be 'nearest', 'bilinear' or 'bicubic'
  196. ax.set_title(label)
  197. cax = ax.imshow(Eabs_data, interpolation = 'nearest', cmap = cm.jet,
  198. origin = 'lower'
  199. , vmin = min_tick, vmax = max_tick
  200. , extent = (min(scale_x), max(scale_x), min(scale_z), max(scale_z))
  201. #,norm = LogNorm()
  202. )
  203. ax.axis("image")
  204. # Add colorbar
  205. cbar = fig.colorbar(cax, ticks = [a for a in scale_ticks])
  206. cbar.ax.set_yticklabels(['%5.3g' % (a) for a in scale_ticks]) # vertically oriented colorbar
  207. pos = list(cbar.ax.get_position().bounds)
  208. #fig.text(pos[0] - 0.02, 0.925, '|E|/|E$_0$|', fontsize = 14)
  209. if crossplane=='XZ':
  210. plt.xlabel('Z, '+WL_units)
  211. plt.ylabel('X, '+WL_units)
  212. elif crossplane=='YZ':
  213. plt.xlabel('Z, '+WL_units)
  214. plt.ylabel('Y, '+WL_units)
  215. elif crossplane=='XY':
  216. plt.xlabel('Y, '+WL_units)
  217. plt.ylabel('X, '+WL_units)
  218. # # This part draws the nanoshell
  219. from matplotlib import patches
  220. from matplotlib.path import Path
  221. for xx in x:
  222. r= xx*WL/2.0/np.pi
  223. s1 = patches.Arc((0, 0), 2.0*r, 2.0*r, angle=0.0, zorder=1.8,
  224. theta1=0.0, theta2=360.0, linewidth=1, color='black')
  225. ax.add_patch(s1)
  226. #
  227. # for flow in range(0,flow_total):
  228. # flow_x, flow_z = GetFlow(scale_x, scale_z, Ec, Hc,
  229. # min(scale_x)+flow*(scale_x[-1]-scale_x[0])/(flow_total-1),
  230. # min(scale_z),
  231. # #0.0,
  232. # npts*16)
  233. # verts = np.vstack((flow_z, flow_x)).transpose().tolist()
  234. # #codes = [Path.CURVE4]*len(verts)
  235. # codes = [Path.LINETO]*len(verts)
  236. # codes[0] = Path.MOVETO
  237. # path = Path(verts, codes)
  238. # patch = patches.PathPatch(path, facecolor='none', lw=1, edgecolor='yellow')
  239. # ax.add_patch(patch)
  240. if (crossplane=='XZ' or crossplane=='YZ') and flow_total>0:
  241. from matplotlib.path import Path
  242. scanSP = np.linspace(-factor*x[-1], factor*x[-1], npts)
  243. min_SP = -factor*x[-1]
  244. step_SP = 2.0*factor*x[-1]/(flow_total-1)
  245. x0, y0, z0 = 0, 0, 0
  246. max_length=factor*x[-1]*8
  247. #max_length=factor*x[-1]*5
  248. max_angle = np.pi/160
  249. for flow in range(0,flow_total*2+1):
  250. #for flow in range(0,flow_total):
  251. if crossplane=='XZ':
  252. x0 = min_SP*2 + flow*step_SP
  253. #x0 = min_SP + flow*step_SP
  254. z0 = min_SP
  255. #y0 = x[-1]/20
  256. elif crossplane=='YZ':
  257. y0 = min_SP*2 + flow*step_SP
  258. #y0 = min_SP + flow*step_SP
  259. z0 = min_SP
  260. #x0 = x[-1]/20
  261. flow_xSP, flow_ySP, flow_zSP = GetFlow3D(x0, y0, z0, max_length, max_angle, x, m)
  262. if crossplane=='XZ':
  263. flow_z_plot = flow_zSP*WL/2.0/np.pi
  264. flow_f_plot = flow_xSP*WL/2.0/np.pi
  265. elif crossplane=='YZ':
  266. flow_z_plot = flow_zSP*WL/2.0/np.pi
  267. flow_f_plot = flow_ySP*WL/2.0/np.pi
  268. verts = np.vstack((flow_z_plot, flow_f_plot)).transpose().tolist()
  269. codes = [Path.LINETO]*len(verts)
  270. codes[0] = Path.MOVETO
  271. path = Path(verts, codes)
  272. #patch = patches.PathPatch(path, facecolor='none', lw=0.2, edgecolor='white',zorder = 2.7)
  273. patch = patches.PathPatch(path, facecolor='none', lw=0.7, edgecolor='white',zorder = 1.9)
  274. ax.add_patch(patch)
  275. #ax.plot(flow_z_plot, flow_f_plot, 'x',ms=2, mew=0.1, linewidth=0.5, color='k', fillstyle='none')
  276. plt.savefig(comment+"-R"+str(int(round(x[-1]*WL/2.0/np.pi)))+"-"+crossplane+"-"
  277. +field_to_plot+".pdf")
  278. plt.draw()
  279. # plt.show()
  280. plt.clf()
  281. plt.close()
  282. finally:
  283. terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(np.array([x]),
  284. np.array([m]))
  285. print("Qabs = "+str(Qabs));
  286. #