d=2
x1=3
V(x,y,z)=1/sqrt(x^2+y^2+(z-d)^2)-1/sqrt(x^2+y^2+(z+d)^2)-1/sqrt((x1-x)^2+y^2+(z-d)^2)+1/sqrt((x1-x)^2+y^2+(z+d)^2)
def LOF(nx, nz, step=.3):
L=[(nx,nz)]
grad=V.gradient()
for i in range(70):
vec=step/grad(nx,0,nz).norm()*grad(nx,0,nz)
nx=nx+vec[0]
nz=nz+vec[2]
T=(nx,nz)
L.append(T)
return line(L)
n=12
cp1=contour_plot(V(x,0,z),(x,-3,6),(z,0,6),contours=(1,2,3,.7,.8,.9,.5,.6,.3,.2,.1), fill=false)
cp2=contour_plot(-V(x,0,z),(x,-3,6),(z,0,6),contours=(1,2,3,.7,.8,.9,.5,.6,.3,.2,.1), fill=false)
LOFs1=[LOF(float(0+.01*cos(t*2*pi/n)),float(2+.01*sin(t*2*pi/n)),-.05) for t in range(n)]
LOFs2=[LOF(float(3+.01*cos(t*2*pi/n)),float(2+.01*sin(t*2*pi/n)),.05) for t in range(n)]
for i in range(n):
cp1=cp1+LOFs1[i]
cp2=cp2+LOFs2[i]
show(cp1+cp2 , xmin=-3, xmax=6, ymin=0, ymax=6,aspect_ratio=1)
