# Revision history [back]

Figured it out.

x,y = var('x y')
C = circle((0,0),1)

# useful functions

def norm2((x,y)):
"""
Norm squared of (x,y)
"""
return x^2 + y^2

def norm((x,y)):
"""
Norm of (x,y)
"""
return norm2((x,y)).sqrt()

def dy(x,y):
r = norm((x,y))
if r <= 1:
return 0
else:
if x == 0:
theta = 0
else:
theta = arctan(y/x)
dr = cos(theta)*(1-1/r^2)
dtheta = -sin(theta)*(1+1/r^2)
return dr*cos(theta)-r*sin(theta)*dtheta

def dx(x,y):
r = norm((x,y))
if r <= 1:
return 0
else:
if x == 0:
theta = 0
else:
theta = arctan(y/x)
dr = cos(theta)*(1-1/r^2)
dtheta = -sin(theta)*(1+1/r^2)
return dr*sin(theta) + r*cos(theta)*dtheta

ep = 0.3
def arrowplot(x,y):
n = norm((x - dx(x,y), y - dy(x,y)))*(1/ep)
if n == 0:
return ((x,y))
else:
return ((x+dx(x,y)/n), (y+dy(x,y)/n))

A = []
scale = 10
for i in range(scale+1):
x = -2 + i*4/scale
for j in range(scale+1):
y = -2 + j*4/scale
if (-2.1 <= arrowplot(x,y)[0] <= 2.1) and (-2.1 <= arrowplot(x,y)[1] <= 2.1):
a = arrow((x,y),arrowplot(x,y),arrowsize=1, color='black') #hue(norm((x - dx(x,y), y - dy(x,y)))*(1/ep)))
A.append(a)

show(sum(A)+C,aspect_ratio=1)