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)