The result is :
This page shows the method to draw electric field line around a point charge adjecent to a grounded sphere using method of image charges.
This code is based on following web sites:
- LECTURE NOTES 6, The Method of Image Charges - UIUC Physics 435 EM Fields & Sources, Prof. Steven Errede -
- Method of image charges, Reflection in a conducting sphere - Wikipedia -
- Visualizing streamlines - Number Crunch -
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import ode as ode
from matplotlib import cm
from itertools import product
from matplotlib.patches import Circle
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
a1 = 1
q1 = 1
r2 = 0.5
b2 = (r2**2)/a1
q2 = -r2/a1*q1
charges=[ charge(q1, [a1, 0]), charge(q2,[b2,0])]
# calculate field lines
x0, x1=-1, 3
y0, y1=-2, 2
R=0.01
# loop over all charges
xs,ys = [],[]
C = charges[0]
# calculate field lines only from point charge outside of the sphere
dt=0.8*R
for alpha in np.linspace(0+2*np.pi/256, 2*np.pi*127/128+2*np.pi/256, 128):
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)
while r.successful():
r.integrate(r.t+dt)
x.append(r.y[0])
y.append(r.y[1])
hit_charge=False
# check if field line left drwaing area or 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 or (not (5*x0<r.y[0] and r.y[0]<5*x1)) or \
(not (5*y0<r.y[1] and r.y[1]<5*y1)):
break
xs.append(x)
ys.append(y)
In [4]:
# calculate electric potential
vvs = []
xxs = []
yys = []
numcalcv = 300
for x,y in product(np.linspace(x0,x1,numcalcv),np.linspace(y0,y1,numcalcv)):
xxs.append(x)
yys.append(y)
vvs.append(V_total(x,y,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.3)
# plot point charges
for C in charges:
if C.q<0:
plt.plot(C.pos[0], C.pos[1], 'bo', ms=8*np.sqrt(-C.q))
if C.q>0:
plt.plot(C.pos[0], C.pos[1], 'ro', ms=8*np.sqrt(C.q))
# plot electric potential
clim0,clim1 = 0,4
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,10,colors="0.3")
plt.tricontourf(xxs,yys,vvs,100,cmap=cm.hot_r)
cbar = plt.colorbar()
cbar.set_clim(clim0,clim1)
cbar.set_ticks(np.linspace(clim0,clim1,11))
cbar.set_label("Electric Potential")
# plot grounded sphere
plt.axes().add_patch(Circle((0, 0), radius=r2, edgecolor="k",facecolor='w', linewidth=1,zorder=10))
# adjustment
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.axes().set_aspect("equal","datalim")
plt.xlim(x0, x1)
plt.ylim(y0, y1)
plt.savefig("field_line_outside_a_sphere.png",bbox_inches="tight",pad_inches=0.02)
plt.show()