####################################################
##  Matematično modeliranje okoljskih procesov    ##
##  3. vaja: 2D resevanje enacbe transporta toka  ##
##                                                ##
##  Izračun in izris grafov                       ##
##                                                ##
##  Izdelal: Lovrenc Pavlin                       ##
##  Datum: 7.4.2015                               ##
####################################################

##  Podatki  ##
b <- 59 #sirina kanala
Q <- 29 #pretok
h <- 1 #globina
q0 <- 0.29 #dotok polutanta
c0 <- 100 #zacetna koncentracija
K <- 0.001 #mesalni koeficient

##  Izracun hitrosti toka  ##
u <- Q/(b*h)

## Izracun koncentracij
#x <- c(50,100,500,1000,5000,10000,15000,20000)
x <- seq(50,20000,10)
#y <- c(0,1,2,3,4,5,10,15,20)
y <- seq(-20,20,0.1)

z <- matrix(0,nrow=length(x),ncol=length(y))
for (j in 1:length(y)){
  for (i in 1:length(x)){
    z[i,j] <- q0*c0/(2*h*u*sqrt(pi*x[i]/u*K))*exp(-u*y[j]**2/(4*K*x[i]))
  }
}

#############
##  Grafi  ##
#############

# Graf celotnega opazovanega obmocja
png(filename="konc_celotno.png",width=1000,height=740,units="px")
filled.contour(x,y,z,
               nlevels=40, 
               color.palette=rainbow, 
               plot.title = title(main = "Koncentracija polutanta v kanalu",
                            xlab = "stacionaža [m]", ylab = "širina [m]"),
               key.title = title(main = "konc.\n[%]"))
dev.off()


# Graf omejenega območja x[50,200],y=[-5,5]
png(filename="konc_delno.png",width=1000,height=740,units="px")
filled.contour(x,y,z,
               nlevels=40, 
               color.palette=rainbow,
               xlim=range(50:200,finite=TRUE),
               ylim=range(-5:5,finite=TRUE),
               plot.title = title(main = "Koncentracija polutanta v kanalu",
                                  xlab = "stacionaža [m]", ylab = "širina [m]"),
               key.title = title(main = "konc.\n[%]"))
dev.off()


# Graf po prečnih prerezih
png(filename="konc_precni.png",width=1000,height=740,units="px")
cl=rainbow(8) #paleta barv
plot(y,z[1,],type="l", col=cl[1],lwd=2, main="Koncentracija po prečnih presekih", xlab="y [m]", ylab="c [%]")
lines(y,z[4,],type="l", col=cl[2],lwd=2)
lines(y,z[30,],type="l", col=cl[3],lwd=2)
lines(y,z[100,],type="l", col=cl[4],lwd=2)
lines(y,z[300,],type="l", col=cl[5],lwd=2)
lines(y,z[500,],type="l", col=cl[6],lwd=2)
lines(y,z[700,],type="l", col=cl[7],lwd=2)
lines(y,z[900,],type="l", col=cl[8],lwd=2)
legend("topleft", lwd=2, legend = c(paste("X =",toString(x[1])),paste("X=",toString(x[4])),paste("X=",toString(x[30])),paste("X=",toString(x[100])),paste("X=",toString(x[300])),paste("X=",toString(x[500])),paste("X=",toString(x[700])),paste("X=",toString(x[900]))), col=cl)
dev.off()


# Graf po vzdolžnih prerezih
png(filename="konc_vzdolzi.png",width=1000,height=740,units="px")
cl2=rainbow(6) #paleta barv
Y=c(201,211,221,231,261,301)
plot(x,z[,Y[1]],type="l", col=cl2[1],lwd=2, main="Koncentracija po vzdolžnih presekih", xlab="x [m]", ylab="c [%]")
lines(x,z[,Y[2]],type="l", col=cl2[2],lwd=2)
lines(x,z[,Y[3]],type="l", col=cl2[3],lwd=2)
lines(x,z[,Y[4]],type="l", col=cl2[4],lwd=2)
lines(x,z[,Y[5]],type="l", col=cl2[5],lwd=2)
lines(x,z[,Y[6]],type="l", col=cl2[6],lwd=2)
legend("topright",
       lwd=2,
       legend = c(paste("y =",toString(y[Y[1]])),
                  paste("y =",toString(y[Y[2]])),
                  paste("y =",toString(y[Y[3]])),
                  paste("y =",toString(y[Y[4]])),
                  paste("y =",toString(y[Y[5]])),
                  paste("y =",toString(y[Y[6]]))),
       col=cl2)
dev.off()