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)