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 & Animating 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).

To make a particle, provide the constructor with some initial conditions:

x0 = rand(); y0 = rand();
φ0 = 2π*rand()
p = Particle(x0, y0, φ0)
Particle{Float64}
position: [0.8351881435134845, 0.8192446867530516]
velocity: [0.39923985891529484, 0.916846516628217]

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.8351881435134845, 0.8192446867530516]
velocity: [0.39923985891529484, 0.916846516628217]
ang. velocity: 0.5

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.

You can initialize several particles with the same direction but slightly different position is the following function:

DynamicalBilliards.particlebeamFunction
particlebeam(x0, y0, φ, N, dx, ω = nothing, T = eltype(x0)) → ps

Make N particles, all with direction φ, starting at x0, y0. The particles don't all have the same position, but are instead spread by up to dx in the direction normal to φ.

The particle element type is T.

source
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, we have the following function:

DynamicalBilliards.randominsideFunction
randominside(bd::Billiard [, ω]) → p

Return a particle p 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 ω.

WARNING : randominside works for any convex billiard but it does not work for non-convex billiards. (this is because it uses distance)

source

For example:

p = randominside(bd)
Particle{Float64}
position: [1.1175298699463314, 0.9116367144893559]
velocity: [-0.8420827488356939, 0.5393483513586759]

and

mp = randominside(bd, ω)
MagneticParticle{Float64}
position: [0.3013316474465286, 0.39466733078013116]
velocity: [-0.5816498723205017, -0.8134392577381202]
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  [1.1175298699463314, 0.9116367144893559]  [-0.8420827488356939, 0.5393483513586759]
0.164  [0.9795685724444261, 1.0]  [-0.8420827488356939, -0.5393483513586759]
1.583  [-0.35336906842675897, 0.14626238328500563]  [0.5384702613023018, 0.842644514426476]
1.013  [0.1921899540146541, 1.0]  [0.5384702613023018, -0.842644514426476]
1.187  [0.8312141831924794, 1.1102230246251565e-16]  [0.5384702613023018, 0.842644514426476]

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.3013316474465286, 0.39466733078013116]  [-0.5816498723205017, -0.8134392577381202]
0.453  [0.08180711581674549, 0.0]  [-0.3843162069304361, 0.9232015235530118]
0.909  [-0.44276600449202685, 0.7322891845656324]  [0.9703118281340256, -0.24185730541624018]
1.623  [1.1159312554457164, 0.986374283870745]  [0.5101447815498561, -0.8600885430334776]
0.607  [1.4996481438592184, 0.5187545284674471]  [-0.6917671909604265, -0.7221205948529107]
0.607  [1.1520662464176716, 0.0236851286171727]  [-0.8811896484766224, 0.4727629463247375]
1.579  [-0.37886291334868244, 0.1737134803750716]  [0.43126722147378904, 0.9022242424598654]
0.82  [-0.18459237148335833, 0.964678013671994]  [0.7117987890916097, -0.702383430789564]
1.829  [1.491223220885588, 0.40672756429261414]  [-0.872142951220924, 0.48925113452669433]
2.071  [-0.4885786796155667, 0.393741476459018]  [0.9953908376225877, -0.09590140967161626]

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)

Evolve the given particle p inside the billiard bd for the condition t and return the x, y, vx, vy timeseries and the time vector. If t is of type AbstractFloat, then evolve for as much time as t. If however t is of type Int, evolve for as many collisions as t. Otherwise, t can be any function, that takes as an input t(n, τ, i, p) and returns true when the evolution should terminate. Here n is the amount of obstacles collided with so far, τ the amount time evolved so far, i the obstacle just collided with and p the particle (so you can access e.g. p.pos).

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.

Internally uses DynamicalBilliards.extrapolate.

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 Matrix{Float64}:
  1.11753   0.911637     -0.842083   0.539348  0.0
  0.979569  1.0          -0.842083  -0.539348  0.163833
 -0.353369  0.146262      0.53847    0.842645  1.74674
  0.19219   1.0           0.53847   -0.842645  2.7599
  0.831214  1.11022e-16   0.53847    0.842645  3.94664

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 Matrix{Float64}:
 0.301332  0.394667  -0.58165   -0.813439  0.0
 0.295536  0.386518  -0.577575  -0.816337  0.01
 0.28978   0.378341  -0.573487  -0.819215  0.02
 0.284066  0.370134  -0.569383  -0.822072  0.03
 0.278393  0.361899  -0.565266  -0.824909  0.04

Sometimes we may need information about which obstacles a particle visited, in which sequence, and when. For this we have the following function:

DynamicalBilliards.visited_obstacles!Function
visited_obstacles!([p::AbstractParticle,] bd::Billiard, t)

Evolve the given particle p inside the billiard bd exactly like evolve!. However return only:

  • ts::Vector{T} : Vector of time points of when each collision occured.
  • obst::Vector{Int} : Vector of obstacle indices in bd that the particle collided with at the time points in ts.

The first entries are 0.0 and 0. Similarly with evolve! the function does not record collisions with periodic walls.

Currently does not support raysplitting. Returns empty arrays for pinned particles.

source
Type of `t`

Remember that the behavior of time evolution depends on the type of the t 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.psosFunction
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, CairoMakie
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 ? "#6D44D0" : "#DA5210"

fig = Figure(); ax = Axis(fig[1,1]; xlabel = L"y", ylabel=L"v_y")

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]
        scatter!(y, vy; color = c(y), markersize = 2)
    end
end
fig
Example block output
`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.escapetimeFunction
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

To create a "door" simply use FiniteWall.

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, 23 escaped
Mean escape time was 6.806621169997773

Of course, escapetime works with MagneticParticle as well

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

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.meancollisiontimeFunction
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.455705047712392

Parallelization

DynamicalBilliards.parallelizeFunction
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:

source

Here are some examples

bd = billiard_stadium()
particles = [randominside(bd) for i in 1:1000]
parallelize(meancollisiontime, bd, 1000, particles)
1000-element Vector{Float64}:
 1.0911931050960788
 1.1143498305587527
 1.0804095934233868
 1.079245473144379
 1.1391768177241426
 1.0375691058152232
 1.119758441109758
 1.074341893002397
 1.1033671241036371
 1.0455108393335315
 ⋮
 1.0950349703130853
 1.081499611829668
 1.065281645515326
 1.0496785538443612
 1.139837366656638
 1.1106213457693743
 1.0759472191909856
 1.0920698019624837
 1.0866983153129903
parallelize(lyapunovspectrum, bd, 1000, particles)
1000-element Vector{Float64}:
 0.8968067178536954
 0.8910371722057758
 0.8415404680237077
 0.8526916351093605
 0.8788544706420249
 0.8338437536736311
 0.8320015812975936
 0.8698092829912604
 0.8621559697161493
 0.5978318497811886
 ⋮
 0.8654841146346096
 0.6990919498610223
 0.906593330845827
 0.9603219703301735
 0.8927064649799904
 0.853310929554065
 0.8263737759593881
 0.8101951096526985
 0.8277969242434021

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_rectangleFunction
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_sinaiFunction
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_bunimovichFunction
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_mushroomFunction
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_polygonFunction
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_verticesFunction
billiard_vertices(v, type = FiniteWall)

Construct a polygon billiard that connects the given vertices v (vector of 2-vectors). The vertices should construct a billiard in a counter-clockwise orientation (i.e. the normal vector always points to the left of v[i+1] - v[i].).

type decides what kind of walls to use. Ths function assumes that the first entry of v should be connected with the last.

source
DynamicalBilliards.billiard_irisFunction
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
DynamicalBilliards.billiard_hexagonal_sinaiFunction
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_showcaseFunction
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
Function
billiard_logo(;h=1.0, α=0.8, r=0.18, off=0.25, T = Float64) -> bd, ray

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

source

Particle types

DynamicalBilliards.ParticleType
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.MagneticParticleType
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