# Two-dimensional collisions

This small Python project is a physical simulation of two-dimensional physics. The animation is carried out using Matplotlib's FuncAnimation method and is implemented by the class Simulation. Each "particle" of the simulation is represented by an instance of the Particle class and depicted as a circle with a fixed radius which undergoes elastic collisions with other particles.

A Particle's position and velocity vectors are represented by NumPy arrays r and v, but their components can be retrieved and set by using the shortcuts Particle.x, Particle.y etc., implemented as class properties. Particles move within the domain $0 \le x < 1$, $0 \le y < 1$ and are reflected (elastically) off the edges (walls) of this domain. The elastic collision between two particles occurs when their centres are separated by the sum of their radii and the collision event changes their velocities to:

$$\boldsymbol{u}_1 = \boldsymbol{v}_1 - \frac{2m_2}{M}\frac{(\boldsymbol{v}_1 - \boldsymbol{v}_2)\cdot(\boldsymbol{r}_1 - \boldsymbol{r}_2)}{|\boldsymbol{r}_1 - \boldsymbol{r}_2|^2}(\boldsymbol{r}_1 - \boldsymbol{r}_2)$$

and

$$\boldsymbol{u}_2 = \boldsymbol{v}_2 - \frac{2m_1}{M}\frac{(\boldsymbol{v}_2 - \boldsymbol{v}_1)\cdot(\boldsymbol{r}_2 - \boldsymbol{r}_1)}{|\boldsymbol{r}_1 - \boldsymbol{r}_2|^2}(\boldsymbol{r}_2 - \boldsymbol{r}_1),$$

as dictated by the requirements of conservation of momentum and kinetic energy. The Wikipedia article on elastic collisions is a good place to learn more about this. Here, $M=m_1+m_2$ is the total mass of the colliding pair, and the particle's mass is considered to be just the square of its radius (i.e. its density is $\pi^{-1}$).

The most up-to-date code is on GitHub but as initially committed is given below.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
from matplotlib import animation
from itertools import combinations

class Particle:
"""A class representing a two-dimensional particle."""

def __init__(self, x, y, vx, vy, radius=0.01, styles=None):
"""Initialize the particle's position, velocity, and radius.

Any key-value pairs passed in the styles dictionary will be passed
as arguments to Matplotlib's Circle patch constructor.

"""

self.r = np.array((x, y))
self.v = np.array((vx, vy))

self.styles = styles
if not self.styles:
# Default circle styles
self.styles = {'edgecolor': 'b', 'fill': False}

# For convenience, map the components of the particle's position and
# velocity vector onto the attributes x, y, vx and vy.
@property
def x(self):
return self.r[0]
@x.setter
def x(self, value):
self.r[0] = value
@property
def y(self):
return self.r[1]
@y.setter
def y(self, value):
self.r[1] = value
@property
def vx(self):
return self.v[0]
@vx.setter
def vx(self, value):
self.v[0] = value
@property
def vy(self):
return self.v[1]
@vy.setter
def vy(self, value):
self.v[1] = value

def overlaps(self, other):
"""Does the circle of this Particle overlap that of other?"""

def draw(self, ax):
"""Add this Particle's Circle patch to the Matplotlib Axes ax."""

return circle

"""Advance the Particle's position forward in time by dt."""

self.r += self.v * dt

# Make the Particles bounce off the walls
if self.x - self.radius < 0:
self.vx = -self.vx
if self.x + self.radius > 1:
self.vx = -self.vx
if self.y - self.radius < 0:
self.vy = -self.vy
if self.y + self.radius > 1:
self.vy = -self.vy

class Simulation:
"""A class for a simple hard-circle molecular dynamics simulation.

The simulation is carried out on a square domain: 0 <= x < 1, 0 <= y < 1.

"""

radius can be a single value or a sequence with n values.

Any key-value pairs passed in the styles dictionary will be passed
as arguments to Matplotlib's Circle patch constructor when drawing
the Particles.

"""

"""Initialize the n Particles of the simulation.

Positions and velocities are chosen randomly; radius can be a single
value or a sequence with n values.

"""

try:
except TypeError:
# r isn't iterable: turn it into a generator that returns the
# same value n times.
for i in range(n):

self.n = n
self.particles = []
# Try to find a random initial position for this particle.
while True:
# Choose x, y so that the Particle is entirely inside the
# domain of the simulation.
# Choose a random velocity (within some reasonable range of
# values) for the Particle.
vr = 0.1 * np.random.random() + 0.05
vphi = 2*np.pi * np.random.random()
vx, vy = vr * np.cos(vphi), vr * np.sin(vphi)
particle = Particle(x, y, vx, vy, rad, styles)
# Check that the Particle doesn't overlap one that's already
# been placed.
for p2 in self.particles:
if p2.overlaps(particle):
break
else:
self.particles.append(particle)
break

def handle_collisions(self):
"""Detect and handle any collisions between the Particles.

When two Particles collide, they do so elastically: their velocities
change such that both energy and momentum are conserved.

"""

def change_velocities(p1, p2):
"""
Particles p1 and p2 have collided elastically: update their
velocities.

"""

M = m1 + m2
r1, r2 = p1.r, p2.r
d = np.linalg.norm(r1 - r2)**2
v1, v2 = p1.v, p2.v
u1 = v1 - 2*m2 / M * np.dot(v1-v2, r1-r2) / d * (r1 - r2)
u2 = v2 - 2*m1 / M * np.dot(v2-v1, r2-r1) / d * (r2 - r1)
p1.v = u1
p2.v = u2

# We're going to need a sequence of all of the pairs of particles when
# we are detecting collisions. combinations generates pairs of indexes
# into the self.particles list of Particles on the fly.
pairs = combinations(range(self.n), 2)
for i,j in pairs:
if self.particles[i].overlaps(self.particles[j]):
change_velocities(self.particles[i], self.particles[j])

"""Advance the animation by dt, returning the updated Circles list."""

for i, p in enumerate(self.particles):
self.circles[i].center = p.r
self.handle_collisions()
return self.circles

for i, p in enumerate(self.particles):
self.handle_collisions()

def init(self):
"""Initialize the Matplotlib animation."""

self.circles = []
for particle in self.particles:
self.circles.append(particle.draw(self.ax))
return self.circles

def animate(self, i):
"""The function passed to Matplotlib's FuncAnimation routine."""

return self.circles

def do_animation(self, save=False):
"""Set up and carry out the animation of the molecular dynamics.

To save the animation as a MP4 movie, set save=True.
"""

fig, self.ax = plt.subplots()
for s in ['top','bottom','left','right']:
self.ax.spines[s].set_linewidth(2)
self.ax.set_aspect('equal', 'box')
self.ax.set_xlim(0, 1)
self.ax.set_ylim(0, 1)
self.ax.xaxis.set_ticks([])
self.ax.yaxis.set_ticks([])
anim = animation.FuncAnimation(fig, self.animate, init_func=self.init,
frames=800, interval=2, blit=True)
if save:
Writer = animation.writers['ffmpeg']
writer = Writer(fps=100, bitrate=1800)
anim.save('collision.mp4', writer=writer)
else:
plt.show()

if __name__ == '__main__':
nparticles = 50
styles = {'edgecolor': 'C0', 'linewidth': 2, 'fill': None}
sim.do_animation(save=False)

Current rating: 5