The result is:
When you search "python electric field line" in Google, you would often see this image:
This image is generated using Mayavi module in python. Although Mayavi is powerful module to visualize 3D or 2D image, it works only in python 2.7, the legacy version that will retire in 2020. By using scipy.ode module, you can plot the electric field line without using Mayavi module.
This page shows alternative ways to plot the electric field lines around the point charges arranged in the same places as above figure.
This code is based on the following web sites:
- Visualizing streamlines - Number Crunch
https://www.numbercrunch.de/blog/2013/05/visualizing-streamlines/ - Using Mayavi to visualize electric fields - Elektromagnetisme -
http://elektromagnetisme.no/2010/09/25/using-mayavi-to-visualize-electric-fields/#more-141
In [1]:
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import ode as ode
from matplotlib import cm
from itertools import product
In [2]:
class charge:
def __init__(self, q, pos):
self.q = q
self.pos = pos
def E_point_charge(q, a, x, y):
return q*(x-a[0])/((x-a[0])**2+(y-a[1])**2)**(1.5), \
q*(y-a[1])/((x-a[0])**2+(y-a[1])**2)**(1.5)
def E_total(x, y, charges):
Ex, Ey = 0, 0
for C in charges:
E = E_point_charge(C.q, C.pos, x, y)
Ex = Ex+E[0]
Ey = Ey+E[1]
return [ Ex, Ey ]
def E_dir(t, y, charges):
Ex, Ey = E_total(y[0], y[1], charges)
n = np.sqrt(Ex**2+Ey*Ey)
return [Ex/n, Ey/n]
def V_point_charge(q, a, x, y):
return q/((x-a[0])**2+(y-a[1])**2)**(0.5)
def V_total(x, y, charges):
V = 0
for C in charges:
Vp = V_point_charge(C.q, C.pos, x, y)
V = V+Vp
return V
In [3]:
# charges and positions
charges = [charge( 1, [0.56, 0.56]),
charge(-1, [0.26, 0.76]),
charge( 1, [0.66, 0.16]),
charge(-1, [0.66, 0.86]) ]
# calculate field lines
R = 0.01
# loop over all charges
xs, ys = [], []
for C in charges:
# plot field lines starting in current charge
dt = 0.5*R
if C.q < 0:
# because the electric field lines start only from positive charge,
# skip the process when the current charges is negative.
continue
# loop over field lines starting in different directions
# around current charge
for alpha in np.linspace(0, 2*np.pi*31/32, 32):
r = ode(E_dir)
r.set_integrator('vode')
r.set_f_params(charges)
x = [ C.pos[0] + np.cos(alpha)*R ]
y = [ C.pos[1] + np.sin(alpha)*R ]
r.set_initial_value([x[0], y[0]], 0)
cnt = 0
while r.successful():
Enorm = E_total(r.y[0],r.y[1],charges)
Enorm = (Enorm[0]**2 + Enorm[1]**2)**0.5
a = 0.1
dt2 = R*a*Enorm**(-1)
#if cnt % 1000 == 0:
# print(r.y[0],r.y[1],Enorm,dt2)
#cnt += 1
r.integrate(r.t+dt)
x.append(r.y[0])
y.append(r.y[1])
hit_charge=False
# check if field line ends in some charge
for C2 in charges:
if np.sqrt((r.y[0]-C2.pos[0])**2+(r.y[1]-C2.pos[1])**2)<R:
hit_charge = True
if hit_charge:
break
xs.append(x)
ys.append(y)
In [4]:
# calculate electric potential
vvs = []
xxs = []
yys = []
numcalcv = 300
for xx,yy in product(np.linspace(0,1,numcalcv),np.linspace(0,1,numcalcv)):
xxs.append(xx)
yys.append(yy)
vvs.append(V_total(xx,yy,charges))
xxs = np.array(xxs)
yys = np.array(yys)
vvs = np.array(vvs)
In [5]:
plt.figure(figsize=(5.5, 4.5),facecolor="w")
# plot field line
for x, y in zip(xs,ys):
plt.plot(x, y, color="k",lw=0.8)
# plot point charges
for C in charges:
if C.q>0:
plt.plot(C.pos[0], C.pos[1], 'ro', ms=8*np.sqrt(C.q))
if C.q<0:
plt.plot(C.pos[0], C.pos[1], 'bo', ms=8*np.sqrt(-C.q))
# plot electric potential
clim0,clim1 = -20,20
vvs[np.where(vvs<clim0)] = clim0*0.999999 # to avoid error
vvs[np.where(vvs>clim1)] = clim1*0.999999 # to avoid error
plt.tricontour(xxs,yys,vvs,20,colors="0.3")
plt.tricontourf(xxs,yys,vvs,100,cmap=cm.jet)
cbar = plt.colorbar()
cbar.set_clim(clim0,clim1)
cbar.set_ticks(np.linspace(clim0,clim1,9))
cbar.set_label("Electric Potential")
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.axes().set_aspect('equal','datalim')
plt.savefig('electric_field_line_wo_mayavi.png',dpi=250,bbox_inches="tight",pad_inches=0.02)
plt.show()
