Brownian Motion

globals [
  ;; time to next collision
  tick-delta
  ;; distance from the center to the walls
  box-edge
  ;; first and second particle colliding.
  particle1 particle2
  ;; list of collisions times [[time particle1 particle2] ...]
  collisions
  big-particle-size big-particle-mass
  small-particle-size small-particle-mass
]

breed [particles particle]

particles-own [
  speed
  mass
  energy
]
;;;;;;;;;;;;;;;;;;;;;; SETUP PROCEDURES ;;;;;;;;;;;;;;;;;;;;;;;;;

to setup
  clear-all
  set-default-shape particles "circle"
  set box-edge max-pxcor
  ask patches with [(abs pxcor = box-edge or abs pycor = box-edge)]
    [ set pcolor yellow ]
  set small-particle-mass 2 * initial-total-energy / number-of-particles
  set big-particle-mass mass-ratio * small-particle-mass
  set small-particle-size 1
  set big-particle-size 5
  make-particles
  set particle1 nobody
  set particle2 nobody
  set collisions []
  reset-ticks
  ask particles [ check-for-wall-collision ]
  ask particles [ check-for-particle-collision ]
end 

to make-particles
  ;;
  ;; Create one big particle at the center
  create-particles 1 [
    set speed 0
    set size big-particle-size
    set mass big-particle-mass
    set energy (mass * speed ^ 2) / 2
    set color yellow
    setxy 0 0
  ]
  create-particles number-of-particles [
    set speed 1
    set size 1
    set mass small-particle-mass
    set energy (mass * speed ^ 2) / 2
    set color green
    while [overlapping?] [position-randomly]
  ]
end 

to-report overlapping?
  report any? other particles in-radius small-particle-size or
  distance patch 0 0 <= (small-particle-size + big-particle-size) / 2
end 

to position-randomly ;; particle procedure
  ;; place particle at random inside the box
  setxy one-of [1 -1] * random-float (box-edge - 0.5 - size / 2)
  one-of [1 -1] * random-float (box-edge - 0.5 - size / 2)
end 

;;;;;;;;;;;;;;;;;; GO PROCEDURES ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to go
  ;; advance time to next collision tick-delta
  ;; and set the colliding particles particle1 and particle2
  choose-next-collision
  ask particles [ jump (speed * tick-delta) ]
  perform-next-collision
  tick-advance tick-delta
  recalculate-particles-that-just-collided
end 

to choose-next-collision
  ;;
  if collisions = [] [ stop ]
  ;; pick the shortes time of projected collision
  let this-next first collisions
  foreach collisions [a -> if first a < first this-next [set this-next a] ]
  let dt item 0 this-next
  ;; if the next collision is more than 1 tick only advance 1 tick
  set tick-delta (dt - ticks)
  if tick-delta > 1 [
    set tick-delta 1
    set particle1 nobody
    set particle2 nobody
    stop
  ]
  set particle1 item 1 this-next
  set particle2 item 2 this-next
end 

to perform-next-collision
  ;;
  if particle1 = nobody [ stop ]
  if is-string? particle2 [
    if particle2 = "left wall" or particle2 = "right wall" [
      ask particle1 [ set heading (- heading) ]
      stop
    ]
    if particle2 = "top wall" or particle2 = "bottom wall" [
      ask particle1 [ set heading (180 - heading) ]
      stop
    ]
  ]
  ask particle1 [ collide-with particle2 ]
end 

to collide-with [ other-particle ]
  ;;
  let mass2 [mass] of other-particle
  let speed2 [speed] of other-particle
  let heading2 [heading] of other-particle
  let theta towards other-particle
  ;; convert velocity to components along theta and perpendicular to theta
  let v1-tangent (speed * cos (theta - heading))
  let v1-perpendicular (speed * sin (theta - heading))
  let v2-tangent (speed2 * cos (theta - heading2))
  let v2-perpendicular (speed2 * sin (theta - heading2))
  ;; compute the velocity of the center of mass along theta
  let vcm (((mass * v1-tangent) + (mass2 * v2-tangent)) / (mass + mass2))
  ;; new velocities after the collision along direction of theta
  set v1-tangent (2 * vcm - v1-tangent)
  set v2-tangent (2 * vcm - v2-tangent)
  ;; convert back to normal speed and heading
  set speed sqrt (v1-tangent ^ 2 + v1-perpendicular ^ 2)
  set energy kinetic-energy
  if v1-perpendicular != 0 or v1-tangent != 0 [
    set heading (theta - (atan v1-perpendicular v1-tangent))
  ]
  ask other-particle [
    set speed sqrt(v2-tangent ^ 2 + v2-perpendicular ^ 2)
    set energy kinetic-energy
    if v2-perpendicular != 0 or v2-tangent != 0 [
      set heading (theta - (atan v2-perpendicular v2-tangent))
    ]
  ]
end 

to recalculate-particles-that-just-collided
  ;; Remove from the list collisions of particles that we're done
  ifelse is-turtle? particle2 [
    set collisions filter [ a ->
      item 1 a != particle1 and
      item 2 a != particle1 and
      item 1 a != particle2 and
      item 2 a != particle2 ] collisions
    ask particle2 [ check-for-wall-collision ]
    ask particle2 [ check-for-particle-collision ]
  ][
    ;; particle2 is a wall. Filter out particle1 of the list
    set collisions filter [ a ->
      item 1 a != particle1 and
      item 2 a != particle1 ] collisions
  ]
  if particle1 != nobody [
    ask particle1 [ check-for-wall-collision ]
    ask particle1 [ check-for-particle-collision ]
  ]
  ;; make sure we fileter out particles1 and particle2
  set collisions filter [ a ->
    item 1 a != particle1 or
    item 2 a != particle2 ] collisions
  ;; done
  set particle1 nobody
  set particle2 nobody
end 

to check-for-wall-collision
  ;; determines when a particle will hit a wall
  ;; first along the horizontal direction
  let vx (speed * dx)
  if vx != 0 [
    ;; how long to hit the right wall
    let t-right (box-edge - 0.5 - xcor - (size / 2)) / vx
    if t-right > 0 [assign-colliding-wall t-right "right wall"]
    ;; how long to hit the left wall
    let t-left ((- box-edge) + 0.5 - xcor + (size / 2)) / vx
    if t-left > 0 [assign-colliding-wall t-left "left wall"]
  ]
  ;; Now along the vertical direction
  let vy (speed * dy)
  if vy != 0 [
    ;; how long to hit the top wall
    let t-top (box-edge - 0.5 - ycor - (size / 2)) / vy
    if t-top > 0 [assign-colliding-wall t-top "top wall"]
    ;; how long to hit the bottom wall
    let t-bottom ((- box-edge) + 0.5 - ycor + (size / 2)) / vy
    if t-bottom > 0 [assign-colliding-wall t-bottom "bottom wall"]
  ]
end 

to assign-colliding-wall [time-to-collision wall]
  ;; store the time to collision, self (the particle) and the wall
  ;; add clock to time-to-collision
  let colliding-pair (list (time-to-collision + ticks) self wall)
  set collisions fput colliding-pair collisions
end 

to check-for-particle-collision
  ;;
  let my-x xcor
  let my-y ycor
  let my-particle-size size
  let my-x-speed speed * dx
  let my-y-speed speed * dy
  ask other particles [
    let dpx (xcor - my-x) ;; relative distance along x direction
    let dpy (ycor - my-y)
    let x-speed (speed * dx)
    let y-speed (speed * dy)
    let dvx (x-speed - my-x-speed) ;; relative speed difference
    let dvy (y-speed - my-y-speed)
    let sum-radii (((my-particle-size) / 2) + (([size] of self) / 2))
    let p-squared (dpx ^ 2 + dpy ^ 2 - sum-radii ^ 2)
    let pv ( 2 * ((dpx * dvx) + (dpy * dvy)))
    let v-squared ( dvx ^ 2 + dvy ^ 2)
    let determinant pv ^ 2 - (4 * v-squared * p-squared)
    let time-to-collision -1
    if determinant > 0 [
      set time-to-collision ( - pv - sqrt determinant) / (2 * v-squared)
    ]
    if time-to-collision > 0 [
      set collisions fput (list (time-to-collision + ticks) self myself)
      collisions
    ]
  ]
end 

to-report kinetic-energy
  report (mass * speed ^ 2) / 2
end