GasLab Moving Piston

globals [
  fast average slow     ;; current counts
  avg-speed avg-energy  ;; current averages
  vsplit vclock         ;; clock variables
  raw-width raw-height  ;; box size variables
  piston-position       ;; position of the piston at any given time
  volume area           ;; because this is 2D, area is the perimeter and volume is the area
  total-pressure        ;; pressure per unit area
  pressure-history      ;; list of 'scale' previous total-pressures
  avg-pressure          ;; mean of the pressure-history
  initspeed initmass    ;; initial speed and initial mass of the all the particles-particle mass always stays the same, but the speed changes.
  going-down?           ;; flag for whether or not the piston is moving up or down
]

turtles-own [
  speed mass energy new-speed            ;; Turtle Info
  v1t v1l tmp-turtle                     ;; Turtle 1 -- Collide
  heading2 mass2 speed2 v2t v2l turtle2  ;; Turtle 2 -- Collide
  theta                                  ;; Turtles 1 + 2 -- Collide
  pressure                               ;; pressure variable
]

;; procedure that setup up initial variables

to setup
  clear-all
  set going-down? true
  set pressure-history []
  set initspeed 10.0
  set initmass 1.0
  set raw-width  round (0.01 * box-width  * max-pxcor)
  set raw-height round (0.01 * box-height * max-pycor)
  set piston-position 0.75 * raw-height
  set area ((4 * raw-width) + (2 * (piston-position + raw-height)))
  set volume ((2 * raw-width) * (piston-position + raw-height))
  make-box
  draw-piston
  ;;set vclock 0
  ;; create the gas molecules
  crt number [
     set new-speed initspeed
     set mass initmass
     random-position
     set shape "circle"
     recolor
     set pressure 0
  ]
  update-variables
  reset-ticks
  setup-histograms
end 

;; updates variables after every time tick

to update-variables
  ;; Gas Variables
  ask turtles
    [ set speed new-speed
      set energy (0.5 * mass * speed * speed) ]
  set average count turtles with [ color = green ]
  set slow    count turtles with [ color = blue  ]
  set fast    count turtles with [ color = red   ]
  set avg-speed  mean [ speed  ] of turtles
  set avg-energy mean [ energy ] of turtles

  ;; System Variables
  calculate-pressure
  set vsplit (round (max [speed] of turtles * 1.2))
end 

;; procedure that runs the model

to go
  ask turtles [bounce]
  ask turtles [move]
  ask turtles [check-for-collision]
  ;; control the piston's motion
  if piston-position < (-0.75 * raw-height)
  [ set going-down? false ]
  if piston-position > (0.75 * raw-height)
  [ set going-down? true ]
  ifelse going-down?
  [ piston-down piston-speed / vsplit ]
  [ piston-up piston-speed / vsplit ]

  set vclock vclock + 1
  ifelse (vclock = vsplit)
  [
    tick
    set vclock 0
    update-variables
    do-plotting
    do-histograms
  ]
  [ display ]
end 


;; turtle procedure for bouncing off of the walls

to bounce
  ; if we're not about to hit a wall (yellow patch)
  ; or the piston (gray+2 patch),
  ; we don't need to do any further checks
  if ([pcolor] of patch-ahead 1 != yellow) and
     ([pcolor] of patch-ahead 1 != gray + 2) [ stop ]
  ; get the coordinates of the patch we'll be on if we go forward 1
  let new-px [pxcor] of patch-ahead 1
  let new-py [pycor] of patch-ahead 1
  ; check: hitting left or right wall?
  if (abs new-px = raw-width)
    ; if so, reflect heading around x axis
    [ set heading (- heading)
      set pressure pressure + abs (dx * mass * speed)
    ]
  ; check: hitting piston or bottom wall?
  if (abs new-py = raw-height) or (new-py = round piston-position)
    ; if so, reflect heading around y axis
    [ set heading (180 - heading)
      set pressure pressure + abs (dy * mass * speed)
    ]
end 

;; turtle procedure that moves all the particles

to move
  jump (speed / vsplit)
end 

;; turtle procedure to check to see if two particles collide

to check-for-collision
  if count other turtles-here = 1
    [ set tmp-turtle one-of other turtles-here
      if ((who > [who] of tmp-turtle) and (turtle2 != tmp-turtle))
        [ collide ]
    ]
end 

;; turtle procedure for when two particles collide

to collide
  get-turtle2-info
  calculate-velocity-components
  set-new-speed-and-headings
end 

;; turtle gets mass and speed info from turtle it is colliding with

to get-turtle2-info
  set turtle2 tmp-turtle
  set mass2 [mass] of turtle2
  set speed2 [new-speed] of turtle2
  set heading2 [heading] of turtle2
end 

;; calculates new turtle velocity after the collision

to calculate-velocity-components
  set theta (random-float 360)
  set v1l (new-speed * sin (theta - heading))
  set v1t (new-speed * cos (theta - heading))
  set v2l (speed2 * sin (theta - heading2))
  set v2t (speed2 * cos (theta - heading2))
  let vcm (((mass * v1t) + (mass2 * v2t)) / (mass + mass2))
  set v1t (vcm + vcm - v1t)
  set v2t (vcm + vcm - v2t)
end 

;; set new speed and headings of each turtles that has had a collision

to set-new-speed-and-headings
  set new-speed sqrt ((v1t * v1t) + (v1l * v1l))
  set heading (theta - (atan v1l v1t))

  let new-speed2 sqrt ((v2t * v2t) + (v2l * v2l))
  let new-heading (theta - (atan v2l v2t))
  ask turtle2 [
    set new-speed new-speed2
    set heading new-heading
  ]

  recolor
  ask turtle2 [ recolor ]
end 

to recolor  ;; turtle procedure
  ifelse new-speed < (0.5 * initspeed)
    [ set color blue ]
    [ ifelse new-speed > (1.5 * initspeed)
        [ set color red ]
        [ set color green ] ]
end 

;; patch procedure to make a box

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 

;; turtle procedure to give turtles a random position within the confined area

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


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

to piston-up [dist]
  if (dist > 0)
  [ ifelse ((piston-position + dist) < raw-height - 1)
    [ undraw-piston
      set piston-position (piston-position + dist)
      draw-piston ]
    [ undraw-piston
      set piston-position (raw-height - 1)
      draw-piston ]
    set volume ((2 * raw-width) * (piston-position + raw-height))
    set area ((4 * raw-width) + (2 * (piston-position + raw-height)))
  ]
end 

to piston-down [dist]
  if (dist > 0)
  [ ifelse (piston-position - dist) > (2 - raw-height)
    [ undraw-piston
      set piston-position (piston-position - dist)
      ask turtles
      [ if (ycor >= (piston-position - 1))
        [ bounce-off-piston ] ]
      draw-piston ]
    [ undraw-piston
      set piston-position (3 - raw-height)
      ask turtles
      [ if (pycor >= 3 - raw-height)
        [ bounce-off-piston ] ]
      draw-piston ]
    set area ((4 * raw-width) + (2 * (piston-position + raw-height)))
    set volume ((2 * raw-width) * (piston-position + raw-height))
  ]
end 

to draw-piston
  ask patches with [ ((pycor = (round piston-position)) and ((abs pxcor) < raw-width)) ]
    [ set pcolor gray + 2 ]
end 

to undraw-piston
  ask patches with [ (pycor = round piston-position) and ((abs pxcor) < raw-width) ]
    [ set pcolor black ]
end 

to bounce-off-piston  ;; Turtles procedure particle bounces off piston
  ifelse ((((2 * piston-position) - (ycor + 2)) < (1 - raw-height)) or
          (((2 * piston-position) - (ycor + 2)) > (piston-position - 2)))
   [ set ycor ((random (raw-height + piston-position - 2)) - (raw-height - 1)) ]
   [ set ycor ((2 * piston-position) - (ycor + 2)) ]
end 

to calculate-pressure  ;; Observer procedure
  set total-pressure 100 * (sum [pressure] of turtles) / area
  ifelse (length pressure-history < scale)
  [ set pressure-history fput total-pressure pressure-history ]
  [ set pressure-history fput total-pressure but-last pressure-history ]
  set avg-pressure mean pressure-history
  ;; rezero pressures in preparation for the next cycle
  ask turtles [ set pressure 0 ]
end 

;;; plotting procedures

to setup-histograms
  ;; Speed Histogram
  set-current-plot "Speed histogram"
  set-plot-x-range 0 (initspeed * 2)
  set-plot-y-range 0 ceiling (number / 6)
  set-current-plot-pen "average"
  set-histogram-num-bars 45
  set-current-plot-pen "fast"
  set-histogram-num-bars 45
  set-current-plot-pen "slow"
  set-histogram-num-bars 45

  ;; Energy histogram
  set-current-plot "Energy histogram"
  set-plot-x-range 0 (0.5 * (initspeed * 2) * (initspeed * 2) * initmass)
  set-plot-y-range 0 ceiling (number / 6)
  set-current-plot-pen "average"
  set-histogram-num-bars 45
  set-current-plot-pen "fast"
  set-histogram-num-bars 45
  set-current-plot-pen "slow"
  set-histogram-num-bars 45
end 

;; does actual plotting (called in Go)

to do-plotting
  set-current-plot "Volume"
  plot volume
  set-current-plot "Pressure"
  plot avg-pressure
  set-current-plot "Temperature"
  plot avg-energy
  set-current-plot "Pressure vs. Volume"
  plotxy volume avg-pressure
  set-current-plot "Pressure * Volume"
  plot avg-pressure * volume / 1000
end 

;; does actual histograms plotting (called in Go)

to do-histograms
  if (histogram?)
    [ histo-energy
      histo-speed ]
end 

;; draw energy histogram

to histo-energy
  set-current-plot "Energy histogram"
  set-current-plot-pen "average"
  histogram [ energy ] of turtles with [ color = green ]
  set-current-plot-pen "slow"
  histogram [ energy ] of turtles with [ color = blue ]
  set-current-plot-pen "fast"
  histogram [ energy ] of turtles with [ color = red ]
  set-current-plot-pen "avg-energy"
  plot-pen-reset
  draw-vert-line avg-energy
end 

;; draw speed histogram

to histo-speed
  set-current-plot "Speed histogram"
  set-current-plot-pen "average"
  histogram [ speed ] of turtles with [ color = green ]
  set-current-plot-pen "slow"
  histogram [ speed ] of turtles with [ color = blue ]
  set-current-plot-pen "fast"
  histogram [ speed ] of turtles with [ color = red ]
  set-current-plot-pen "avg-speed"
  plot-pen-reset
  draw-vert-line avg-speed
end 

; draws a vertical line at xval on the current-plot with the current plot-pen

to draw-vert-line [xval]
  plotxy xval plot-y-min
  plot-pen-down
  plotxy xval plot-y-max
  plot-pen-up
end 


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