# 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()
