# Restricting vector field to a subset of a rectangular region

I'm trying to restrict a vector field plot to a specific region of the plane. I found the code for the plot_vector_field3d function online, and just modified that accordingly, but I'm having trouble with the 2d case. plot_vector_field always picks points in a rectangular region, but there's a part of the region I want to skip. Here's what I have so far.

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))
theta = arctan(y/x)
dr = cos(theta)*(1-1/r^2)
dtheta = -sin(theta)*(1+1/r^2)
if r<= 1:
return 0
else:
return dr*cos(theta)-r*sin(theta)*dtheta

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

VF = plot_vector_field((dx(x,y),dy(x,y)),(x,-2,2),(y,-2,2))

show(VF + C,aspect_ratio=1)


I'm not able to upload an image file, but if you run this, I'm trying to get rid of all of the vectors within the circle.

For comparison, here's the working code in 3d. This generates a flow on a torus, with the vectors only coming from points on or in the torus, instead of the usual plot_vector_field3d, which puts the points on an entire cube. dx, dy, and dz are the functions that give the correct vector components.

twopi = 2 * RDF.pi()
DT = [] # discretized tori

n = 20 # points in the boundary
m = 5 #layers in the torus
#Discretizes T0
for k in xrange(m):
for i in xrange(n):
DT.append((2*cos(i*twopi/(n)), 2*sin(i*twopi/(n)),0))
for j in xrange(n/2):
x0 = (2 + cos(j*twopi/(n/2)))*cos(i*twopi/(n))
y0 = (2 + cos(j*twopi/(n/2)))*sin(i*twopi/(n))
z0 = ((k+1)/m)*sin(j*twopi/(n/2))
DT.append((x0,y0,z0))
points = [vector(DT[i]) for i in range(len(DT))]
vectors = [vector((dx(*point), dy(*point), dz(*point))) for point in points]

max_len = max(v.norm() for v in vectors)
scaled_vectors = [v/max_len for v in vectors]

VF = sum([plot(v,color=hue(v.norm())).translate(p) for v,p  in zip(scaled_vectors,points)])

edit retag close merge delete

Sort by » oldest newest most voted

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)

more