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"