Connected Chemistry 3 Circular Particles

globals
[
  tick-delta                ;; how much we advance the tick counter this time through
  max-tick-delta            ;; the largest tick-delta is allowed to be
  box-edge                  ;; distance of box edge from axes
  colliding-particles
  sorted-colliding-particles
  colliding-particle-1
  colliding-particle-2
  colliding-pair
  original-tick-delta
  last-view-update
  manage-view-updates?
  view-update-rate          ;; specifies the minimum amount of simulation time that must
                            ;; pass before the view is updated
  winners
  smallest-particle-size
  largest-particle-size
  kinetic-energy-of-one
  kinetic-energy-of-two
  pressure
  total-particle-number
  total-energy
]

breed [ particles particle ]

particles-own
[
  speed
  mass
]

;;;;;;;;;;;;;;;;

to setup
  clear-all reset-ticks
  set total-particle-number initial-number
  set-default-shape particles "circle"
  set manage-view-updates? true
  set view-update-rate 0.2
  set box-edge (max-pxcor - 1)
  set smallest-particle-size 4
  set largest-particle-size 4
  set kinetic-energy-of-one 0
  set kinetic-energy-of-two 0
  make-box
  make-particles

  ;; set variable tick delta based on fastest particle.   If the fastest particle has a speed of 1,
  ;; then tick-delta is 1.  If the fastest particles has a speed of 10, then tick-delta is 1/10.
  set tick-delta (1 / (ceiling max [speed] of particles))
  set original-tick-delta tick-delta
  set colliding-particles []
  ask particles [check-for-wall-collision]
  ask particles [check-for-particle-collision ]
  set colliding-particle-1 nobody
  set colliding-particle-2 nobody
        ifelse labels?
    [ask turtles [ set label who set label-color orange + 3 ]]
    [ask turtles [ set label ""]]
end 

to go
  ;;Since only collisions involving the particles that collided most recently could be affected,
  ;;we filter those out of colliding-particles.  Then we recalculate all possible collisions for the
  ;;particles that collided last.  The ifelse statement is necessary because colliding-particle-2
  ;;can be either a particle or a string representing a wall.  If it is a wall, we don't want to
  ;;invalidate all collisions involving that wall (because the wall's position wasn't affected, those
  ;;collisions are still valid.
  ifelse is-turtle? colliding-particle-2
  [
    set colliding-particles filter [item 1 ? != colliding-particle-1 and
                                    item 2 ? != colliding-particle-1 and
                                    item 1 ? != colliding-particle-2 and
                                    item 2 ? != colliding-particle-2]
                              colliding-particles
    ask colliding-particle-2 [check-for-wall-collision]
    ask colliding-particle-2 [check-for-particle-collision]
  ]
  [
    set colliding-particles filter [item 1 ? != colliding-particle-1 and
                                    item 2 ? != colliding-particle-1]
                              colliding-particles
  ]
  if colliding-particle-1 != nobody [ask colliding-particle-1 [check-for-wall-collision]]
  if colliding-particle-1 != nobody [ask colliding-particle-1 [check-for-particle-collision]]

  sort-collisions
  calculate-energies

  ask particles [ jump speed * tick-delta ]
  if winners != []
      [collide-winners]
  tick-advance tick-delta
  do-plots

    ask particles [
    if show-speed-as-color? = "red-green-blue" [ recolor-banded ]
    if show-speed-as-color? = "violet shades" [ recolor-shaded ]
    if show-speed-as-color?  = "one color" [ recolor-none ]
    if show-speed-as-color?  = "custom color" [ ]
    ]
     ifelse labels?
    [ask turtles [ set label who set label-color orange + 3 ]]
    [ask turtles [ set label ""]]

  ;; flag that updates display only after enough simulation time has passed.
  ;; the display-update-rate sets the minimum simulation time that must pass
  ;; before updating the display.  This avoids many redisplays of the view for
  ;; a series of small time steps in the simulation (which would make the view show
  ;; what looks like particles slowing down as they get near multiple collision points)
  if manage-view-updates? [
    if (ticks - last-view-update) > view-update-rate
    [ display
      set last-view-update ticks ]
      ]
end 

to calculate-energies
  ask particle 0 [
     set kinetic-energy-of-one ( 0.5 * mass * (speed * speed))

  ]
  ask particle 1 [
     set kinetic-energy-of-two ( 0.5 * mass *  (speed * speed))
  ]
end 

to-report convert-heading-x [heading-angle]
  report sin heading-angle
end 

to-report convert-heading-y [heading-angle]
  report cos heading-angle
end 

to-report future-time? [time]
  ;;This reporter is necessary because of slight discrepancies in the floating point math.
  ;;Sometimes particles that are colliding at that instant are reported as colliding again
  ;;imperceptibly in the future; this causes the model to hang.
  ;;This function ensures that all expected collisions are legitimate collisions, and not just
  ;;ghosts of floating point errors.
  report time > .0000000001
end 

to check-for-particle-collision

;; check-for-particle-collision is a particle procedure that determines the time it takes to the collision between
;; two particles (if one exists).  It solves for the time by representing the equations of motion for
;; distance, velocity, and time in a quadratic equation of the vector components of the relative velocities
;; and changes in position between the two particles and solves for the time until the next collision

  let my-x xcor
  let my-y ycor
  let my-particle-size size
  let my-x-speed (speed * convert-heading-x heading)
  let my-y-speed (speed * convert-heading-y heading)

  ask particles with [self != myself]
  [
         let dpx (xcor - my-x)   ;; relative distance between particles in the x direction
         let dpy (ycor - my-y)    ;; relative distance between particles in the y direction
         let x-speed (speed * convert-heading-x heading) ;; speed of other particle in the x direction
         let y-speed (speed * convert-heading-y heading) ;; speed of other particle in the x direction
         let dvx (x-speed - my-x-speed) ;; relative speed difference between particles in the x direction
         let dvy (y-speed - my-y-speed) ;; relative speed difference between particles in the y direction
         let sum-r (((my-particle-size) / 2 ) + (([size] of self) / 2 )) ;; sum of both particle radii



        ;; To figure out what the difference in position (P1) between two particles at a future time (t) would be,
        ;; one would need to know the current difference in position (P0) between the two particles
        ;; and the current difference in the velocity (V0) between of the two particles.

        ;; The equation that represents the relationship would be:   P1 = P0 + t * V0

        ;; we want find when in time (t), P1 would be equal to the sum of both the particle's radii (sum-r).
        ;; When P1 is equal to is equal to sum-r, the particles will just be touching each other at
        ;; their edges  (a single point of contact).

        ;; Therefore we are looking for when:   sum-r =  P0 + t * V0

        ;; This equation is not a simple linear equation, since P0 and V0 should both have x and y components
        ;;  in their two dimensional vector representation (calculated as dpx, dpy, and dvx, dvy).


        ;; By squaring both sides of the equation, we get:     (sum-r) * (sum-r) =  (P0 + t * V0) * (P0 + t * V0)

        ;;  When expanded gives:   (sum-r ^ 2) = (P0 ^ 2) + (t * PO * V0) + (t * PO * V0) + (t ^ 2 * VO ^ 2)

        ;;  Which can be simplified to:    0 = (P0 ^ 2) - (sum-r ^ 2) + (2 * PO * V0) * t + (VO ^ 2) * t ^ 2

        ;;  Below, we will let p-squared represent:   (P0 ^ 2) - (sum-r ^ 2)
        ;;  and pv represent: (2 * PO * V0)
        ;;  and v-squared represent: (VO ^ 2)

        ;;  then the equation will simplify to:     0 = p-squared + pv * t + v-squared * t^2


         let p-squared   ((dpx * dpx) + (dpy * dpy)) - (sum-r ^ 2)   ;; p-squared represents difference of the
                                                                     ;; square of the radii and the square
                                                                     ;; of the initial positions

         let pv  (2 * ((dpx * dvx) + (dpy * dvy)))  ;;the vector product of the position times the velocity
         let v-squared  ((dvx * dvx) + (dvy * dvy)) ;; the square of the difference in speeds
                                                    ;; represented as the sum of the squares of the x-component
                                                    ;; and y-component of relative speeds between the two particles


         ;; p-squared, pv, and v-squared are coefficients in the quadratic equation shown above that
         ;; represents how distance between the particles and relative velocity are related to the time,
         ;; t, at which they will next collide (or when their edges will just be touching)

         ;; Any quadratic equation that is the function of time (t), can represented in a general form as:
         ;;   a*t*t + b*t + c = 0,
         ;; where a, b, and c are the coefficients of the three different terms, and has solutions for t
         ;; that can be found by using the quadratic formula.  The quadratic formula states that if a is not 0,
         ;; then there are two solutions for t, either real or complex.

         ;; t is equal to (b +/- sqrt (b^2 - 4*a*c)) / 2*a

         ;; the portion of this equation that is under a square root is referred to here
         ;; as the determinant, D1.   D1 is equal to (b^2 - 4*a*c)
         ;; and:   a = v-squared, b = pv, and c = p-squared.


         let D1 pv ^ 2 -  (4 * v-squared * p-squared)


         ;; the next line next line tells us that a collision will happen in the future if
         ;; the determinant, D1 is >= 0,  since a positive determinant tells us that there is a
         ;; real solution for the quadratic equation.  Quadratic equations can have solutions
         ;; that are not real (they are square roots of negative numbers).  These are referred
         ;; to as imaginary numbers and for many real world systems that the equations represent
         ;; are not real world states the system can actually end up in.

         ;; Once we determine that a real solution exists, we want to take only one of the two
         ;; possible solutions to the quadratic equation, namely the smaller of the two the solutions:

         ;;  (b - sqrt (b^2 - 4*a*c)) / 2*a
         ;;  which is a solution that represents when the particles first touching on their edges.

         ;;  instead of (b + sqrt (b^2 - 4*a*c)) / 2*a
         ;;  which is a solution that represents a time after the particles have penetrated
         ;;  and are coming back out of each other and when they are just touching on their edges.


         let time-to-collision  -1

         if D1 >= 0 and v-squared > 0
            [set time-to-collision (- pv - sqrt D1) / (2 * v-squared) ]        ;;solution for time step


         ;; if time-to-collision is still -1 there is no collision in the future - no valid solution
         ;; note:  negative values for time-to-collision represent where particles would collide
         ;; if allowed to move backward in time.
         ;; if time-to-collision is greater than 1, then we continue to advance the motion
         ;; of the particles along their current trajectories.  They do not collide yet.

         if future-time? time-to-collision
             [
             ;; time-to-collision is relative (ie, a collision will occur one second from now)
             ;; We need to store the absolute time (ie, a collision will occur at time 48.5 seconds.
             ;; So, we add ticks to time-to-collision when we store it.
              set colliding-pair (list (time-to-collision + ticks) self myself) ;; sets a three element list of
                                                        ;; time to collision and the colliding pair
              set colliding-particles lput colliding-pair colliding-particles  ;; adds above list to collection
                                                                               ;; of colliding pairs and time
                                                                               ;; steps
             ]
  ]
end 

to check-for-wall-collision ;; particle procedure for determining if a particle will hit one of the
                            ;; four walls of the box

  let x-speed (speed * convert-heading-x heading)
  let y-speed (speed * convert-heading-y heading)
  let xpos-plane (box-edge - 0.5)      ;;inside boundary of right side of the box
  let xneg-plane (- box-edge + 0.5)    ;;inside boundary of left side of the box
  let ypos-plane (box-edge - 0.5)      ;;inside boundary of top side of the box
  let yneg-plane (- box-edge + 0.5)    ;;inside boundary of bottom side of the box

  ;; find point of contact on edge of circle
  ;; points of contact located at 1 radius above, below, to the left, and to the right
  ;; of the center of the particle

  let contact-point-xpos (xcor + (size / 2))
  let contact-point-xneg (xcor - (size / 2))
  let contact-point-ypos (ycor + (size / 2))
  let contact-point-yneg (ycor - (size / 2))

  ;; find difference in position between plane location and edge of circle

  let dpxpos (xpos-plane - contact-point-xpos)
  let dpxneg (xneg-plane - contact-point-xneg)
  let dpypos (ypos-plane - contact-point-ypos)
  let dpyneg (yneg-plane - contact-point-yneg)

  let t-plane-xpos 0

  ;; solve for the time it will take the particle to reach the wall by taking
  ;; the distance to the wall and dividing it by the speed in the direction to the wall

  ifelse  x-speed != 0 [set t-plane-xpos (dpxpos / x-speed)] [set t-plane-xpos 0]
   if future-time? t-plane-xpos
      [
       assign-colliding-wall t-plane-xpos "plane-xpos"
      ]

  let t-plane-xneg 0
  ifelse  x-speed != 0 [set t-plane-xneg (dpxneg / x-speed)] [set t-plane-xneg 0]
   if future-time? t-plane-xneg
      [
       assign-colliding-wall t-plane-xneg "plane-xneg"
      ]
  let t-plane-ypos 0
  ifelse  y-speed != 0 [set t-plane-ypos (dpypos / y-speed)] [set t-plane-ypos 0]
   if future-time? t-plane-ypos
      [
       assign-colliding-wall t-plane-ypos "plane-ypos"
      ]

  let t-plane-yneg 0
  ifelse  y-speed != 0 [set t-plane-yneg (dpyneg / y-speed)] [set t-plane-yneg 0]
  if future-time? t-plane-yneg
      [
       assign-colliding-wall t-plane-yneg "plane-yneg"
      ]
end 

to assign-colliding-wall [time-to-collision wall]
  ;; this procedure is used by the check-for-wall-collision procedure
  ;; to assemble the correct particle-wall pair
  ;; time-to-collision is relative (ie, a collision will occur one second from now)
  ;; We need to store the absolute time (ie, a collision will occur at time 48.5 seconds.
  ;; So, we add ticks to time-to-collision when we store it.

  set colliding-pair (list (time-to-collision + ticks) self wall)
  set colliding-particles lput colliding-pair colliding-particles
end 

to sort-collisions
  ;; Sort the list of projected collisions between all the particles into an ordered list.
  ;; Take the smallest time-step from the list (which represents the next collision that will
  ;; happen in time).  Use this time step as the tick-delta for all the particles to move through

  ifelse colliding-particles != []
   [
    set winners first colliding-particles
    foreach colliding-particles [if first ? < first winners [set winners ?]]
    ;;winners is now the collision that will occur next
    let dt first winners
    set colliding-particle-1 item 1 winners
    set colliding-particle-2 item 2 winners
    if dt > 0
    [
      ;;If the next collision is more than 1 in the future,
      ;;clear the winners and advance the simulation one tick.
      ;;This helps smooth the model on smaller particle counts.
      ifelse dt - ticks <= 1
      ;;We have to subtract ticks back out because now we want the relative time until collision,
      ;;not the absolute time the collision will occur.
      [set tick-delta dt - ticks]
      ;;Since there are no collisions in the next second, we will set winners to [] to keep from
      ;;mistakenly colliding any particles that shouldn't collide yet.
      [set tick-delta 1 set winners []]
    ]

   ]
   ;; When there are no collisions for the next time step,
   ;; tick-delta goes back to the value of original-tick-delta
   [set tick-delta original-tick-delta]
end 

to collide-winners  ;; deal with 3 possible cases of collisions:
                    ;; particle and one wall, particle and two walls, and two particles

    ;; deal with a case where the next collision in time is between a particle and a wall

    if colliding-particle-2 = "plane-xpos" or colliding-particle-2 = "plane-xneg"
         [ask colliding-particle-1 [set heading (- heading)]
          stop]
    if colliding-particle-2 = "plane-ypos" or colliding-particle-2 = "plane-yneg"
         [ask colliding-particle-1 [set heading (180 - heading)]
          stop]

    ;; deal with the remaining case of the next collision in time being between two particles.

    ask colliding-particle-1 [collide-with colliding-particle-2]
end 

to collide-with [ other-particle ] ;; particle procedure

  ;;; PHASE 1: initial setup

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

  ;;modified so that theta is heading toward other particle
  let theta towards other-particle

  ;;; 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
  let v1t (speed * cos (theta - heading))
  let v1l (speed * sin (theta - heading))

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

  ;;; PHASE 3: manipulate vectors to implement collision

  ;; compute the velocity of the system's center of mass along theta
  let 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)) ]
  ]

  ;; PHASE 5: final updates

  ;; now recolor, since color is based on quantities that may have changed
  recolor ask other-particle [ recolor ]
end 

to recolor
    if show-speed-as-color? = "red-green-blue" [ recolor-banded ]
    if show-speed-as-color? = "violet shades" [ recolor-shaded ]
    if show-speed-as-color?  = "one color" [ recolor-none ]
end 

to recolor-banded  ;; particle procedure
  let avg-speed 1
  ;; avg-speed is assumed to be 0.5, since particles are assigned a random speed between 0 and 1
  ;; particle coloring procedures for visualizing speed with a color palette,
  ;; red are fast particles, blue slow, and green in between.

  ifelse speed < (0.5 * avg-speed) ;; at lower than 50% the average speed
  [
    set color blue       ;; slow particles colored blue
  ]
  [
    ifelse speed > (1.5 * avg-speed) ;; above 50% higher the average speed
      [ set color red ]        ;; fast particles colored blue
      [ set color green ]      ;; medium speed particles colored green
  ]
end 

to recolor-shaded
  let avg-speed 1
 ;; avg-speed is assumed to be 0.5, since particles are assigned a random speed between 0 and 1
 ;; a particle shading gradient is applied to all particles less than speed 1.5,
 ;; the uppermost threshold speed to apply the shading gradient to.

  ifelse speed < (3 * avg-speed)
  [ set color (violet - 3.001) + (8 * speed / (3 * avg-speed)) ]
  [ set color (violet + 4.999)]
end 

to recolor-none
  set color green - 1
end 


;;;
;;; drawing procedures
;;;

to make-box
  ask patches with [ ((abs pxcor = box-edge) and (abs pycor <= box-edge)) or
                     ((abs pycor = box-edge) and (abs pxcor <= box-edge)) ]
    [ set pcolor yellow ]
end 

;; creates some particles

to make-particles
  create-particles initial-number
  [

    set speed 1
    set size smallest-particle-size + random-float (largest-particle-size - smallest-particle-size)
    set mass (size * size) ;; set the mass proportional to the area of the particle
    random-position
    recolor
  ]
end 

;; place particle at random location inside the box.

to random-position ;; particle procedure
  let count-loops 0
  let radius (size / 2)
  let largest-radius (largest-particle-size / 2)
  let both-radii (radius + largest-radius)

  ;;this loop only works with < 25 particles
  while [(count particles in-radius both-radii) > 1 and count-loops < 5000]
  [
    setxy ((radius - (box-edge - size )) + random-float ((2 * (box-edge - size )) - (2 * radius)))
          ((radius - (box-edge - size )) + random-float ((2 * (box-edge - size )) - (2 * radius)))
          set count-loops count-loops + 1
  ]
end 

to do-plots
  set-current-plot "Particle 0 Heading vs. Time"
  plotxy ticks ([heading] of particle 0)

  set-current-plot "Particle 0 Speed vs. Time"
  plotxy ticks ([speed] of particle 0)

  set-current-plot "Kinetic Energy vs. Time"
  set-current-plot-pen "particle 0"
  plotxy ticks kinetic-energy-of-one
  set-current-plot-pen "particle 1"
  plotxy ticks kinetic-energy-of-two
end 


; Copyright 2005 Uri Wilensky.
; See Info tab for full copyright and license.