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}

1.1s

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)))

0.0s

Then, we define the obstacle structs. 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

0.0s

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

0.4s

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

3.0s

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)

0.0s

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))

0.2s

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

0.2s

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)

0.2s

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

0.3s

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

0.5s

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[1]]; yt = [p.pos[2]]; 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

0.1s

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

0.3s

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

0.8s

and also initialize a particle

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

1.1s

and evolve it

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

2.0s

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)
    PyPlot.gca()[:add_artist](circle1)
end
function plot(w::Wall)
    PyPlot.plot([w.sp[1],w.ep[1]],[w.sp[2],w.ep[2]];
        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()

20.7s

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()

0.9s

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()

1.1s

11.
Showcase 1: Performance & Metaprogramming

using BenchmarkTools

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

6.6s

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)

4.5s

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

0.0s

The methods we need to extend are only these:

normalvec
collision

0.7s

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[2]/p.vel[1]
    ψ = pc[2] - μ*pc[1]

    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

0.2s

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))

0.7s

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)
    PyPlot.gca()[:add_artist](ellipse)
end

figure(); plot(bd2)
gcf()

0.5s

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()

1.1s

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()

0.9s

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

0.1s

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)), ω)

0.4s

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

0.2s

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[1] + a*(p1[1] - pc[1])/d + h*(p1[2] - pc[2])/d,
        pc[2] + a*(p1[2] - pc[2])/d - h*(p1[1] - pc[1])/d
    )
    I2 = SV(
        pc[1] + a*(p1[1] - pc[1])/d - h*(p1[2] - pc[2])/d,
        pc[2] + a*(p1[2] - pc[2])/d + h*(p1[1] - pc[1])/d
    )

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

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

0.2s

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[2], -p.vel[1]), abs(1/p.ω))

0.4s

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[1]*PC[2] - PI[2]*PC[1])*ω

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

0.4s

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[2], p.vel[1])
    p.pos = pos
    p.vel = SV(cossin(p.ω*t + φ0))
    return
end

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

0.4s

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[2], prevvel[1])
    s0, c0 = sincos(φ0)
    x0 = prevpos[1]; y0 = prevpos[2]
    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

0.2s

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()

1.2s

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()

0.8s

Simulating billiards in Julia

1. Author

2. This tutorial

3. Features of Julia highlighted

4. The system and algorithm

5. The Algorithm

6. The basic structs we will need

7. Collisions

8. Evolving a particle in a billiard

9. Running the code

10. Plotting

11. Showcase 1: Performance & Metaprogramming

12. Showcase 2: Extendability

13. Showcase 3: Extendability, again.

13.1. The type

13.2. Extending collision

13.3. Propagation & extrapolation

13.4. Evolve the magnetic particle

1.
Author

2.
This tutorial

3.
Features of Julia highlighted

4.
The system and algorithm

5.
The Algorithm

6.
The basic `struct`s we will need

7.
Collisions

8.
Evolving a particle in a billiard

9.
Running the code

10.
Plotting

11.
Showcase 1: Performance & Metaprogramming

12.
Showcase 2: Extendability

13.
Showcase 3: Extendability, again.

13.1.
The type

13.2.
Extending `collision`

13.3.
Propagation & extrapolation

13.4.
Evolve the magnetic particle