field-nanoshell.py 4.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134
  1. #!/usr/bin/env python
  2. # -*- coding: UTF-8 -*-
  3. #
  4. # Copyright (C) 2009-2015 Ovidio Peña Rodríguez <ovidio@bytesfall.com>
  5. # Copyright (C) 2013-2015 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 the following reference:
  22. # [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by
  23. # a multilayered sphere," Computer Physics Communications,
  24. # vol. 180, Nov. 2009, pp. 2348-2354.
  25. #
  26. # You should have received a copy of the GNU General Public License
  27. # along with this program. If not, see <http://www.gnu.org/licenses/>.
  28. # This test case calculates the electric field in the
  29. # XY plane, for a silver nanoshell embedded in water.
  30. # Refractive index values correspond to the wavelength
  31. # where maximum of the surface plasmon resonance (and,
  32. # hence, of electric field) is expected.
  33. from scattnlay import fieldnlay
  34. import numpy as np
  35. import time
  36. n1 = 1.53413
  37. n2 = 0.565838 + 7.23262j
  38. nm = 1.3205
  39. r1 = 0.05
  40. r2 = 0.06
  41. x = 2.0*np.pi*nm*np.array([[r1, r2]], dtype = np.float64)/1.064
  42. m = np.array((n1/nm, n2/nm), dtype = np.complex128)
  43. print("x =", x)
  44. print("m =", m)
  45. npts = 501
  46. scan = np.linspace(-3.0*x[0, 0], 3.0*x[0, 0], npts)
  47. coordX, coordY = np.meshgrid(scan, scan)
  48. coordX.resize(npts*npts)
  49. coordY.resize(npts*npts)
  50. coordZ = np.zeros(npts*npts, dtype = np.float64)
  51. start_time = time.time()
  52. terms, E, H = fieldnlay(x, m, coordX, coordY, coordZ)
  53. elapsed_time = time.time() - start_time
  54. print("Time: ", elapsed_time)
  55. Er = np.absolute(E)
  56. # |E|/|Eo|
  57. Eh = np.sqrt(Er[0, :, 0]**2 + Er[0, :, 1]**2 + Er[0, :, 2]**2)
  58. result = np.vstack((coordX, coordY, coordZ, Eh)).transpose()
  59. try:
  60. import matplotlib.pyplot as plt
  61. from matplotlib import cm
  62. from matplotlib.colors import LogNorm
  63. min_tick = 0.0
  64. max_tick = 5.0
  65. edata = np.resize(Eh, (npts, npts))
  66. fig = plt.figure()
  67. ax = fig.add_subplot(111)
  68. # Rescale to better show the axes
  69. scale_x = 1000.0*np.linspace(min(coordX)*1.064/2.0/np.pi/nm, max(coordX)*1.064/2.0/np.pi/nm, npts)
  70. scale_y = 1000.0*np.linspace(min(coordY)*1.064/2.0/np.pi/nm, max(coordY)*1.064/2.0/np.pi/nm, npts)
  71. # Define scale ticks
  72. min_tick = min(min_tick, np.amin(edata))
  73. max_tick = max(max_tick, np.amax(edata))
  74. #scale_ticks = np.power(10.0, np.linspace(np.log10(min_tick), np.log10(max_tick), 6))
  75. scale_ticks = np.linspace(min_tick, max_tick, 6)
  76. # Interpolation can be 'nearest', 'bilinear' or 'bicubic'
  77. cax = ax.imshow(edata, interpolation = 'nearest', cmap = cm.jet,
  78. origin = 'lower', vmin = min_tick, vmax = max_tick,
  79. extent = (min(scale_x), max(scale_x), min(scale_y), max(scale_y)))
  80. # Add colorbar
  81. cbar = fig.colorbar(cax, ticks = [a for a in scale_ticks])
  82. cbar.ax.set_yticklabels(['%4.2g' % (a) for a in scale_ticks]) # vertically oriented colorbar
  83. pos = list(cbar.ax.get_position().bounds)
  84. fig.text(pos[0] - 0.02, 0.925, '|E|/|E$_0$|', fontsize = 14)
  85. plt.xlabel('X ( nm )')
  86. plt.ylabel('Y ( nm )')
  87. # This part draws the nanoshell
  88. from matplotlib import patches
  89. s1 = patches.Arc((0, 0), 2000.0*r1, 2000.0*r1, angle=0.0, zorder=2,
  90. theta1=0.0, theta2=360.0, linewidth=1, color='#00fa9a')
  91. ax.add_patch(s1)
  92. s2 = patches.Arc((0, 0), 2000.0*r2, 2000.0*r2, angle=0.0, zorder=2,
  93. theta1=0.0, theta2=360.0, linewidth=1, color='#00fa9a')
  94. ax.add_patch(s2)
  95. # End of drawing
  96. plt.draw()
  97. plt.show()
  98. plt.clf()
  99. plt.close()
  100. finally:
  101. np.savetxt("field-nanoshell.txt", result, fmt = "%.5f")
  102. print(result)