GasLab Adiabatic Piston

globals
[
  tick-delta                      ;; how much we advance the tick counter this time through
  max-tick-delta                  ;; the largest tick-delta is allowed to be
  raw-width raw-height            ;; box size variables
  gravity                         ;; acceleration of the piston
  piston-energy total-energy      ;; current energies
  piston-height                   ;; piston variables
  piston-vel                      ;; piston speed
  pressure
  pressure-history
  length-horizontal-surface       ;; the size of the wall surfaces that run horizontally - the top and bottom of the box
  length-vertical-surface         ;; the size of the wall surfaces that run vertically - the left and right of the box
  init-avg-speed init-avg-energy  ;; initial averages
  avg-speed avg-energy            ;; current averages
  tot-particle-energy             ;; sum of the energy of all particles taken together
  piston-kinetic-energy           ;; piston's kinetic energy
  piston-potential-energy         ;; piston's potential energy
  fast medium slow                ;; current counts
]

breed [ particles particle ]
breed [ flashes flash ]
breed [ pistons piston]


flashes-own [birthday]

particles-own
[
  speed mass energy          ;; particle info
  momentum-difference        ;; used to calculate pressure from wall hits
  last-collision             ;; used to prevent particles from colliding multiple times

]

pistons-own
[
  speed mass energy
]

to setup
  clear-all
  set-default-shape particles "circle"
  set-default-shape flashes "plane"
  set max-tick-delta 0.1073
  ;; box has constant size.
  set raw-width  round (0.01 * box-width  * max-pxcor)
  set raw-height round (0.01 * box-height * max-pycor)
  set piston-height raw-height
  make-box
  make-piston
  set piston-vel 0
  set gravity 0.125
  ;;; the length of the horizontal or vertical surface of
  ;;; the inside of the box must exclude the two patches
  ;;; that are the where the perpendicular walls join it,
  ;;; but must also add in the axes as an additional patch
  ;;; example:  a box with an box-edge of 10, is drawn with
  ;;; 19 patches of wall space on the inside of the box
  set length-horizontal-surface  ( 2 * (raw-height - 1) + 1)
  set length-vertical-surface  (piston-height)
  make-box
  make-particles
  set pressure-history [0 0 0]  ;; plotted pressure will be averaged over the past 3 entries
  update-variables
  set init-avg-speed avg-speed
  set init-avg-energy avg-energy
  reset-ticks
end 

to go
  if piston-height < 3
  [ user-message "The piston reached the bottom of the chamber. The simulation will stop."
    stop
  ]
  if piston-height >= 2 * raw-height - 1
  [ user-message "The piston reached the top of the chamber. The simulation will stop."
    stop
  ]
  ask particles [ bounce ]
  ask particles [ move ]
  ask particles [ check-for-collision ]
  move-piston
  tick-advance tick-delta
  if floor ticks > floor (ticks - tick-delta)
  [
    calculate-pressure
    update-variables
    update-plots
  ]
  calculate-tick-delta
  ;; we check for pcolor = black to make sure flashes that are left behind by the piston die
  ask flashes with [ticks - birthday > 0.4 or pcolor = black]
    [ die ]
  display
end 

to update-variables
  ;; particle variables
  set medium count particles with [color = green]
  set slow   count particles with [color = blue]
  set fast   count particles with [color = red]
  set avg-speed  mean [speed] of particles
  set avg-energy mean [energy] of particles
  set tot-particle-energy sum [energy] of particles

  ;; piston Variables
  set piston-kinetic-energy (0.5 * piston-mass * (piston-vel ^ 2))
  set piston-potential-energy (piston-mass * gravity * piston-height)
  set piston-energy (piston-kinetic-energy + piston-potential-energy)

  ;; system Variables
  set total-energy (tot-particle-energy + piston-energy)
  set length-vertical-surface (piston-height)
  calculate-pressure
end 

to calculate-tick-delta
  ;; tick-delta is calculated in such way that even the fastest
  ;; particle (or the piston) will jump at most 1 patch length in a
  ;; tick. As particles jump (speed * tick-delta) at every
  ;; tick, making tick length the inverse of the speed of the
  ;; fastest particle (1/max speed) assures that. Having each particle
  ;; advance at most one patch-length is necessary for them not to
  ;; "jump over" a wall or the piston.
  ifelse any? particles with [speed > 0]
    [ set tick-delta min list
                            (1 / (ceiling max (sentence ([speed] of particles) ([speed] of one-of pistons))))
                            max-tick-delta ]
    [ set tick-delta max-tick-delta ]
end 

;;; Pressure is defined as the force per unit area.  In this context,
;;; that means the total momentum per unit time transferred to the walls
;;; by particle hits, divided by the surface area of the walls.  (Here
;;; we're in a two dimensional world, so the "surface area" of the walls
;;; is just their length.)  Each wall contributes a different amount
;;; to the total pressure in the box, based on the number of collisions, the
;;; direction of each collision, and the length of the wall.  Conservation of momentum
;;; in hits ensures that the difference in momentum for the particles is equal to and
;;; opposite to that for the wall.  The force on each wall is the rate of change in
;;; momentum imparted to the wall, or the sum of change in momentum for each particle:
;;; F = SUM  [d(mv)/dt] = SUM [m(dv/dt)] = SUM [ ma ], in a direction perpendicular to
;;; the wall surface.  The pressure (P) on a given wall is the force (F) applied to that
;;; wall over its surface area.  The total pressure in the box is sum of each wall's
;;; pressure contribution.

to calculate-pressure
  ;; by summing the momentum change for each particle,
  ;; the wall's total momentum change is calculated
  set pressure 15 * sum [momentum-difference] of particles
  set pressure-history lput pressure but-first pressure-history
  ask particles
    [ set momentum-difference 0 ]  ;; once the contribution to momentum has been calculated
                                   ;; this value is reset to zero till the next wall hit
end 

to bounce  ;; particles procedure
  ;; get the coordinates of the patch we'll be on if we go forward 1
  let new-patch patch-ahead 1
  let new-px [pxcor] of new-patch
  let new-py [pycor] of new-patch
  ; if we're not about to hit a wall (yellow patch) or piston (orange patch)
  ; we don't need to do any further checks
  if ([pcolor] of new-patch != yellow and [pcolor] of new-patch != orange)
    [ stop ]

  ; check: hitting left or right wall?
  if (abs new-px = raw-width)
    ; if so, reflect heading around x axis
    [
      ;;  if the particle is hitting a vertical wall, only the horizontal component of the speed
      ;;  vector can change.  The change in velocity for this component is 2 * the speed of the particle,
      ;; due to the reversing of direction of travel from the collision with the wall
      set momentum-difference momentum-difference + (abs (dx * 2 * mass * speed) / length-vertical-surface)
      set heading (- heading) ]

  ; check: hitting top or bottom wall? (Should never hit top, but this would handle it.)
  if (abs new-py = raw-height)
  [
    ;;  if the particle is hitting a horizontal wall, only the vertical component of the speed
    ;;  vector can change.  The change in velocity for this component is 2 * the speed of the particle,
    ;; due to the reversing of direction of travel from the collision with the wall
    set momentum-difference momentum-difference + (abs (dy * 2 * mass * speed) / length-horizontal-surface)
    set heading (180 - heading)
  ]

  ; check: hitting piston?
  if (new-py = [pycor] of one-of pistons and (speed * dy) > piston-vel)
  [
    ;;  if the particle is hitting the piston, only the vertical component of the speed
    ;;  vector can change.  The change in velocity for this component is 2 * the speed of the particle,
    ;; due to the reversing of direction of travel from the collision with the wall
    ;; make sure that each particle finishes exchanging energy before any others can
    set momentum-difference momentum-difference + (abs (dy * 2 * mass * speed) / length-horizontal-surface)
    exchange-energy-with-piston
  ]

  ask patch new-px new-py
    [ sprout-flashes 1 [
      set color pcolor - 2
      set birthday ticks
      set heading 0
    ]
  ]
end 

to move  ;; particle procedure
  if patch-ahead (speed * tick-delta) != patch-here
    [ set last-collision nobody ]
  jump (speed * tick-delta)
end 

to check-for-collision  ;; particle procedure
  ;; Here we impose a rule that collisions only take place when there
  ;; are exactly two particles per patch.

  if count other particles-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 particles-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 particles 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 particles colliding are self and other-particle, and while the
;; collision is performed from the point of view of self, both particles 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
  ;;; 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

  ;; since particles are modeled as zero-size points, theta isn't meaningfully
  ;; defined. we can assign it randomly without affecting the model's outcome.
  let 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
  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 ^ 2) + (v1l ^ 2))
  set energy (0.5 * mass * (speed ^ 2))
  ;; 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))
    set energy (0.5 * mass * (speed ^ 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  ;; particle procedure
  ifelse speed < (0.5 * 10)
  [
    set color blue
  ]
  [
    ifelse speed > (1.5 * 10)
      [ set color red ]
      [ set color green ]
  ]
end 

;;;
;;; drawing procedures
;;;

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

;; creates initial particles

to make-particles
  create-particles number-of-particles
  [
    setup-particle
    random-position
    recolor
  ]
  calculate-tick-delta
end 

to setup-particle  ;; particle procedure
  set speed init-particle-speed
  set mass particle-mass
  set energy (0.5 * mass * (speed  ^ 2))
  set last-collision nobody
  set momentum-difference 0
end 

;; place particle at random location inside the box.

to random-position  ;; particle procedure
  setxy ((1 - raw-width)  + random-float (2 * raw-width - 2))
        ((1 - raw-height) + random-float (raw-height - 2))
end 


;; ------ Piston ----------

to make-piston
  ask patches with [pycor = 0 and (abs pxcor < raw-width)]
  [ sprout-pistons 1
    [ set color orange
      set heading 0
      set pcolor color
      ht
    ]
  ]
end 

to move-piston
  let old-piston-vel piston-vel
  set piston-vel (old-piston-vel - gravity * tick-delta)    ;;apply gravity
  let movement-amount ((old-piston-vel * tick-delta) - (gravity * (0.5 * (tick-delta  ^ 2))))
  ;; Setting the pcolor makes the piston look like a wall to the particles.
  ask pistons
  [ set pcolor black
    while [(piston-vel * tick-delta) >= 1.0]
      [ calculate-tick-delta ]
    ifelse piston-height + movement-amount <= 2 * raw-height - 1
    [
      fd movement-amount
      set piston-height (raw-height + [ycor] of one-of pistons)
    ]
    [
      set ycor raw-height - 1
      set piston-height 2 * raw-height - 1
      if piston-vel > 0
        [ set piston-vel 0 ]
    ]
    set speed piston-vel ;; just used for tick-delta calculations
    set pcolor color
    if (piston-vel < 0) ;; piston can't hit particles when moving upwards
    [ if (any? particles-here with [(speed * dy) > piston-vel])
      [ ;; only bounce particles that are moving down slower than the piston
        ;; faster ones should outrun it
        ask particles-here with [(speed * dy) > piston-vel]
        [
          ;;  if the particle is hitting the piston, only the vertical component of the speed
          ;;  vector can change.  The change in velocity for this component is 2 * the speed of the particle,
          ;; due to the reversing of direction of travel from the collision with the wall
          ;; make sure that each particle finishes exchanging energy before any others can
          set momentum-difference momentum-difference + (abs (dy * 2 * mass * speed) / length-horizontal-surface)
          exchange-energy-with-piston
        ]
      ]
    ]
  ]
end 

to exchange-energy-with-piston  ;; particle procedure -- piston and particle exchange energy
  let vx (speed * dx)         ;;only along x-axis
  let vy (speed * dy)         ;;only along y-axis
  let old-vy vy
  let old-piston-vel piston-vel
  set piston-vel ((((piston-mass - mass) / (piston-mass + mass)) * old-piston-vel) +
                  (((2 * mass) / (piston-mass + mass)) * old-vy))
  set vy ((((2 * piston-mass) / (piston-mass + mass)) * old-piston-vel) -
         (((piston-mass - mass) / (piston-mass + mass)) * old-vy))
  set speed (sqrt ((vx ^ 2) + (vy ^ 2)))
  set energy (0.5 * mass * (speed  ^ 2))
  set heading atan vx vy
end 


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