Knowing the velocity field, how can I find the trajectory of a "massless" particle inside that velocity field?

Question

x,y = np.meshgrid(np.linspace(-8,8,30),np.linspace(-8,8,30))
q=3
w=10
freq=2
wavelenght=0.6
r=x**2+y**2
u=np.zeros((len(x),len(y)))
v=np.zeros((len(x),len(y)))
for i in range(0,len(x)):
    for j in range (0,len(y)):
        if (r[i,j]<=q**(3/4)):
            x[i,j]=0
            y[i,j]=0
        if (r[i,j]>q**(3/4)):
            u[i,j]=freq*wavelenght 

This is my velocity field and this is how it looks velocity field

I tried some tips I found on other questions similar to mine however I was getting blank graphs or lines that do not make any sense. I guess it's partly due to the zeros in the middle of the graph. What I would appreciate is a method to send massless particles from an initial point with an initial direction and to see how it moves in this field.

Thank you !

Answer 1

It depends what you really need. I can offer something I wrote, but there might be better python packages or functions that can do better. Since I do not really have your vector field, I just tested the code with a mesh-grid vector field coming from the restriction of some polynomial vector field onto the vertices of the mesh-grid. Feel free to use your own mesh-grid vector field

import numpy as np
import matplotlib.pyplot as plt

# left and right ends of the interval on the x-axis:
x_left = -8 
x_right = 8
# bottom and top ends of the interval of the y-axis:
y_bottom = -8
y_top = 8
# number of points to be included in the mesh-grid along the x-axis
N_points_x = 200
# number of points to be included in the mesh-grid along the y-axis
N_points_y = 200

# generate the mesh-gird
x, y = np.meshgrid(np.linspace(x_left, x_right, N_points_x),np.linspace(y_bottom, y_top, N_points_y))


# a function that calculates the indices of the 
# closest point from the mesh-grid to a given arbitrary point p
def indxs(p, x, y):
    i = np.int(x.shape[1]*(p[0] - x[0,0])/(x[0,-1] - x[0,0]))
    j = np.int(y.shape[0]*(p[1] - y[0,0])/(y[-1,0] - y[0,0]))
    return i, j


# a simple quadratic interpolation of the mesh-grid vector filed
# to a quadratically interpolated vector field at a point p inside 
# mesh-grid the square in which p is located
def VF(p, x, y, V):
    i, j = indxs(p, x, y)
    if 0 < i and i < x.shape[1]-1 and 0 < j and j < y.shape[0]-1: 
        a = ( p[0] - x[0, i] ) / (x[0, i+1] - x[0, i]) 
        b = ( p[1] - y[j, 0] ) / (y[j+1, 0] - y[j, 0])
        W = (1-a) * (1-b) * V[:, j, i]  +  (1-a) * b * V[:, j, i+1]
        W = W  +  a * (1-b) * V[:, j+1, i]  +  a * b * V[:, j+1, i+1]   
        return W #/ np.linalg.norm(W) # you can also normalize the vector field to get only the trajecotry, without accounting for parametrization
    else:
        return np.array([0.0, 0.0])


# integrating the v#ector field one time step 
# starting from a given point p
# uses Runge-Kutta 4 integrations, which 
# allows you to sample the vector fields at four different near-by points and 
# 'average' the result 
def VF_flow(p, x, y, V, t_step):
    k1 = VF(p, x, y, V)
    k2 = VF(p + t_step*k1/2, x, y, V)
    k3 = VF(p + t_step*k2/2, x, y, V)
    k4 = VF(p + t_step*k3, x, y, V)
    return p + t_step * (k1 + 2*k2 + 2*k3 + k4) / 6


def VF_trajectory(p, x, y, V, t_step, n_iter):
    traj = np.empty( (2, n_iter), dtype = float)
    traj[:, 0] = p
    #m = 0
    #while m < n_iter and (x[0,0] < traj[0,m]     traj[1,m]): 
    for m in range(n_iter-1):
        traj[:, m+1] = VF_flow(traj[:, m], x, y, V, t_step) 
        m = m+1
    return traj

'''
This part is just to set up a test example, to run the code on. 
When needed, generate your own vector field
'''
# initialize an empty array for the values of the vector filed at each point from the mesh-grid
V = np.empty((2, x.shape[1], y.shape[0]), dtype=float)


# a function that turns a planar polynomial vector field into a vector field on the mesh-grid
def U(x, y):
    return -x + 0.3*x*y + 0.7*y**2, -0.5*y + 0.2*x*y - 0.1*x**2


# turns a planar polynomial vector field into a vector field on the mesh-grid
V[0,:,:], V[1,:,:] = U(x, y)
'''
End of the test example setup
'''

n_iterations = 700

# initial point
p0 = np.array([ -7.3, 3.1])
#p0 = np.array([-3, 5])

t_step = 0.005


trajectory = VF_trajectory(p0 , x, y, V, t_step, n_iterations)

# this function actually allows you to generate a plot of the flow-trajectories of 
# a mesh-grid vector field
plt.streamplot(x, y, V[0,:,:], V[1,:,:])
# here I am plotting the trajectory generated by my functions on top of the plot of the 
# flow-portrait of the mesh-grid vector field
plt.plot(trajectory[0, :], trajectory[1, :], '-r')
plt.axis('square')
plt.axis([x_left, x_right, y_bottom, y_top])
plt.show()