# paste this code at the end of VectorFieldPlot 3.2
# https://commons.wikimedia.org/wiki/User:Geek3/VectorFieldPlot
doc = FieldplotDocument('VFPt_superconductor_cylinder_E-field_potential+contour',
    width=600, height=600, commons=True)
unit = 100.
E0 = [0.0, -1.0]
sphere = {'p':sc.array([0., 0.]), 'r':1.2}

field_outside = Field([ ['homogeneous', {'Fx':E0[0], 'Fy':E0[1]}],
    ['dipole2d', {'x':sphere['p'][0], 'y':sphere['p'][1],
     'px':2*pi*sphere['r']**2 * E0[0],
     'py':2*pi*sphere['r']**2 * E0[1]}] ])

# E-field inside is actually 0. Set to finite value, so that lines continue,
# i.e. assume susceptibility chi slightly below 1.
Efield_inside = Field([ ['homogeneous', {'Fx':1e-8*E0[0], 'Fy':1e-8*E0[1]}] ])

def sphere_Efield(xy):
    if vabs(xy - sphere['p']) < sphere['r']:
        return Efield_inside.F(xy)
    else:
        return field_outside.F(xy)

def sphere_potential(xy):
    if vabs(xy - sphere['p']) < sphere['r']:
        return Efield_inside.V(xy)
    else:
        return field_outside.V(xy)

field = Field([ ['custom', {'F':sphere_Efield, 'V':sphere_potential}] ])

U0 = field.V([3, 3])
doc.draw_scalar_field(func=field.V, cmap=doc.cmap_AqYlFs, vmin=-U0, vmax=U0)
doc.draw_contours(func=field.V, levels=sc.linspace(-3, 3, 11),
    linewidth=1, linecolor='#444444')

# draw the superconducting cylinder
cylinder = doc.draw_object('g', {'id':'metal_cylinder'})

def triangle_path(phi1, phi2, r):
    x1, y1 = r * cos(radians(phi1)), r * sin(radians(phi1))
    x2, y2 = r * cos(radians(phi2)), r * sin(radians(phi2))
    d = 'M {:.4f},{:.4f}'.format(x1, y1)
    d += ' A {:.4f},{:.4f} 0 0 1 {:.4f},{:.4f}'.format(r, r, x2, y2)
    d += ' L {:.4f},{:.4f}'.format(-x2, -y2)
    d += ' A {:.4f},{:.4f} 0 0 0 {:.4f},{:.4f}'.format(r, r, -x1, -y1)
    d += ' L {:.4f},{:.4f} Z'.format(x1, y1)
    return d

def grey(bright):
    return '#' + 3 * ('%02x' % int(256. * bright - 0.5))

doc.draw_object('circle',
    {'cx':sphere['p'][0], 'cy':sphere['p'][1], 'r':'{:.4f}'.format(sphere['r']),
    'style':'fill:' + grey(0.75) + '; stroke:none'}, group=cylinder)

for phi0 in [0]:
    ncolors = 25
    for a in sc.linspace(.5 / ncolors, 1 - .5 / ncolors, ncolors):
        bright = 0.75 + 0.15 * a
        phi1 = phi0 - 60 * (acos(2. * a - 1) / pi)**1.5
        phi2 = phi0 + 60 * (acos(2. * a - 1) / pi)**1.5
        d = triangle_path(phi1, phi2, sphere['r'])
        doc.draw_object('path', {'d':d,
            'style':'fill:' + grey(bright) + '; stroke:none'}, group=cylinder)

doc.draw_object('circle',
    {'cx':sphere['p'][0], 'cy':sphere['p'][1], 'r':'{:.4f}'.format(sphere['r']),
    'style':'fill:none; stroke:black; stroke-width:0.02'}, group=cylinder)
cylinder_charges = doc.draw_object('g', {'style':'stroke-width:0.02; stroke-linecap:square'}, group=cylinder)

nlines = 24
for iline in range(nlines):
    a = -3.6 + 7.2 * (0.5 + iline) / nlines
    line = FieldLine(field, [a, 200], maxr=210, directions='forward', pass_dipoles=1)
    doc.draw_line(line, linewidth=2.4, arrows_style={'at_potentials':[-2.1, 2.1]})

    # draw little charge signs near the surface
    path_minus = 'M {0:.5f},0 h {1:.5f}'.format(-2./unit, 4./unit)
    path_plus = 'M {0:.5f},0 h {1:.5f} M 0,{0:.5f} v {1:.5f}'.format(-2./unit, 4./unit)
    
    # check if fieldline crosses sphere surface
    tlist = sc.linspace(0., 1., 1001)
    for i in range(1, len(tlist)):
        in0 = vabs(line.get_position(tlist[i-1]) - sphere['p']) <= sphere['r']
        in1 = vabs(line.get_position(tlist[i]) - sphere['p']) <= sphere['r']
        if in0 != in1:
            # find the point where the field line cuts the surface
            t = optimize.brentq(lambda t: vabs(line.get_position(t)
                - sphere['p']) - sphere['r'], tlist[i-1], tlist[i])
            pr = line.get_position(t) - sphere['p']
            cpos = 0.92 * sphere['r'] * pr / vabs(pr)
            if in1:
                path_d = path_minus
            else:
                path_d = path_plus
            doc.draw_object('path', {'stroke':'black', 'd':path_d,
                'transform':'translate({:.5f},{:.5f})'.format(
                    round(unit*cpos[0])/unit, round(unit*cpos[1])/unit)},
                    group=cylinder_charges)

doc.write()