# paste this code at the end of VectorFieldPlot 1.10
# https://commons.wikimedia.org/wiki/User:Geek3/VectorFieldPlot
doc = FieldplotDocument('VFPt_flat_magnets_gap', commons=True,
    width=800, height=800)
field = Field({'coils':[[0, 1, pi/2, 2, 0.5 ,1],
                        [0, -1, pi/2, 2, 0.5 ,1]]})

doc.draw_magnets(field)

# distribute field lines proportionally to the field strength at y=0
x0 = op.brentq(lambda x: field.F([x, 0.])[1], 0, 3)
Fsum = ig.quad(lambda x: field.F([x, 0.])[1], -x0, x0)[0]
nlines = 22
xstart = [op.brentq(lambda xs: ig.quad(lambda x: field.F([x, 0.])[1],
    -x0, xs)[0] / Fsum - (0.2+i)/(nlines-0.6), -x0, x0) for i in range(nlines)]

for i, x in enumerate(xstart):
    line = FieldLine(field, [x, 0.], directions='both', maxr=12)
    if i < 2 or nlines - i <= 2:
        min_arrows = 4
    else:
        min_arrows = 3
    doc.draw_line(line, arrows_style={'dist':3, 'min_arrows':min_arrows,
        'offsets':[1, 0.5, 0.5, 0], 'scale':1.3})

for x0, y0 in ((-1, -1), (-1, 1), (1, -1), (1, 1)):
    line = FieldLine(field, [2.3 * x0, 1. * y0], directions='both', maxr=5)
    doc.draw_line(line, arrows_style={'dist':2,
        'offsets':[1, 0.5, 0.5, 0], 'scale':1.3})

doc.write()