

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

5. If yes, perform the ray-splitting algorithm (see the ray-splitting page).

6. If not, then DynamicalBilliards.resolvecollision! of the particle with the obstacle (specular or periodic conditions).

7. 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.distance — Function.

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 this means that the particle is outside the billiard.

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.collision — Function.

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.

Returns 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_collision — Function.

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

Obstacle related

# DynamicalBilliards.normalvec — Function.

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.cellsize — Function.

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

# DynamicalBilliards.periodicity! — Function.

periodicity!(p::AbstractParticle, w::PeriodicWall)

Perform periodicity conditions of w on p.

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!