rm(list=ls())
#-------------------------------------------------------------------------------
# Auswertung eines Durchbruchsexperiments durch inverse Modellierung
#-------------------------------------------------------------------------------
library(VGAM)
#
erfc <- function(x) {erfc <- 1-erf(x)}
#
# Umsetzung der analytischen Lösung für eine Durchbruchsfunktion in R
#
C.btc.rel <- function(D,R,pv) {
  C.btc.rel <- 0.5*erfc((R-pv)*(sqrt((v*L/D)/(4*R*pv))))
}
# siehe hierzu:
# http://de.wikipedia.org/wiki/Fehlerfunktion
#
# --- Messdaten ----
#
daten <- read.table("BTC_meas.txt",header=T)
pv.lab   <- daten[,1]
conc.lab <- daten[,2]
plot(pv.lab,conc.lab)
#
# --- Versuchsbedingungen ---
#
L <- 15                     # Länge der Säule
q <- 10                    # Flussdichte
theta <- 0.4               # volumetrischer Wassergehalt
v <- q/theta               # Abstandsgeschwindigkeit des Wassers
#
# --- Startschätzungen ---
#
R.ini <- 2.                    # Retardationskoeffizient
D.ini <- 3.                    # Diffusions-/Dispersionskoeffizient
#
# --- Parameteroptimierung ---
#
y <- conc.lab; x <- pv.lab
model <- nls(y~C.btc.rel(D,R,x),model=TRUE,
             start=list(D=D.ini,R=R.ini))
cat("\n Zusammenfassung:\n") 
cat(" ----------------\n")
print(summary(model))
cat("\n Details: \n")
cat(" --------\n")
print(model)
par <- model$m$getAllPars()
D.est <- par[1]
R.est <- par[2]
cat("\n Parameter: \n")
cat(" -----------\n")
cat("\n optimierter Parameter D: ",D.est,"\n" )
cat("\n optimierter Parameter R: ",R.est,"\n" )
C.rel.est <- C.btc.rel(D.est,R.est,x)
plot(conc.lab,C.rel.est)
#
pv   <- seq(0,max(pv.lab),0.05)         # Zeit-array für Berechnungen
Cest <- C.btc.rel(D.est,R.est,pv)        # Berechnung einer dimensionslosen Zeit
                           # was ist pv physikalisch bzw. anschaulich?
plot(pv,Cest,type="l",lwd=2,col="blue",ylim=c(0,1.2),
     main="Messdaten und Modell")
points(pv.lab,conc.lab,pch=15,col="darkgreen",cex=1.2)
#
# --- Ende ---