123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237 |
- import scattnlay
- from scattnlay import fieldnlay
- from scattnlay import scattnlay
- import numpy as np
- import cmath
- def get_index(array,value):
- idx = (np.abs(array-value)).argmin()
- return idx
- def GetFlow(scale_x, scale_z, Ec, Hc, a, b, nmax):
-
- flow_x = [a]
- flow_z = [b]
- for n in range(0, nmax):
-
-
- x_pos = flow_x[-1]
- z_pos = flow_z[-1]
- x_idx = get_index(scale_x, x_pos)
- z_idx = get_index(scale_z, z_pos)
- S=np.cross(Ec[npts*z_idx+x_idx], np.conjugate(Hc[npts*z_idx+x_idx]) ).real
- Snorm=S/np.linalg.norm(S)
-
- dpos = abs(scale_z[0]-scale_z[1])/4.0
- dx = dpos*Snorm[0]
- dz = dpos*Snorm[2]
- x_pos = x_pos+dx
- z_pos = z_pos+dz
-
- flow_x.append(x_pos)
- flow_z.append(z_pos)
- return flow_x, flow_z
- epsilon_Si = 13.64 + 0.047j
- epsilon_Ag = -28.05 + 1.525j
- index_Si = np.sqrt(epsilon_Si)
- index_Ag = np.sqrt(epsilon_Ag)
- print(index_Si)
- print(index_Ag)
- WL=800
- core_width = 17.74
- inner_width = 23.31
- outer_width = 22.95
- core_r = core_width
- inner_r = core_r+inner_width
- outer_r = inner_r+outer_width
- nm = 1.0
- x = np.ones((1, 3), dtype = np.float64)
- x[0, 0] = 2.0*np.pi*core_r/WL
- x[0, 1] = 2.0*np.pi*inner_r/WL
- x[0, 2] = 2.0*np.pi*outer_r/WL
- m = np.ones((1, 3), dtype = np.complex128)
- m[0, 0] = index_Si/nm
- m[0, 1] = index_Ag/nm
- m[0, 2] = index_Si/nm
- print "x =", x
- print "m =", m
- npts = 241
- factor=2.2
- scan = np.linspace(-factor*x[0, 2], factor*x[0, 2], npts)
- coordX, coordZ = np.meshgrid(scan, scan)
- coordX.resize(npts*npts)
- coordZ.resize(npts*npts)
- coordY = np.zeros(npts*npts, dtype = np.float64)
- coord = np.vstack((coordX, coordY, coordZ)).transpose()
- terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(x, m)
- terms, E, H = fieldnlay(x, m, coord)
- Er = np.absolute(E)
- Hr = np.absolute(H)
- Eabs = np.sqrt(Er[0, :, 0]**2 + Er[0, :, 1]**2 + Er[0, :, 2]**2)
- Ec = E[0, :, :]
- Hc = H[0, :, :]
- Eangle = np.angle(E[0, :, 0])/np.pi*180
- P=[]
- for n in range(0, len(E[0])):
- P.append(np.linalg.norm( np.cross(Ec[n], np.conjugate(Hc[n]) ).real/2 ))
- Habs= np.sqrt(Hr[0, :, 0]**2 + Hr[0, :, 1]**2 + Hr[0, :, 2]**2)
- Hangle = np.angle(H[0, :, 1])/np.pi*180
- try:
- import matplotlib.pyplot as plt
- from matplotlib import cm
- from matplotlib.colors import LogNorm
- min_tick = 0.0
- max_tick = 1.0
-
- Eabs_data = np.resize(Eabs, (npts, npts)).T
-
-
-
- fig, ax = plt.subplots(1,1)
-
-
- scale_x = np.linspace(min(coordX)*WL/2.0/np.pi/nm, max(coordX)*WL/2.0/np.pi/nm, npts)
- scale_z = np.linspace(min(coordZ)*WL/2.0/np.pi/nm, max(coordZ)*WL/2.0/np.pi/nm, npts)
-
- min_tick = min(min_tick, np.amin(Eabs_data))
- max_tick = max(max_tick, np.amax(Eabs_data))
-
-
- scale_ticks = np.linspace(min_tick, max_tick, 11)
-
- ax.set_title('Eabs')
- cax = ax.imshow(Eabs_data, interpolation = 'nearest', cmap = cm.jet,
- origin = 'lower'
- , vmin = min_tick, vmax = max_tick
- , extent = (min(scale_x), max(scale_x), min(scale_z), max(scale_z))
-
- )
- ax.axis("image")
-
- cbar = fig.colorbar(cax, ticks = [a for a in scale_ticks])
- cbar.ax.set_yticklabels(['%5.3g' % (a) for a in scale_ticks])
- pos = list(cbar.ax.get_position().bounds)
- fig.text(pos[0] - 0.02, 0.925, '|E|/|E$_0$|', fontsize = 14)
- plt.xlabel('Z, nm')
- plt.ylabel('X, nm')
-
- from matplotlib import patches
- s1 = patches.Arc((0, 0), 2.0*core_r, 2.0*core_r, angle=0.0, zorder=2,
- theta1=0.0, theta2=360.0, linewidth=1, color='black')
- s2 = patches.Arc((0, 0), 2.0*inner_r, 2.0*inner_r, angle=0.0, zorder=2,
- theta1=0.0, theta2=360.0, linewidth=1, color='black')
- s3 = patches.Arc((0, 0), 2.0*outer_r, 2.0*outer_r, angle=0.0, zorder=2,
- theta1=0.0, theta2=360.0, linewidth=1, color='black')
- ax.add_patch(s1)
- ax.add_patch(s2)
- ax.add_patch(s3)
- from matplotlib.path import Path
-
- flow_total = 21
- for flow in range(0,flow_total):
- flow_x, flow_z = GetFlow(scale_x, scale_z, Ec, Hc,
- min(scale_x)+flow*(scale_x[-1]-scale_x[0])/(flow_total-1),
- min(scale_z), npts*12)
- verts = np.vstack((flow_z, flow_x)).transpose().tolist()
- codes = [Path.CURVE4]*len(verts)
-
- codes[0] = Path.MOVETO
- path = Path(verts, codes)
- patch = patches.PathPatch(path, facecolor='none', lw=1, edgecolor='white')
- ax.add_patch(patch)
-
- plt.savefig("SiAgSi-flow.png")
- plt.draw()
- plt.show()
- plt.clf()
- plt.close()
- finally:
- print("Qabs = "+str(Qabs));
|