Child of Osmotic Pressure

globals [
  left-side                                 ;; left side of the membrane
  right-side                                ;; right side of the membrane
  split                                     ;; location of membrane
  pre-split                                 ;; a value for calculating the change in membrane location
  equilibrium                               ;; the difference in number of particles moving each direction across membrane
  tick-delta                                ;; how much we advance the tick counter this time through
  collisions        ;; list used to keep track of future collisions
  particle1         ;; first particle currently colliding
  particle2         ;; second particle currently colliding
]

breed [particles particle]
breed [membranes membrane]

turtles-own [
  speed
  mass
  energy
  old-pos
  new-pos
  part-type
]

to setup
  ;; (for this model to work with NetLogo's new plotting features,
  ;; __clear-all-and-reset-ticks should be replaced with clear-all at
  ;; the beginning of your setup procedure and reset-ticks at the end
  ;; of the procedure.)
  clear-all
  reset-ticks
  set-default-shape particles "circle"
  create-container
  spawn-particles
  set split 0
  calculate-wall
  set particle1 nobody
  set particle2 nobody
  set collisions []
  ask particles [ check-for-wall-collision ]
  ask particles [ check-for-particle-collision ]
  my-update-plots
end 

to create-container
  ask patches with [pycor = max-pycor or pycor = min-pycor] [              ;; sets the walls of the container
    set pcolor blue
  ]
  ask patches with [pxcor = max-pxcor or pxcor = min-pxcor] [
    set pcolor blue
  ]
  set left-side patches with [pxcor < 0]                                   ;; defines the left and right sides of the container equally
  set right-side patches with [pxcor > 0]
  set equilibrium 0                                                        ;; sets equilibrium to starting value (the middle of the container)
  set split 0                                                              ;; sets the starting location of the membrane
  ask patches with [pxcor = 0 and pycor != max-pycor and pycor != min-pycor and abs pycor mod 2 != 0] [
    sprout-membranes 1 [                                                   ;; create the membrane
      set shape "square"
      set color red
      set size 1.3
      set heading 90
    ]
  ]
end 

to setup-particles                                                      ;; set variable values for each new particle
  set speed 0.5
  if part-type = "solvent" [
    set size 1
  ]
  if part-type = "solute" [
    set size 2
  ]
  set mass (size * size)
  set energy (0.5 * mass * speed * speed)
end 

to-report overlapping?  ;; particle procedure
  ;; here, we use IN-RADIUS just for improved speed; the real testing
  ;; is done by DISTANCE
  report any? other turtles in-radius ((size + 2) / 2)
                              with [distance myself < (size + [size] of myself) / 2]
end 

to-report wall-lapping?  ;; particle procedure
  report any? patches in-radius ((size + 2) / 2) with [pcolor = blue] or
  any? membranes in-radius ((size + 2) / 2)
end 

to start-loc
  move-to one-of patches with [pcolor = black]                        ;; randomly distribute them throughout the world
    while [overlapping?] [
      move-to one-of patches with [pcolor = black]
    ]
    while [wall-lapping?] [
      move-to one-of patches with [pcolor = black]
    ]
end 

to start-right-loc
  move-to one-of right-side with [pcolor = black]                        ;; randomly distribute them throughout the world
    while [overlapping?] [
      move-to one-of right-side with [pcolor = black]
    ]
    while [wall-lapping?] [
      move-to one-of right-side with [pcolor = black]
    ]
end 

to start-left-loc
  move-to one-of left-side with [pcolor = black]                        ;; randomly distribute them throughout the world
    while [overlapping?] [
      move-to one-of left-side with [pcolor = black]
    ]
    while [wall-lapping?] [
      move-to one-of left-side with [pcolor = black]
    ]
end 

to spawn-particles
  if solute = "Sugar" [                                                ;; checks to see the identity of the solute based on the chooser
    output-type "C12H22O11"                                            ;; write the chemical formula in the output area
    create-particles solute-left [                                     ;; create solute particles based on slider value
      set color white
      set shape "circle"
      set part-type "solute"                                        ;; define these as solute particles
      setup-particles
      start-left-loc
    ]
    create-particles solute-right [                                    ;; create solute particles based on slider value
      set color white
      set shape "circle"
      set part-type "solute"
      setup-particles
      start-right-loc
    ]
  ]

  if solute = "Sodium Chloride" [
    output-type "NaCl"
    create-particles solute-left [
      set color white
      set shape "circle"
      set part-type "solute"
      setup-particles
      start-left-loc
      ask patch-here [
        sprout-particles 1 [                                   ;; solutes that dissociate in water split into their respective ions
          set color white                                      ;; since NaCl splits into two ions, one new turtle is sprouted for each solute molecule
          set shape "circle"
          set part-type "solute"
          setup-particles
          fd 0.2
          ]
      ]
    ]
    create-particles solute-right [
      set color white
      set shape "circle"
      set part-type "solute"
      setup-particles
      start-right-loc
      ask patch-here [
        sprout-particles 1 [
          set color white
          set shape "circle"
          set part-type "solute"
          setup-particles
          fd 0.2
          ]
      ]
    ]
  ]

  if solute = "Magnesium Chloride" [
    output-type "MgCl2"
    create-particles solute-left [
      set color white
      set shape "circle"
      set part-type "solute"
      setup-particles
      start-left-loc
      ask patch-here [
        sprout-particles 2 [                                        ;; splits into 3 ions..so 2 new turtles are sprouted
          set color white
          set shape "circle"
          set part-type "solute"
          setup-particles
          fd 0.2
          ]
      ]
    ]
    create-particles solute-right [
      set color white
      set shape "circle"
      set part-type "solute"
      setup-particles
      start-right-loc
      ask patch-here [
        sprout-particles 2 [
          set color white
          set shape "circle"
          set part-type "solute"
          setup-particles
          fd 0.2
          ]
      ]
    ]
  ]

  if solute = "Aluminum Chloride" [
    output-type "AlCl3"
    create-particles solute-left [
      set color white
      set shape "circle"
      set part-type "solute"
      setup-particles
      start-left-loc
      ask patch-here [
        sprout-particles 3 [                                     ;; splits into 4 ions...so 3 new turtles are sprouted
          set color white
          set shape "circle"
          set part-type "solute"
          setup-particles
          fd 0.2
          ]
      ]
    ]
    create-particles solute-right [
      set color white
      set shape "circle"
      set part-type "solute"
      setup-particles
      start-right-loc
      ask patch-here [
        sprout-particles 3 [
          set color white
          set shape "circle"
          set part-type "solute"
          setup-particles
          fd 0.2
          ]
      ]
    ]
  ]
  create-particles (100 - count particles with [part-type = "solute" and xcor < split]) [                                               ;; creates 1000 water molecules
    set color blue + 2
    set part-type "solvent"                                             ;; define these as solvents
    setup-particles
    start-left-loc
  ]
  create-particles (100 - count particles with [part-type = "solute" and xcor > split]) [                                               ;; creates 1000 water molecules
    set color blue + 2
    set part-type "solvent"                                             ;; define these as solvents
    setup-particles
    start-right-loc
  ]
end 

to particle-jump
  ask particles with [xcor >= pre-split and xcor <= split and part-type = "solutes"] [     ;; check for solutes in the way of a membrane jump
    set xcor (split - pre-split) + 1
;    recalculate-particles-that-just-collided
  ]
  ask particles with [xcor <= pre-split and xcor >= split and part-type = "solutes"] [
    set xcor (pre-split - split) + 1
;    recalculate-particles-that-just-collided
  ]
end 

to calculate-wall
  set pre-split split                                                  ;; save location of the membrane before it moves
  let nudge ((equilibrium * -1) / 10)                                  ;; calculates the amount to move the membrane
  set split split + nudge                                              ;; set split to be the new location of the membrane
  particle-jump                                                        ;; move solutes in the way of the membrane jump
  ask membranes [
    set xcor xcor + nudge                                              ;; move membrane turtles according to nudge value
  ]
  
  set left-side patches with [pxcor < split]                           ;; redefine right and left sides
  set right-side patches with [pxcor > split]
  set equilibrium 0                               ;; reset equilibrium so we can keep track of what happens during the next tick
end 

to go
  if count turtles-on right-side = 0 or count turtles-on left-side = 0 [                    ;; stops model if membrane reaches either edge
    user-message "The membrane has burst! Make sure you have some solute on both sides!"
    stop
  ]
  ask particles [ set old-pos xcor ]
  choose-next-collision

  ask particles with [ part-type = "solute" ] [      ;; bouncing for solutes hitting membrane
    if any? membranes in-cone 2 180 [
      set heading (- heading)
    ]
  ]

  ask particles [ jump speed * tick-delta ]
  perform-next-collision

  ask particles with [ part-type = "solvent" ] [
    if old-pos < split and xcor >= split [                      ;; if a solvent moves from the left of the membrane to the right of the membrane
      set equilibrium equilibrium + 1                           ;; add one to equilibrium
    ]
    if old-pos > split and xcor <= split [                      ;; if a solvent moves from the right of the membrane to the left of the membrane
      set equilibrium equilibrium - 1                           ;; subract one from equilibrium
    ]
  ]

  calculate-wall                                                ;; recalculate the new location of the wall
  tick-advance tick-delta
  recalculate-particles-that-just-collided
  my-update-plots
end 

to recalculate-particles-that-just-collided
  ;; Since only collisions involving the particles that collided most recently could be affected,
  ;; we filter those out of collisions.  Then we recalculate all possible collisions for
  ;; the particles that collided last.  The ifelse statement is necessary because
  ;; particle2 can be either a particle or a string representing a wall.  If it is a
  ;; wall, we don't want to invalidate all collisions involving that wall (because the wall's
  ;; position wasn't affected, those collisions are still valid.
  ifelse is-turtle? particle2
    [
      set collisions filter [item 1 ? != particle1 and
                             item 2 ? != particle1 and
                             item 1 ? != particle2 and
                             item 2 ? != particle2]
                            collisions
      ask particle2 [ check-for-wall-collision ]
      ask particle2 [ check-for-particle-collision ]
    ]
    [
      set collisions filter [item 1 ? != particle1 and
                             item 2 ? != particle1]
                            collisions
    ]
  if particle1 != nobody [ ask particle1 [ check-for-wall-collision ] ]
  if particle1 != nobody [ ask particle1 [ check-for-particle-collision ] ]
  ;; Slight errors in floating point math can cause a collision that just
  ;; happened to be calculated as happening again a very tiny amount of
  ;; time into the future, so we remove any collisions that involves
  ;; the same two particles (or particle and wall) as last time.
  set collisions filter [item 1 ? != particle1 or
                         item 2 ? != particle2]
                        collisions
  ;; All done.
  set particle1 nobody
  set particle2 nobody
end 

;; check-for-particle-collision 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

to check-for-particle-collision
  let my-x xcor
  let my-y ycor
  let my-particle-size size
  let my-x-speed speed * sin heading
  let my-y-speed speed * cos heading
  ask other particles
  [
    let dpx (xcor - my-x)   ;; relative distance between particles in the x direction
    let dpy (ycor - my-y)    ;; relative distance between particles in the y direction
    let x-speed (speed * sin heading) ;; speed of other particle in the x direction
    let y-speed (speed * cos heading) ;; speed of other particle in the x direction
    let dvx (x-speed - my-x-speed) ;; relative speed difference between particles in x direction
    let dvy (y-speed - my-y-speed) ;; relative speed difference between particles in y 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) will 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 the two particles.
    ;;
    ;; The equation that represents the relationship is:
    ;;   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)) - (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)))  ;; vector product of the position times the velocity
    let v-squared  ((dvx * dvx) + (dvy * dvy)) ;; 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 a function of time (t) can be represented 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 test 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.

    if time-to-collision > 0
    [
      ;; time-to-collision is relative (ie, a collision will occur one second from now)
      ;; We need to store the absolute time (ie, a collision will occur at time 48.5 seconds.
      ;; So, we add clock to time-to-collision when we store it.
      ;; The entry we add is a three element list of the time to collision and the colliding pair.
      set collisions fput (list (time-to-collision + ticks) self myself)
                          collisions
    ]
  ]
end 


;; determines when a particle will hit any of the four walls or membrane

to check-for-wall-collision  ;; particle procedure
  ;; right & left walls
  let x-speed (speed * sin heading)
  if x-speed != 0
    [ ;; solve for how long it will take particle to reach right wall
      let right-interval (max-pxcor - 0.5 - xcor - size / 2) / x-speed
      if right-interval > 0
        [ assign-colliding-wall right-interval "right wall" ]
      ;; solve for time it will take particle to reach left wall
      let left-interval ((- max-pxcor) + 0.5 - xcor + size / 2) / x-speed
      if left-interval > 0
        [ assign-colliding-wall left-interval "left wall" ] ]

;  ;; membrane
;  if x-speed != 0 and part-type = "solute" and xcor < split
;    [ ;; solve for how long it will take particle to reach the membrane from the left
;      let left-interval2 (split - xcor - size / 2) / x-speed
;      if left-interval2 > 0
;        [ assign-colliding-wall left-interval2 "membrane-left" ] ]
;   if x-speed != 0 and part-type = "solute" and xcor < split
;    [ ;; solve for time it will take particle to reach the membrane from the right
;      let right-interval2 (split - xcor + size / 2) / x-speed
;      if right-interval2 > 0
;        [ assign-colliding-wall right-interval2 "membrane-right" ] ]

  ;; top & bottom walls
  let y-speed (speed * cos heading)
  if y-speed != 0
    [ ;; solve for time it will take particle to reach top wall
      let top-interval (max-pycor - 0.5 - ycor - size / 2) / y-speed
      if top-interval > 0
        [ assign-colliding-wall top-interval "top wall" ]
      ;; solve for time it will take particle to reach bottom wall
      let bottom-interval ((- max-pycor) + 0.5 - ycor + size / 2) / y-speed
      if bottom-interval > 0
        [ assign-colliding-wall bottom-interval "bottom wall" ] ]
end 

to assign-colliding-wall [time-to-collision wall]  ;; particle procedure
  ;; this procedure is used by the check-for-wall-collision procedure
  ;; to assemble the correct particle-wall pair
  ;; time-to-collision is relative (ie, a collision will occur one second from now)
  ;; We need to store the absolute time (ie, a collision will occur at time 48.5 seconds.
  ;; So, we add clock to time-to-collision when we store it.
  let colliding-pair (list (time-to-collision + ticks) self wall)
  set collisions fput colliding-pair collisions
end 

to choose-next-collision
  if collisions = [] [ stop ]
  ;; Sort the list of projected collisions between all the particles into an ordered list.
  ;; Take the smallest time-step from the list (which represents the next collision that will
  ;; happen in time).  Use this time step as the tick-delta for all the particles to move through
  let winner first collisions
  foreach collisions [ if first ? < first winner [ set winner ? ] ]
  ;; winner is now the collision that will occur next
  let dt item 0 winner
  ;; If the next collision is more than 1 in the future,
  ;; only advance the simulation one tick, for smoother animation.
  set tick-delta dt - ticks
  if tick-delta > 1
    [ set tick-delta 1
      set particle1 nobody
      set particle2 nobody
      stop ]
  set particle1 item 1 winner
  set particle2 item 2 winner
end 

to perform-next-collision
  ;; deal with 3 possible cases:
  ;; 1) no collision at all
  if particle1 = nobody [ stop ]
  ;; 2) particle meets wall
  if is-string? particle2
    [ if particle2 = "left wall" or particle2 = "right wall" ;or particle2 = "membrane-left" or particle2 = "membrane-right"
        [ ask particle1 [ set heading (- heading) ]
          stop ]
      if particle2 = "top wall" or particle2 = "bottom wall"
        [ ask particle1 [ set heading 180 - heading ]
          stop ] ]
  ;; 3) particle meets particle
  ask particle1 [ collide-with particle2 ]
end 

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
  ;; modified so that theta is heading toward other particle
  let theta towards other-particle

  ;;; 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 * v1t) + (v1l * v1l))
  ;; 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.
  set energy (0.5 * mass * speed ^ 2)

  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)) ]
  ]
end 

to my-update-plots
  set-current-plot "Water#"
  set-current-plot-pen "left"
  plot count particles with [pxcor < split and part-type = "solvent"]
  set-current-plot-pen "right"
  plot count particles with [pxcor > split and part-type = "solvent"]
end