billiards

A 2D physics engine for simulating dynamical billiards

billiards is a python library that implements a very simple physics engine: It simulates the movement and elastic collisions of hard, disk-shaped particles in a two-dimensional world.

Features

• Collisions are found and resolved exactly. No reliance on time steps, no tunneling of high-speed bullets!
• Quick state updates thanks to numpy, especially if there are no collisions between the given start and end times.
• Static obstacles to construct a proper billiard table.
• Balls with zero radii behave like point particles, useful for simulating dynamical billiards (although this library is not optimized for point particles).
• Optional features: plotting and animation with matplotlib, interaction with pyglet.
• Free and open source software under the MIT license.

Installation

billiards is a library for Python 3. It only depends on numpy.

Billiard systems can be visualized with matplotlib (and tqdm to display progress in visualize_matplotlib.animate). Interaction with the simulation is possible via pyglet. These visualization features are optional.

Clone the repository from GitHub and install the package:

git clone https://github.com/markus-ebke/python-billiards.git
cd python-billiards/
pip install .[visualize]

Quickstart

All important classes (the billiard simulation and obstacles) are accessible from the top-level module. The visualization modules must be imported separately and will load matplotlib or pyglet. For the following examples we will use matplotlib visualizations.

>>> import billiards  # access to Billiard, Disk and InfiniteWall
>>> import billiards.visualize_matplotlib as visualize  # for plot and animate
>>> import matplotlib.pyplot as plt  # for plt.show()

Example: Computing π with pool

Let's compute the first few digits of π using a billiard simulation following the setup of Gregory Galperin. We need a billiard table with a vertical wall and two balls:

>>> obstacles = [billiards.obstacles.InfiniteWall((0, -1), (0, 1), blocked="right")]
>>> bld = billiards.Billiard(obstacles)
>>> bld.add_ball((3, 0), (0, 0), radius=0.2, mass=1)  # returns index of new ball
0
>>> bld.add_ball((6, 0), (-1, 0), radius=1, mass=100**5)
1

Using the visualize module, let's see how this initial state looks:

>>> visualize.plot(bld)
(<Figure size 800x600 with 1 Axes>, <Axes: >)
>>> plt.show()

The Billiard.evolve method simulates our billiard system from bld.time until a given end time. It returns a list of collisions (ball-ball and ball-obstacle collisions).

>>> bld.next_collision  # (time, ball index, ball index or obstacle)-triplet
(1.8000000000000005, 0, 1)
>>> total_collisions = 0
>>> for i in [1, 2, 3, 4, 5]:
...     total_collisions += sum(bld.evolve(i))
...     print(f"Until t = {bld.time}: {total_collisions} collisions")
...
Until t = 1: 0 collisions
Until t = 2: 1 collisions
Until t = 3: 1 collisions
Until t = 4: 4 collisions
Until t = 5: 314152 collisions

The first collision happened at time t = 1.8. Until t = 4 there were only 4 collisions, but then between t = 4 and t = 5 there were several thousands. Let's see how the situation looks now:

>>> bld.time  # current time
5
>>> visualize.plot(bld)
(<Figure size 800x600 with 1 Axes>, <Axes: >)
>>> plt.show()

Let's advance the simulation to t = 16. As we can check, there won't be any other collisions after this time:

>>> total_collisions += sum(bld.evolve(16))
>>> bld.balls_velocity  # nx2 numpy array where n is the number of balls
array([[0.73463055, 0.        ],
       [1.        , 0.        ]])
>>> bld.next_ball_ball_collision
(inf, -1, 0)
>>> bld.next_ball_obstacle_collision
(inf, 0, None)
>>> visualize.plot(bld)
(<Figure size 800x600 with 1 Axes>, <Axes: >)
>>> plt.show()

Both balls are moving towards infinity, the smaller ball to slow to catch the larger one. What is the total number of collisions?

>>> total_collisions
314159
>>> import math
>>> math.pi
3.141592653589793

The first six digits match! For an explanation why this happens, see Galperin's paper Playing pool with π (the number π from a billiard point of view) or the series of youtube videos by 3Blue1Brown starting with The most unexpected answer to a counting puzzle.

Lastly, I want to point out that all collisions were elastic, i.e. they conserved the kinetic energy (within floating point accuracy):

>>> 100**5 * (-1) ** 2 / 2  # kinetic energy = m v^2 / 2 at the beginning
5000000000.0
>>> v_squared = (bld.balls_velocity**2).sum(axis=1)
>>> (bld.balls_mass * v_squared).sum() / 2  # kinetic energy now
np.float64(4999999999.989935)

The video examples/pi_with_pool.mp4 replays the whole billiard simulation (it was created using visualize.animate).

More Examples

Setup:

>>> import matplotlib.pyplot as plt
>>> import billiards
>>> import billiards.visualize_matplotlib as visualize

First shot in Pool (no friction)

Construct the billiard table:

>>> width, length = 112, 224
>>> bounds = [
...     billiards.InfiniteWall((0, 0), (length, 0)),  # bottom side
...     billiards.InfiniteWall((length, 0), (length, width)),  # right side
...     billiards.InfiniteWall((length, width), (0, width)),  # top side
...     billiards.InfiniteWall((0, width), (0, 0)),  # left side
... ]
>>> bld = billiards.Billiard(obstacles=bounds)

Arrange the balls in a pyramid shape:

>>> from math import sqrt
>>> radius = 2.85
>>> for i in range(5):
...     for j in range(i + 1):
...         x = 0.75 * length + radius * sqrt(3) * i
...         y = width / 2 + radius * (2 * j - i)
...         bld.add_ball((x, y), (0, 0), radius)
...

Add the white ball and give it a push, then view the animation:

>>> bld.add_ball((0.25 * length, width / 2), (length / 3, 0), radius)
>>> anim, fig, ax = visualize.animate(bld, end_time=10, figsize=(10, 5.5))
>>> plt.show()

See examples/pool.mp4

Brownian motion

The billiard table is a square box:

>>> obs = [
...     billiards.InfiniteWall((-1, -1), (1, -1)),  # bottom side
...     billiards.InfiniteWall((1, -1), (1, 1)),  # right side
...     billiards.InfiniteWall((1, 1), (-1, 1)),  # top side
...     billiards.InfiniteWall((-1, 1), (-1, -1)),  # left side
... ]
>>> bld = billiards.Billiard(obstacles=obs)

Distribute small particles (atoms) uniformly in the square, moving in random directions but with the same speed:

>>> from math import cos, pi, sin
>>> from random import uniform
>>> for i in range(250):
...     pos = [uniform(-1, 1), uniform(-1, 1)]
...     angle = uniform(0, 2 * pi)
...     vel = [cos(angle), sin(angle)]
...     bld.add_ball(pos, vel, radius=0.01, mass=1)
...

Add a bigger ball (like a dust particle)

>>> idx = bld.add_ball((0, 0), (0, 0), radius=0.1, mass=10)

and simulate until t = 50, recording the position of the bigger ball at each collision (this will take some time)

>>> poslist = [bld.balls_position[idx].copy()]  # record initial position
>>> def record(t, p, u, v, i_o):
...     poslist.append(p)
...
>>> bld.evolve(50, ball_callbacks={idx: record})
>>> poslist.append(bld.balls_position[idx].copy())  # record last position

Plot the billiard and overlay the path of the particle

>>> fig, ax = visualize.plot(bld, arrow_size=0, figsize=(7, 7))
>>> poslist = np.asarray(poslist)
>>> ax.plot(poslist[:, 0], poslist[:, 1], color="red")
[<matplotlib.lines.Line2D object at 0x...>]
>>> plt.show()

Authors

• Markus Ebke - https://github.com/markus-ebke