Revision history [back]

I figured something out, using a manual grid, and scaling the drawn arrows appropriately. I also modified it so any level set between z = -1 and z = 1 will work. Thank you for the help with cleaning up the code!

z = 0 #choose a z-value for the level set

x,y = var('x y')

# Use "real double field" floating-point values for efficiency
zero = RDF.zero()
one = RDF.one()
four = 4 * one
twopi = 2 * RDF.pi()

r1 = 2 + cos(arcsin(z)) #sets radii for the inner and outer circles
r2 = 2 + cos(arcsin(z) + pi)

# Circles as graphic objects
Cr1 = circle((0,0), r1,color='black')
Cr2 = circle((0,0), r2,color='black')
C4 = circle((0,0), 4,color='black')
C = Cr1 + Cr2 + C4

# 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()

# h0 and h1

def h0(x,y):
"""
    Return piecewise affine function of `r = norm(x,y)`

    With r = norm((x,y)) = sqrt(x^2+y^2), h0(x,y) is:
    * r inside C_r1
    * 1 between C_r1 and C_r2
    * 4 - r between C_r2 and C_4
    * 0 outside C_4
"""
    r = norm((x,y))
    if r < r1:
        return r
    elif r1 <= r <= r2:
        return one
    elif r2 <= r <= 4:
        return four - r
    else:
        return zero

def h1(x,y):
"""
    Return piecewise affine function of `r = norm(x,y)`

    With r = norm((x,y)) = sqrt(x^2+y^2), h0(x,y) is:
    * r inside C_r1
    * 0 between C_r1 and C_r2
    * 4 - r between C_r2 and C_4
    * 0 outside C_4
"""
    r = norm((x,y))
    if r < r1:
        return r
    elif r1 <= r <= r2:
        return zero
    elif r2 <= r <= 4:
        return four - r
    else:
        return zero

def dx(x,y):
    return y * h0(x,y)

def dy(x,y):
    return -x * h0(x,y) + -z * h1(x,y)

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

A = []
scale = 15
for i in range(scale+1):
    x = -4 + i*8/scale
    for j in range(scale+1):
        y = -4 + j*8/scale
        a = arrow((x,y),arrowplot(x,y),arrowsize=3, color=hue(norm((x - dx(x,y), y - dy(x,y)))*(1/ep)))
        A.append(a)

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