Smog

; This model simulates the formation of photochemical smog over time. 

globals [
  earth-top      ;; y-coordinate that splits earth from atmosphere
  timeofday      ;; time of day in hours (0-24)
  sunlight       ;; level of sunlight (0-100)
  p-rate         ;; adjuster for emission rate of pollutants (0-1)
  NO2-rate       ;; rate of NO2-splitting reaction
  O3-rate        ;; rate of O3-formation reaction
  sunrise        ;; time at which sun rises
  sunset ]       ;; time at which sun sets

patches-own [
  patch-type ]   ;; 0 = earth, 1 = buildings, 2 = sky

breed [VOC]   ;; volatile organic chemicals
breed [NO2]   ;; nitrogen dioxide
breed [NO]    ;; nitrogen oxide
breed [O1]    ;; oxygen radicals
breed [O2]    ;; oxygen
breed [O3]    ;; ozone
breed [Ox]    ;; photochemical oxidants


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;; SETUP PROCEDURES ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to Setup
  clear-all
  set-default-shape VOC "dot"   ;; set default shapes for all turtle breeds
  set-default-shape NO2 "dot"
  set-default-shape O3 "dot"
  set-default-shape NO "dot"
  set-default-shape O1 "dot"
  set-default-shape O2 "dot"  
  set-default-shape Ox "dot"  
  setup-world
  set sunrise 420
  set sunset 1140
  set timeofday 0                ;; reset clock 
  reset-ticks                    ;; reset tick counter
end 

to setup-world
  set earth-top -12
  set sunlight sunlight-intensity
  ask patches with [             ;; create grass
    pycor < earth-top ] [
    set patch-type 0
    set pcolor 64 ]
  ask patches with [
    pycor >= earth-top ] [        ;; create sky
    set patch-type 2 
    set pcolor scale-color sky sunlight 0 100 ]
  ask patches with [             ;; create buildings
    pycor >= earth-top and
    pycor <= earth-top + 1 and
    pxcor >= -8 and
    pxcor <= 8 ]
    [ set patch-type 1 
    set pcolor 3 ]
  ask patches with [
    pycor >= earth-top and
    pycor <= earth-top + 4 and
    pxcor >= 2 and
    pxcor <= 3 ]
    [ set patch-type 1 
    set pcolor 3 ]
  ask patches with [
    pycor >= earth-top and
    pycor <= earth-top + 6 and
    pxcor >= -1 and
    pxcor <= 1 ]
    [ set patch-type 1 
    set pcolor 3 ]    
  ask patches with [
    pycor >= earth-top and
    pycor <= earth-top + 2 and
    pxcor >= -5 and
    pxcor <= -3 ]
    [ set patch-type 1 
    set pcolor 3 ]
  ask patches with [
    pycor >= earth-top and
    pycor <= earth-top + 4 and
    pxcor > -4 and
    pxcor <= -3 ]
    [ set patch-type 1 
    set pcolor 3 ]
end 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;; RUNTIME PROCEDURES ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to Go
  ifelse timeofday < 1440 [        ;; time progresses, one minute per tick
    set timeofday timeofday + 1 ][
    set timeofday 0]  
  update-sunlight
  ifelse pollutants-cycle [        ;; add primary pollutants to atmosphere
    cycle-pollutants]
    [add-pollutants]  
  ifelse atmosphere-inversion [    ;; move the pollutants around 
    run-pollutants-inversion]
    [run-pollutants]             
  run-reactions                    ;; perform chemical reactions on primary pollutants 
  tick
end 

to update-sunlight
  ifelse day-night-cycle [    
    if timeofday < sunrise or timeofday >= sunset [      ;; night
      set sunlight 0]
    if timeofday >= sunrise and timeofday < sunset [     ;; day
      set sunlight sunlight-intensity ]
  ]
    [set sunlight sunlight-intensity]
  ask patches with [
    patch-type = 2 ] [ 
    set pcolor scale-color sky sunlight 0 100 ]             ;; color sky according to sunlight level  
end 

to cycle-pollutants            ;;;;;;;;;; Add pollutants to the atmosphere in a cycle ;;;;;;;;;
  if timeofday < sunrise or timeofday >= sunset [            ;; night (no pollution)
    set p-rate 0]
  if timeofday >= sunrise and timeofday < sunset  [           ;; day (max pollution rate)
    set p-rate 1 ]
  create-NO2 p-rate * nitrogen-dioxide [
    set color brown
    setxy random-xcor earth-top + random-float max-pycor]
  create-VOC p-rate * volatile-organic-compounds [
    set color gray
    setxy random-xcor earth-top + random-float max-pycor]   
end 

to add-pollutants             ;;;;;;;;;; Add pollutants to the atmosphere (no cycle) ;;;;;;;;;
  create-NO2 nitrogen-dioxide [
    set color brown
    setxy 0 -6
    set heading -90 + random-float 180 ]
  create-VOC volatile-organic-compounds [
    set color gray
    setxy 0 -6 
    set heading -90 + random-float 180 ] 
end 

to run-pollutants-inversion    ;;;;;;;;;; Make pollutants move around, with atmospheric inversion ;;;;;;;;;
  ask turtles [
    let dist 0.5 * random-float 1
    ifelse can-move? dist
      [ fd dist ]
      [ set heading heading + 180 
        fd dist ]                     ;; if turtle hits the edge of the world, it turns around
    if ycor <= earth-top 
      [ set heading random-float 80 ] ;; if turtle hits ground, return it to atmosphere
  ]
end 

to run-pollutants              ;;;;;;;;;; Make pollutants move around (no inversion) ;;;;;;;;;
  ask turtles [
    let dist 0.5 * random-float 1
    ifelse can-move? dist
      [ fd dist ]
      [ die ]                  ;; if turtle hits the edge of the world, it dies
    if ycor <= earth-top 
      [ set heading random 80 ] ;; if turtle hits ground, return it to atmosphere
  ]
end 

to run-reactions                 ;;;;;;;;;;; Chemical reactions WITHOUT O2 turtles ;;;;;;;;;;
  set NO2-rate sunlight                   ;; set reaction rate for splitting NO2
  set O3-rate sunlight                    ;; set reaction rate for forming O3
  ask NO2 [                               ;; (1) SPLIT NO2 
    if random-float 200 < NO2-rate [        ;; roll dice to see if there'e enough energy for reaction
      set breed O1                            ;; if there is, convert the NO2 turtle to O1
      set color yellow
      hatch-NO 1 [                            ;; and create a new turtle that is NO
        set color orange 
        set heading random-float 360 ]]      
  ]
  ask O1 [                                ;; (2) CREATE OZONE
    if random-float 200 < O3-rate [         ;; roll dice to see if there's enough energy for reaction
        set breed O3                          ;; if there is, convert O1 turtle to O3
        set color violet
        set heading random-float 360 ]       
  ]
  ask NO [                                ;; (3) DESTROY OZONE
    if any? O3-here [                       ;; see if there are any O3 in the patch
      set breed NO2                           ;; if there are, convert NO turtle to NO2
      set color brown
      set heading random-float 360 ]  
    let reactant one-of O3-here               ;; destroy used up O3 turtle 
    if reactant != nobody [
      ask reactant [ die ]]
  ]  
  ask VOC [                               ;; (4) CREATE OXIDANTS
    if any? NO-here [                       ;; see if any NO are in the patch
      set breed Ox                            ;; if there are, convert VOC turtle to Ox
      set color red
      set heading random-float 360 ]     
    let reactant one-of NO-here               ;; destroy used up NO turtle
    if reactant != nobody [
      ask reactant [die]]
  ] 
end