How to create 'billiard ball' reflection boundary condition in python?

Question

According to Erwin Schrodinger (in What is Life?), diffusion can be explained entirely by the random motion of particles. I want to test this myself by create a program the creates a time-step visualization of the diffusion of "gas molecules" in a closed container. The initial conditions would have two partitions, one with low and one with high concentration. After t0 the partition is removed and the gas is allowed to diffuse. The only mechanism I want to use is adding displacement random vectors to each molecule. The initial conditions would look like this.

The part of the problem that I con't figure out is how to create a simple billiards type reflection when the molecule hits the bounding surfaces. I am assume simple symmetrical reflections (angle in = angle out at boundaries). I haven't started the code at all because I don't know how to deal with this part, while I know how to do the rest of it. I know this is more of a math question, but how can I create these boundary conditions in python? Ideally, I would like to have to program this functionality myself so that I can understand it, rather than using a pre-built package that could do this. This is what I am looking for, for any given molecule.

Ultimately, all I really need is this: given the initial location (x1,y2), a vector magnitude v, an angle theta, and the box size and location, what is the final resting position of the molecule (x2,y2).

Answer 1

You don't need to calculate the reflection angle, just decompose the problem in two: one for x and one for y. In both cases, you need the particle to "go back" when it exceeds the boundary.

I've done this time ago for an excercise studying particle density in fluids. The easiest thing is to consider a (0, 1) boundary in both directions. The following code should do it (tip: a proper use of abs will create the equivalent of a reflection):

x0 = [.1, .9]
delta = [-0.2, 0.3]
x1 = [(1-abs(abs(xi + di)-1)) for xi, di in zip(x0, delta)]
print(x1)
# 0.1, 0.8
#or using numpy:
x1 = 1-np.abs(np.abs(np.asarray(x0) + np.asarray(delta))-1)
print(x1)
>> [0.09999999999999998, 0.8]
   array([0.1, 0.8])

I'm assuming from your question that you are neglecting particle-particle collision and particle-particle "non-superposition"

Answer 2

Here is a simple implementation. I change the movement vectors only every tenth step, that way one can visually check boundary reflections. Particles flash red when the movement vectors are updated.

The trick as described ħere is to "unroll" the bounding box. Instead we let the particles move unconstrained and then fold space into the bounding box.

import numpy as np
import pylab
from matplotlib.animation import FuncAnimation

xy = np.random.uniform(-1, 1, (2, 200))
xy[0, :160] = np.abs(xy[0, :160])
xy[0, 160:] = -np.abs(xy[0, 160:])
xy += 1

f, a = pylab.subplots()
pxy, = pylab.plot(*xy, 'o')

def init():
    a.set_xlim(0, 2)
    a.set_ylim(0, 2)
    return pxy,

def update(frame):
    global inc, xy
    if frame % 1 < 0.01:
        inc = np.random.normal(0, 0.01, xy.shape)
        pxy.set_markerfacecolor('red')
    elif frame % 1 < 0.11:
        pxy.set_markerfacecolor('blue')        
    xy += inc
    fxy = np.abs((xy+2)%4-2)
    pxy.set_data(*fxy)
    return pxy,

anim = FuncAnimation(f, update, frames=np.arange(1200) / 10,
                     init_func=init, blit=True)

pylab.show()

Answer 3

So a couple things to keep in mind:

1. You need a friction component or else the particle will keep moving forever (conservation of energy). Friction happens as a function of velocity in this case, and friction will also happen on bounces.

2. If it's just a single particle, you can calculate it by defining bounding boxes like x is between 0 and 5, y is between 0 and 3 for example. Then you can calculate the intercept with a wall by plugging in the x = 5 value and then solving for y in the equation for the line.

For one particle, you don't have to do it parametrically with increments of t_0, you can calculate the intercepts and basically zoom it on over there. For multiple, you would have to calculate inter-molecular diffusion and collision forces...which is a much harder problem that should be done parametrically.

You would have to calculate collisions, which is when the centres of two molecules are 2*radius away from each other, then do a collision that conserves momentum.