1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
| from matplotlib.tri import ( Triangulation, UniformTriRefiner, CubicTriInterpolator) import matplotlib.pyplot as plt import matplotlib.cm as cm import numpy as np
def dipole_potential(x, y): """ The electric dipole potential V """ r_sq = x**2 + y**2 theta = np.arctan2(y, x) z = np.cos(theta)/r_sq return (np.max(z) - z) / (np.max(z) - np.min(z))
n_angles = 30 n_radii = 10 min_radius = 0.2 radii = np.linspace(min_radius, 0.95, n_radii)
angles = np.linspace(0, 2 * np.pi, n_angles, endpoint=False) angles = np.repeat(angles[..., np.newaxis], n_radii, axis=1) angles[:, 1::2] += np.pi / n_angles
x = (radii*np.cos(angles)).flatten() y = (radii*np.sin(angles)).flatten() V = dipole_potential(x, y)
triang = Triangulation(x, y)
triang.set_mask(np.hypot(x[triang.triangles].mean(axis=1), y[triang.triangles].mean(axis=1)) < min_radius)
refiner = UniformTriRefiner(triang) tri_refi, z_test_refi = refiner.refine_field(V, subdiv=3)
tci = CubicTriInterpolator(triang, -V)
(Ex, Ey) = tci.gradient(triang.x, triang.y) E_norm = np.sqrt(Ex**2 + Ey**2)
fig, ax = plt.subplots() ax.set_aspect('equal')
ax.use_sticky_edges = False ax.margins(0.07)
ax.triplot(triang, color='0.8')
levels = np.arange(0., 1., 0.01) cmap = cm.get_cmap(name='hot', lut=None) ax.tricontour(tri_refi, z_test_refi, levels=levels, cmap=cmap, linewidths=[2.0, 1.0, 1.0, 1.0])
ax.quiver(triang.x, triang.y, Ex/E_norm, Ey/E_norm, units='xy', scale=10., zorder=3, color='blue', width=0.007, headwidth=3., headlength=4.)
ax.set_title('Gradient plot: an electrical dipole') plt.show()
|