Fix elastic collision

src/spaces/continuous.jl

@@ -321,10 +321,11 @@ function elastic_collision!(a, b, f = nothing)
     # https://en.wikipedia.org/wiki/Elastic_collision#Two-dimensional_collision_with_two_moving_objects
     v1, v2, x1, x2 = a.vel, b.vel, a.pos, b.pos
     length(v1) ≠ 2 && error("This function works only for two dimensions.")
-    r1 = x1 .- x2
+    r1 = x1 .- x2 # B to A
     n = norm(r1)^2
     n == 0 && return false # do nothing if they are at the same position
-    r2 = x2 .- x1
+    dv = a.vel .- b.vel
+    r2 = x2 .- x1 # A to B
     m1, m2 = f === nothing ? (1.0, 1.0) : (getfield(a, f), getfield(b, f))
     # mass weights
     m1 == m2 == Inf && return false
@@ -339,14 +340,13 @@ function elastic_collision!(a, b, f = nothing)
         v2 = ntuple(x -> zero(eltype(v1)), length(v1))
         f1, f2 = 2.0, 0.0
     else
-        # Check if disks face each other, to avoid double collisions
-        !(dot(r2, v1) > 0 && dot(r2, v1) > 0) && return false
+        # Check if disks face or overtake each other, to avoid double collisions
+        dot(dv, r2) ≤ 0 && return false
         f1 = (2m2 / (m1 + m2))
         f2 = (2m1 / (m1 + m2))
     end
-    dv = a.vel .- b.vel
     a.vel = v1 .- f1 .* (dot(dv, r1) / n) .* (r1)
-    b.vel = v2 .- f2 .* (dot(dv, r2) / n) .* (r2)
+    b.vel = v2 .+ f2 .* (dot(dv, r2) / n) .* (r2)
     return true
 end
 

test/continuous_space_tests.jl

@@ -297,7 +297,9 @@ using StableRNGs
         end
 
         function kinetic(model)
+            # Kinetic enrgy
             K = sum(sum(abs2.(a.vel)) for a in allagents(model))
+            # Momentum
             p = (0.0, 0.0)
             for a in allagents(model)
                  p = p .+ a.vel
@@ -322,12 +324,10 @@ using StableRNGs
         @test y == 0
 
         # x, which is the changed velocities
-        # should be at least the amount of collisions happened divided by 2
-        # because half the agents are unmovable (in a collision at most one agent
-        # must change velocity)
-        @test x > 0
+        # should be at most half the agents
+        # because half the agents are unmovable
+        @test 0 < x ≤ 50
         @test model.c > 0
-        @test x ≥ model.c/2
         K1, p1 = kinetic(model)
         @test p1 != p0
         # TODO: This test fails but I do not know why. Must be fixed later.