Splitter escapetime (#203)

* Add FiniteSplitterWall to obstacles.

* Add FiniteSplitterWall to exports.

* Fix FiniteSplitterWall typo.

* Make FiniteSplitterWall mutable.

* Fix FiniteSplitterWall color.

* Make FiniteSplitterWalls render dashed.

* Fix FiniteSplitterWall normal function.

* Impement distance_init for FiniteSplitterWalls.

* Add finite collisions to FiniteSplitterWalls.

* Add FiniteSplitterWall test.

* Update Project.toml and CHANGELOG.md.

* Add FiniteSplitterWall to docs.

* Add raysplitting to the escapetime! function.

* Fix typo.

* Add ray-splitting escapetime test.

* Increment patch version.

Co-authored-by: owend <[email protected]>

Author: omdowley

Project.toml

@@ -1,7 +1,7 @@
 name = "DynamicalBilliards"
 uuid = "4986ee89-4ee5-5cef-b6b8-e49ba721d7a5"
 repo = "https://github.com/JuliaDynamics/DynamicalBilliards.jl.git"
-version = "3.12.0"
+version = "3.12.1"
 
 [deps]
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"

src/timeevolution/highleveltimes.jl

@@ -23,7 +23,7 @@ escapetime(p, bd, t; warning = false) =
 escapetime(bd::Billiard, t; kwargs...) =
     escapetime(randominside(bd), bd, t; kwargs...)
 
-function escapetime!(p::AbstractParticle{T}, bd::Billiard{T}, t;
+function escapetime!(p::AbstractParticle{T}, bd::Billiard{T}, t, raysplitters = nothing;
         warning::Bool=false)::T where {T<:AbstractFloat}
 
     ei = escapeind(bd)
@@ -32,10 +32,13 @@ function escapetime!(p::AbstractParticle{T}, bd::Billiard{T}, t;
     end
 
     totalt = zero(T); count = zero(t); t_to_write = zero(T)
+    isray = !isa(raysplitters, Nothing)
+    isray && acceptable_raysplitter(raysplitters, bd)
+    raysidx = raysplit_indices(bd, raysplitters)
 
     while count < t
 
-        i, tmin, pos, vel = bounce!(p, bd)
+        i, tmin, pos, vel = bounce!(p, bd, raysidx, raysplitters)
         t_to_write += tmin
 
         if isperiodic(i, bd)

test/raysplt_tests.jl

@@ -126,36 +126,67 @@ end
 
 billiards_testset("RAY Stays inside", identity; caller = inside_antidot)
 
-@testset "FiniteSplitterWall" begin
-bdr = billiard_rectangle(1.0, 1.0)
-splitter = FiniteSplitterWall([0.5,0.0], [0.5,1.0], [-1.0,0.0])
-bd = Billiard(splitter, bdr...)
-
-refraction(ϕ, pflag, ω) = pflag ? 0.5ϕ : 2.0ϕ
-transmission(ϕ, pflag, ω) = begin
+test_refraction(ϕ, pflag, ω) = pflag ? 0.5ϕ : 2.0ϕ
+test_transmission(ϕ, pflag, ω) = begin
     if pflag
         0.5*exp(-(ϕ)^2/2(π/8)^2)
     else
         abs(ϕ) < π/4 ? 0.5*exp(-(ϕ)^2/2(π/4)^2) : 0.0
     end
 end
 
-N = 500
-rs = (RaySplitter([1], transmission, refraction),)
-ps = particlebeam(0.01, 0.5, 0, N, 0)
-positions = []
+@testset "FiniteSplitterWall" begin
+    bdr = billiard_rectangle(1.0, 1.0)
+    splitter = FiniteSplitterWall([0.5,0.0], [0.5,1.0], [-1.0,0.0])
+    bd = Billiard(splitter, bdr...)
+
+    N = 500
+    rs = (RaySplitter([1], test_transmission, test_refraction),)
+    ps = particlebeam(0.01, 0.5, 0, N, 0)
+    positions = []
+
+    for p in ps
+        ct, pos, vel = evolve!(p, bd, 2, rs)
+        push!(positions, pos[end][1])
+        reset_billiard!(bd)
+    end
+
+    particles_on_left = length(filter(x ->x < 0.4, positions))
+    particles_on_right = length(filter(x -> x > 0.6, positions))
+    particles_in_middle = length(filter(x -> 0.4 <= x <= 0.6, positions))
 
-for p in ps
-    ct, pos, vel = evolve!(p, bd, 2, rs)
-    push!(positions, pos[end][1])
-    reset_billiard!(bd)
+    @test 0.3N < particles_on_left < 0.7N
+    @test 0.3N < particles_on_right < 0.7N
+    @test particles_in_middle == 0
 end
 
-particles_on_left = length(filter(x ->x < 0.4, positions))
-particles_on_right = length(filter(x -> x > 0.6, positions))
-particles_in_middle = length(filter(x -> 0.4 <= x <= 0.6, positions))
+@testset "Ray-splitting escapetime" begin
+    verts = [[0.0,0.0], [0.0,1.0], [1.0,1.0], [1.0,0.0]]
+    bdr = []
+
+    # create rect billiard where vertical walls are doors
+    for i in eachindex(verts)
+        s = verts[i]
+        e = verts[mod1(i+1, length(verts))]
+        w = s - e
+        n = [-w[2], w[1]]
+        push!(bdr, FiniteWall(s, e, n, mod1(i, 2)==1))
+    end
+
+    splitter = FiniteSplitterWall([0.75,0.0], [0.75,1.0], [-1.0,0.0])
+    bd = Billiard(splitter, bdr...)
+
+    N = 500
+    rs = (RaySplitter([1], test_transmission, test_refraction),)
+    ps = particlebeam(0.01, 0.5, 0, N, 0)
+    ts = []
+
+    for p in ps
+        push!(ts, escapetime!(p, bd, 3, rs))
+        reset_billiard!(bd)
+    end
 
-@test 0.3N < particles_on_left < 0.7N
-@test 0.3N < particles_on_right < 0.7N
-@test particles_in_middle == 0
+    mean_t = +(ts...)/N
+    expected_mean_t = 0.5 * 1 + 0.5 * 0.75 * 2
+    @test abs(expected_mean_t - mean_t) < 0.2
 end