Demonic Particles

globals [ color-range ]
turtles-own [ speed mass last-collision neighbor-speeds ]

to setup

  clear-all

  draw-walls

  create-turtles population [

    setxy random-xcor random-ycor
    avoid-walls

    set shape "circle"
    set size 0.5
    set mass 20

    ;; make sure the speed is never 0
    set speed random-float 1
    if speed = 0 [ set speed 0.1 ]

    ;; set the particle's color according to its speed. Red is fast, Sky-blue is
    ;; slow and white-ish is average
    ifelse (speed >= 0.5)
    [ set color scale-color red speed 1.5 0.5 ]
    [ set color scale-color sky speed -0.5 0.5 ]

    ;; keeps a record of the speeds of the other turtles in
    ;; its room
    set neighbor-speeds ( list speed )
  ]

  reset-ticks
end 

to draw-walls

  ;; draw membrane
  ask patches with [ pxcor = 0 ] [
    set pcolor yellow
  ]

  ;; draw the surrounding walls
  ask patches with [
      pxcor = max-pxcor or
      pxcor = min-pxcor or
      pycor = max-pycor or
      pycor = min-pycor ] [ set pcolor brown ]
end 

;; reposition turtles that are in a wall patch

to avoid-walls

  if pxcor = 0 [
      set xcor (xcor + one-of [1 -1])
    ]
    if pxcor = max-pxcor [
      set xcor xcor - 2
    ]
    if pxcor = min-pxcor [
      set xcor xcor + 2
    ]
    if pycor = max-pycor [
      set ycor ycor - 2
    ]
    if pycor = min-pycor [
      set ycor ycor + 2
    ]
end 

to go

  ask turtles [
    check-wall-collision
    update-neighbor-speeds
    if collide? [ check-for-particlecollision ]

    let max-speed max [ speed ] of turtles

    ifelse (speed >= 0.5)
    [ set color scale-color red ( speed / max-speed ) 1.5 0.5 ]
    [ set color scale-color sky ( speed / max-speed ) -0.5 0.5 ]

    ; normalize the speed
    forward speed / max-speed
  ]
  tick
end 

;; bounce when collides with any wall

to check-wall-collision

  ; check collision with middle wall
  if [ pxcor ] of patch-ahead 1 = 0 [
    ifelse [ pcolor ] of patch-ahead 1 = yellow
      [ decide-traspase-door ]
      [ set heading (- heading) ]
  ]

  ; check collision with borders
  if abs ( [ pxcor ] of patch-ahead 1 ) = max-pxcor [
    set heading (- heading)
  ]
  if abs [ pycor ] of patch-ahead 1 = max-pycor [
    set heading (180 - heading)
  ]
end 

to decide-traspase-door

  ;; Calculate the difference between a particle speed and its log of other particles' speed
  let speed-dif abs ( mean neighbor-speeds - speed )

  ifelse speed-dif < traspasing-threshold
    [ set heading (- heading) ]
    [
      jump 2
      set neighbor-speeds ( list speed )
    ]
end 

;; update the mean speed that each particle perceives according to
;; their encounters with other turtles

to update-neighbor-speeds

  set neighbor-speeds sentence neighbor-speeds [ speed ] of turtles-on neighbors

  if length neighbor-speeds > memory [
    set neighbor-speeds sublist neighbor-speeds memory length neighbor-speeds
  ]
end 


;; From now on, collision is implemented. This code comes from the Connected Chemistry 1
;; (Bike Tire) model, which can be found in the models library. Copyright to Uri Wilensky,
;; licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 3.0 License.

to check-for-particlecollision ;; particle procedure
  if count other turtles-here  >= 1
  [
    ;; the following conditions are imposed on collision candidates:
    ;;   1. they must have a lower who number than my own, because collision
    ;;      code is asymmetrical: it must always happen from the point of view
    ;;      of just one particle.
    ;;   2. they must not be the same particle that we last collided with on
    ;;      this patch, so that we have a chance to leave the patch after we've
    ;;      collided with someone.
    let candidate one-of other turtles-here with
      [who < [who] of myself and myself != last-collision]
    ;; we also only collide if one of us has non-zero speed. It's useless
    ;; (and incorrect, actually) for two turtles with zero speed to collide.
    if (candidate != nobody) and (speed > 0 or [speed] of candidate > 0)
    [
      collide-with candidate
      set last-collision candidate
      ask candidate [ set last-collision myself ]
    ]
  ]
end 

;; implements a collision with another particle.
;;
;; THIS IS THE HEART OF THE PARTICLE SIMULATION, AND YOU ARE STRONGLY ADVISED
;; NOT TO CHANGE IT UNLESS YOU REALLY UNDERSTAND WHAT YOU'RE DOING!
;;
;; The two turtles colliding are self and other-particle, and while the
;; collision is performed from the point of view of self, both turtles are
;; modified to reflect its effects. This is somewhat complicated, so I'll
;; give a general outline here:
;;   1. Do initial setup, and determine the heading between particle centers
;;      (call it theta).
;;   2. Convert the representation of the velocity of each particle from
;;      speed/heading to a theta-based vector whose first component is the
;;      particle's speed along theta, and whose second component is the speed
;;      perpendicular to theta.
;;   3. Modify the velocity vectors to reflect the effects of the collision.
;;      This involves:
;;        a. computing the velocity of the center of mass of the whole system
;;           along direction theta
;;        b. updating the along-theta components of the two velocity vectors.
;;   4. Convert from the theta-based vector representation of velocity back to
;;      the usual speed/heading representation for each particle.
;;   5. Perform final cleanup and update derived quantities.

to collide-with [ other-particle ] ;; particle procedure
  let mass2 0
  let speed2 0
  let heading2 0
  let theta 0
  let v1t 0
  let v1l 0
  let v2t 0
  let v2l 0
  let vcm 0



  ;;; PHASE 1: initial setup

  ;; for convenience, grab some quantities from other-particle
  set mass2 [mass] of other-particle
  set speed2 [speed] of other-particle
  set heading2 [heading] of other-particle

  ;; since turtles are modeled as zero-size points, theta isn't meaningfully
  ;; defined. we can assign it randomly without affecting the model's outcome.
  set theta (random-float 360)


  ;;; PHASE 2: convert velocities to theta-based vector representation

  ;; now convert my velocity from speed/heading representation to components
  ;; along theta and perpendicular to theta
  set v1t (speed * cos (theta - heading))
  set v1l (speed * sin (theta - heading))

  ;; do the same for other-particle
  set v2t (speed2 * cos (theta - heading2))
  set v2l (speed2 * sin (theta - heading2))


  ;;; PHASE 3: manipulate vectors to implement collision

  ;; compute the velocity of the system's center of mass along theta
  set vcm (((mass * v1t) + (mass2 * v2t)) / (mass + mass2) )

  ;; now compute the new velocity for each particle along direction theta.
  ;; velocity perpendicular to theta is unaffected by a collision along theta,
  ;; so the next two lines actually implement the collision itself, in the
  ;; sense that the effects of the collision are exactly the following changes
  ;; in particle velocity.
  set v1t (2 * vcm - v1t)
  set v2t (2 * vcm - v2t)


  ;;; PHASE 4: convert back to normal speed/heading

  ;; now convert my velocity vector into my new speed and heading
  set speed sqrt ((v1t * v1t) + (v1l * v1l))
  ;; if the magnitude of the velocity vector is 0, atan is undefined. but
  ;; speed will be 0, so heading is irrelevant anyway. therefore, in that
  ;; case we'll just leave it unmodified.
  if v1l != 0 or v1t != 0
    [ set heading (theta - (atan v1l v1t)) ]

  ;; and do the same for other-particle
  ask other-particle [
    set speed sqrt ((v2t ^ 2) + (v2l ^ 2))
    if v2l != 0 or v2t != 0
      [ set heading (theta - (atan v2l v2t)) ]
  ]
end