GasLab Two Gas 3D

globals
[
  tick-delta                         ;; how much we advance the tick counter this time through
  min-tick-delta                     ;; the smallest tick-delta is allowed to be
  init-avg-speed init-avg-energy     ;; initial averages
  collision-times                    ;; a list that of times of pending collisions

  ;; averages from the end of the last tick
  avg-speed avg-energy
  avg-speed-magenta
  avg-speed-cyan
  avg-energy-magenta
  avg-energy-cyan

  open?                              ;; is there an opening in the divider?
]

breed [ dividers divider ]
breed [ walls wall ]

breed [ particles particle ]
particles-own
[
  vx vy vz                   ;; velocities rel axes
  speed mass energy          ;; particle info
  collision-time             ;; to determine when collision is
  collision-with             ;; to determine who the collision is with
  last-collision             ;; so they don't collide with one another many times
]

to setup
  clear-all
  ;; the wall shape is a custom 3D shape contained in "wall.txt"
  load-shapes-3d "wall.txt"
  set-default-shape particles "circle"
  set-default-shape walls "wall"
  set tick-delta .01
  set min-tick-delta .0000001
  make-box
  make-particles
  update-variables
  set init-avg-speed avg-speed
  set init-avg-energy avg-energy
  set open? false
  set avg-speed-cyan mean [speed] of turtles with [color = cyan]
  set avg-speed-magenta mean [speed] of turtles with [color = magenta]
  set avg-energy-cyan mean [energy] of turtles with [color = cyan]
  set avg-energy-magenta mean [ energy ] of turtles with [color = magenta]
  do-plotting
end 

to go
  set collision-times [] ;; empty this out for new input
  ask particles
  [
    set collision-time tick-delta
    set collision-with nobody
      if collide? [
        detect-collisions
        detect-wall-collisions
      ]
  ]
  set collision-times sort collision-times

  ifelse first collision-times < tick-delta   ;; if something will collide before the tick
  [
    ask particles [ jump speed * first collision-times ] ; most particles to first collision
    tick-advance first collision-times ;; now, collide all the particles that are ready
    ask particles with [ collision-time = first collision-times ]
    [
      ifelse is-particle? collision-with
      [
        if collision-with > self [ ;; so that we don't collide the same particles twice
          collide collision-with
          set last-collision collision-with
          ask collision-with [ set last-collision myself ]
        ]
      ]
      [ wall-collide collision-with ]
    ]
  ]
  [
    ask particles [ jump speed * tick-delta ]
    tick-advance tick-delta
  ]

  ask particles [
    if last-collision != nobody and is-particle? last-collision
    [
      if distance last-collision > ( ( ( [size] of last-collision ) / 2 ) + ( size / 2 ) ) * 1.1
        [ set last-collision nobody ]
    ]
  ]

  if floor ticks > floor (ticks - tick-delta)
    [ update-variables ]

  set avg-speed-cyan mean [speed] of turtles with [color = cyan]
  set avg-speed-magenta mean [speed] of turtles with [color = magenta]
  set avg-energy-cyan mean [energy] of turtles with [color = cyan]
  set avg-energy-magenta mean [energy] of turtles with [color = magenta]

  do-plotting
  display
end 

to update-variables
  set avg-speed  mean [speed] of particles
  set avg-energy  mean [energy] of particles
end 

;;;
;;; distance and collision procedures
;;;

to detect-collisions ;; particle procedure

;; detect-collisions 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-speed x-velocity heading pitch speed
  let my-y-speed y-velocity heading pitch speed
  let my-z-speed z-velocity pitch speed

  ask other particles with [self != [last-collision] of myself]
  [
    let dpx 0
    let dpy 0
    let dpz 0

   ;; since our world is wrapped, we can't just use calcs like xcor - my-x. Instead, we take the smallest
   ;; of either the wrapped or unwrapped distance for each dimension

    set dpx xcor - [xcor] of myself
    set dpy ycor - [ycor] of myself
    set dpz zcor - [zcor] of myself

    let x-speed x-velocity heading pitch speed  ;; speed of other particle in the x direction
    let y-speed y-velocity heading pitch speed  ;; speed of other particle in the y direction
    let z-speed z-velocity pitch speed          ;; speed of other particle in the z 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 dvz z-speed - my-z-speed ;; relative speed difference between particles in the z direction

    let sum-r ([size] of myself / 2) + (size / 2) ;; sum of both particle radii

    let p-squared   ((dpx * dpx) + (dpy * dpy) + (dpz * dpz)) - (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) + (dpz * dvz))  ;;the vector product of the position times the velocity
    let v-squared (dvx * dvx) + (dvy * dvy) + (dvz * dvz) ;; 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)

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

    let time-to-collision  -1

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

    if time-to-collision < tick-delta and time-to-collision > min-tick-delta
    [
      set collision-with myself
      set collision-time time-to-collision
      set collision-times lput time-to-collision collision-times
    ]
    if time-to-collision < min-tick-delta and time-to-collision > 0
    [
      set collision-with myself
      set collision-time min-tick-delta
      set collision-times lput min-tick-delta collision-times
    ]
  ]
end 

to detect-wall-collisions ;; particle procedure
  update-component-vectors
  let my-vx vx * tick-delta
  let my-vy vy * tick-delta
  let my-vz vz * tick-delta

  detect-wall-collision "xy"
                        ( wall-max-pzcor - ( size / 2 ) )
                        ( wall-min-pzcor + ( size / 2 ) )
                        zcor
                        my-vz
  detect-wall-collision "yz"
                        ( wall-max-pxcor - ( size / 2 ) )
                        ( wall-min-pxcor + ( size / 2 ) )
                        xcor
                        my-vx
  detect-wall-collision "xz"
                        ( wall-max-pycor - ( size / 2 ) )
                        ( wall-min-pycor + ( size / 2 ) )
                        ycor
                        my-vy
  detect-divider-collision

  if collision-time < min-tick-delta [
    set collision-time min-tick-delta
  ]
  set collision-times lput collision-time collision-times
end 

;; detect-wall-collision plane of wall, wall cor, wall cor, cor of particle, speed of particle

to detect-wall-collision [ plane max-wall min-wall cor cor-speed ] ;; particle procedure
  if ( cor + cor-speed ) > max-wall or
     ( cor + cor-speed ) < min-wall
  [
    let distance-to-wall abs( max-wall - cor )
    let time-fraction ( distance-to-wall / cor-speed )
    if ( time-fraction * tick-delta ) < collision-time and last-collision != plane
    [
      set collision-time time-fraction * tick-delta
      set collision-with plane
    ]
  ]
end 

to detect-divider-collision ;; particle procedure
  let my-vx vx * tick-delta
  if xcor > 0 and my-vx < 0 and ( xcor + my-vx ) < (size / 2)
  [
    let distance-to-wall xcor - (size / 2)
    let time-fraction distance-to-wall / my-vx
    if ( time-fraction * tick-delta ) < collision-time and last-collision != "divider"
    [
      ;; where particle will be once it hits divider
      let future-ycor ycor + ( vy * time-fraction * tick-delta )
      let future-zcor zcor + ( vz * time-fraction * tick-delta )
      if not ( ( abs( future-ycor ) <= ( wall-max-pycor * opening-size / 100 ) and
                 abs( future-zcor ) <= ( wall-max-pzcor * opening-size / 100 ) ) and
               open? )
      [
        set collision-time time-fraction * tick-delta
        set collision-with "divider"
      ]
    ]
  ]
  if xcor < 0 and my-vx > 0 and ( xcor + my-vx ) > (- size / 2)
  [
    let distance-to-wall abs( xcor + (size / 2) )
    let time-fraction distance-to-wall / my-vx
    if ( time-fraction * tick-delta ) < collision-time and last-collision != "divider"
    [
      let future-ycor ycor + ( vy * time-fraction * tick-delta ) ;; where particle will be once it hits divider
      let future-zcor zcor + ( vz * time-fraction * tick-delta )
      if not ( ( abs( future-ycor ) <= ( wall-max-pycor * opening-size / 100 ) and
                 abs( future-zcor ) <= ( wall-max-pzcor * opening-size / 100 ) ) and
               open? )
      [
        set collision-time time-fraction * tick-delta
        set collision-with "divider"
      ]
    ]
  ]
end 

to collide [ particle2 ] ;; particle procedure
  update-component-vectors
  ask particle2 [ update-component-vectors ]

  ;; find heading and pitch from the center of particle1 to the center of particle2
  let theading towards particle2
  let tpitch towards-pitch particle2

  ;; use these to determine the x, y, z components of theta
  let tx x-velocity theading tpitch 1
  let ty y-velocity theading tpitch 1
  let tz z-velocity tpitch 1

  ;; find the speed of particle1 in the direction of n
  let particle1-to-theta orth-projection vx vy vz tx ty tz

  ;; express particle1's movement along theta in terms of xyz
  let x1-to-theta particle1-to-theta * tx
  let y1-to-theta particle1-to-theta * ty
  let z1-to-theta particle1-to-theta * tz

  ;; now we can find the x, y and z components of the particle's velocity that
  ;; aren't in the direction of theta by subtracting the x, y, and z
  ;; components of the velocity in the direction of theta from the components
  ;; of the overall velocity of the particle
  let x1-opp-theta vx - x1-to-theta
  let y1-opp-theta vy - y1-to-theta
  let z1-opp-theta vz - z1-to-theta

  ;; do the same for particle2
  let particle2-to-theta orth-projection [vx] of particle2 [vy] of particle2 [vz] of particle2 tx ty tz

  let x2-to-theta particle2-to-theta * tx
  let y2-to-theta particle2-to-theta * ty
  let z2-to-theta particle2-to-theta * tz

  let x2-opp-theta [vx] of particle2 - x2-to-theta
  let y2-opp-theta [vy] of particle2 - y2-to-theta
  let z2-opp-theta [vz] of particle2 - z2-to-theta

  ;; calculate the velocity of the center of mass along theta
  let vcm ( ( mass * particle1-to-theta ) + ( [mass] of particle2 * particle2-to-theta ) )
            / ( mass + [mass] of particle2 )

  ;; switch momentums along theta
  set particle1-to-theta 2 * vcm - particle1-to-theta
  set particle2-to-theta 2 * vcm - particle2-to-theta

  ;; determine the x, y, z components of each particle's new velocities
  ;; in the direction of theta
  set x1-to-theta particle1-to-theta * tx
  set y1-to-theta particle1-to-theta * ty
  set z1-to-theta particle1-to-theta * tz

  set x2-to-theta particle2-to-theta * tx
  set y2-to-theta particle2-to-theta * ty
  set z2-to-theta particle2-to-theta * tz

  ;; now, we add the new velocities along theta to the unchanged velocities
  ;; opposite theta to determine the new heading, pitch, and speed of each particle
  set vx x1-to-theta + x1-opp-theta
  set vy y1-to-theta + y1-opp-theta
  set vz z1-to-theta + z1-opp-theta
  set heading vheading vx vy vz
  set pitch vpitch vx vy vz
  set speed vspeed vx vy vz
  set energy 0.5 * mass * speed ^ 2

  ask particle2 [
    set vx x2-to-theta + x2-opp-theta
    set vy y2-to-theta + y2-opp-theta
    set vz z2-to-theta + z2-opp-theta
    set heading vheading vx vy vz
    set pitch vpitch vx vy vz
    set speed vspeed vx vy vz
    set energy 0.5 * mass * speed ^ 2
  ]
end 

to open-middle
  if opening-size = 20 [
    ask dividers [ set shape "opening20" ]
  ]
  if opening-size = 40 [
    ask dividers [ set shape "opening40" ]
  ]
  if opening-size = 60 [
    ask dividers [ set shape "opening60" ]
  ]
  if opening-size = 80 [
    ask dividers [ set shape "opening80" ]
  ]
  if opening-size = 100 [
    ask dividers [ die ]
  ]
  set open? true
end 

to close-middle
  ask dividers [ set shape "flash" ]
  set open? false
end 

to wall-collide [ collision-wall ] ;; particle procedure
  update-component-vectors

  ifelse collision-wall = "yz"
  [
    set heading vheading (- vx ) vy vz
  ][
  ifelse collision-wall = "xz"
  [
    set heading vheading vx (- vy ) vz
  ][
  ifelse collision-wall = "xy"
  [
    set pitch vpitch vx vy ( - vz )
   ]
   [
     set heading vheading (- vx ) vy vz
   ] ] ]
end 

;;;
;;; drawing procedures
;;;

;; creates box

to make-box
  create-walls 1 [ ;; bottom wall
    set heading 0
    set zcor wall-min-pzcor
  ]
  create-walls 1 [ ;; top wall
    set heading 0
    set pitch 180
    set zcor wall-max-pzcor
  ]
  create-walls 1 [ ;; upper wall
    set heading 0
    set pitch 90
    set ycor wall-max-pycor
  ]
  create-walls 1 [ ;; lower wall
    set heading 0
    set pitch -90
    set ycor wall-min-pycor
  ]
  create-walls 1 [ ;; right wall
    set heading 0
    set roll -90
    set xcor wall-max-pxcor
  ]
  create-walls 1 [ ;; left wall
    set heading 0
    set roll 90
    set xcor wall-min-pxcor
  ]
  create-dividers 1 [ ;; center wall
    set heading 0
    set color grey
    set roll 90
    set xcor 0
    set shape "flash"
    set size box-width
  ]
  ask walls [
    set heading 0
    set color grey + random 3
    set size box-width
    set shape "wall"
  ]
end 

;; creates initial particles

to make-particles
  create-particles num-magentas
  [
    setup-particle magenta-init-speed magenta-mass magenta
    random-position "left-half"
  ]
  create-particles num-cyans
  [
    setup-particle cyan-init-speed cyan-mass cyan
    random-position "right-half"
  ]
  check-initial-positions 0
  check-center-divider
end 

to setup-particle [ my-speed my-mass my-color ] ;; particle procedure
  set speed my-speed
  set mass my-mass
  set energy 0.5 * mass * (speed ^ 2)
  set color my-color
  set size mass ^ 0.33
end 

;; makes sure particles aren't overlapped at setup

to check-initial-positions [iterations]
  let check-again? false
  ask particles
  [
    if particle-overlap?
    [
      ifelse color = cyan
      [ random-position "right-half" ]
      [ random-position "left-half" ]
      set check-again? true
    ]
  ]
  ifelse iterations < 50
  [
    if check-again?
      [ check-initial-positions iterations + 1 ]
  ]
  [
    beep
    user-message "Not enough room for all these particles!"
  ]
end 

;; makes sure particles don't go through center wall at setup

to check-center-divider
  let check-again? false
  ask particles
  [
    if abs( xcor ) < ( size / 2 )
    [
      ifelse color = cyan
      [ random-position "right-half" ]
      [ random-position "left-half" ]
      set check-again? true
    ]
  ]
  if check-again?
  [
    check-initial-positions 0
    check-center-divider
  ]
end 

to-report particle-overlap? ;; particle procedure
  report any? other particles with [ distance myself <= ((size + [size] of myself) / 2 ) ]
end 

;; place particle at random location inside the box.

to random-position [ side ] ;; particle procedure
  ifelse side = "left-half" [
    setxyz ( random-float ( ( box-width / 2 ) - 2 ) +
             ( wall-min-pxcor + 1 ) )
           ( random-float ( box-height - 2 ) +
             ( wall-min-pycor + 1 ) )
           ( random-float ( box-depth - 2 ) +
             ( wall-min-pzcor + 1 ) )
  ] [
    setxyz ( random-float ( ( box-width / 2 ) - 2 ) +
             ( 1 ) )
           ( random-float ( box-height - 2 ) +
             ( wall-min-pycor + 1 ) )
           ( random-float ( box-depth - 2 ) +
             ( wall-min-pzcor + 1 ) )
  ]
  tilt-up asin (1.0 - random-float 2.0)
  roll-right random-float 360
  update-component-vectors
  set heading vheading vx vy vz
  set pitch vpitch vx vy vz
end 

;;;
;;; math procedures
;;;

;; consider the desired box-size

to-report box-width
  report ( world-width - 1 ) * ( box-size / 100 )
end 

to-report box-height
  report ( world-height - 1 ) * ( box-size / 100 )
end 

to-report box-depth
  report ( world-depth - 1 ) * ( box-size / 100 )
end 

to-report wall-max-pxcor
  report max-pxcor * ( box-size / 100 )
end 

to-report wall-max-pycor
  report max-pycor * ( box-size / 100 )
end 

to-report wall-max-pzcor
  report max-pzcor * ( box-size / 100 )
end 

to-report wall-min-pxcor
  report min-pxcor * ( box-size / 100 )
end 

to-report wall-min-pycor
  report min-pycor * ( box-size / 100 )
end 

to-report wall-min-pzcor
  report min-pzcor * ( box-size / 100 )
end 

;; makes sure that the values stored in vx, vy, vz actually reflect
;; the appropriate heading, pitch, speed

to update-component-vectors ;; particle procedure
  set vx x-velocity heading pitch speed
  set vy y-velocity heading pitch speed
  set vz z-velocity pitch speed
end 

;; reports velocity of a vector at a given angle and pitch
;; in the direction of x.

to-report x-velocity [ vector-angle vector-pitch vector-speed ]
  report sin( vector-angle ) * abs( cos( vector-pitch ) ) * vector-speed
end 

;; reports velocity of a vector at a given angle and pitch
;; in the direction of y.

to-report y-velocity [ vector-angle vector-pitch vector-speed ]
  report cos( vector-angle ) * abs( cos( vector-pitch ) ) * vector-speed
end 

;; reports velocity of a vector at a given angle and pitch
;; in the direction of z.

to-report z-velocity [ vector-pitch vector-speed ]
  report sin( vector-pitch ) * vector-speed
end 

;; reports speed of a vector given xyz coords

to-report vspeed [ x y z ]
  report sqrt( x ^ 2 + y ^ 2 + z ^ 2 )
end 

;; reports xt heading of a vector given xyz coords

to-report vheading [ x y z ]
  report atan x y
end 

;; reports pitch of a vector given xyz coords

to-report vpitch [ x y z ]
  report asin ( z / ( vspeed x y z ) )
end 

;; called by orthprojection

to-report dot-product [ x1 y1 z1 x2 y2 z2 ]
  report ( x1 * x2 ) + ( y1 * y2 ) + ( z1 * z2 )
end 

;; component of 1 in the direction of 2 (Note order)

to-report orth-projection [ x1 y1 z1 x2 y2 z2 ]
  let dproduct dot-product x1 y1 z1 x2 y2 z2
  let speed-of-2 vspeed x2 y2 z2
  ;; if speed is 0 then there's no projection anyway
  ifelse speed-of-2 > 0
  [ report dproduct / speed-of-2 ]
  [ report 0 ]
end 

;;;
;;; plotting procedures
;;;

to do-plotting
  set-current-plot "Average Speeds"
  set-current-plot-pen "cyan"
  plotxy ticks avg-speed-cyan
  set-current-plot-pen "magenta"
  plotxy ticks avg-speed-magenta

  set-current-plot "Average Energies"
  set-current-plot-pen "cyan"
  plotxy ticks avg-energy-cyan
  set-current-plot-pen "magenta"
  plotxy ticks avg-energy-magenta
end 


; Copyright 2007 Uri Wilensky. All rights reserved.
; The full copyright notice is in the Information tab.