rm(list=ls())
library(VGAM)
erfc <- function(x) {erfc <- 1-erf(x)}
#
# Umsetzung der analytischen Lösung für eine Stufenfunktion in R
#
C.btc.rel <- function(D,R,L,pv,v) {
  C.btc.rel <- 0.5*erfc((R-pv)*(sqrt((v*L/D)/(4*R*pv))))
}
#
# http://de.wikipedia.org/wiki/Fehlerfunktion
#
R <- 2.8                   # 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,3,0.05)         # Zeit-array für Berechnungen
pv <- q*t/(theta*L)        # Berechnung einer dimensionslosen Zeit
                           # was ist pv physikalisch bzw. anschaulich?
C.rel <- C.btc.rel(D,R,L,pv,v)
print(C.rel)
plot(pv,C.rel)
result <- data.frame(round(pv,3), round(C.rel+rnorm(length(pv),0.0,0.02),3)) 
# Ausgabe wird mit rnorm() etwas verrauscht (warum?)
names(result) <- c("pv","Conc")
write.table(result,"BTC_meas.txt",sep="\t",row.names=F)
plot(result[,1],result[,2])
#
# --- Ende ---