Internals

This page is not part of the public API defined by DynamicalBilliards. Consider it something like a developer's guide.

Implementation

Before talking about the low level methods that enable everything to work nicely together, let's talk about how this package works.

Firstly one defines a Billiard and optionally some RaySplitter instances. Then one creates a particle inside the defined billiard. The algorithm for the propagation of a particle is the following:

  1. Calculate the collision of the particle with all obstacles in the billiard.
  2. Find the collision that happens first (in time), and the obstacle corresponding to that.
  3. DynamicalBilliards.relocate! the particle, and ensure that it is inside the billiard. This means that DynamicalBilliards.distance between particle and obstacle is either positive or close to machine precision.
  4. (Optionally) check if there is transmission for ray-splitting: T(φ) > rand()
  1. Continue this loop for a given amount of time.

Notice that the DynamicalBilliards.relocate! step is very important because it takes care that all particles remain inside the billiard.

The exposed bounce! function bundles steps 1-4 together.

Where is "inside"?

If for some reason (finite numeric precision) a particle goes outside a billiard, then it will escape to infinity. But what is inside?

"Inside" is defined on obstacle level by the function distance:

DynamicalBilliards.distanceFunction
distance(p::AbstractParticle, o::Obstacle)

Return the signed distance between particle p and obstacle o, based on p.pos. Positive distance corresponds to the particle being on the allowed region of the Obstacle. E.g. for a Disk, the distance is positive when the particle is outside of the disk, negative otherwise.

distance(p::AbstractParticle, bd::Billiard)

Return minimum distance(p, obst) for all obst in bd. If the distance(p, bd) is negative and bd is convex, this means that the particle is outside the billiard.

WARNING : distance(p, bd) may give negative values for non-convex billiards, or billiards that are composed of several connected sub-billiards.

All distance functions can also be given a position (vector) instead of a particle.

source

Notice that for very small negative values of distance, collision takes care of finite precision issues and does not return wrong collisions.

Numerical Precision

All core types of DynamicalBilliards are parametrically constructed, with parameter T <: AbstractFloat. This means that the fields of all particles and obstacles contain numbers strictly of type T. You will understand why this choice happened as you continue reading this paragraph.

The main concerns during evolution in a billiard table are:

  1. The particle must never leak out of the billiard table. This is simply translated to the distance function being positive after any collision and that collision takes care of extreme cases with very small (but negative) distance.
  2. The collision time is never infinite, besides the cases of Pinned Particles in a magnetic billiard.

These are solved with two ways:

  1. After the next collision is computed, relocate! brings the particle to that point and calculates the distance with the colliding obstacle. If it is negative, it translates the particle's position by this distance, along the normal vector.
  2. collision takes care of cases where the distance between particle and obstacle is less than accuracy(::T). (This is necessary only for magnetic propagation, as for straight propagation checking the velocity direction with respect to the normal is always enough).

Adjusting the global precision of DynamicalBilliards is easy and can be done by choosing the floating precision you would like. This is done by initializing your billiard table with parametric type T, e.g. bd = billiard_sinai(Float16(0.3)). This choice will propagate to the entire bd, all particles resulting from randominside, as well as the entire evolution process.

BigFloats

Evolution with BigFloat in DynamicalBilliards is on average 3 to 4 orders of magnitude slower than with Float64.


Collision Times

DynamicalBilliards.collisionFunction
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.

Return Inf, SV(0, 0) if the collision is not possible or if the collision happens backwards in time.

It is the duty of collision to avoid incorrect collisions when the particle is on top of the obstacle (or very close).

source
DynamicalBilliards.next_collisionFunction
next_collision(p::AbstractParticle, bd::Billiard) -> i, tmin, cp

Compute the collision across all obstacles in bd and find the minimum one. Return the index of colliding obstacle, the time and the collision point.

source

Non-Exported Internals

DynamicalBilliards.normalvecFunction
normalvec(obst::Obstacle, position)

Return the vector normal to the obstacle's boundary at the given position (which is assumed to be very close to the obstacle's boundary).

source
DynamicalBilliards.cellsizeFunction
cellsize(bd)

Return the delimiters xmin, ymin, xmax, ymax of the given obstacle/billiard.

Used in randominside(), error checking and plotting.

source

Propagation

DynamicalBilliards.propagate!Function
propagate!(p::AbstractParticle, t)

Propagate the particle p for given time t, changing appropriately the the p.pos and p.vel fields.

propagate!(p, position, t)

Do the same, but take advantage of the already calculated position that the particle should end up at.

source
DynamicalBilliards.resolvecollision!Function
resolvecollision!(p::AbstractParticle, o::Obstacle)

Resolve the collision between particle p and obstacle o, depending on the type of o (do specular! or periodicity!).

resolvecollision!(p, o, T::Function, θ::Function, new_ω::Function)

This is the ray-splitting implementation. The three functions given are drawn from the ray-splitting dictionary that is passed directly to evolve!(). For a calculated incidence angle φ, if T(φ) > rand(), ray-splitting occurs.

source
DynamicalBilliards.relocate!Function
relocate!(p::AbstractParticle, o::Obstacle, t, cp)

Propagate the particle to cp and propagate velocities for time t. Check if it is on the correct side of the obstacle. If not, change the particle position by distance along the normalvec of the obstacle.

source
DynamicalBilliards.specular!Function
specular!(p::AbstractParticle, o::Obstacle)

Perform specular reflection based on the normal vector of the Obstacle.

In the case where the given obstacle is a RandomObstacle, the specular reflection randomizes the velocity instead (within -π/2+ε to π/2-ε of the normal vector).

source

Timeseries

DynamicalBilliards.extrapolateFunction
extrapolate(particle, prevpos, prevvel, ct, dt[, ω]) → x, y, vx, vy, t

Create the timeseries that connect a particle's previous position and velocity prevpos, prevvel with the particle's current position and velocity, provided that the collision time between previous and current state is ct.

dt is the sampling time and if the particle is MagneticParticle then you should provide ω, the angular velocity that governed the free flight.

Here is how this function is used (for example)

prevpos, prevvel = p.pos + p.current_cell, p.vel
for _ in 1:5
    i, ct, pos, vel = bounce!(p, bd)
    x, y, x, vy, t = extrapolate(p, prevpos, prevvel, ct, 0.1)
    # append x, y, ... to other vectors or do whatever with them
    prevpos, prevvel = p.pos + p.current_cell, p.vel
end
source

Cyclotron center is a field of `MagneticParticle`

For almost all operations involving a MagneticParticle, the center of the cyclotron is required. In order to compute this center only when it physically changes, we have made it a field of the struct.

This means that after changing the position or velocity of the particle, this center must be changed by doing mp.center = DynamicalBilliards.find_cyclotron(mp). The bounce! function takes care of that in the most opportune moment, but if you want to write your own specific low level function, do not forget this point!