#########################################################
# Ben Hyman
# Pollution Diffusion Framework For ETS Project (J-PAL)
# SCREEN Model (Gaussian Plume) - See Associated PDF
# Last Edited: July 1, 2013
#########################################################

#This code runs under the assumption that each firm emits 150g/s (15g/s * 10 stacks)
#Stack heights are randomized across firms (but not within firms) between 5 - 30 meters
#Real data on source strength, # of stacks, and stack height, will require alterations.

setwd("C:/Users/bhyman/Documents/GIS/GIS Tutorial Materials")
library("foreign", lib.loc="C:/Program Files/R/R-3.0.1/library")
y <- read.csv("distancematrix.csv",header=T,sep= ",")             #load distance matrix
sub_districts <- read.dbf("Split_Subdistricts_GU.dbf")            #load sub_district data

#Set-Up and Assumptions
length <- 488                     #sets width of distance matrix (# of sub-districts)
width <- 404                      #sets length of distance matrix (# of ETS firms)
x <- matrix(0,length,width)       #ground level concentration (micro-g/m^3)
Q <- .015                         #assume source strength = 15mg/s per stack (.015g)
s <- 10                           #assume 10 stacks per firm all at the same height
alpha <- .008                     #assumes neutral horizontal dispersion
beta <- .006                      #assumes neutral vertical dispersion
u <- 10                           #assumes constant wind speed in all directions
PM_before <- 0                    #sets initial PM Level Vector Before Intervention
PM_after <- 0                     #sets initial PM Level Vector After Intervention
PM_data <- 0                      #Sets initial value for final data
set.seed(1)                       #saves specific randomized values on next line
H <- sample(5:30,404,replace=T)   #random stack heights between 5 and 30 meters

#Gaussian Plume Model
for (j in 1:width) {
  for (i in 1:length) {
    x[i,j] <- (Q*s / (pi*(y[i,j]*alpha/((1+.0001*y[i,j])^.5)*(y[i,j]*beta/(1+.00015*y[i,j])^.5))*u))*
  (exp(-.5*(H[j]/(y[i,j]*beta/(1+.00015*y[i,j])^.5))^2))
  }
}

PM_before <- rowSums(x)*1000 #sum over all industry contributions to PM, convert from g to mg
PM_after <- .6 * PM_before #assume the treatment reduces PM to 60% original concentration
PM_data <- cbind(sub_districts,PM_before,PM_after)
write.dbf(PM_data, "Split_Subdistricts_GU.dbf")
#Note, it may not be obsevable in R but the new vairables bind, as can be seen when opening
#the .dbf file in GIS