This is a suggestion for making #47 possible.

At the moment, raysplitting automatically changes the pflag field of the collided object to its opposite value.

Maybe we can introduce one more function argument to the "ray-splitter", something like an affect! function, that does some in-place operation to the entire billiard table, e.g. inverting the pflag of all obstacles.

(it must be ensured that the collided obstacle changes its pflag though!)

the distance function measures the signed distance between particle and obstacle. The collisiontime function calculates the time until collision with a given obstacle. The function propagate! moves a particle forwards in time.

The basic way the package works is that it first calculates the collision time and then propagates a particle for that amount of time. After that other processes happen (e.g. specular reflection).

The problem is that after the particle is propagated, it should be really close (if not exactly on top) of an obstacle. However, measuring the distance between the particle and the obstacle, it can be as large as of order 1e-9. I think that this is really inaccurate and I cannot understand the reason why. Both the functions collisiontime and distance are supposed to result in exact answers in the "theoretical sense". How come that the numeric implementation is so far from it?

@dpsanders what do you think? Is it reasonable that after propagation the particle can still be so far away from an obstacle's borders?

Also decide what to do with #54 and #53

I suggest to consider using Plots.jl, rather than straight PyPlot for plotting.

This has a few benefits:

• Allows to use various backends
• Has @gif for automatically generating GIFs.

(The current plotting looks very nice, though!)

The link to the docs in the README is missing.

Issue Description

In Ray-splitting, magnetic particles cannot be turned into normal particles, e.g. without magnetic field. Using an extremely large magnetic field does not help, because it leads to extreme inaccuracies in the internals

Solution

We must somehow either overhaul or re-think how we define magnetic particles so that this is possible.

https://github.com/yuyichao/FunctionWrappers.jl

Julia does not replace this:

@inline timeprec(::Type{T}) where {T} = eps(T)^(4/5)

with a constant value, but actually evaluates stuff
We can use @eval and a loop over the possible {T} and simply "manually" inline the function.

Or we can even use @generated as well

And as a result, the current implementation that is called "poincaresos" should be renamed to boundarymap.

A DoorWall will simply be a FiniteWall, where when the particle collides with it, it "escapes" the billiard.

This allows for computation of escape times.

I am guessing that the DoorWall will be super easy to implement, but one would need to also create an
escapetime function, such that when the index of collision becomes the index of the DoorWall, the escape function stops and returns the time of escape!

raysplit magnetic is IMPOSSIBLY AND EXTREMELY slow sometimes. There are some
problematic orbits for sure (because sometimes its fast).
Or maybe there is some really slow regression happening during relocate.
One really needs to debug this. (That "one" is me).

A function arclength that reduces the billiard table to a single variable (the arclength) will be helpful.

The function can either get a position (assumed on the obstacle) and spit out s.

However, one can associate the function with a billiard table. Given the table and an angle φ, it can again spit out s.

issue already in progress by Lukas

I also have an idea, but not sure yet if possible:

When calling functions like psos one can use a precomputed arc length array with some kind of resolution ds ... as to not compute arclengths at every collision but instead extrapolate...? will be tricky, because I dont see how to extrapolate from given position.

I cannot replicate that issue in my computer. No matter how many times I run the tests they always pass.

If somebody can replicate this please investigate.
The function to run is raysplit_magnetic() at the /test/raysplit.jl file. If this function freezes for you I can tell you what to look for.

Plotting in DynamicalSystems uses the Requires module. E.G.:

using Requires
@require PyPlot include("visualizations.jl")

and then noone has to ever do the nonsense enableplotting. The plotting functions become available the moment you do using PyPlot.

If this is done, one should also be careful to update the documentation page.

BigFloat and raysplitting does not work for magnetic propagation!

pretty printing works by overloading show. e;g;

#pretty print:
show(io::IO, w::Wall{T}) where {T} = print(io, "$(w.name) {$T}\n",
"start point: $(w.sp)\nend point: $(w.ep)\nnormal vector: $(w.normal)")

the printing of billiard table should list all obstacle names line by line

I am pretty sure that the code is type stable at almost all function calls but it's better to be 100% sure than "pretty sure"

With the new billiard table, we should split particles and obstacles.

I also believe that may method declarations could be reduced. In general many functions are copied and pasted all over the place. This is clearly not a good approach and I am sure we can do better.

• split code into more files
• randominside 1 function
• relocate (without raysplitting) one function
• escape time 1 function
• reduce source code of evolve! functions

(Related with #53 )

Should we create a dedicated struct for a "raysplitter", and make the ray-splitting dictionary a tuple of ray splitters???

something like:

struct RaySplitter{A,B,C,D}
  transmission::A
  refraction::B
  angular::C
  affect!::D
end

with the default implementation of affect! discussed in #53, and also angular defaulting to nothing?

Changing the billiard table to a tuple completely removes any allocations in the case of straight propagation:

bt = billiard_buminovich()
p = randominside()
bt2 = (bt...)
@btime next_collision($p, $bt);

  1.237 μs (8 allocations: 128 bytes)
@btime next_collision($p, $bt2);

  36.266 ns (0 allocations: 0 bytes)

however, it does not in the case of magnetic fields. We all know why, because of creating the vector of intersection points before finding which one is the point of first collision.

Inverse magnetic billiards are systems where there is magnetic field outside a region, but not inside, making the particle exit the boundary and return again inside the billiard. An excellent example is this one:

where you have inverse Buminovich, taken from 10.1103/PhysRevE.67.065202

Making this possible would allow all kind of craziness in inverse periodic billiards, SUPER COOL STUFF.

Heavily relies on somehow solving #46 first.

A standard evolve call has at least 100,000 allocations even though everything is done through SVectors.

I do not know how or why this happens and also how to solve it.

Just do

using DynamicalBilliards, BenchmarkTools
bt = billiard_sinai()
p = randominside(bt)
evolve!(p, bt, 1.0) # compile
@btime evolve!(p, bt, 1000.0)

to see what I mean.

, which will represent the billiard table from now on.

The plan is to incorporate this with #30, and have a tuple instead of a Vector of Obstacles.

A second field will be a vector of indices, that sorts the obstacles to an increasing arclength order.

It is probably necessary to parameterize the table to the entire tuple O, but maybe it is possible to parameterize only on T, the number type (only if there are no performance hits).

I noticed that the most current release of DynamicalBilliards is version 1.5.1, while the documentation seems to be geared towards v1.6. In case this is intended: Maybe it would be nice to keep two versions of the documentation, like latest/stable? (In that way the reputation of a nice package doesn't get tarnished by failing examples in the documentation 😉 )

I am quite sure it will be much faster to use static vectors and not mutable static vectors all over lyapunovs.

Also, there is no need for mutable matrix, because you can simply vcat or hcat static vectors and it is super efficient, and you simply use static matrices.

Also, accessing static matrices with ranges like 1:2 is bad for performance. They should be accessed with static vectors of indices, like SVector(1, 2) instead.

The function cyclotron is computed at all collisiontime functions. This means that the same computation is done like 5 times. We should make both center and radius fields of MagneticParticle and instead have the functions resolvecollision to change the fields when they e.g. reflect the particle.

Is it setup? There's no badge.

tmin, i = next_collision(p, bt)

(0.388…, 2)
@btime relocate!($p, $bt[i], $tmin)

  59.723 ns (2 allocations: 32 bytes)

I think there is no reason for relocate! to allocate, besides bad programming.

The function will propagate a particle in a billiard table with a set of initial conditions (either number->random i.c. or vector of vectors).

The will compute the poincare sos for when the j variable crosses some value either from positive to negative or vise versa.

The crossing is calculated via interpolation, and same technique applies to magnetic and standard:

1. propagate (either for full time or step by step).
2. Start iterating through j variable timeseries. When sign*(a[j] - a[j-1]) > crossing point then you have crossing of psos.
3. interpolate between a[j] and a[j-1] linearly.
4. profit

We need a function that evolves the particle with a single step, for only a single collision, without doing any unnecessary stuff like collecting vectors or whatever. This function should be non-allocating! A good name inspired by DiffEq is step!. But I really feel like we could come up with something even better.

Then, the high level function evolve! or evolve (that first deepcopies the particle) will call this new step! function as necessary.

This will make many other aspects of code clearer and fasterer and betterer. Especially the process used for plotting!

I have made major changes in the way this package work. I don't remember how fast it was before, but it certainly is fast now. And very precise as well, with the user being able to control the precision!

In the section https://juliadynamics.github.io/DynamicalBilliards.jl/latest/physics/#numerical-precision I am describing how this process works.

As far as the lyapunov exponents are concerned, the functions that are introduced/changed are:

1. timeprec(T) at https://github.com/JuliaDynamics/DynamicalBilliards.jl/blob/master/src/Propagation.jl#L11
2. specular at https://github.com/JuliaDynamics/DynamicalBilliards.jl/blob/master/src/Propagation.jl#L37 (which is made to depend on parametric types)
3. resolvecollision at https://github.com/JuliaDynamics/DynamicalBilliards.jl/blob/master/src/Propagation.jl#L77 now simply decides whether to do specular reflection or not
4. Most importantly, the function relocate! at https://github.com/JuliaDynamics/DynamicalBilliards.jl/blob/master/src/Propagation.jl#L87 has been fundamentally changed.
5. Due to the above change, the way propagate! works and is used in evolve!() is different now. See https://github.com/JuliaDynamics/DynamicalBilliards.jl/blob/master/src/Propagation.jl#L162 and https://github.com/JuliaDynamics/DynamicalBilliards.jl/blob/master/src/Propagation.jl#L321

@dapias Since you have contributed the Lyapunov exponents, would you also like to update your functions so that they work with this new and major change to the package? If you do that, the lyapunovspectrum will be included in the 1.5.0 release version!

Also, notice that the tests have been updated. I have created an entirely new function that takes care of the lyapunov exponents: https://github.com/JuliaDynamics/DynamicalBilliards.jl/blob/master/test/various.jl#L61 and you don't have to do anything about that (when the functions in the LyapunovSpectrum.jl file work properly).

Hello. I have found that this piece of code does not work adequately

using Dynamical Billiards
a = 1.0
b = 3.0
bt = billiard_polygon(6, b; setting = "periodic")
disc = Disk([0., 0.], a)
push!(bt, disc)
p = randominside(bt)
t = 5
ts, pos, vel = evolve!(p, bt, t)

And the reason is that when we create the normal vector as

normal = inr*normalize([-w[2], w[1]])

defined in the function billiard_polygon (line 151 of StandardBilliards.jl), the periodicity! function updates the position of the particle wrongly as w.normal is normalized. What I have tried to fix the problem is to change that line by

normal = normalize([-w[2], w[1]])*sqrt(3)*r

As $sqrt(3) \times r$ is the distance between any two parallel sides of a hexagon. The new problem is that this just works for a triangular billiard (hexagonal unit cell). As the function billiard_polygon works only for $4$ or $6$ sides (by mathematical reasons) and the square case works correctly, I propose to implement a specific function for the triangular lattice.

Additionally, a specific function for plotting this kind of unit cell needs to be implemented as the existing function does not translate correctly the scatterers when the periodicity is along non-coordinate vectors. I have made some code for the case discussed here but we could consider to open another issue after solving this one.

Some kind of interface would allow specific indexes of the billiard table to be excluded from search during minimum collision time estimation.

See https://datseris.github.io/DynamicalBilliards.jl/latest/physics/#implementation

By reading the paper, we know that in actuality only a method that propagates the offset vectors in magnetic field needs to be added, all the reflection methods are the same!

Since the exponents are known to be λ1 = -λ4 and λ2=λ3=0 , Only return λ1?

Although having all 4 serves as a good test case. Still, we can do the computation for all 4 (which will happen anyway) and only return the largest one.

The function randominside takes about 2.3µs per particle in some billiards (see benchmark report attached below). Maybe there is a way to make it faster, which would be great for applications that require lots of particles.

Benchmark report for randominside

This code

using DynamicalBilliards, BenchmarkTools
bt = billiard_mushroom()
@benchmark randominside($bt)

produced this output

BenchmarkTools.Trial: 
  memory estimate:  1.05 KiB
  allocs estimate:  30
  --------------
  minimum time:     716.871 ns (0.00% GC)
  median time:      1.574 μs (0.00% GC)
  mean time:        2.284 μs (28.37% GC)
  maximum time:     136.644 μs (97.92% GC)
  --------------
  samples:          10000
  evals/sample:     124

This could possibly increase performance.

With combination of the last answers of this discourse post: https://discourse.julialang.org/t/is-there-a-way-to-forward-getindex-for-tuples-in-a-type-stable-way/2889/20
it could lead to major improvements.

One needs to do proper benchmarking and profiling during this change!

(this issue should be solved at the same time with #51 )

The function sincos (currently on Julia master) returns the sine and cosine of an angle more efficiently than calling the individual functions.

1. Rename FiniteWall to InfiniteWall
2. Create a FiniteWall.
3. The collision time will return Inf if the collision happens at a point which does not belong to the actual length of the wall.
4. This will be only a small amount of additional operations, and will allow the creation of non-convex billiards, like the mushroom one.

most notably all evolve like functions