I am attempting to write a Python code to simulate many particles in a confined box. These particles behave in such a way that they move in the box in straight lines with a slight angular noise (small changes in the direction of the particle path). They should interact by acknowledging the other particle and ‘shuffle/squeeze’ past each other and continue on their intended path, much like humans on a busy street. Eventually, the particles should cluster together when the density of particles (or packing fraction) reaches a certain value, but I haven’t got to this stage yet.

The code currently has particle interactions and I have attempted the angular noise but without success so far.

However, I have a feeling there are parts of the code that are inefficient or which could be either sped up or written more conveniently.

If anyone has any improvements for the code speed or ideas which may help with the interactions and/or angular noise that would be much appreciated. I will also leave an example of an animation which is my aim: https://warwick.ac.uk/fac/sci/physics/staff/research/cwhitfield/abpsimulations

The above link shows the animation I am looking for, although I don’t need the sliders, just the box, and moving particles. The whole code is shown below:

```
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
def set_initial_coordinates():
x_co = (np.random.uniform(0, 2) for i in range(n_particles))
y_co = (np.random.uniform(0, 2) for i in range(n_particles))
return x_co, y_co
def set_initial_velocities():
x_vel = np.array((np.random.uniform(-1, 1) for i in range(n_particles)))
y_vel = np.array((np.random.uniform(-1, 1) for i in range(n_particles)))
return x_vel, y_vel
def init():
ax.set_xlim(-0.05, 2.05)
ax.set_ylim(-0.07, 2.07)
return ln,
def update(dt):
xdata = initialx + vx * dt
ydata = initialy + vy * dt
fx = np.abs((xdata + 2) % 4 - 2)
fy = np.abs((ydata + 2) % 4 - 2)
for i in range(n_particles):
for j in range(n_particles):
if i == j:
continue
dx = fx(j) - fx(i) # distance in x direction
dy = fy(j) - fy(i) # distance in y direction
dr = np.sqrt((dx ** 2) + (dy ** 2)) # distance between x
if dr <= r:
force = k * ((2 * r) - dr) # size of the force if distance is less than or equal to radius
# Imagine a unit vector going from i to j
x_comp = dx / dr # x component of force
y_comp = dy / dr # y component of force
fx(i) += -x_comp * force # x force
fy(i) += -y_comp * force # y force
ln.set_data(fx, fy)
return ln,
# theta = np.random.uniform(0, 2) for i in range(n_particles)
n_particles = 10
initialx, initialy = set_initial_coordinates()
vx, vy = set_initial_velocities()
fig, ax = plt.subplots()
x_co, y_co = (), ()
ln, = plt.plot((), (), 'bo', markersize=15) # radius 0.05
plt.xlim(0, 2)
plt.ylim(0, 2)
k = 1
r = 0.1
t = np.linspace(0, 10, 1000)
ani = FuncAnimation(fig, update, t, init_func=init, blit=True, repeat=False)
plt.show()
```
```