GasLab Second Law

globals [
  fast average slow               ;; current counts
  avg-speed avg-energy            ;; current averages
  avg-speed-init avg-energy-init  ;; initial averages
  vsplit vclock                   ;; clock variables
  left-count right-count          ;; # of particles on each side
  propeller-velocity              ;; current velocity of the propeller
  left-pressure right-pressure    ;; pressure in each chamber
  left-walls right-walls          ;; agentsets of the walls on each side (used when calculating pressure)
  propeller-angle                 ;; current angular position of the propeller
]

turtles-own [
  speed mass energy new-speed           ;; turtle info
  v1t v1l tmp-turtle                    ;; collision info (turtle 1)
  heading2 mass2 speed2 v2t v2l turtle2 ;; collision info (turtle 2)
  theta                                 ;; collision info (both particles)
]

patches-own [
  wall?       ;; is this patch part of the wall?
  pressure    ;; sum of momentums of particles that have bounced here during this time slice
]

to setup [mode]
  clear-all
  make-box
  set vclock  0
  crt number
  [ set new-speed 10.0
    set shape "circle"
    set mass 1.0
    setup-position mode
    recolor ]
  update-variables
  set avg-speed-init avg-speed
  set avg-energy-init avg-energy
  rotate-propeller
  reset-ticks
end 

to setup-position [mode]  ;; turtle procedure
  if mode = "corner"
    [ setxy (- (max-pxcor - 1)) (- (max-pycor - 1))
      set heading random-float 90
      fd random-float 8 ]
  if mode = "one side"
    [ setxy (- (1 + random-float (max-pxcor - 2)))
            (random-float (world-height - 3) + min-pycor + 1) ]
  if mode = "both sides"
    [ setxy (random-float (world-width - 3) + min-pxcor + 1)
            (random-float (world-height - 3) + min-pycor + 1) ]
end 

to update-variables
  ask turtles
    [ set speed new-speed
      set energy (0.5 * speed * speed * mass) ]
  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
  set vsplit (round ((max [speed] of turtles) * 1.2))
  set left-count count turtles with [xcor < 0]
  set right-count count turtles with [xcor >= 0]
  set left-pressure sum [pressure] of left-walls
  ask left-walls [ set pressure 0 ]
  set right-pressure sum [pressure] of right-walls
  ask right-walls [ set pressure 0 ]
end 

to go
  ask turtles [ bounce ]
  ask turtles [ move ]
  ask turtles [ check-for-collision ]
  set vclock (vclock + 1)
  rotate-propeller
  ifelse (vclock = vsplit)
  [
    tick
    set vclock 0
    update-variables
  ]
  [ display ]
end 

to rotate-propeller
  set-current-plot "Propeller"
  plot-pen-reset
  set propeller-angle propeller-angle - propeller-velocity / 2
  plot-pen-up
  plotxy cos propeller-angle sin propeller-angle
  plot-pen-down
  plotxy 0 - cos propeller-angle 0 - sin propeller-angle
  plot-pen-up
  plotxy cos (propeller-angle + 90) sin (propeller-angle + 90)
  plot-pen-down
  plotxy 0 - cos (propeller-angle + 90) 0 - sin (propeller-angle + 90)
  ;; slow down the propeller due to friction
  set propeller-velocity propeller-velocity * 0.999
end 

to bounce ;; turtle procedure
  ; if we're not about to hit a wall (yellow patch),
  ; we don't need to do any further checks
  if [pcolor] of patch-ahead 1 != yellow [ 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
  let pressure-increment mass * speed
  ask patch new-px new-py [ set pressure pressure + pressure-increment ]
  ; check: hitting left, right, or middle wall?
  if (abs new-px = max-pxcor) or (pxcor != 0 and new-px = 0)
    ; if so, reflect heading around x axis
    [ set heading (- heading) ]
  ; check: hitting top or bottom wall?
  if (abs new-py = max-pycor) or (pxcor = 0)
    ; if so, reflect heading around y axis
    [ set heading (180 - heading) ]
end 

to move  ;; turtle procedure
  let old-xcor xcor
  jump (speed / vsplit)
  if (old-xcor < 0) and (xcor >= 0)
    [ set propeller-velocity propeller-velocity + 0.03 * speed ]
  if (old-xcor > 0) and (xcor <= 0)
    [ set propeller-velocity propeller-velocity - 0.03 * speed ]
end 

to check-for-collision  ;; turtle procedure
  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 

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

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

to calculate-velocity-components  ;; turtle procedure
  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))
  ;; CM vel. along dir. theta
  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 < 5.0
    [ set color blue ]
    [ ifelse new-speed > 15.0
        [ set color red ]
        [ set color green ] ]
end 

to make-box
  ask patches
    [ set pressure 0
      set wall? false
      if count neighbors != 8 or
         ((pxcor = 0) and (abs pycor > propeller-radius))
        [ set pcolor yellow
          set wall? true ]
      if (pxcor = 0) and (abs pycor <= propeller-radius)
        [ set pcolor gray ] ]
  set left-walls patches with [wall? and (pxcor < 0)]
  set right-walls patches with [wall? and (pxcor > 0)]
end 

to-report calculate-order
  let x-patches-per-grid-cell (world-width / 5)
  let y-patches-per-grid-cell (world-height / 5)
  let counts-list []
  let gridx -2
  repeat 5
    [ let gridy -2
      repeat 5
        [ let gridcount count turtles with [int (pxcor / x-patches-per-grid-cell) = gridx
                                            and int (pycor / y-patches-per-grid-cell) = gridy]
          set counts-list lput gridcount counts-list
          set gridy gridy + 1 ]
      set gridx gridx + 1 ]
   ;; show counts-list
  report variance counts-list
end 


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