Finite splitter wall (#202)

* 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.

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

Author: omdowley

CHANGELOG.md

@@ -1,3 +1,5 @@
+# 3.12.0
+* New `FiniteSplitterWall` structure.
 # 3.11.0
 * `particlebeam` and `randominside_xyφ` functions.
 # 3.10.0

Project.toml

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

docs/src/ray-splitting.md

@@ -6,7 +6,7 @@ them using a simple example in this documentation page.
 ## 1. Ray-Splitting Obstacles
 The first step is that an [`Obstacle`](@ref) that supports ray-splitting is required to be present in your billiard table. The only new feature these obstacles have is an additional Boolean field called `pflag` (propagation flag). This field notes on which side of the obstacle the particle is currently propagating.
 
-The normal vector as well as the distance from boundary change sign depending on the value of `pflag`. The obstacles [`Antidot`](@ref) and [`SplitterWall`](@ref) are the equivalents of disk and wall for ray-splitting.
+The normal vector as well as the distance from boundary change sign depending on the value of `pflag`. The obstacles [`Antidot`](@ref), [`SplitterWall`](@ref) and [`FiniteSplitterWall`](@ref) are the equivalents of disk, wall and finite-wall for ray-splitting.
 
 Let's create a billiard with a bunch of ray-splitting obstacles!
 ```@example ray

docs/src/tutorials/billiard_table.md

@@ -101,7 +101,7 @@ evolution:
      means that during collision time estimation, if the collision point that was
      calculated lies *outside* of the boundaries of the `FiniteWall`, then the
      returned collision time is `Inf` (no collision). `FiniteWall` is slower
-     than `InfiniteWall` for that reason.
+     than `InfiniteWall` for that reason. [`FiniteSplitterWall`](@ref) behaves like a `FiniteWall` during evolution.
 
 
 If you wish to create a billiard table that you know will be convex, you should
@@ -133,6 +133,7 @@ RandomWall
 PeriodicWall
 SplitterWall
 FiniteWall
+FiniteSplitterWall
 ```
 
 ---

src/billiards/obstacles.jl

@@ -1,6 +1,6 @@
 export Obstacle, Disk, Antidot, RandomDisk, Wall, Circular,
 InfiniteWall, PeriodicWall, RandomWall, SplitterWall, FiniteWall,
-Semicircle, Ellipse
+FiniteSplitterWall, Semicircle, Ellipse
 export translate
 
 using InteractiveUtils
@@ -314,6 +314,52 @@ function SplitterWall(sp::AbstractVector, ep::AbstractVector,
 end
 SplitterWall(sp, ep, n, name::String = "Splitter wall") =
 SplitterWall(sp, ep, n, true, name)
+
+"""
+    FiniteSplitterWall{T<:AbstractFloat} <: Wall{T}
+Finite wall obstacle allowing fow ray-splitting (mutable type).
+### Fields:
+* `sp::SVector{2,T}` : Starting point of the Wall.
+* `ep::SVector{2,T}` : Ending point of the Wall.
+* `normal::SVector{2,T}` : Normal vector to the wall, pointing to where the
+  particle *will come from before a collision* (pointing towards the inside of the
+  billiard). The size of the vector is irrelevant
+  since it is internally normalized.
+* `isdoor::Bool` : Flag of whether this `FiniteSplitterWall` instance is a "Door".
+  Defaults to `false`.
+* `pflag::Bool` : Flag that keeps track of where the particle is currently
+  propagating (`pflag` = propagation flag).
+  `true` is associated with the `normal` vector the wall is
+  instantiated with. Defaults to `true`.
+* `name::String` : Name of the obstacle, given for user convenience.
+  Defaults to "Finite Splitter Wall".
+"""
+mutable struct FiniteSplitterWall{T<:AbstractFloat} <: Wall{T}
+    sp::SVector{2,T}
+    ep::SVector{2,T}
+    normal::SVector{2,T}
+    width::T
+    center::SVector{2,T}
+    isdoor::Bool
+    pflag::Bool
+    name::String
+end
+function FiniteSplitterWall(sp::AbstractVector, ep::AbstractVector,
+    n::AbstractVector, isdoor::Bool = false, pflag::Bool = true, name::String = "Finite Splitter Wall")
+    T = eltype(sp)
+    n = normalize(n)
+    d = dot(n, ep-sp)
+    if abs(d) > 10eps(T)
+        error("Normal vector is not actually normal to the wall: dot = $d")
+    end
+    T = eltype(sp) <: Integer ? Float64 : eltype(sp)
+    w = norm(ep - sp)
+    center = @. (ep+sp)/2
+    return FiniteSplitterWall{T}(SVector{2,T}(sp), SVector{2,T}(ep), SVector{2,T}(n),
+    w, SVector{2,T}(center), isdoor, pflag, name)
+end
+FiniteSplitterWall(a, b, c, n::String) = FiniteSplitterWall(a, b, c, false, true, n)
+
 #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)")
@@ -410,6 +456,7 @@ assumed to be very close to the obstacle's boundary).
 @inline normalvec(wall::Wall, pos) = wall.normal
 @inline normalvec(w::PeriodicWall, pos) = normalize(w.normal)
 @inline normalvec(w::SplitterWall, pos) = w.pflag ? w.normal : -w.normal
+@inline normalvec(w::FiniteSplitterWall, pos) = w.pflag ? w.normal : -w.normal
 @inline normalvec(disk::Circular, pos) = normalize(pos - disk.c)
 @inline normalvec(a::Antidot, pos) =
     a.pflag ? normalize(pos - a.c) : -normalize(pos - a.c)
@@ -494,11 +541,11 @@ function distance(pos::SV, e::Ellipse{T})::T where {T}
 end
 
 # The entire functionality of `distance_init` is necessary only for
-# FiniteWall !!!
+# FiniteWall and FiniteSplitterWall !!!
 distance_init(p::AbstractParticle, a::Obstacle) = distance_init(p.pos, a)
 distance_init(pos::SVector, a::Obstacle) = distance(pos, a)
 
-function distance_init(pos::SVector{2,T}, w::FiniteWall{T})::T where {T}
+function distance_init(pos::SVector{2,T}, w::Union{FiniteWall{T}, FiniteSplitterWall{T}})::T where {T}
 
     n = normalvec(w, pos)
     posdot = dot(w.sp .- pos, n)

src/plotting/obstacles.jl

@@ -1,11 +1,11 @@
 obcolor(::Obstacle) = (0,0.6,0)
 obcolor(::Union{RandomWall, RandomDisk}) = (149/255, 88/255, 178/255)
-obcolor(::Union{SplitterWall, Antidot, Ellipse}) = (0.8,0.0,0)
+obcolor(::Union{FiniteSplitterWall, SplitterWall, Antidot, Ellipse}) = (0.8,0.0,0)
 obcolor(::PeriodicWall) = (0.8,0.8,0)
 obalpha(::Obstacle) = 0.5
 obalpha(::Union{Antidot, Ellipse}) = 0.1
 obls(::Obstacle) = "solid"
-obls(::Union{SplitterWall, Antidot, Ellipse}) = "dashed"
+obls(::Union{FiniteSplitterWall, SplitterWall, Antidot, Ellipse}) = "dashed"
 obls(::PeriodicWall) = "dotted"
 
 """

src/timeevolution/collisions.jl

@@ -40,7 +40,7 @@ function collision(p::Particle{T}, w::Wall{T}) where {T}
     end
 end
 
-function collision(p::Particle{T}, w::FiniteWall{T}) where {T}
+function collision(p::Particle{T}, w::Union{FiniteWall{T}, FiniteSplitterWall{T}}) where {T}
     n = normalvec(w, p.pos)
     denom = dot(p.vel, n)
     # case of velocity pointing away of wall:

test/raysplt_tests.jl

@@ -125,3 +125,37 @@ function inside_antidot(args...)
 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
+    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 = []
+
+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))
+
+@test 0.3N < particles_on_left < 0.7N
+@test 0.3N < particles_on_right < 0.7N
+@test particles_in_middle == 0
+end