# paste this code at the end of VectorFieldPlot 1.10 # https://commons.wikimedia.org/wiki/User:Geek3/VectorFieldPlot u = 100.0 doc = FieldplotDocument('VFPt_metal_balls_largesmall_transparent',     commons=True, width=800, height=600, center=[400, 300], unit=u)  # define two spheres with position, radius and charge s1 = {'p':sc.array([-1.0, 0.]), 'r':1.5} s2 = {'p':sc.array([2.0, 0.]), 'r':0.5}  # 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 = (s2['p'] - s1['p']) / d  # compute series of charges https://dx.doi.org/10.2174/1874183500902010032 charges = [[s1['p'][0], s1['p'][1], s1['q']], [s2['p'][0], s2['p'][1], 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([p1[0], p1[1], q1])     charges.append([p2[0], p2[1], q2])     if max(fabs(q1), fabs(q2)) < 1e-3 * q0:         break  field = Field({'monopoles':charges})  # draw symbols for c in charges:     doc.draw_charges(Field({'monopoles':[c]}), scale=0.6*sqrt(fabs(c[2])))  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)  ball_charges = [] for ib in range(2):     ball = doc.draw_object('g', {'id':'metal_ball{:}'.format(ib+1),         'transform':'translate({:.3f},{:.3f})'.format(*([s1, s2][ib]['p'])),         'style':'fill:none; stroke:#000;stroke-linecap:square', 'opacity':0.5})          # draw rods     if ib == 0:         x1, x2 = -4.1 - s1['p'][0], -0.9 * s1['r']     else:         x1, x2 = 0.9 * s2['r'], 4.1 - s2['p'][0]     doc.draw_object('rect', {'x':x1, 'width':x2-x1,         'y':-0.1/1.2+0.01, 'height':0.2/1.2-0.02,         'style':'fill:url(#rod_shade); stroke-width:0.02'}, group=ball)          # draw metal balls     doc.draw_object('circle', {'cx':0, 'cy':0, 'r':[s1, s2][ib]['r'],         'style':'fill:url(#metal_spot); stroke-width:0.02'}, group=ball)     ball_charges.append(doc.draw_object('g',         {'style':'stroke-width:0.02'}, group=ball))  # find well-distributed start positions of field lines def get_startpoint_function(startpath, field):     '''     Given a vector function startpath(t), this will return a new     function such that the scalar parameter t in [0,1] progresses     indirectly proportional to the orthogonal field strength.     '''     def dstartpath(t):         return (startpath(t+1e-6) - startpath(t-1e-6)) / 2e-6     def FieldSum(t0, t1):         return ig.quad(lambda t: sc.absolute(sc.cross(             field.F(startpath(t)), dstartpath(t))), t0, t1)[0]     Ftotal = FieldSum(0, 1)     def startpos(s):         t = op.brentq(lambda t: FieldSum(0, t) / Ftotal - s, 0, 1)         return startpath(t)     return startpos  startp = [] def startpath1(t):     phi = 2. * pi * t     return (sc.array(s2['p']) + 1.5 * sc.array([cos(phi), sin(phi)])) start_func1 = get_startpoint_function(startpath1, field) nlines1 = 16 for i in range(nlines1):     startp.append(start_func1((0.5 + i) / nlines1))  def startpath2(t):     phi = 2. * pi * (0.195 + 0.61 * t)     return (sc.array(s1['p']) + 1.5 * sc.array([cos(phi), -sin(phi)])) start_func2 = get_startpoint_function(startpath2, field) nlines2 = 14 for i in range(nlines2):     startp.append(start_func2((0.5 + i) / nlines2))  # draw the field lines for p0 in startp:     line = FieldLine(field, p0, directions='both', maxr=7.)          arrow_d = 2.0     of = [0.5 + s1['r'] / arrow_d, 0.5, 0.5, 0.5 + s2['r'] / arrow_d]     doc.draw_line(line, arrows_style={'dist':arrow_d, 'offsets':of}) doc.write()