

High Level API

DynamicalBilliards was created with ease-of-use as its main cornerstone. With 3 simple steps, the user can get the output of the propagation of a particle in a billiard.

In general, the workflow of DynamicalBilliards follows these simple steps:

1. Create a billiard.
2. Create particles inside that billiard.
3. Get the output you want by using one of the high level functions.

Adding more complexity in your billiard does not add complexity in your code. For example, to implement a ray-splitting billiard you only need to define one additional variable, a RaySplitter and pass it to the high level functions.

After reading through this page, you will be able to use almost all aspects of DynamicalBilliards with minimal effort.

Visualizations

Visualizing the billiards, particles and their motion is one of the most important parts of the DynamicalBilliards. It is not discussed in this page however, but rather in the Visualizing page.

---

Billiard

A Billiard is simply a collection of Obstacle subtypes. Particles are propagating inside a Billiard, bouncing from obstacle to obstacle while having constant velocity in-between.

There is a tutorial on how to create your own billiard. In addition, there are many pre-defined billiards that can be found in the Standard Billiards Library section. That is why knowing how to construct a Billiard is not important at this point.

In this page we will be using the Bunimovich billiard as an example:

using DynamicalBilliards
bd = billiard_bunimovich()

Billiard{Float64} with 4 obstacles:
  Bottom wall
  Right semicircle
  Top wall
  Left semicircle

Particles

A "particle" is that thingy that moves around in the billiard. It always moves with velocity of measure 1, by convention.

Currently there are two types of particles:

• Particle, which propagates as a straight line.
• MagneticParticle, which propagates as a circle instead of a line (similar to electrons in a perpendicular magnetic field).

There are two ways to create a particle. The first one is to provide the constructor with some initial conditions:

x0 = rand(); y0 = rand();
φ0 = 2π*rand()
p = Particle(x0, y0, φ0)

Particle{Float64}
position: [0.132251, 0.444436]
velocity: [0.746217, -0.665703]

To create a MagneticParticle simply provide the constructor with one more number, the angular velocity:

ω = 0.5
mp = MagneticParticle(x0, y0, φ0, ω)

MagneticParticle{Float64}
position: [0.132251, 0.444436]
velocity: [0.746217, -0.665703]
ang. velocity: 0.5

Why the {Float64} ?

When creating a billiard or a particle, the object is printed with {Float64} at the end. This shows what type of numbers are used for all numerical operations. If you are curious you can learn more about it in the numerical precision page.

Particles must be inside the Billiard!

Keep in mind that the particle must be initialized inside a billiard for any functionality to work properly and make sense. If you are not sure what we mean by that, then you should check out the Internals page.

Random initial conditions

If you have a Billiard which is not a rectangle, creating many random initial conditions inside it can be a pain. Fortunately, the second way to create a particle is to use the following function:

# DynamicalBilliards.randominside — Function.

randominside(bd::Billiard [, ω])

Return a particle with random allowed initial conditions inside the given billiard. If supplied with a second argument the type of the returned particle is MagneticParticle, with angular velocity ω.

source

---

For example:

p = randominside(bd)

Particle{Float64}
position: [-0.358445, 0.23933]
velocity: [-0.944648, -0.328084]

and

mp = randominside(bd, ω)

MagneticParticle{Float64}
position: [0.735833, 0.00240983]
velocity: [-0.690733, 0.72311]
ang. velocity: 0.5

randominside always creates particles with same number type as the billiard.

evolve & timeseries

Now that we have created a billiard and a particle inside, we want to evolve it! There is a simple function for that, called evolve! (or evolve if you don't want to mutate the particle), which returns the time, position and velocities at the collision points:

# DynamicalBilliards.evolve! — Function.

evolve!([p::AbstractParticle,] bd::Billiard, t)

Evolve the given particle p inside the billiard bd. If t is of type AbstractFloat, evolve for as much time as t. If however t is of type Int, evolve for as many collisions as t. Return the states of the particle between collisions.

This function mutates the particle, use evolve otherwise. If a particle is not given, a random one is picked through randominside.

Return

• ct::Vector{T} : Collision times.
• poss::Vector{SVector{2,T}} : Positions at the collisions.
• vels::Vector{SVector{2,T}}) : Velocities exactly after the collisions.
• ω, either T or Vector{T} : Angular velocity/ies (returned only for magnetic particles).

The time ct[i+1] is the time necessary to reach state poss[i+1], vels[i+1] starting from the state poss[i], vels[i]. That is why ct[1] is always 0 since poss[1], vels[1] are the initial conditions. The angular velocity ω[i] is the one the particle has while propagating from state poss[i], vels[i] to i+1.

Notice that at any point, the velocity vector vels[i] is the one obdained after the specular reflection of the i-1th collision.

Ray-splitting billiards

evolve!(p, bd, t, raysplitters)

To implement ray-splitting, the evolve! function is supplemented with a fourth argument, raysplitters which is a tuple of RaySplitter instances. Notice that evolve always mutates the billiard if ray-splitting is used! For more information and instructions on using ray-splitting please visit the official documentation.

source

---

Forget the ray-splitting part for now (see Ray-Splitting).

Let's see an example:

ct, poss, vels = evolve(p, bd, 100)
for i in 1:5
  println(round(ct[i], digits=3), "  ", poss[i], "  ", vels[i])
end

0.0  [-0.358445, 0.23933]  [-0.944648, -0.328084]
0.059  [-0.414219, 0.219959]  [0.65645, 0.75437]
1.034  [0.26457, 1.0]  [0.65645, -0.75437]
1.306  [1.12169, 0.0150333]  [0.222545, 0.974922]
0.891  [1.32007, 0.884127]  [-0.91877, -0.394794]

Similarly, for magnetic propagation

ct, poss, vels, ω = evolve(mp, bd, 100)
for i in 1:10
  println(round(ct[i], digits=3), "  ", poss[i], "  ", vels[i])
end

0.0  [0.735833, 0.00240983]  [-0.690733, 0.72311]
1.394  [-0.488075, 0.608549]  [0.947236, -0.320537]
1.903  [1.32345, 0.881284]  [-0.445214, -0.895424]
0.9  [1.11443, 0.01327]  [-0.455572, 0.890199]
1.55  [-0.0319696, 0.998977]  [-0.900216, -0.435443]
0.488  [-0.441251, 0.735155]  [-0.102913, -0.99469]
0.488  [-0.432027, 0.2483]  [0.791998, -0.610523]
0.488  [-0.0130471, 0.000170256]  [0.935738, 0.352695]
1.512  [1.0786, 0.993783]  [0.138244, -0.990398]
0.877  [1.38335, 0.179001]  [-0.920846, 0.389926]

The above return types are helpful in some applications. In other applications however one prefers to have the time series of the individual variables. For this, the timeseries function is used:

# DynamicalBilliards.timeseries! — Function.

timeseries!([p::AbstractParticle,] bd::Billiard, t; dt, warning)

Evolves the given particle p inside the billiard bd. If t is of type AbstractFloat, evolve for as much time as t. If however t is of type Int, evolve for as many collisions as t. Returns the time series for position and velocity as well as the time vector.

This function mutates the particle, use timeseries otherwise. If a particle is not given, a random one is picked through randominside.

The keyword argument dt is the time step used for interpolating the time series in between collisions. dt is capped by the collision time, as the interpolation always stops at collisions. For straight propagation dt = Inf, while for magnetic dt = 0.01.

For pinned magnetic particles, timeseries! issues a warning and returns the trajectory of the particle. If t is integer, the trajectory is evolved for one full circle only

Return:

• x position time-series
• y position time-series
• x velocity time-series
• y velocity time-series
• time vector

Ray-splitting billiards

timeseries!(p, bd, t, raysplitters; ...)

To implement ray-splitting, the timeseries! function is supplemented with a fourth argument, raysplitters which is a tuple of RaySplitter instances. Notice that timeseries always mutates the billiard if ray-splitting is used! For more information and instructions on using ray-splitting please visit the official documentation.

source

---

For example:

xt, yt, vxt, vyt, t = timeseries(p, bd, 100)

# print as a matrix:
hcat(xt, yt, vxt, vyt, t)[1:5, :]

5×5 Array{Float64,2}:
 -0.358445  0.23933    -0.944648  -0.328084  0.0
 -0.414219  0.219959    0.65645    0.75437   0.0590417
  0.26457   1.0         0.65645   -0.75437   1.03403
  1.12169   0.0150333   0.222545   0.974922  1.30568
  1.32007   0.884127   -0.91877   -0.394794  0.891449

Same story for magnetic particles:

# evolve the magnetic particle instead:
xt, yt, vxt, vyt, t = timeseries(mp, bd, 100)

# print as a matrix:
hcat(xt, yt, vxt, vyt, t)[1:5, :]

5×5 Array{Float64,2}:
 0.735833  0.00240983  -0.690733  0.72311   0.0
 0.728907  0.00962363  -0.69434   0.719647  0.01
 0.721946  0.0168027   -0.69793   0.716166  0.02
 0.714949  0.0239469   -0.701502  0.712668  0.03
 0.707916  0.031056    -0.705056  0.709151  0.04

Type of t

Remember that the behavior of evolve! depends on the type of the third argument, which represents "total amount". If it is AbstractFloat, it represents total amount of time, but if it is Int it represents total number of collisions.

Poincaré Sections

# DynamicalBilliards.psos — Function.

psos(bd::Billiard, plane::InfiniteWall, t, particles)

Compute the Poincaré section of the particles with the given plane, by evolving each one for time t (either integer or float) inside bd.

The plane can be an InfiniteWall of any orientation, however only crossings of the plane such that dot(velocity, normal) < 0 are allowed, with normal the normal unit vector of the plane.

particles can be:

• A single particle.
• A Vector of particles.
• An integer n optionally followed by an angular velocity ω.

Return the positions poss and velocities vels at the instances of crossing the plane. If given more than one particle, the result is a vector of vectors of vectors.

Notice - This function can handle pinned particles. If a pinned particle can intersect with the plane, then an intersection is returned. If however it can't then empty vectors are returned.

source

---

For example, the surface of section in the periodic Sinai billiard with magnetic field reveals the mixed nature of the phase-space:

using DynamicalBilliards, PyPlot
t = 100; r = 0.15
bd = billiard_sinai(r, setting = "periodic")

# the direction of the normal vector is IMPORTANT!!!
# (always keep in mind that ω > 0  means counter-clockwise rotation!)
plane = InfiniteWall([0.5, 0.0], [0.5, 1.0], [-1.0, 0.0])

posvector, velvector = psos(bd, plane, t, 1000, 2.0)
c(a) = length(a) == 1 ? "C1" : "C0"

figure()
for i in 1:length(posvector)
    poss = posvector[i] # vector of positions
    vels = velvector[i] # vector of velocities at the section
    L = length(poss)
    if L > 0
        #plot y vs vy
        y = [a[2] for a in poss]
        vy = [a[2] for a in vels]

        plot(y, vy, ls = "None", color = c(y), ms = 2.0, alpha = 0.75, marker = "o")
    end
end
xlabel("\$y\$"); ylabel("\$v_y\$")

psos operates on the unit cell

The psos function always calculates the crossings within the unit cell of a periodic billiard. This means that no information about the "actual" position of the particle is stored, everything is modulo the unit cell.

This can be seen very well in the above example, where there aren't any entries in the region 0.5 - r ≤ y ≤ 0.5 + r.

Of course it is very easy to "re-normalize" the result such that it is coherent. The only change we have to do is simply replace the line y = [a[2] for a in poss] with

y = [a[2] < 0.5 ? a[2] + 1 : a[2]  for a in poss]

Escape Times

It is very easy to create your own function that calculates an "escape time": the time until the particle leaves the billiard by meeting a specified condition. There is also a high-level function for this though:

# DynamicalBilliards.escapetime — Function.

escapetime([p,] bd, t; warning = false)

Calculate the escape time of a particle p in the billiard bd, which is the time until colliding with any "door" in bd. As a "door" is considered any FiniteWall with field isdoor = true.

If the particle evolves for more than t (integer or float) without colliding with the Door (i.e. escaping) the returned result is Inf.

A warning can be thrown if the result is Inf. Enable this using the keyword warning = true.

If a particle is not given, a random one is picked through randominside. See parallelize for a parallelized version.

source

Creating a door

To create a "door" simply visit the library page to learn more about the individual obstacle types (specifically FiniteWall). To be able to combine them into a Billiard you should also check out the tutorial on defining your own billiard.

For example, the default implementation of the mushroom billiard has a "door" at the bottom of the stem. Thus,

using Statistics
bd = billiard_mushroom()
et = zeros(100)
for i ∈ 1:100
    particle = randominside(bd)
    et[i] = escapetime(particle, bd, 10000)
end
println("Out of 100 particles, $(count(x-> x != Inf, et)) escaped")
println("Mean escape time was $(mean(et[et .!= Inf]))")

Out of 100 particles, 31 escaped
Mean escape time was 4.393499313417378

Of course, escapetime works with MagneticParticle as well

escapetime(randominside(bd, 1.0), bd, 10000)

15.722571714230952

Mean Collision Times

Because the computation of a mean collision time (average time between collisions in a billiard) is often a useful quantity, the following function computes it

# DynamicalBilliards.meancollisiontime — Function.

meancollisiontime([p,] bd, t) → κ

Compute the mean collision time κ of the particle p in the billiard bd by evolving for total amount t (either float for time or integer for collision number).

Collision times are counted only between obstacles that are not PeriodicWall.

If a particle is not given, a random one is picked through randominside. See parallelize for a parallelized version.

source

For example,

bd = billiard_sinai()
meancollisiontime(randominside(bd), bd, 10000.0)

0.45209982482588307

Parallelization

# DynamicalBilliards.parallelize — Function.

parallelize(f, bd::Billiard, t, particles; partype = :threads)

Parallelize function f across the available particles. The parallelization type can be :threads or :pmap, which use threads or a worker pool initialized with addprocs before using DynamicalBilliards.

particles can be:

• A Vector of particles.
• An integer n optionally followed by an angular velocity ω. This uses randominside.

The functions usable here are:

• meancollisiontime
• escapetime
• lyapunovspectrum (returns only the maximal exponents)
• boundarymap (returns vector of vectors of 2-vectors and arcintervals)

source

---

Here are some examples

bd = billiard_stadium()
particles = [randominside(bd) for i in 1:1000]
parallelize(meancollisiontime, bd, 1000, particles)

1000-element Array{Float64,1}:
 1.0866577108051312
 1.1318093170582553
 1.104723060808063
 1.0693283722486462
 1.0887855833108986
 1.1104116928030154
 1.1180584910261002
 1.123470428953834
 1.1261543227129203
 1.086978801368143
 ⋮
 1.1066030472561408
 1.1287157765332596
 1.096005106667699
 1.0571619057556725
 1.0974733124153602
 1.0116652573929137
 1.068258345000985
 1.0731689707286707
 1.1084525354439871

parallelize(lyapunovspectrum, bd, 1000, particles)

1000-element Array{Float64,1}:
 0.7860980007045091
 0.8921790633011215
 0.9054480972191328
 0.919862143145693
 0.8796603015697422
 0.9036069882418831
 0.8826732135341996
 0.8731006583458032
 0.8546467687243233
 0.8999362331175333
 ⋮
 0.9047251857760783
 0.8920523659932186
 0.8634399465537317
 0.8467094836924418
 0.7463858607107943
 0.7931415672873441
 0.8519561989567872
 0.8427253972768933
 0.8811405537402082

It's all about bounce!

The main propagation algorithm used by DynamicalBilliards is bundled in the following well-behaving function:

# DynamicalBilliards.bounce! — Function.

bounce!(p::AbstractParticle, bd::Billiard) → i, t, pos, vel

"Bounce" the particle (advance for one collision) in the billiard. Takes care of finite-precision issues.

Return:

• index of the obstacle that the particle just collided with
• the time from the previous collision until the current collision t
• position and velocity of the particle at the current collision (after the collision has been resolved!). The position is given in the unit cell of periodic billiards. Do pos += p.current_cell for the position in real space.

bounce!(p, bd, raysplit) → i, t, pos, vel

Ray-splitting version of bounce!.

source

---

bounce! is the function used internally by all high-level functions, like evolve!, boundarymap, escapetime, etc.

This is the function a user should use if they want to calculate other things besides what is already available in the high level API.

Standard Billiards Library

You can also use keywords!

All standard billiards have a function version that accepts keyword arguments instead of positional arguments, for ease of use.

# DynamicalBilliards.billiard_rectangle — Function.

billiard_rectangle(x=1.0, y=1.0; setting = "standard")

Return a vector of obstacles that defines a rectangle billiard of size (x, y).

Settings

• "standard" : Specular reflection occurs during collision.
• "periodic" : The walls are PeriodicWall type, enforcing periodicity at the boundaries
• "random" : The velocity is randomized upon collision.
• "ray-splitting" : All obstacles in the billiard allow for ray-splitting.

source

# DynamicalBilliards.billiard_sinai — Function.

billiard_sinai(r=0.25, x=1.0, y=1.0; setting = "standard")

Return a vector of obstacles that defines a Sinai billiard of size (x, y) with a disk in its center, of radius r.

In the periodic case, the system is also known as "Lorentz Gas".

Settings

• "standard" : Specular reflection occurs during collision.
• "periodic" : The walls are PeriodicWall type, enforcing periodicity at the boundaries
• "random" : The velocity is randomized upon collision.
• "ray-splitting" : All obstacles in the billiard allow for ray-splitting.

source

# DynamicalBilliards.billiard_bunimovich — Function.

billiard_bunimovich(l=1.0, w=1.0)

Return a vector of Obstacles that define a Buminovich billiard, also called a stadium. The length is considered without the attached semicircles, meaning that the full length of the billiard is l + w. The left and right edges of the stadium are Semicircles.

billiard_stadium is an alias of billiard_bunimovich.

source

# DynamicalBilliards.billiard_mushroom — Function.

billiard_mushroom(sl = 1.0, sw = 0.2, cr = 1.0, sloc = 0.0; door = true)

Create a mushroom billiard with stem length sl, stem width sw and cap radius cr. The center of the cap (which is Semicircle) is always at [0, sl]. The center of the stem is located at sloc.

Optionally, the bottom-most Wall is a Door (see escapetime).

source

# DynamicalBilliards.billiard_polygon — Function.

billiard_polygon(n::Int, R, center = [0,0]; setting = "standard")

Return a vector of obstacles that defines a regular-polygonal billiard with n sides, radius r and given center.

Note: R denotes the so-called outer radius, not the inner one.

Settings

• "standard" : Specular reflection occurs during collision.
• "periodic" : The walls are PeriodicWall type, enforcing periodicity at the boundaries. Only available for n=4 or n=6.
• "random" : The velocity is randomized upon collision.

source

# DynamicalBilliards.billiard_hexagonal_sinai — Function.

billiard_hexagonal_sinai(r, R, center = [0,0]; setting = "standard")

Create a sinai-like billiard, which is a hexagon of outer radius R, containing at its center (given by center) a disk of radius r. The setting keyword is passed to billiard_polygon.

source

# DynamicalBilliards.billiard_raysplitting_showcase — Function.

billiard_raysplitting_showcase(x=2.0, y=1.0, r1=0.3, r2=0.2) -> bd, rayspl

Showcase example billiard for ray-splitting processes. A rectangle (x,y) with a SplitterWall at x/2 and two disks at each side, with respective radii r1, r2.

Notice: This function returns a billiard bd as well as a rayspl dictionary!

source

# DynamicalBilliards.billiard_logo — Function.

billiard_logo(;h=1.0, α=0.8, r=0.18, off=0.25) -> bd, ray

Create the billiard used as logo of DynamicalBilliards and return it along with the tuple of raysplitters.

source

# DynamicalBilliards.billiard_iris — Function.

billiard_iris(a=0.2, b=0.4, w=1.0; setting = "standard")

Return a billiard that is a square of side w enclosing at its center an ellipse with semi axes a, b.

source

Particle types

# DynamicalBilliards.Particle — Type.

Particle(ic::Vector{T}) #where ic = [x0, y0, φ0]
Particle(x0, y0, φ0)
Particle(pos::SVector, vel::SVector)

Create a particle with initial conditions x0, y0, φ0. It propagates as a straight line.

The field current_cell shows at which cell of a periodic billiard is the particle currently located.

source

# DynamicalBilliards.MagneticParticle — Type.

MagneticParticle(ic::AbstractVector{T}, ω::Real) # where ic = [x0, y0, φ0]
MagneticParticle(x0, y0, φ0, ω)
MagneticParticle(pos::SVector, vel::SVector, ω)
MagneticParticle(p::AbstractParticle, ω)

Create a magnetic particle with initial conditions x0, y0, φ0 and angular velocity ω. It propagates as a circle instead of a line, with radius 1/abs(ω).

The field current_cell shows at which cell of a periodic billiard is the particle currently located.

source