# paste this code at the end of VectorFieldPlot 1.10
# https://commons.wikimedia.org/wiki/User:Geek3/VectorFieldPlot
doc = FieldplotDocument('VFPt_four_charges_potential',
    width=600, height=600, commons=True)
charges = [[-1,0,-1], [0,1,1], [1,0,-1], [0,-1,1]]
field = Field({'monopoles':charges})
doc.draw_charges(field)

# draw potential in background
p_array = sc.array([c[:2] for c in charges])
q_array = sc.array([c[2] for c in charges])
def potential(xy):
    return sc.dot(q_array, 1. / sc.linalg.norm(xy - p_array, axis=1))

from matplotlib import colors
# colormap from aqua through yellow to fuchsia
cmap = colors.ListedColormap([sc.clip((2*x, 2*(1-x), 4*(x-0.5)**2), 0, 1)
    for x in sc.linspace(0., 1., 2048)])

U0 = potential(sc.array([0, 1.15]))
doc.draw_scalar_field(func=potential, cmap=cmap, vmin=-U0, vmax=U0)

n = 8
Fsum = ig.quad(lambda t: sc.dot(field.F([.5*sin(t), 1-0.5*cos(t)]),
    [sin(t), -cos(t)]), 0, pi)[0]
for i in range(n):
    a = (i + 0.5) / n
    w = op.brentq(lambda w: ig.quad(lambda t: sc.dot(field.F([0.5*sin(t),
    	1-0.5*cos(t)]), [sin(t), -cos(t)]), 0, w)[0] - a * Fsum, 0, pi)
    for w1 in range(4):
    	line = FieldLine(field, rot([.5*sin(w), 1.0-0.5*cos(w)], w1 * pi/2),
    	    directions='both')
    	doc.draw_line(line, arrows_style={'dist':1.5})
doc.write()