GasLab Single Collision 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
  collision-times     ;; a list that of times of pending collisions

  ;; coordinates of the center of mass
  x-center
  y-center
  z-center
]

breed [ centers-of-mass center-of-mass ]

breed [ particles particle ]
particles-own
[
  ;; properties specific to particles
  speed
  energy
  mass

  collision-time ;; time in the future of the next collision
  collision-with ;; particle to collide with
  last-collision ;; particle we last collided with so we don't collide with the same particle many times

  ;; component vectors of the speed and direction of the particle
  x-vector
  y-vector
  z-vector
]

to setup
  clear-all
  create-particles 1 [
    set color pink
    set speed init-pink-speed
    setxyz 0 0 0
    set heading 0
    set pitch 0
    bk speed
    if wiggle? [ wiggle3d 5 ]
    set mass mass-of-pink
    set energy 0.5 * mass * speed ^ 2
    pen-down
    set last-collision nobody
  ]
  create-particles 1 [
    set color blue
    set speed init-blue-speed
    setxyz 0 0 0
    set heading collision-angle-xy
    set pitch collision-angle-z
    bk speed
    if wiggle? [ wiggle3d 5 ]
    set mass mass-of-blue
    set energy 0.5 * mass * speed ^ 2
    pen-down
    set last-collision nobody
  ]
  ask particles [
    set shape "circle"
    set size sqrt mass
  ]

  let total-mass sum [mass] of particles
  set x-center ( sum [ xcor * mass ] of particles ) / total-mass
  set y-center ( sum [ ycor * mass ] of particles ) / total-mass
  set z-center ( sum [ zcor * mass ] of particles ) / total-mass

  create-centers-of-mass 1 [
    set shape "x"
    set size 1
    set color gray
    setxyz x-center
           y-center
           z-center
    ifelse show-center-of-mass?
    [ pen-down ]
    [ hide-turtle ]
  ]

  set tick-delta .01
  set min-tick-delta 0
end 

to go-mode
  if run-mode = "one-collision" [ go-once ]
  if run-mode = "all-collision-angles" [ go-all-angles ]
  if run-mode = "all-collision-pitches" [ go-all-angles-and-pitches ]
end 

to go-once
  while [ go ]
    []
end 

to go-all-angles
  let angle 15
  set wiggle? false
  repeat 23 [
    set collision-angle-xy angle
    setup
    while [ go ]
      []
    set angle angle + 15
  ]
end 

to go-all-angles-and-pitches
  let zangle -90
  set wiggle? false
  repeat 11 [
    set collision-angle-z zangle
    setup
    while [ go ] [ ]
    set zangle zangle + 15
  ]
end 

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

  set x-center ( sum [ xcor * mass ] of particles ) / ( sum [mass] of particles )
  set y-center ( sum [ ycor * mass ] of particles ) / ( sum [mass] of particles )
  set z-center ( sum [ zcor * mass ] of particles ) / ( sum [mass] of particles )

  ask centers-of-mass
  [
    ifelse show-center-of-mass?  ;; marks a gray path along the particles' center of mass
      [ show-turtle
        pen-down ]
      [ hide-turtle
        pen-up ]
    setxyz x-center y-center z-center
  ]

  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
  [
    tick-advance first collision-times
    ;; precheck that all the particles can move before we move any
    ;; this is so that the center of mass doesn't get slightly
    ;; off path because one particle moves and the other doesn't
    if any? particles with [ not can-move-cube? (speed * first collision-times) ]
      [ report false ]
    ask particles
    [
      jump speed * first collision-times
    ] ; most particles to first collision
    ask particle 0 [
      collide particle 1
      set last-collision collision-with
      ask particle 1 [ set last-collision myself ]
    ]
  ][
    ;; precheck that all the particles can move before we move any
    ;; this is so that the center of mass doesn't get slightly
    ;; off path because one particle moves and the other doesn't
    if any? particles with [ not can-move-cube? (speed * tick-delta) ]
      [ report false ]
    ask particles
    [
       jump speed * tick-delta
    ]
    tick-advance tick-delta
  ]

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

  do-plotting
  display
  report true
end 

;; we don't yet have different topologies in 3D so we have
;; to prevent the particles from wrapping around the world
;; manually (just like in pre-3.1 NetLogo).  We don't want the
;; particles to wrap around the world because it is confusing
;; when the center-of-mass is turned on.

to-report can-move-cube? [dist]
  let x xcor + (dx * dist)
  let y ycor + (dy * dist)
  let z zcor + (dz * dist)
  report not (x < min-pxcor + 0.5 or x > max-pxcor - 0.5 or
              y < min-pycor + 0.5 or y > max-pxcor - 0.5 or
              z < min-pzcor + 0.5 or z > max-pzcor - 0.5)
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-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

    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 ([size] of myself / 2 ) + (size / 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 > 0 [
      set collision-with myself
      set collision-time time-to-collision
      set collision-times lput time-to-collision 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 particle1-to-theta orth-projection x-vector y-vector z-vector 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 ( ( x-vector ) - ( x1-to-theta ) )
  let y1-opp-theta ( ( y-vector ) - ( y1-to-theta ) )
  let z1-opp-theta ( ( z-vector ) - ( z1-to-theta ) )

  ;; do the same for particle2
  let particle2-to-theta orth-projection [x-vector] of particle2 [y-vector] of particle2 [z-vector] 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 ( ( [x-vector] of particle2 ) - ( x2-to-theta ) )
  let y2-opp-theta ( ( [y-vector] of particle2 ) - ( y2-to-theta ) )
  let z2-opp-theta ( ( [z-vector] 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 x-vector x1-to-theta + x1-opp-theta
  set y-vector y1-to-theta + y1-opp-theta
  set z-vector z1-to-theta + z1-opp-theta
  set heading vheading x-vector y-vector z-vector
  set pitch vpitch x-vector y-vector z-vector
  set speed vspeed x-vector y-vector z-vector
  set energy ( 0.5 * mass * speed ^ 2 )
  set shape "circle 2"
  stamp
  set shape "circle"

  ask particle2 [
    set x-vector x2-to-theta + x2-opp-theta
    set y-vector y2-to-theta + y2-opp-theta
    set z-vector z2-to-theta + z2-opp-theta
    set heading vheading x-vector y-vector z-vector
    set pitch vpitch x-vector y-vector z-vector
    set speed vspeed x-vector y-vector z-vector
    set energy ( 0.5 * mass * speed ^ 2 )
    set shape "circle 2"
    stamp
    set shape "circle"
  ]
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 x-vector speed * sin heading * cos pitch
  set y-vector speed * cos heading * cos pitch
  set z-vector 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 ]
  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 round 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 

;; wiggle up to angle in random direction

to wiggle3d [ angle ]
  roll-right random-float 360
  tilt-up acos (1 - random-float ( 1 - cos angle ) )
end 

;;; plotting procedure

to do-plotting
  set-current-plot "Speeds"

  ifelse [speed] of turtle 0 = [speed] of turtle 1
  [
    set-current-plot-pen "both"
    plotxy ticks [speed] of turtle 0
  ]
  [
    set-current-plot-pen "pink"
    plotxy ticks [speed] of turtle 0
    set-current-plot-pen "blue"
    plotxy ticks [speed] of turtle 1
  ]
end 


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