MolecularThermodynamicsInBoxes

globals
[
  tick-delta                          ;; how much we advance the tick counter this time through
  max-tick-delta                      ;; the largest tick-delta is allowed to be
  avg-speed-init avg-energy-init      ;; initial averages
  avg-speed avg-energy                ;; current averages
  fast medium slow lost-particles     ;; current counts
  percent-medium percent-slow         ;; percentage of current counts
  percent-fast percent-lost-particles ;; percentage of current counts
  box-edge1
  box-edge2
  leftbox-edge
  rightbox
  floor1
  floor2
  entropy2
  R
  enthalpy1
  Eafloor
  deltaG
  deltaG0
  deltaG0e
  deltaS
  deltaH
  K
  percent1
  yboltz
  tracer
  punta

  min-xcor
  max-xcor
  min-ycor
  max-ycor
  p/r

  ; entropy1
  ; enthalpy2
  ; scale
  ; initialstate ; particles in initially in state 1
  ; temperature
]

breed [invisibles invisible]
breed [ particles particle ]


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

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

to setup
  clear-all
  set-default-shape particles "circle"
  set max-tick-delta 0.1073
  set R 1.9858775 / 1000 ;constante de gases en kcal/mol/K
  set entropy2 1 - entropy1 ; entropy expresed as probability
  set enthalpy1 0
  ifelse enthalpy1 <= enthalpy2
   [set floor1 0
     set floor2 (enthalpy2 * scale)] ; scale 1kcal= scale y coordenate
   [set floor2 0
      set floor1 (-1 * enthalpy2 * scale)] ; scale 1kcal=scale y coordinate
  set Eafloor max (list floor1 floor2) + Ea * scale
  set min-xcor min-pxcor
  set max-xcor max-pxcor
  set min-ycor min-pycor
  set max-ycor max-pycor
  set box-edge2 max-xcor
  set box-edge1 min-xcor + (entropy1 * (max-xcor - min-xcor))
  set punta patch box-edge1 Eafloor
  set deltaH enthalpy2 - enthalpy1
  set deltaS  R * ln (entropy2 / entropy1) ;calculo entropía
  set deltaG0 deltaH - temperature * deltaS ;calculo delta G cero teórica
  ;; make floor
  make-box
  make-separation
  make-particles
  make-invisibles
  update-variables
  set avg-speed-init avg-speed
  set avg-energy-init avg-energy
  reset-ticks
  set tracer one-of particles
end 

to go



  if not any? particles [stop]  ;; particles can die when they float too high
  ask particles
  [
    if collide? [check-for-collision]
    change-position
    ]
  if tracer = nobody
  [ set tracer one-of particles ]
  ifelse trace?
  [ ask tracer [ pen-down ] ]
  [ ask tracer [ pen-up ] ]
  tick-advance tick-delta
  if floor ticks > floor (ticks - tick-delta)
  [
    ask invisibles [ bounce ]
    update-variables
    update-plots
  ]
  calculate-tick-delta
  display
end 

to update-variables
  set K count invisibles with [state = 2] / (count invisibles with [state = 1] + 1E-20)
  set deltaG0e (- R * temperature * ln (K + 1E-20))
  set deltaG0 deltaH - temperature * deltaS ;calculo delta G cero teórica
  set deltaG R * temperature * ln (K + 1E-20) + deltaG0
  set avg-speed  mean [speed] of particles
  set avg-energy  mean [energy] of particles
  set percent1 count invisibles with [state = 1] / (number-of-particles) * 100
  set p/r count invisibles with [state = 2] / (count invisibles with [state = 1] + 0.000001)
end 

to bounce  ;; particle procedure

 ifelse state = 1
     [
     set yboltz floor1 - scale * ln (random-float 1) * R * temperature
     if yboltz >= Eafloor
     [set xcor min-pxcor + random (max-pxcor - min-pxcor)]
     ]
     [
     set yboltz floor2 - scale * ln (random-float 1) * R * temperature
     if yboltz >= Eafloor
     [set xcor min-pxcor + random (max-pxcor - min-pxcor)]
     ]

  ifelse yboltz > max-pycor
       [set ycor  max-pycor]
       [set ycor yboltz]

  ifelse (xcor >= min-pxcor and xcor < box-edge1)
   [set state 1
    set color blue]
  [set state 2
    set color yellow]
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 change-position
   let oldx xcor
   let oldy ycor
   move
   if xcor = oldx and ycor = oldy
   [print "no se mueve"]
end 

to move  ;; particle procedure

  let new-patch patch-ahead 1
  ;; if we're  about to hit a wall, bounce and then move
  ;if new-patch = nobody or [pcolor] of new-patch != 0
 ;[bounce]

  let xdiff dx * speed * tick-delta
  let ydiff dy * speed * tick-delta - gravity-acceleration * (0.5 * (tick-delta ^ 2))
  let new-x xcor + xdiff
  let new-y ycor + ydiff
  if new-x <= min-xcor
  [
    set new-x min-xcor
    set heading (- heading)
        boltzman-n
    ]
  if new-x >= max-xcor
  [
    set new-x max-xcor
    set heading (- heading)
        boltzman-n
    ]
    if state = 1 and ycor <= Eafloor and new-x > box-edge1
  [
    set new-x box-edge1 - 1
    set heading (- heading)
        boltzman-n
    ]
  if state = 2 and ycor <= Eafloor and new-x <= box-edge1
  [
    set new-x box-edge1 + 1
    set heading (- heading)
    boltzman-n

    ]
  if new-y >= max-ycor
  [set new-y max-ycor
    set heading (180 - heading)
        boltzman-n
    ]
  if state = 1 and new-y <= floor1
  [
    set new-y floor1
    boltzman-n
    set heading (180 - heading)
    ]
  if state = 2 and new-y <= floor2
  [
    set new-y floor2
    boltzman-n
    set heading (180 - heading)
    ]
  setxy new-x  new-y
  ifelse (xcor >= min-pxcor and xcor < box-edge1)
  [set state 1
   set color green]
  [set state 2
    set color red]
  factor-gravity
end 

to boltzman-n


      set yboltz -1 * scale * (ln (random-float 1)) * R * temperature; 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))
  set energy (0.5 * mass * (speed ^ 2))
  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.

  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
  ;; local copies of other-particle's relevant quantities:
  ;; mass2 speed2 heading2
  ;;
  ;; quantities used in the collision itself
  ;; theta   ;; heading of vector from my center to the center of other-particle.
  ;; v1t     ;; velocity of self along direction theta
  ;; v1l     ;; velocity of self perpendicular to theta
  ;; v2t v2l ;; velocity of other-particle, represented in the same way
  ;; vcm     ;; velocity of the center of mass of the colliding particles,
             ;;   along direction theta

  ;;; 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 ^ 2))
  ;; 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
end 



;;;
;;; drawing procedures

to make-box

     ; white the part of the  box that is inactive
     ifelse floor1 <= floor2
     [
       ask patches with [ (pxcor > box-edge1) and (pycor < floor2) ]
     [ set pcolor 8 ]
     ]
     [
       ask patches with [ (pxcor < box-edge1) and (pycor < floor1) ]
     [ set pcolor 8 ]
     ]
end 
;;;

to make-separation
  create-particles 1
  [
    setxy box-edge1 0
    set pen-size world-width / 50
    set color 45
    pen-down
    move-to patch box-edge1 Eafloor
    die
  ]
end 

;; creates initial particles

to make-particles
  create-particles number-of-particles
  [
    setup-particle
    ifelse random-float 1 < initialstate
    [
      set xcor min-xcor + random (box-edge1 - min-xcor - 1)
      set yboltz floor1  - scale * ln (random-float 1) * R * temperature
      ifelse yboltz > max-ycor
       [set ycor  max-ycor]
       [set ycor yboltz]
    ]

    [
      set xcor box-edge1 + 1 + random (max-xcor - box-edge1 - 1)
      set yboltz floor2 - scale * ln (random-float 1) * R * temperature
      ifelse yboltz > max-ycor
       [set ycor  max-ycor]
       [set ycor yboltz]
    ]

   ifelse (xcor < box-edge1)
  [set state 1
    set color green]
  [set state 2
    set color red]
    ;;random-position

  ]
 ; set K count particles with [state = 2] / (count particles with [state = 1] + 1E-20)
 ;set deltaG0e (- R * temperature * ln (K + 1E-20))
 ;set deltaG R * temperature * ln (K + 1E-20) + deltaG0

  calculate-tick-delta
end 

to make-invisibles
  create-invisibles number-of-particles
  [
      set size 5
      set hidden? true
    ifelse random-float 1 < initialstate
    [
      set xcor min-xcor + random (box-edge1 - min-xcor - 1)
      set yboltz floor1  - scale * ln (random-float 1) * R * temperature
      ifelse yboltz > max-ycor
       [set ycor  max-ycor]
       [set ycor yboltz]
    ]

    [
      set xcor box-edge1 + 1 + random (max-xcor - box-edge1 - 1)
      set yboltz floor2 - scale * ln (random-float 1) * R * temperature
      ifelse yboltz > max-ycor
       [set ycor  max-ycor]
       [set ycor yboltz]
    ]

   ifelse (xcor < box-edge1)
  [set state 1
    set color blue]
  [set state 2
    set color yellow]
    ;;random-position

  ]
  set K count invisibles with [state = 2] / (count invisibles with [state = 1] + 1E-20)
  set deltaG0e (- R * temperature * ln (K + 1E-20))
  set deltaG R * temperature * ln (K + 1E-20) + deltaG0
end 

to setup-particle  ;; particle procedure
  set size 10
  set speed init-particle-speed
  set heading 90 - random 180
  set mass particle-mass
  set energy (0.5 * mass * (speed ^ 2))
  set last-collision nobody
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 1997 Uri Wilensky.
; See Info tab for full copyright and license.