rm(list=ls())
par(mfrow=c(1,1))
#-------------------------------------------------------------------------------
#  Änderung des Grundwasserspiegel bei nach einem Infiltrationsereignis 
#  Annahme: hydraulisches Gleichgewicht (q(z) = 0) vor und nach dem Ereignis
#-------------------------------------------------------------------------------
# Aufgaben: a) Variieren Sie die Infiltrationsmenge  
#           b) Variieren Sie die Parameter des MVG-Modells
#
# --- Parametrisierung der Wasserspannungskurve ---
#
theta_VG <- function(alpha,n,theta_s,theta_r,psi) {
    m <- 1-1/n
    theta_VG <- theta_r + (theta_s-theta_r)/(1+(alpha*psi)^n)^m
}
#
# --- Mualemmodell für die relative hydraulische Leitfähigkeit ---
#
K_VGrel <- function(alpha,n,theta_s,theta_r,psi) {
    m <- 1-1/n
    K_VGrel <- (1-(alpha*psi)^(n-1)*(1+(alpha*psi)^n)^(-m))^2/
                 (1+(alpha*psi)^n)^(m/2)
}
#
# --- Modellparameter ---       
#
theta_s <- 0.55   # Sättigungswassergehalt
theta_r <- 0.15   # Residualwassergehalt
K_s <- 100        # Wasserleitfähigkeit in cm/d
alpha <- 0.02     # emp. Formparameter
nvg <- 3          # emp. Formparameter
#
#inf   <- 8         # Infiltrationsmenge
# 
# --- Wassermenge im Profil ---
#
gwl.1 <- 100 
dth   <- 0.1                                                                      
z   <- seq(0,gwl.1,dth) # Tiefe
gwl.2 <- 75.
#
theta <- function(gwl){
  n   <- length(z)
  th <- 1:n
  for (i in 1:length(z)){
    if(z[i] >= gwl) {
      th[i] <- theta_s
    } else {
      psi <- z[i] - gwl
      #print(z);print(psi)
      th[i] <- theta_VG(alpha,nvg,theta_s,theta_r,abs(psi))
    }
  }
  th
}
storage <- function(gwl){
  th1 <- theta(gwl.1)
  th2 <- theta(gwl)
  stg <- sum((th2-th1)*dth)
  list(th1,th2,stg)
}  
b <- storage(gwl.2); 
cat("\n\n Ergebnisse: \n\n")
cat(" Infiltrationsmenge: ",10*b[[3]]," Liter/qm \n")
cat(" GW-Spiegel vorher:  ",gwl.1," cm unter Bodenoberfläche \n")  
cat(" GW-Spiegel nachher: ",gwl.2," cm unter Bodenoberfläche \n\n")  
#
par(mfrow=c(1,2))
th1 <- b[[1]]
th2 <- b[[2]]
plot(th1,z,xlim=c(0,1.2*max(th2)),ylim=c(gwl.1,0),type="l",lwd=2,col="blue",main="Vorher")
plot(th2,z,xlim=c(0,1.2*max(th2)),ylim=c(gwl.1,0),type="l",lwd=2,col="blue",main="Nachher")
stop("++++ hier ++++")     
# Skriptende