rm(list=ls())
par(mfrow=c(1,1),mar=c(5, 6, 3, 3) + 0.1)
#library(VGAM) # enthält eine Fehlerfunktion erf
#
# Eine selbst geschriebene Fehlerfunktion
erf  <- function(x) {
  fun <- function(x){exp(-x^2)}
  z <- x
  intgr <- integrate(fun,0,z)
  erf.def <- (2/sqrt(pi))*intgr$value
}
# komplementäre Fehlerfunktion
erfc <- function(x) {erfc <- 1-erf(x)}
#
# Umsetzung der analytischen Lösung für eine Stufenfunktion in R
#
C.btc.rel <- function(D,R,z,t,v) {
  C.btc.rel <- 0.5*erfc((R*z-v*t)/(2*sqrt(D*R*t)))
}
#
# http://de.wikipedia.org/wiki/Fehlerfunktion
#
R <- 2.                    # Retardationskoeffizient
D <- 3.                    # Diffusions-/Dispersionskoeffizient
L <- 15                    # Länge der Säule
q <- 10                    # Flussdichte
theta <- 0.4               # volumetrischer Wassergehalt
v <- q/theta               # Abstandsgeschwindigkeit des Wassers
t <- seq(0.1,1,0.05)       # Zeit-Array für Berechnungen
z <- seq(0.0,L,0.5)        # Tiefen-Array für Berechnungen
datmat <- outer(z,t)       # erzeugt eine Matrix für die Konzentrationen
C.rel <- datmat
for(i in 1:(length(z))){
  for(j in 1:length(t)){
      C.rel[i,j] <- C.btc.rel(D,R,z[i],t[j],v)
  }
}
# C.rel <- C.btc.rel()  # dimensionslose Konzentrationen: 0 <= C <=1
#
#-------------------------------------------------------------------------------
# Graphiken
#-------------------------------------------------------------------------------
#filled.contour(z,t,C.rel,color = topo.colors,nlevels=50
filled.contour(t,z,t(C.rel),color = topo.colors,nlevels=50, ylim=c(max(z),0)
   , main="breakthrough curve",xlab="Zeit (d)",ylab="\n\n Tiefe (cm) "
   ,plot.axes={axis(1);axis(2)
   ; contour(t,z,t(C.rel),nlevels=10,add=T,lwd=1,col="black", axes=F)
})
#
persp(z,t,C.rel,theta = 30,phi = 20, expand = 0.7, 
     col = "lightblue",ltheta = 120, shade = 0.75, ticktype = "detailed",
     xlab = "Tiefe", ylab = "\n Porenvolumina", zlab = "Konzentration"
     , main="breakthrough curve",r=4,border="lightblue")