particles in a box Master Equation

globals
[
  box-edge1
  box-edge2
  leftbox-edge
  rightbox
  floor1
  floor2
  entropy1
 ; entropy2
 ; temperature
   prev-temp
  R
  yboltz
  enthalpy1
  ; enthalpy2
 ; scale
  Eafloor
 initialstate ; particles in initially in state 1
  ΔG
  ΔG0
  ΔG0e
  -ΔG0e
  ΔS
  ΔH
  Kin1
  Kin2
  K
  tracer
  a
    P1 ; probabilidad estado 1
  P2 ; probabilidad estado 2
  alfa ; probabilidad de saltar a estado 2 estando en 1
  beta ; probabilidad de saltar a estado 1 estando en 2
  rr ; factor probabilidades
  pinicial; probabilidad inicial de estar en estado 1
  averageEnergy
  time
]

breed [ particles particle ]
breed [ flashes flash ]

flashes-own [birthday]

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

to setup
  clear-all
  set-default-shape particles "circle"
  set-default-shape flashes "plane"
  set R 1.9858775 / 1000 ;constante de gases en kcal/mol/K
 ; set entropy1 0.3
   set initialState 1 - avance
  set prev-temp temperature
  set entropy1 1 - entropy2
  set enthalpy1 0
  ifelse enthalpy1 <= enthalpy2
   [set floor1 0
     set floor2 enthalpy2 * scale] ; scale 1kcal= scale y coordenate
   [set floor2 0
     set floor1 -1 * scale * enthalpy2] ; scale 1kcal=scale y coordinate
  set Eafloor max (list floor1 floor2) + Ea * scale
 ; set initialstate 0.1
  set box-edge2 (max-pxcor)
  set box-edge1 min-pxcor + (entropy1 * max-pxcor * 2)
  set ΔH enthalpy2 - enthalpy1
  set ΔS  R * ln (entropy2 / entropy1) ;calculo entropía
  set ΔG0 ΔH - temperature * ΔS ;calculo Δ G cero teórica
  set pinicial initialState
  make-box
  set initialstate 1 - avance
  make-particles
  set prev-temp 0
  update-variables
  reset-ticks
  set a 1
    set a 1
end 

to go
  ask particles [ bounce ]
    update-variables
    update-plots
    display
    set a a + 1
  if single = true [stop]
;  if a > 1 [stop]
end 

to update-variables
  set ΔG0 ΔH - temperature * ΔS ;calculo Δ G cero teórica
    if temperature != prev-temp [
    masterEquation
  set prev-temp temperature
  set time 0
  ]
  ; ecuación maestra
  set P1 1 / (alfa + beta) * (beta + rr * exp(-(alfa + beta) * time ))
  set P2 1 / (alfa + beta) * (alfa - rr * exp(-(alfa + beta) * time ))
  set averageEnergy R * temperature * ln (P2 / (P1 + 1E-10)) + ΔG0

  set Kin1 count turtles with [state = 1]
  set Kin2 count turtles with [state = 2]
  set K count turtles with [state = 2] / (count turtles with [state = 1] + 0.0000000000000001)
  set ΔG0e (- R * temperature * ln (K + 1E-10))
  set -ΔG0e -1 * ΔG0e
  set ΔG R * temperature * ln (K + 1E-10) + ΔG0
end 

to masterEquation
   ; establece cantiadades para la ecuación maestra
  ifelse enthalpy1 <= enthalpy2
  ; si deltaH es positivo
  [set alfa entropy2 * exp(- (Ea + ΔH) / R / temperature);si deltaH es positivo, para ir a rojo tiene que saltar
  set beta entropy1 * exp(- (Ea) / R / temperature)
  ]
  ;si deltaH es negativo
  [set alfa entropy2 * exp(- (Ea) / R / temperature)
  set beta entropy1 * exp(- (Ea - ΔH) / R / temperature); si deltaH es negativo para ir a verde tiene que saltar
  ]
  set Kin1 count turtles with [state = 1]
  set Kin2 count turtles with [state = 2]
  set pinicial Kin1 / (Kin1 + Kin2)
  set rr alfa * pinicial - beta * (1 - pinicial)
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 green]
  [set state 2
    set color red]
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 white ]
     ]
     [
       ask patches with [ (pxcor < box-edge1) and (pycor < floor1) ]
     [ set pcolor white ]
     ]

       ; limite superior
    ask patches with [ pycor > max-pycor - 5 ]
    [ set pcolor blue ]

    ; limite inferior
    ask patches with [ (pycor < floor1 + 5 and pycor > floor1 ) and (pxcor <= box-edge1) ]
    [ set pcolor yellow ]

    ask patches with [ (pycor < floor2 + 5 and pycor > floor2) and (pxcor >= box-edge1) ]
    [ set pcolor yellow ]
    ; limite izquierdo
    ask patches with [ (pxcor < min-pxcor + 5) and (pycor >= floor1) ]
    [ set pcolor yellow ]

    ; limite derecho
    ask patches with [ (pxcor > max-pxcor - 5) and (pycor >= floor2) ]
    [ set pcolor yellow ]

    ; limite central
    ask patches with [ (pxcor < box-edge1 + 3 and pxcor > box-edge1 - 3) and (pycor <= Eafloor) ]
    [ set pcolor yellow ]
end 
;;;


;; creates initial particles

to make-particles
  create-particles number-of-particles
  [

    set size 6
    ifelse random-float 1 < initialstate
    [
      set xcor min-pxcor + random (box-edge1 - min-pxcor)
      set yboltz floor1 - scale * ln (random-float 1) * R * temperature
      ifelse yboltz > max-pycor
       [set ycor  max-pycor]
       [set ycor yboltz]
    ]

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


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

  ]
   set K count turtles with [state = 2] / (count turtles with [state = 1] + 0.00000000000000000000001)
  set ΔG0e (- R * temperature * ln (K + 0.0000000000001))
  set ΔG R * temperature * ln (K + 0.0000000000001) + ΔG0
end 




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