# paste this code at the end of VectorFieldPlot 3.3
# https://commons.wikimedia.org/wiki/User:Geek3/VectorFieldPlot
u = 100.0
doc = FieldplotDocument('VFPt_metal_balls_plusminus_potential+contour',
    commons=True, width=800, height=600, center=[400, 300], unit=u)

# define two spheres with position, radius and charge
s1 = {'p':array([-1.5, 0.]), 'r':1.0}
s2 = {'p':array([1.5, 0.]), 'r':1.0}
spheres = s1, s2

# make charge proportional to capacitance, which is proportional to radius.
s1['Q'] = s1['r']
s2['Q'] = -s2['r']
d = vabs(s2['p'] - s1['p'])
v12 = vnorm(s2['p'] - s1['p'])

# compute series of charges https://dx.doi.org/10.2174/1874183500902010032
charges = [{'x':s1['p'][0], 'y':s1['p'][1], 'Q':s1['Q']},
           {'x':s2['p'][0], 'y':s2['p'][1], 'Q':s2['Q']}]
r1 = r2 = 0.
q1, q2 = s1['Q'], s2['Q']
q0 = max(fabs(q1), fabs(q2))
for i in range(10):
    q1, q2 = -s1['r'] * q2 / (d - r2), -s2['r'] * q1 / (d - r1), 
    r1, r2 = s1['r']**2 / (d - r2), s2['r']**2 / (d - r1)
    p1, p2 = s1['p'] + r1 * v12, s2['p'] - r2 * v12
    charges.append({'x':p1[0], 'y':p1[1], 'Q':q1})
    charges.append({'x':p2[0], 'y':p2[1], 'Q':q2})
    if max(fabs(q1), fabs(q2)) < 1e-3 * q0:
        break

field = Field([ ['monopole', c] for c in charges ])

U1 = field.V(s1['p'] - s1['r'] * array([1., 0.]))
U2 = field.V(s2['p'] + s2['r'] * array([1., 0.]))
doc.draw_scalar_field(func=field.V, cmap=doc.cmap_AqYlFs, vmin=U2, vmax=U1)
doc.draw_contours(func=field.V, linewidth=1, linecolor='#111111',
    levels=sc.linspace(U2, U1, 11)[1:-1], attributes={'opacity':'0.7'})

# draw symbols
for c in charges:
    doc.draw_charges(Field([ ['monopole', c] ]), scale=0.6*sqrt(fabs(c['Q'])))

gradr = doc.draw_object('linearGradient', {'id':'rod_shade', 'x1':0, 'x2':0,
    'y1':0, 'y2':1, 'gradientUnits':'objectBoundingBox'}, group=doc.defs)
for col, of in (('#666', 0), ('#ddd', 0.6), ('#fff', 0.7), ('#ccc', 0.75),
    ('#888', 1)):
    doc.draw_object('stop', {'offset':of, 'stop-color':col}, group=gradr)
gradb = doc.draw_object('radialGradient', {'id':'metal_spot', 'cx':'0.53',
    'cy':'0.54', 'r':'0.55', 'fx':'0.65', 'fy':'0.7',
    'gradientUnits':'objectBoundingBox'}, group=doc.defs)
for col, of in (('#fff', 0), ('#e7e7e7', 0.15), ('#ddd', 0.25),
    ('#aaa', 0.7), ('#888', 0.9), ('#666', 1)):
    doc.draw_object('stop', {'offset':of, 'stop-color':col}, group=gradb)

for si, s in enumerate(spheres):
    ball = doc.draw_object('g', {'id':'metal_ball{:}'.format(si+1),
        'transform':'translate({:.3f},{:.3f})'.format(*(s['p'])),
        'style':'fill:none; stroke:#000;stroke-linecap:square', 'opacity':1})
    
    # draw metal balls
    doc.draw_object('circle', {'cx':0, 'cy':0, 'r':s['r'],
        'style':'fill:url(#metal_spot); stroke-width:0.02'}, group=ball)
    s['cgroup'] = doc.draw_object('g', {'style':'stroke-width:0.02'}, group=ball)

# find start positions of field lines
def start(t):
    phi = 2. * pi * t
    return array(s1['p']) + 1.5 * array([cos(phi), sin(phi)])

# draw the field lines
nlines = 24
p0_list = Startpath(field, start).npoints(nlines)
for iline in range(nlines):
    line = FieldLine(field, p0_list[iline], directions='both', maxr=1e4)
    
    # draw little charge signs near the surface
    path_minus = 'M {0:.5f},0 h {1:.5f}'.format(-2 / u, 4 / u)
    path_plus = 'M {0:.5f},0 h {1:.5f} M 0,{0:.5f} v {1:.5f}'.format(-2 / u, 4 / u)
    for si, s in enumerate(spheres):
        # check if fieldline ends inside the sphere
        for ci in range(2):
            if vabs(line.get_position(ci) - s['p']) < s['r']:
                # find the point where the field line cuts the surface
                t = optimize.brentq(lambda t: vabs(line.get_position(t)
                    - s['p']) - s['r'], 0., 1.)
                pr = line.get_position(t) - s['p']
                cpos = 0.9 * s['r'] * vnorm(pr)
                doc.draw_object('path', {'stroke':'black', 'd':
                    [path_plus, path_minus][ci],
                    'transform':'translate({:.5f},{:.5f})'.format(
                        round(u*cpos[0])/u, round(u*cpos[1])/u)},
                        group=s['cgroup'])
    
    pot = []
    if iline != 6 and iline != nlines - 1 - 6:
        pot.append(0.)
    if iline >= 6 and iline < nlines - 6:
        pot += [0.27 * U1 + 0.73 * U2, 0.73 * U1 + 0.27 * U2]
    
    doc.draw_line(line, linewidth=2.4, arrows_style={'at_potentials':pot})

doc.write()