GasLab Free Gas 3D

globals
[
  tick-delta                         ;; how much we advance the tick counter this time through
  max-tick-delta                     ;; the largest tick-delta is allowed to be
  min-tick-delta                     ;; the smallest tick-delta is allowed to be
  init-avg-speed init-avg-energy     ;; initial averages
  avg-speed avg-energy               ;; current averages
  fast medium slow                   ;; current counts
  percent-fast percent-medium        ;; percentage of the counts
  percent-slow                       ;; percentage of the counts

  collision-times                    ;; a list that of times of pending collisions
]

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
  ca
  set-default-shape particles "circle"
  set tick-delta .01
  set min-tick-delta .0000001
  make-particles
  check-initial-positions
  update-variables
  set init-avg-speed avg-speed
  set init-avg-energy avg-energy
  setup-histograms
  do-plotting
  if trace?
  [ ask one-of particles [ pd ] ]
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 ]
  ]

  ifelse( empty? collision-times )
  [ set collision-times lput ( tick-delta ) collision-times ]
  [ 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 ) ] [
        if collision-with > who [ ;; so that we don't collide the same particles twice
            collide ( turtle collision-with )
            set last-collision collision-with
            ask turtle collision-with [ set last-collision [who] of myself ]
        ]
      ]
    ] [
      ask particles [ jump speed * tick-delta ]
      tick-advance tick-delta
    ]

  ask particles [
    if last-collision != nobody
    [
      if distance turtle last-collision > ( ( ( [size] of turtle last-collision ) / 2 ) + ( size / 2 ) ) * 1.5
      [ set last-collision nobody ]
    ]
  ]

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

  display
end 

to update-variables
  set medium count particles with [ speed < ( 1.5 * 10 ) and speed > ( 0.5 * 10 ) ]
  set slow count particles with [ speed < ( 0.5 * 10 ) ]
  set fast count particles with [ speed > ( 1.5 * 10 ) ]
  set percent-medium (medium / count particles) * 100
  set percent-slow (slow / count particles) * 100
  set percent-fast (fast / count particles) * 100
  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 xcor
  let my-y ycor
  let my-z zcor
  let my-particle-size size
  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 )
  let my-last-collision last-collision

  ask other particles with [who != my-last-collision]
  [
    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

    ifelse ( abs ( xcor - my-x ) <= abs ( ( xcor - my-x ) - world-width ) )
      [ set dpx (xcor - my-x) ]
      [ set dpx (xcor - my-x) - world-width ]  ;; relative distance between particles in the x direction
    ifelse ( abs ( ycor - my-y ) <= abs ( ( ycor - my-y ) - world-height ) )
      [ set dpy (ycor - my-y) ]
      [ set dpy (ycor - my-y) - world-height ]    ;; relative distance between particles in the y direction
    ifelse ( abs ( zcor - my-z ) <= abs ( ( zcor - my-z ) - world-depth ) )
      [ set dpz (zcor - my-z) ]
      [ set dpz (zcor - my-z) - world-depth ]       ;; relative distance between particles in the z direction

    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 (((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) + (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)

    ;; 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
      [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.
    ;; to keep the model from slowing down too much, if the particles are going to collide
    ;; at a time before min-tick-delta, just collide them a min-tick-delta instead

    if ( time-to-collision < tick-delta and time-to-collision > min-tick-delta ) [
      set collision-with ( [who] of 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 ( [who] of myself )
      set collision-time ( min-tick-delta )
      set collision-times ( lput ( min-tick-delta ) collision-times )
    ]
  ]
end 

to collide [ particle2 ] ;; turtle 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 particle1totheta ortho-projection vx vy vz tx ty tz

  ;; express particle1's movement along theta in terms of xyz
  let x1totheta particle1totheta * tx
  let y1totheta particle1totheta * ty
  let z1totheta particle1totheta * 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 x1opptheta ( ( vx ) - ( x1totheta ) )
  let y1opptheta ( ( vy ) - ( y1totheta ) )
  let z1opptheta ( ( vz ) - ( z1totheta ) )

  ;; do the same for particle2
  let particle2totheta ortho-projection [vx] of particle2 [vy] of particle2 [vz] of particle2 tx ty tz

  let x2totheta particle2totheta * tx
  let y2totheta particle2totheta * ty
  let z2totheta particle2totheta * tz

  let x2opptheta ( ( [vx] of particle2 ) - ( x2totheta ) )
  let y2opptheta ( ( [vy] of particle2 ) - ( y2totheta ) )
  let z2opptheta ( ( [vz] of particle2 ) - ( z2totheta ) )

  ;; calculate the velocity of the center of mass along theta
  let vcm ( ( ( mass * particle1totheta ) + ( [mass] of particle2 * particle2totheta ) )
      / ( mass + [mass] of particle2 ) )

  ;; switch momentums along theta
  set particle1totheta (2 * vcm - particle1totheta)
  set particle2totheta (2 * vcm - particle2totheta)

  ;; determine the x, y, z components of each particle's new velocities
  ;; in the direction of theta
  set x1totheta particle1totheta * tx
  set y1totheta particle1totheta * ty
  set z1totheta particle1totheta * tz

  set x2totheta particle2totheta * tx
  set y2totheta particle2totheta * ty
  set z2totheta particle2totheta * 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 x1totheta + x1opptheta
  set vy y1totheta + y1opptheta
  set vz z1totheta + z1opptheta
  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 )
  if particle-color = "red-green-blue" [ recolor ]
  if particle-color = "purple shades" [ recolorshade ]
  if particle-color = "one color" [ recolornone ]

  ask particle2 [
    set vx x2totheta + x2opptheta
    set vy y2totheta + y2opptheta
    set vz z2totheta + z2opptheta
    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 )
    if particle-color = "red-green-blue" [ recolor ]
    if particle-color = "purple shades" [ recolorshade ]
    if particle-color = "one color" [ recolornone ]
  ]
end 


;;;
;;; drawing procedures
;;;

;; creates initial particles

to make-particles
  create-particles number-of-particles
  [
    setup-particle
    random-position
    if particle-color = "red-green-blue" [ recolor ]
    if particle-color = "purple shades" [ recolorshade ]
    if particle-color = "one color" [ recolornone ]
  ]
end 

to setup-particle  ;; particle procedure
  set speed init-particle-speed
  set mass particle-mass
  set energy (0.5 * mass * (speed ^ 2))
end 

;; makes sure particles aren't overlapped at setup

to check-initial-positions
  let check-again? false
  ask particles [
    if particle-overlap? [
      random-position
      set check-again? true
    ]
  ]
  if check-again? [ check-initial-positions ]
end 

to-report particle-overlap? ; particle procedure
  let me ( self )
  let overlap false
  ask other particles [
    if distance ( me ) < ( ( size / 2 ) + ( [size] of me / 2 ) + .1 ) [
      set overlap true
    ]
  ]
  report overlap
end 


;; place particle at random location inside the box.

to random-position ;; particle procedure
  setxyz ((1 + min-pxcor) + random-float ((2 * max-pxcor) - 2))
         ((1 + min-pycor) + random-float ((2 * max-pycor) - 2))
         ((1 + min-pzcor) + random-float ((2 * max-pzcor) - 2)) ;; added for 3d
  set heading random-float 360
  set pitch random-float 360
end 

to move  ;; particle procedure
  jump (speed * tick-delta)
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 

to recolorshade ;; particle procedure
  ifelse speed < 27
  [ set color 111 + speed / 3 ]
  [ set color 119.999 ]
end 

to recolornone ;; particle procedure
  set color blue + 1
end 


;;;
;;; math procedures
;;;

;; 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 ( speed * sin ( heading ) * cos ( pitch ) )
  set vy ( speed * cos ( heading ) * cos ( pitch ) )
  set vz ( speed * sin ( pitch ) )
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 ]
  let xvel sin( vector-angle ) * abs( cos( vector-pitch ) ) * vector-speed
  report xvel
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 ]
  let yvel cos( vector-angle ) * abs( cos( vector-pitch ) ) * vector-speed
  report yvel
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 ]
  let zvel ( sin( vector-pitch ) * vector-speed )
  report zvel
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 round asin ( z / ( vspeed x y z ) )
end 

;; called by ortho-projection

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 ortho-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 setup-histograms
  set-current-plot "Speed Histogram"
  set-plot-x-range 0 (init-particle-speed * 2)
  set-plot-y-range 0 ceiling (number-of-particles / 6)
  set-current-plot-pen "medium"
  set-histogram-num-bars 40
  set-current-plot-pen "slow"
  set-histogram-num-bars 40
  set-current-plot-pen "fast"
  set-histogram-num-bars 40
  set-current-plot-pen "init-avg-speed"
  draw-vert-line init-avg-speed

  set-current-plot "Energy Histogram"
  set-plot-x-range 0 (0.5 * (init-particle-speed * 2) * (init-particle-speed * 2) * particle-mass)
  set-plot-y-range 0 ceiling (number-of-particles / 6)
  set-current-plot-pen "medium"
  set-histogram-num-bars 40
  set-current-plot-pen "slow"
  set-histogram-num-bars 40
  set-current-plot-pen "fast"
  set-histogram-num-bars 40
  set-current-plot-pen "init-avg-energy"
  draw-vert-line init-avg-energy
end 

to do-plotting
  set-current-plot "Speed Counts"
  set-current-plot-pen "fast"
  plotxy ticks percent-fast
  set-current-plot-pen "medium"
  plotxy ticks percent-medium
  set-current-plot-pen "slow"
  plotxy ticks percent-slow

  plot-histograms
end 

to plot-histograms
  set-current-plot "Energy histogram"
  set-current-plot-pen "fast"
  histogram [ energy ] of particles with [ speed > ( 1.5 * 10 ) ]
  set-current-plot-pen "medium"
  histogram [ energy ] of particles with [ speed < ( 1.5 * 10 ) and speed > ( 0.5 * 10 ) ]
  set-current-plot-pen "slow"
  histogram [ energy ] of particles with [ speed < ( 0.5 * 10 ) ]
  set-current-plot-pen "avg-energy"
  plot-pen-reset
  draw-vert-line avg-energy

  set-current-plot "Speed histogram"
  set-current-plot-pen "fast"
  histogram [ speed ] of particles with [ speed > ( 1.5 * 10 ) ]
  set-current-plot-pen "medium"
  histogram [ speed ] of particles with [ speed < ( 1.5 * 10 ) and speed > ( 0.5 * 10 ) ]
  set-current-plot-pen "slow"
  histogram [ speed ] of particles with [ speed < ( 0.5 * 10 ) ]
  set-current-plot-pen "avg-speed"
  plot-pen-reset
  draw-vert-line avg-speed
end 

;; histogram procedure

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


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