# paste this code at the end of VectorFieldPlot 2.4
doc = FieldplotDocument('VFPt_superconductor_cylinder_B-field_potential+contour',
    width=600, height=600, commons=True)

B0 = [0.0, -1.0]
sphere = {'p':sc.array([0., 0.]), 'r':1.2}

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

Hfield_inside = Field([ ['homogeneous', {'Fx':2.*B0[0], 'Fy':2.*B0[1]}] ])

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

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

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

U0 = field.V([0, 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.6, 3.6, 13),
    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':0, 'cy':0, 'r':'{:.4f}'.format(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, r)
        doc.draw_object('path', {'d':d,
            'style':'fill:' + grey(bright) + '; stroke:none'}, group=cylinder)

doc.draw_object('circle', {'cx':0, 'cy':0, 'r':str(r),
    'style':'fill:none; stroke:black; stroke-width:0.02'}, group=cylinder)

n_lines = 20
for i in range(n_lines):
    a = -3 + 6 * (0.5 + i) / n_lines
    line = FieldLine(field, [a, 8], maxr=16, pass_dipoles=1)
    doc.draw_line(line, linewidth=2.4, arrows_style={'at_potentials':[-2.7, 2.7]})

doc.write()