# Simulating billiards in Julia

## 1. Author

George Datseris, contributor of JuliaDynamics and JuliaMusic

## 2. This tutorial

In this tutorial we will simulate a dynamical billiard. This is a very simple system where a point particle is propagating inside a domain with constant speed. When encountering a boundary, the particle undergoes specular reflection.

There is a Julia package to simulate these kind of systems, `DynamicalBilliards`. In this tutorial we are simply creating a simplified version of this package that has less features, less performance and less "safety control". Other than that, the core of how `DynamicalBilliards` works is identical to the present tutorial.

## 3. Features of Julia highlighted

• Intuitive coding of the problem
• High-end Performance without doing anything
• Simple metaprogramming for performance gain
• Extendability with least possible effort
• Extandability applies to all aspects of the system (Multiple Dispatch)

## 4. The system and algorithm

A billiard is composed of three things:

• Particles that propagate inside the billiard
• Obstacles that compose the billiard
• Processes that perform the propagation

It is important that these three things, "particles", "obstacles" and "processes" remain independent. This will allow our code to be extendable!

## 5. The Algorithm

To simulate a billiard one has to follow a very simple algorithm. Assuming that a particle `p` has already been initialized in the billiard, the steps are

• Calculate the collisions of `p` with all obstacles in the billiard
• Find the collision that happens first in time, and the obstacle it corresponds to
• Propagate the particle to this collision point
• Do specular reflection at the boundary of the obstacle to be collided with
• Rinse and repeat!

The most central part of this algorithm is a function that given a particle and obstacle it returns the collision (point and time, if they exists). What is crucial is that this function does not "belong" anywhere. It is a process, independent of what obstacle and particle we give it. This independence allows the code to become highly extendable (see below).

This independence is possible only through multiple dispatch. Without it, one way or the other, this method would be made to "belong" in either an `Obstacle` "class" or a `Particle` "class". This would imped extendability though! Not having to "assing" a process to a specific "class" is a huge benefit of Julia.

## 6. The basic `struct`s we will need

All vectors we'll use are two dimensional. We can therefore take advantage of `StaticArrays`. For a shortcut, I'll define `SV` to be a 2D vector of Floats (all the vectors we will use will be of this type)

```using StaticArrays
const SV = SVector{2,Float64}```

We now define the particle `struct`. For extendability, it is best to define it based on an abstract type.

```abstract type AbstractParticle end

mutable struct Particle <: AbstractParticle
pos::SV
vel::SV
end
Particle(x0, y0, φ0) = Particle(SV(x0, y0), SV(cos(φ0), sin(φ0)))```

Then, we define the obstacle `struct`s. Again, to allow extendability, abstraction and generalization, it is best to define all obstacles as a subtype of an abstract type.

```abstract type Obstacle end

struct Wall <: Obstacle
sp::SV
ep::SV
normal::SV
end

struct Disk <: Obstacle
c::SV
r::Float64
end```

A "billiard" in the following will simply be a `Tuple` of `Obstacle` subtypes. Specifically, the type of the billiard is `NTuple{N, Obstacle} where {N}`, so it is convenient to define:

`const Billiard = NTuple{N, Obstacle} where N`

## 7. Collisions

As said in the introduction, the most important part is a function that finds collisions between particles and obstacles. Here is how this looks for collision of a standard particle with a wall:

```using LinearAlgebra: dot, normalize

"""
collision(p::AbstractParticle, o::Obstacle) → t, cp
Find the collision (if any) between given particle and obstacle.
Return the time until collision and the estimated collision point `cp`.
"""
@inline function collision(p::Particle, w::Wall)
n = normalvec(w, p.pos)
denom = dot(p.vel, n)
if denom ≥ 0.0
return nocollision()
else
t = dot(w.sp - p.pos, n)/denom
return t, p.pos + t * p.vel
end
end

normalvec(w::Wall, pos) = w.normal```

and here is the same function but for collisions with disks instead:

```@inline function collision(p::Particle, d::Disk)
dotp = dot(p.vel, normalvec(d, p.pos))
dotp ≥ 0.0 && return nocollision()

dc = p.pos - d.c
B = dot(p.vel, dc)           #pointing towards circle center: B < 0
C = dot(dc, dc) - d.r*d.r    #being outside of circle: C > 0
Δ = B*B - C

Δ ≤ 0.0 && return nocollision()
sqrtD = sqrt(Δ)
# Closest point:
t = -B - sqrtD
return t, p.pos + t * p.vel
end

normalvec(d::Disk, pos) = normalize(pos - d.c)```

You can see that there are cases where collisions are not possible or they happen backwards in time. By convention, this is the value we return then:

`@inline nocollision() = (Inf, SV(0.0, 0.0))`

`next_collision` is a useful function that finds the "true" next collision. It simply loops over the obstacles in a billiard. It simply checks which obstacle has the least collision time:

```function next_collision(p::AbstractParticle, bd)
j, ct, cp = 0, Inf, SV(0.0, 0.0)
for i in eachindex(bd)
t, c = collision(p, bd[i])
if t < ct
j = i
ct = t
cp = c
end
end
return j, ct, cp
end```

## 8. Evolving a particle in a billiard

We need a simple function to propagate a particle to the found collision point. We will also a give the amount of time the propagation "should" take, for extendability. For standard particles, where the velocity vector is constant while travelling, this does not matter.

`propagate!(p::Particle, pos, t) = (p.pos = pos)`

This is an in-place function (notice `!` at the end)

We also have to define a function that performs specular reflection i.e. changes the velocity of the particle (after collision)

```function resolvecollision!(p::AbstractParticle, o::Obstacle)
n = normalvec(o, p.pos)
p.vel = p.vel - 2*dot(n, p.vel)*n
end```

We are now ready to wrap things up. Let's define a function that takes a particle and evolves in a billiard (tuple of obstacles) and returns the timeseries of the positions of the particle.

For convenience, it is worthwhile to define the following function:

```"""
bounce!(p, bd)
Evolve the particle for one collision (in-place).
"""
@inline function bounce!(p::AbstractParticle, bd)
i::Int, tmin::Float64, cp::SV = next_collision(p, bd)
if tmin != Inf
propagate!(p, cp, tmin)
resolvecollision!(p, bd[i])
end
return i, tmin, p.pos, p.vel
end```

Then, we can use it inside a bigger function that calls `bounce!` until a specified amount of time:

```"""
timeseries!(p::AbstractParticle, bd, n) -> xt, yt, t
Evolve the particle in the billiard `bd` for `n` collisions
and return the position timeseries `xt, yt` along with time vector `t`.
"""
function timeseries!(p::AbstractParticle, bd, n::Int)

t = [0.0]; xt = [p.pos]; yt = [p.pos]; c = 0

while c < n
prevpos = p.pos; prevvel = p.vel
i, ct = bounce!(p, bd)
xs, ys = extrapolate(p, prevpos, prevvel, ct)
push!(t, ct)
append!(xt, xs)
append!(yt, ys)
c += 1
end

return xt, yt, t
end```

`extrapolate` simply creates the position timeseries in between two collisions. For a standard particle there is no "extrapolation" needed, one just uses the final position:

`extrapolate(p::Particle, prevpos, prevvel, ct) = p.pos`

Why does this `extrapolate` function exist? See below when we extend our code for magnetic particles!

## 9. Running the code

let's put this to the test now! We'll create the famous Sinai billiard

```x, y, r = 1.0, 1.0, 0.3
sp = [0.0,y]; ep = [0.0, 0.0]; n = [x,0.0]
leftw = Wall(sp, ep, n)
sp = [x,0.0]; ep = [x, y]; n = [-x,0.0]
rightw = Wall(sp, ep, n)
sp = [x,y]; ep = [0.0, y]; n = [0.0,-y]
topw = Wall(sp, ep, n)
sp = [0.0,0.0]; ep = [x, 0.0]; n = [0.0,y]
botw = Wall(sp, ep, n)
disk = Disk([x/2, y/2], r)

bd = (botw, rightw, topw, leftw, disk)
bd isa Billiard```

and also initialize a particle

`p = Particle(0.1, 0.1, 2π*rand())`

and evolve it

`xt, yt, t = timeseries!(p, bd, 10)`

## 10. Plotting

Let's define some simple methods for plotting and plot the result!

```using PyPlot
import PyPlot: plot

const EDGECOLOR = (0,0.6,0)
function plot(d::Disk)
facecolor = (EDGECOLOR..., 0.5)
circle1 = PyPlot.plt[:Circle](d.c, d.r;
edgecolor = EDGECOLOR, facecolor = facecolor, lw = 2.0)
end
function plot(w::Wall)
PyPlot.plot([w.sp,w.ep],[w.sp,w.ep];
color=EDGECOLOR, lw = 2.0)
end
function plot(bd::Billiard)
for o ∈ bd; plot(o); end
gca()[:set_aspect]("equal")
end

figure(); plot(bd)
gcf()```

Awesome! Now let's see it with the orbit as well

```figure(); plot(bd)
xt, yt, t = timeseries!(p, bd, 10)
plot(xt, yt)
gcf()```

Plot 10 orbits!

```figure(); plot(bd)
p = Particle(0.1, 0.5, 2π*rand())
for j in 1:10
xt, yt = timeseries!(p, bd, 20)
plot(xt, yt, alpha = 0.5)
end
gcf()```

## 11. Showcase 1: Performance & Metaprogramming

```using BenchmarkTools

p = Particle(0.1, 0.1, 2π*rand())
@btime bounce!(\$p, \$bd)```

It is already very fast to propagate a particle for one collision, however there are some allocations (even though the function is in place).

These allocations come from type instability in `next_collision`, since the `Tuple` contains elements of different types. However, using metaprogramming it is easy to solve this type instability because `Tuple` has known size!

What we do in the following definition is using metaprogramming to "unroll" the loop

```@generated function next_collision(p::AbstractParticle, bd::Billiard)
L = length(bd.types) # notice that bd stands for the TYPE of bd here!
out = :(ind = 0; tmin = Inf; cp = SV(0.0, 0.0))
for j=1:L
push!(out.args, quote
let x = bd[\$j]
tcol, pcol = collision(p, x)
# Set minimum time:
if tcol < tmin
tmin = tcol
ind = \$j
cp = pcol
end
end
end
)
end
push!(out.args, :(return ind, tmin, cp))
return out
end

@btime bounce!(\$p, \$bd)```

This number is insane!!! Notice that this code is billiard agnostic! You could pass any tuple of obstacles and it would still be as performant!!! The time of `bounce!` scales linearly with the number of obstacles in the billiard.

## 12. Showcase 2: Extendability

Let's say we want to add one more obstacle to this "billiard package" we are making. Do you we have to re-write everything for it? Nope! In the end we only need to extend two methods ! Only two!

To show this let's create an ellipse as an obstacle, with semi-axes `a, b`

```struct Ellipse <: Obstacle
c::SV
a::Float64
b::Float64
end```

The methods we need to extend are only these:

```normalvec
collision```

Yes!!! Only two! So let's get to it! `normalvec` is pretty easy:

```function normalvec(e::Ellipse, pos)
x₀, y₀ = pos
h, k = e.c
return normalize(SV((x₀-h)/(e.a*e.a), (y₀-k)/(e.b*e.b)))
end

using LinearAlgebra: norm
function collision(p::Particle, e::Ellipse)
dotp = dot(p.vel, normalvec(e, p.pos))
dotp ≥ 0.0 && return nocollision()

a = e.a; b = e.b
pc = p.pos - e.c
μ = p.vel/p.vel
ψ = pc - μ*pc

denomin = a*a*μ*μ + b*b
Δ² = denomin - ψ*ψ
Δ² ≤ 0 && return nocollision()
Δ = sqrt(Δ²); f1 = -a*a*μ*ψ; f2 = b*b*ψ # just factors
I1 = SV(f1 + a*b*Δ, f2 + a*b*μ*Δ)/denomin
I2 = SV(f1 - a*b*Δ, f2 - a*b*μ*Δ)/denomin

d1 = norm(pc - I1); d2 = norm(pc - I2)
return d1 < d2 ? (d1, I1 + e.c) : (d2, I2 + e.c)
end```

Alright so now let's create a billiard with both an ellipse and a disk, for the fun of it

```el = Ellipse([0.4, 0.2 ], 0.3, 0.1)
di = Disk([0.6, 0.7], 0.25)

bd2 = Billiard((bd[1:4]..., el, di))```

and plot it

```function plot(e::Ellipse)
facecolor = (EDGECOLOR..., 0.5)
ellipse = PyPlot.matplotlib[:patches][:Ellipse](e.c, 2e.a, 2e.b;
edgecolor = EDGECOLOR, facecolor = facecolor, lw = 2.0)
end

figure(); plot(bd2)
gcf()```

We are now ready to evolve a particle in this brand new billiard:

```p = Particle(0.1, 0.1, 2π*rand())
xt, yt, t = timeseries!(p, bd2, 20)
figure(); plot(bd2)
plot(xt, yt)
gcf()```

plot a bunch more!

```figure(); plot(bd2)
for j in 1:10
p = Particle(0.1, 0.1, 2π*rand())
xt, yt = timeseries!(p, bd2, 20)
plot(xt, yt, alpha = 0.5)
end
gcf()```

## 13. Showcase 3: Extendability, again.

Alright, so it turned out to be almost trivial to add an extra obstacle to our code. But what about an extra particle?

I am not talking about one more instance of `Particle`. I am talking about a new type of particle, that moves around in a different way.

In this part we will create this new type, `MagneticParticle` that moves around in circles instead of straight lines! But how many functions do we need to define? Provided you have already defined the type `MagneticParticle`, then that many:

```collision # for each obstacle we want to support
propagate!
extrapolate```

and yeap, that's it. It may be hard to believe that it only takes so little, but it's true!!!

### 13.1. The type

```mutable struct MagneticParticle <: AbstractParticle
pos::SV
vel::SV
ω::Float64
end
MagneticParticle(x0, y0, φ0, ω) = MagneticParticle(SV(x0, y0), SV(cos(φ0), sin(φ0)), ω)```

This particle moves in circles with angular velocity `ω`.

### 13.2. Extending `collision`

To extend `collision`, we simply have to find intersections of circle-line and circle-circle, for collisions with `Wall` and `Disk`. I won't go into details of how to do this, and instead I'll copy-paste functions from `DynamicalBilliards`. The versions in `DynamicalBilliards` also have a lot of comments that explain what is going on.

Here is the collision with wall:

```function collision(p::MagneticParticle, w::Wall)
ω = p.ω
pc, pr = cyclotron(p)
P0 = p.pos
P2P1 = w.ep - w.sp
P1P3 = w.sp - pc

a = dot(P2P1, P2P1)
b = 2*dot(P2P1, P1P3)
c = dot(P1P3, P1P3) - pr*pr
Δ = b^2 -4*a*c

Δ ≤ 0.0 && return nocollision()
u1 = (-b - sqrt(Δ))/2a
u2 = (-b + sqrt(Δ))/2a
cond1 = 0.0 ≤ u1 ≤ 1.0
cond2 = 0.0 ≤ u2 ≤ 1.0
θ, I = nocollision()
if cond1 || cond2
dw = w.ep - w.sp
for (u, cond) in ((u1, cond1), (u2, cond2))
Y =  w.sp + u*dw
if cond
φ = realangle(p, w, Y)
φ < θ && (θ = φ; I = Y)
end
end
end
return θ*pr, I
end```

and here is the collision with a disk:

```function collision(p::MagneticParticle, o::Disk)
ω = p.ω
pc, rc = cyclotron(p)
p1 = o.c
r1 = o.r
d = norm(p1-pc)
if (d >= rc + r1) || (d <= abs(rc-r1))
return nocollision()
end

a = (rc^2 - r1^2 + d^2)/2d
h = sqrt(rc^2 - a^2)

I1 = SV(
pc + a*(p1 - pc)/d + h*(p1 - pc)/d,
pc + a*(p1 - pc)/d - h*(p1 - pc)/d
)
I2 = SV(
pc + a*(p1 - pc)/d - h*(p1 - pc)/d,
pc + a*(p1 - pc)/d + h*(p1 - pc)/d
)

θ1 = realangle(p, o, I1)
θ2 = realangle(p, o, I2)

return θ1 < θ2 ? (θ1*rc, I1) : (θ2*rc, I2)
end```

The functions `cyclotron` and `realangle` are helper functions. The first one finds the center and radius of the cyclotron traced by the particle.

`cyclotron(p) = (p.pos - (1/p.ω)*SV(p.vel, -p.vel), abs(1/p.ω))`

`realangle` has a simple purpose: the intersections of a circle with any obstacle are always 2. But which one happens first, from a temporal perspective? `realangle` gives the correct angle until the collision point, in forward time.

```function realangle(p::MagneticParticle, o::Obstacle, i)

pc, pr = cyclotron(p); ω = p.ω
P0 = p.pos
PC = pc - P0
d2 = dot(i-P0,i-P0)

if d2 ≤ 1e-8
dotp = dot(p.vel, normalvec(o,  p.pos))
dotp ≥ 0 && return Inf
end

d2r = (d2/(2pr^2))
d2r > 2 && (d2r = 2.0)
θprime = acos(1.0 - d2r)

PI = i - P0
side = (PI*PC - PI*PC)*ω

side < 0 && (θprime = 2π-θprime)
return θprime
end```

The complexity of the functions `collision` and `realangle` exists solely due to the geometry of intersections between circles. What we want to point out is how few methods we have to extend. How easy is defining these new methods is not relevant, blame math and physics for that! So don't be taken aback because these functions are "long"!

### 13.3. Propagation & extrapolation

`propagate!` for a `MagneticParticle` must evolve it in an arc of a circle, so as you can see we have to change the velocity vector!

```function propagate!(p::MagneticParticle, pos, t)
φ0 = atan(p.vel, p.vel)
p.pos = pos
p.vel = SV(cossin(p.ω*t + φ0))
return
end

cossin(x) = ( (y, z) = sincos(x); (z, y) )```

`extrapolate` should simply create the arc that connects the previous point with the current one

```function extrapolate(p::MagneticParticle, prevpos, prevvel, t)
φ0 = atan(prevvel, prevvel)
s0, c0 = sincos(φ0)
x0 = prevpos; y0 = prevpos
xt = [x0]; yt = [y0]; ω = p.ω
tvec = 0.0:0.01:t
for td in tvec
s, c = sincos(p.ω*td + φ0)
push!(xt, s/ω + x0 - s0/ω)
push!(yt, -c/ω + y0 + c0/ω) #vx0 is cos(φ0)
end
return xt, yt
end```

### 13.4. Evolve the magnetic particle

```p = MagneticParticle(0.1, 0.1, 2π*rand(), 3.0)
xt, yt, t = timeseries!(p, bd, 20)
figure(); plot(bd)
plot(xt, yt)
gcf()```

plot a bunch of these!

```figure(); plot(bd)
for j in 1:4
p = MagneticParticle(0.1, 0.1, 2π*rand(), 2.0)
xt, yt = timeseries!(p, bd, 20)
plot(xt, yt, alpha = 0.5)
end
gcf()```