CellularMetabolismLabyrinth

globals
[
  tick-delta                                  ;; how much we advance the tick counter this time through
  max-tick-delta                              ;; the largest tick-delta is allowed to be
  box-edge                                    ;; distance of box edge from axes
  init-avg-speed init-avg-energy              ;; initial averages
  avg-speed avg-energy                        ;; current averages
  fast medium slow                            ;; current counts
  percent-slow percent-medium percent-fast    ;;  percentage of current counts
  aggregate-list                              ;; list of the sums of the temperature at each height
  i
  z
  yboltz
  R
  temp
]

breed [ particles particle ]
breed [ flashes flash ]

flashes-own [birthday]
patches-own [d]

particles-own
[
  speed mass energy          ;; particle info
  last-collision
]

to setup
  clear-all
  set-default-shape particles "circle"
  set-default-shape flashes "plane"
  set max-tick-delta 0.1073
  set R 1.9858775 / 1000
  set temp 310
  make-particles
  ;; box has constant size...
  set box-edge (max-pxcor)
  ;; make floor
;  ask patches with [ pycor = min-pycor ]
;    [ set pcolor brown ]

  ask n-of 10 patches with [pxcor mod 5 = 0 and pycor mod 5 = 0
  and pxcor > min-pxcor and pxcor < max-pxcor
  and pycor > min-pycor and pycor < max-pycor
  ] [set pcolor blue]
  ask patches with [pcolor = blue]
  [
  build-boxes
  ]
  make-box
  ;update-variables
  ;set init-avg-speed avg-speed
  ;set init-avg-energy avg-energy
  reset-ticks
end 

to update-variables
;  set medium count particles with [color = green]
;  set slow   count particles with [color = blue]
;  set fast   count particles with [color = red]
;  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 

to go
  if random 10 < 1 [pathway]
  make-box
  if(count particles < 10) [
    make-particles]
  ask particles [ bounce ]
  ask particles [ move ]
  ask particles
  [ if collide? [check-for-collision] ]

  ifelse (trace?)
    [ ask particle 0 [ pen-down ] ]
    [ ask particle 0 [ pen-up ] ]

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

  ask flashes with [ticks - birthday > 1]
    [ die ]
  display
end 

to calculate-tick-delta
  ;; tick-delta is calculated in such way that even the fastest
  ;; particle will jump at most 1 patch length in a tick. As
  ;; particles jump (speed * tick-delta) at every tick, making
  ;; tick length the inverse of the speed of the fastest particle
  ;; (1/max speed) assures that. Having each particle advance at most
  ;; one patch-length is necessary for it not to "jump over" a wall
  ;; or another particle.
  ifelse any? particles with [speed > 0]
    [ set tick-delta min list (1 / (ceiling max [speed] of particles)) max-tick-delta ]
    [ set tick-delta max-tick-delta ]
end 

to bounce  ;; particle procedure  ;; get the coordinates of the patch we will be on if we go forward 1
  let new-patch patch-ahead 1
  let new-px [pxcor] of new-patch
  let new-py [pycor] of new-patch
  if [pycor] of new-patch <= min-pycor
  [ask new-patch
    [ sprout-flashes 1 [
      set color  brown; white
      set birthday ticks
      set heading 0
    ]
  ]
  ]
  if new-py <= min-pycor[
    ifelse ([who] of self = 0)
    [
     pen-up
     setxy 0  max-pycor - 2; random-xcor  max-pycor - 2
     set color (random-float 4.9) + 5
     set heading (random 180) + 90
     stop
    ]
    [die]
  ]

  ;; if we're not about to hit a wall, we don't need to do any further checks
  if [pcolor] of new-patch != brown and [pcolor] of new-patch != blue
    [ stop ]
  if [pcolor] of new-patch = blue
    [
      set heading heading + 180
      boltzmann
      stop
  ]

  ;; if hitting the top or bottom, reflect heading around y axis
  if ([pcolor] of new-patch = brown and [d] of new-patch = 1)
    [ set heading (180 - heading)
      boltzmann
      stop
  ]
  if ([pcolor] of new-patch = brown and [d] of new-patch = 0)
  [ set heading (- heading)]

;  ask patches with [[pycor] of new-patch <= min-pycor]
;  [ sprout-flashes 1 [
;      set color white
;      set birthday ticks
;      set heading 0
;    ]
;  ]
end 

to move  ;; particle procedure
  ;; In other GasLab models, we use "jump speed * tick-delta" to move the
  ;; turtle the right distance along its current heading.  In this
  ;; model, though, the particles are affected by gravity as well, so we
  ;; need to offset the turtle vertically by an additional amount.  The
  ;; easiest way to do this is to use "setxy" instead of "jump".
  ;; Trigonometry tells us that "jump speed * tick-delta" is equivalent to:
  ;;   setxy (xcor + dx * speed * tick-delta)
  ;;         (ycor + dy * speed * tick-delta)
  ;; so to take gravity into account we just need to alter ycor
  ;; by an additional amount given by the classical physics equation:
  ;;   y(t) = 0.5*a*t^2 + v*t + y(t-1)
  ;; but taking tick-delta into account, since tick-delta is a multiplier of t.
  setxy (xcor + dx * speed * tick-delta)
        (ycor + dy * speed * tick-delta - gravity-acceleration * (0.5 * tick-delta * tick-delta))
  factor-gravity
end 

to build-boxes
  set i  1
    while [[pcolor] of patch-at 0 i != blue]
  [
      ask patch-at 0 i [set pcolor brown
    set d 0]
      set i i + 1
  ]
   set i  1
        while [[pcolor] of patch-at i 0 != blue]
  [
      ask patch-at i 0 [set pcolor brown
    set d 1]
      set i i + 1
  ]
end 

to pathway

  ask up-to-n-of random 2 patches with [pcolor = blue] [
    set i random 40 + 1
    set  z random 4
 ;   show z
;    if z = 1 or z = 3 [if random 2 < 2 [stop]] para que hayan más agujeros verticales
    while [([pcolor] of patch-at-heading-and-distance (z * 90) i) = brown
    and i < random 200][
      ask patch-at-heading-and-distance (z * 90) i [set pcolor black]
      set i i + 1
 ;     show i
    ]
;    ask self [set pcolor blue]
  ]

    ask up-to-n-of random 10 patches with [pcolor = blue] [
    set i random 40 + 1
    set z random 4
    while [([pcolor] of patch-at-heading-and-distance (z * 90) i) = black
    and i < random 200][
      ask patch-at-heading-and-distance (z * 90) i [set pcolor brown]
      set i i + 1
;      show i
    ]
;    ask self [set pcolor blue]
    ]
end 

to boltzmann
      set yboltz -1 * (ln (random-float 1)) * R * temp; from the potencial energy that should have a boltzmann distribution;
    ;  ifelse state = 1 [set yboltz  yboltz + floor2][set yboltz  yboltz + floor1]
      set speed 1 * (2 * gravity-acceleration * yboltz) ^ 0.5   ; the speed is calculated from the free fall formula.
      ; v^2 = 2gh is the speed of an object falling from an heigth h with a gravity g
      ;The 1.15 factor is a correction factor for ???????????????
  ;    set heading (180 - heading)
end 

to factor-gravity  ;; turtle procedure
  let vx (dx * speed)
  let vy (dy * speed) - (gravity-acceleration * tick-delta)
  set speed sqrt ((vy ^ 2) + (vx ^ 2))
;  recolor
  set heading atan vx vy
end 

to check-for-collision  ;; particle procedure
  ;; Here we impose a rule that collisions only take place when there
  ;; are exactly two particles per patch.  We do this because when the
  ;; student introduces new particles from the side, we want them to
  ;; form a uniform wavefront.
  ;;
  ;; Why do we want a uniform wavefront?  Because it is actually more
  ;; realistic.  (And also because the curriculum uses the uniform
  ;; wavefront to help teach the relationship between particle collisions,
  ;; wall hits, and pressure.)
  ;;
  ;; Why is it realistic to assume a uniform wavefront?  Because in reality,
  ;; whether a collision takes place would depend on the actual headings
  ;; of the particles, not merely on their proximity.  Since the particles
  ;; in the wavefront have identical speeds and near-identical headings,
  ;; in reality they would not collide.  So even though the two-particles
  ;; rule is not itself realistic, it produces a realistic result.  Also,
  ;; unless the number of particles is extremely large, it is very rare
  ;; for three or more particles to land on the same patch (for example,
  ;; with 400 particles it happens less than 1% of the time).  So imposing
  ;; this additional rule should have only a negligible effect on the
  ;; aggregate behavior of the system.
  ;;
  ;; Why does this rule produce a uniform wavefront?  The particles all
  ;; start out on the same patch, which means that without the only-two
  ;; rule, they would all start colliding with each other immediately,
  ;; resulting in much random variation of speeds and headings.  With
  ;; the only-two rule, they are prevented from colliding with each other
  ;; until they have spread out a lot.  (And in fact, if you observe
  ;; the wavefront closely, you will see that it is not completely smooth,
  ;; because some collisions eventually do start occurring when it thins out while fanning.)

  if count other particles-here = 1
  [
    ;; the following conditions are imposed on collision candidates:
    ;;   1. they must have a lower who number than my own, because collision
    ;;      code is asymmetrical: it must always happen from the point of view
    ;;      of just one particle.
    ;;   2. they must not be the same particle that we last collided with on
    ;;      this patch, so that we have a chance to leave the patch after we've
    ;;      collided with someone.
    let candidate one-of other particles-here with
      [who < [who] of myself and myself != last-collision]
    ;; we also only collide if one of us has non-zero speed. It's useless
    ;; (and incorrect, actually) for two particles with zero speed to collide.
    if (candidate != nobody) and (speed > 0 or [speed] of candidate > 0)
    [
      collide-with candidate
      set last-collision candidate
      ask candidate [ set last-collision myself ]
    ]
  ]
end 

;; implements a collision with another particle.
;;
;; THIS IS THE HEART OF THE PARTICLE SIMULATION, AND YOU ARE STRONGLY ADVISED
;; NOT TO CHANGE IT UNLESS YOU REALLY UNDERSTAND WHAT YOU'RE DOING!
;;
;; The two particles colliding are self and other-particle, and while the
;; collision is performed from the point of view of self, both particles are
;; modified to reflect its effects. This is somewhat complicated, so I'll
;; give a general outline here:
;;   1. Do initial setup, and determine the heading between particle centers
;;      (call it theta).
;;   2. Convert the representation of the velocity of each particle from
;;      speed/heading to a theta-based vector whose first component is the
;;      particle's speed along theta, and whose second component is the speed
;;      perpendicular to theta.
;;   3. Modify the velocity vectors to reflect the effects of the collision.
;;      This involves:
;;        a. computing the velocity of the center of mass of the whole system
;;           along direction theta
;;        b. updating the along-theta components of the two velocity vectors.
;;   4. Convert from the theta-based vector representation of velocity back to
;;      the usual speed/heading representation for each particle.
;;   5. Perform final cleanup and update derived quantities.

to collide-with [ other-particle ] ;; particle procedure
  ;;; PHASE 1: initial setup

  ;; for convenience, grab some quantities from other-particle
  let mass2 [mass] of other-particle
  let speed2 [speed] of other-particle
  let heading2 [heading] of other-particle

  ;; since particles are modeled as zero-size points, theta isn't meaningfully
  ;; defined. we can assign it randomly without affecting the model's outcome.
  let theta (random-float 360)



  ;;; PHASE 2: convert velocities to theta-based vector representation

  ;; now convert my velocity from speed/heading representation to components
  ;; along theta and perpendicular to theta
  let v1t (speed * cos (theta - heading))
  let v1l (speed * sin (theta - heading))

  ;; do the same for other-particle
  let v2t (speed2 * cos (theta - heading2))
  let v2l (speed2 * sin (theta - heading2))



  ;;; PHASE 3: manipulate vectors to implement collision

  ;; compute the velocity of the system's center of mass along theta
  let vcm (((mass * v1t) + (mass2 * v2t)) / (mass + mass2) )

  ;; now compute the new velocity for each particle along direction theta.
  ;; velocity perpendicular to theta is unaffected by a collision along theta,
  ;; so the next two lines actually implement the collision itself, in the
  ;; sense that the effects of the collision are exactly the following changes
  ;; in particle velocity.
  set v1t (2 * vcm - v1t)
  set v2t (2 * vcm - v2t)



  ;;; PHASE 4: convert back to normal speed/heading

  ;; now convert my velocity vector into my new speed and heading
  set speed sqrt ((v1t ^ 2) + (v1l ^ 2))
  set energy (0.5 * mass * speed * speed)
  ;; if the magnitude of the velocity vector is 0, atan is undefined. but
  ;; speed will be 0, so heading is irrelevant anyway. therefore, in that
  ;; case we'll just leave it unmodified.
  if v1l != 0 or v1t != 0
    [ set heading (theta - (atan v1l v1t)) ]

  ;; and do the same for other-particle
  ask other-particle [
    set speed sqrt ((v2t ^ 2) + (v2l ^ 2))
    set energy (0.5 * mass * (speed ^ 2))
    if v2l != 0 or v2t != 0
      [ set heading (theta - (atan v2l v2t)) ]
  ]


  ;; PHASE 5: final updates

  ;; now recolor, since color is based on quantities that may have changed
;  recolor
;  ask other-particle
;    [ recolor ]
end 

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

;;;
;;; drawing procedures
;;;

;; draws the box

to make-box
  ask patches with [ pxcor = min-pxcor or pxcor = max-pxcor or pycor = min-pycor or pycor = max-pycor]
    [ set pcolor brown]
end 

;; creates initial particles

to make-particles
  create-particles number-of-particles
  [
    setup-particle
    random-position
;    recolor
  ]
  calculate-tick-delta
end 

to setup-particle  ;; particle procedure
  set speed init-particle-speed
  set mass particle-mass
  set energy (0.5 * mass * speed * speed)
  set last-collision nobody
  set size 2
  set color red
  if ([who] of self = 0)
  [set color white set pen-size 2]
end 

;; place particle at random location inside the box.

to random-position ;; particle procedure
  setxy 0  max-pycor - 2; random-float (world-height - 3) + 1
  set heading random-float 360
end 

to draw-energy-graph
  let n min-pycor
  repeat (max-pycor / 2)
  [ let x particles with [ (pycor >= n) and (pycor < (n + 4)) ]
    ifelse count x = 0
     [ plot 0 ]
     [ let temperature mean [ energy ] of x
       plot temperature
       let list-pos (n + max-pycor) / 4
       set aggregate-list (replace-item list-pos aggregate-list ((item list-pos aggregate-list) + temperature))
     ]
    set n n + 4
  ]
end 

to draw-aggregate-graph [lst]
  foreach lst plot
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 

to-report last-n [n the-list]
  ifelse n >= length the-list
    [ report the-list ]
    [ report last-n n butfirst the-list ]
end 


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