Ask Your Question
1

Restricting vector field to a subset of a rectangular region

asked 2015-09-14 22:01:05 +0100

jford1906 gravatar image

updated 2015-09-16 00:26:20 +0100

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 flag offensive close merge delete

1 Answer

Sort by ยป oldest newest most voted
1

answered 2015-09-17 13:32:50 +0100

jford1906 gravatar image

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)
edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2015-09-14 22:01:05 +0100

Seen: 541 times

Last updated: Sep 17 '15