|
|
@@ -74,9 +74,14 @@ def GetFlow3D(x0, y0, z0, max_length, max_angle, x, m):
|
|
|
Ec, Hc = E[0, 0, :], H[0, 0, :]
|
|
|
S = np.cross(Ec, Hc.conjugate()).real
|
|
|
Snorm_prev = S/np.linalg.norm(S)
|
|
|
+ Sprev = S
|
|
|
length = 0
|
|
|
dpos = step
|
|
|
+ count = 0
|
|
|
while length < max_length:
|
|
|
+ count = count + 1
|
|
|
+ if (count>3000): # Limit length of the absorbed power streamlines
|
|
|
+ break
|
|
|
if step<max_step:
|
|
|
step = step*2.0
|
|
|
r = np.sqrt(flow_x[-1]**2 + flow_y[-1]**2 + flow_z[-1]**2)
|
|
|
@@ -86,7 +91,7 @@ def GetFlow3D(x0, y0, z0, max_length, max_angle, x, m):
|
|
|
dx = dpos*Snorm_prev[0];
|
|
|
dy = dpos*Snorm_prev[1];
|
|
|
dz = dpos*Snorm_prev[2];
|
|
|
- #Test the next position not to turn more than max_angle
|
|
|
+ #Test the next position not to turn\chang size for more than max_angle
|
|
|
coord = np.vstack(([flow_x[-1]+dx], [flow_y[-1]+dy], [flow_z[-1]+dz])).transpose()
|
|
|
terms, E, H = fieldnlay(np.array([x]), np.array([m]), coord)
|
|
|
Ec, Hc = E[0, 0, :], H[0, 0, :]
|
|
|
@@ -99,13 +104,14 @@ def GetFlow3D(x0, y0, z0, max_length, max_angle, x, m):
|
|
|
Hc[i] = 0+0j
|
|
|
S = np.cross(Ec, Hc.conjugate()).real
|
|
|
Snorm = S/np.linalg.norm(S)
|
|
|
- # diff = Snorm-Snorm_prev
|
|
|
- # if np.linalg.norm(diff)<0.05:
|
|
|
- angle = angle_between(Snorm, Snorm_prev)
|
|
|
- if abs(angle) < max_angle:
|
|
|
+ diff = (S-Sprev)/max(np.linalg.norm(S), np.linalg.norm(Sprev))
|
|
|
+ if np.linalg.norm(diff)<max_angle:
|
|
|
+ # angle = angle_between(Snorm, Snorm_prev)
|
|
|
+ # if abs(angle) < max_angle:
|
|
|
break
|
|
|
step = step/2.0
|
|
|
#3. Save result
|
|
|
+ Sprev = S
|
|
|
Snorm_prev = Snorm
|
|
|
dx = dpos*Snorm_prev[0];
|
|
|
dy = dpos*Snorm_prev[1];
|
|
|
@@ -114,7 +120,6 @@ def GetFlow3D(x0, y0, z0, max_length, max_angle, x, m):
|
|
|
flow_x.append(flow_x[-1] + dx)
|
|
|
flow_y.append(flow_y[-1] + dy)
|
|
|
flow_z.append(flow_z[-1] + dz)
|
|
|
-
|
|
|
return np.array(flow_x), np.array(flow_y), np.array(flow_z)
|
|
|
|
|
|
|
|
|
@@ -162,5 +167,144 @@ def GetField(crossplane, npts, factor, x, m):
|
|
|
# for n in range(0, len(E[0])):
|
|
|
# P.append(np.linalg.norm( np.cross(Ec[n], np.conjugate(Hc[n]) ).real/2 ))
|
|
|
return Ec, Hc, P, coordPlot1, coordPlot2
|
|
|
+###############################################################################
|
|
|
+def fieldplot(x,m, WL, comment='', WL_units=' ', crossplane='XZ', field_to_plot='Pabs',npts=101, factor=2.1, flow_total=11):
|
|
|
+ Ec, Hc, P, coordX, coordZ = GetField(crossplane, npts, factor, x, m)
|
|
|
+ Er = np.absolute(Ec)
|
|
|
+ Hr = np.absolute(Hc)
|
|
|
+ try:
|
|
|
+ import matplotlib.pyplot as plt
|
|
|
+ from matplotlib import cm
|
|
|
+ from matplotlib.colors import LogNorm
|
|
|
+ if field_to_plot == 'Pabs':
|
|
|
+ Eabs_data = np.resize(P, (npts, npts)).T
|
|
|
+ label = r'$\operatorname{Re}(E \times H)$'
|
|
|
+ elif field_to_plot == 'Eabs':
|
|
|
+ Eabs = np.sqrt(Er[ :, 0]**2 + Er[ :, 1]**2 + Er[ :, 2]**2)
|
|
|
+ Eabs_data = np.resize(Eabs, (npts, npts)).T
|
|
|
+ label = r'$|E|$'
|
|
|
+ elif field_to_plot == 'Habs':
|
|
|
+ Habs= np.sqrt(Hr[ :, 0]**2 + Hr[ :, 1]**2 + Hr[ :, 2]**2)
|
|
|
+ Eabs_data = np.resize(Habs, (npts, npts)).T
|
|
|
+ label = r'$|H|$'
|
|
|
+ elif field_to_plot == 'angleEx':
|
|
|
+ Eangle = np.angle(Ec[ :, 0])/np.pi*180
|
|
|
+ Eabs_data = np.resize(Eangle, (npts, npts)).T
|
|
|
+ label = r'$arg(E_x)$'
|
|
|
+ elif field_to_plot == 'angleHy':
|
|
|
+ Hangle = np.angle(Hc[ :, 1])/np.pi*180
|
|
|
+ Eabs_data = np.resize(Hangle, (npts, npts)).T
|
|
|
+ label = r'$arg(H_y)$'
|
|
|
+
|
|
|
+ fig, ax = plt.subplots(1,1)
|
|
|
+ # Rescale to better show the axes
|
|
|
+ scale_x = np.linspace(min(coordX)*WL/2.0/np.pi, max(coordX)*WL/2.0/np.pi, npts)
|
|
|
+ scale_z = np.linspace(min(coordZ)*WL/2.0/np.pi, max(coordZ)*WL/2.0/np.pi, npts)
|
|
|
+
|
|
|
+ # Define scale ticks
|
|
|
+ min_tick = np.amin(Eabs_data[~np.isnan(Eabs_data)])
|
|
|
+ max_tick = np.amax(Eabs_data[~np.isnan(Eabs_data)])
|
|
|
+ scale_ticks = np.linspace(min_tick, max_tick, 6)
|
|
|
+
|
|
|
+ # Interpolation can be 'nearest', 'bilinear' or 'bicubic'
|
|
|
+ ax.set_title(label)
|
|
|
+ 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))
|
|
|
+ #,norm = LogNorm()
|
|
|
+ )
|
|
|
+ ax.axis("image")
|
|
|
+
|
|
|
+ # Add colorbar
|
|
|
+ cbar = fig.colorbar(cax, ticks = [a for a in scale_ticks])
|
|
|
+ cbar.ax.set_yticklabels(['%5.3g' % (a) for a in scale_ticks]) # vertically oriented colorbar
|
|
|
+ pos = list(cbar.ax.get_position().bounds)
|
|
|
+ #fig.text(pos[0] - 0.02, 0.925, '|E|/|E$_0$|', fontsize = 14)
|
|
|
+ if crossplane=='XZ':
|
|
|
+ plt.xlabel('Z, '+WL_units)
|
|
|
+ plt.ylabel('X, '+WL_units)
|
|
|
+ elif crossplane=='YZ':
|
|
|
+ plt.xlabel('Z, '+WL_units)
|
|
|
+ plt.ylabel('Y, '+WL_units)
|
|
|
+ elif crossplane=='XY':
|
|
|
+ plt.xlabel('Y, '+WL_units)
|
|
|
+ plt.ylabel('X, '+WL_units)
|
|
|
+
|
|
|
+
|
|
|
+ # # This part draws the nanoshell
|
|
|
+ from matplotlib import patches
|
|
|
+ from matplotlib.path import Path
|
|
|
+ for xx in x:
|
|
|
+ r= xx*WL/2.0/np.pi
|
|
|
+ s1 = patches.Arc((0, 0), 2.0*r, 2.0*r, angle=0.0, zorder=1.8,
|
|
|
+ theta1=0.0, theta2=360.0, linewidth=1, color='black')
|
|
|
+ ax.add_patch(s1)
|
|
|
+ #
|
|
|
+ # 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),
|
|
|
+ # #0.0,
|
|
|
+ # npts*16)
|
|
|
+ # verts = np.vstack((flow_z, flow_x)).transpose().tolist()
|
|
|
+ # #codes = [Path.CURVE4]*len(verts)
|
|
|
+ # codes = [Path.LINETO]*len(verts)
|
|
|
+ # codes[0] = Path.MOVETO
|
|
|
+ # path = Path(verts, codes)
|
|
|
+ # patch = patches.PathPatch(path, facecolor='none', lw=1, edgecolor='yellow')
|
|
|
+ # ax.add_patch(patch)
|
|
|
+ if (crossplane=='XZ' or crossplane=='YZ') and flow_total>0:
|
|
|
+
|
|
|
+ from matplotlib.path import Path
|
|
|
+ scanSP = np.linspace(-factor*x[-1], factor*x[-1], npts)
|
|
|
+ min_SP = -factor*x[-1]
|
|
|
+ step_SP = 2.0*factor*x[-1]/(flow_total-1)
|
|
|
+ x0, y0, z0 = 0, 0, 0
|
|
|
+ max_length=factor*x[-1]*8
|
|
|
+ #max_length=factor*x[-1]*5
|
|
|
+ max_angle = np.pi/160
|
|
|
+ for flow in range(0,flow_total*2+1):
|
|
|
+ #for flow in range(0,flow_total):
|
|
|
+ if crossplane=='XZ':
|
|
|
+ x0 = min_SP*2 + flow*step_SP
|
|
|
+ #x0 = min_SP + flow*step_SP
|
|
|
+ z0 = min_SP
|
|
|
+ #y0 = x[-1]/20
|
|
|
+ elif crossplane=='YZ':
|
|
|
+ y0 = min_SP*2 + flow*step_SP
|
|
|
+ #y0 = min_SP + flow*step_SP
|
|
|
+ z0 = min_SP
|
|
|
+ #x0 = x[-1]/20
|
|
|
+ flow_xSP, flow_ySP, flow_zSP = GetFlow3D(x0, y0, z0, max_length, max_angle, x, m)
|
|
|
+ if crossplane=='XZ':
|
|
|
+ flow_z_plot = flow_zSP*WL/2.0/np.pi
|
|
|
+ flow_f_plot = flow_xSP*WL/2.0/np.pi
|
|
|
+ elif crossplane=='YZ':
|
|
|
+ flow_z_plot = flow_zSP*WL/2.0/np.pi
|
|
|
+ flow_f_plot = flow_ySP*WL/2.0/np.pi
|
|
|
+
|
|
|
+ verts = np.vstack((flow_z_plot, flow_f_plot)).transpose().tolist()
|
|
|
+ codes = [Path.LINETO]*len(verts)
|
|
|
+ codes[0] = Path.MOVETO
|
|
|
+ path = Path(verts, codes)
|
|
|
+ #patch = patches.PathPatch(path, facecolor='none', lw=0.2, edgecolor='white',zorder = 2.7)
|
|
|
+ patch = patches.PathPatch(path, facecolor='none', lw=0.7, edgecolor='white',zorder = 1.9)
|
|
|
+ ax.add_patch(patch)
|
|
|
+ #ax.plot(flow_z_plot, flow_f_plot, 'x',ms=2, mew=0.1, linewidth=0.5, color='k', fillstyle='none')
|
|
|
+
|
|
|
+ plt.savefig(comment+"-R"+str(int(round(x[-1]*WL/2.0/np.pi)))+"-"+crossplane+"-"
|
|
|
+ +field_to_plot+".pdf")
|
|
|
+ plt.draw()
|
|
|
+
|
|
|
+ # plt.show()
|
|
|
+
|
|
|
+ plt.clf()
|
|
|
+ plt.close()
|
|
|
+ finally:
|
|
|
+ terms, Qext, Qsca, Qabs, Qbk, Qpr, g, Albedo, S1, S2 = scattnlay(np.array([x]),
|
|
|
+ np.array([m]))
|
|
|
+ print("Qabs = "+str(Qabs));
|
|
|
+ #
|
|
|
|
|
|
|
