Saturnvoyager
Saturnvoyager

Reputation: 33

Plotting a vector field using quiver

I'm trying to reproduce a 2D vector map with components

 v = 100/a * exp(-1/a^2 * ((x+0.55)^2+y^2))(-y,x) - 100/a * exp(-1/a^2 * ((x-0.55)^2+y^2))(-y,x)

and here are my codes. It did not give the map I want (see attached vector map). Could someone please help me with it?

import numpy as np
import matplotlib.pyplot as plt
import math

grid_resolution = 25
grid_size = 2*grid_resolution+1
a = 0.2
x = np.linspace(-1,1,grid_size)
y = np.linspace(-1,1,grid_size)
X,Y = np.meshgrid(x, y)

vx = np.zeros((grid_size,grid_size))
vy = np.zeros((grid_size,grid_size))

for i in range(0,grid_size):
    for j in range(0,grid_size):
        x0 = x[j]
        y0 = y[i]
        xx = (x0 + 0.55) ** 2 + y0 ** 2
        yy = (x0 - 0.55) ** 2 + y0 ** 2
        expf1 = math.exp(-xx / (a ** 2))
        expf2 = math.exp(-yy / (a ** 2))
        vx[i,j] = 100 / a * (-expf1 + expf2) * y0
        vy[i,j] = 100 / a * (expf1 - expf2) * x0

fig, ax = plt.subplots()
ax.quiver(X, Y, vx, vy)
ax.set_aspect('equal')

plt.show()

Upvotes: 3

Views: 263

Answers (1)

Zephyr
Zephyr

Reputation: 12524

In the last passage, when you compute vx[i,j] and vy[i,j], you are computing vector field components in (x0, y0), while you should compute it in the current point, so (x0 ± 0.55, y0). Moreover, you should change the sign of vx and vy in order to draw a vector field like the one you linked.

import numpy as np
import matplotlib.pyplot as plt
import math

grid_resolution = 25
grid_size = 2*grid_resolution + 1
a = 0.2
x = np.linspace(-1, 1, grid_size)
y = np.linspace(-1, 1, grid_size)
X, Y = np.meshgrid(x, y)

vx = np.zeros((grid_size, grid_size))
vy = np.zeros((grid_size, grid_size))

for i in range(0, grid_size):
    for j in range(0, grid_size):
        x0 = x[j]
        y0 = y[i]
        xx = (x0 + 0.55)**2 + y0**2
        yy = (x0 - 0.55)**2 + y0**2
        expf1 = math.exp(-xx/(a**2))
        expf2 = math.exp(-yy/(a**2))
        vx[i, j] = -100/a*(-expf1 + expf2)*y0
        if x0 > 0:
            vy[i, j] = -100/a*(expf1 - expf2)*(x0 - 0.55)
        else:
            vy[i, j] = -100/a*(expf1 - expf2)*(x0 + 0.55)

fig, ax = plt.subplots()
ax.quiver(X,Y,vx,vy)
ax.set_aspect('equal')

plt.show()

enter image description here

Upvotes: 3

Related Questions