#--------------------------------------------------------------
#  Tiefenverlauf von Wassergehalt, Matrixpotential 
#                und Gravitationspotential
#  
#  Annahme: hydraulisches Gleichgewicht (q(z) = 0
#           2 Horizonte  
#--------------------------------------------------------------
# Aufgaben: a) Variieren Sie die Lage der Horizonte  
#           b) Verändern Sie die Dicke der Horizonte
#           c) Variieren Sie die MVG-Parameter der Horizonte
#
# --- 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 Horizont 1 ---       
#
theta_s.h1 <- 0.55   # Sättigungswassergehalt
theta_r.h1 <- 0.15   # Residualwassergehalt
K_s.h1     <- 100    # Wasserleitfähigkeit in cm/d
alpha.h1   <- 0.1    # emp. Formparameter
n.h1       <- 2.5    # emp. Formparameter
#
# --- Modellparameter Horizont 2 ---       
#
theta_s.h2 <- 0.50   # Sättigungswassergehalt
theta_r.h2 <- 0.25   # Residualwassergehalt
K_s.h2     <-  50    # Wasserleitfähigkeit in cm/d
alpha.h2   <- 0.04    # emp. Formparameter
n.h2       <- 2.0    # emp. Formparameter
#
# --- Profil- und Geometriedaten ---
#
Hor.u <- 80       # Untergrenze des Horizonts
Hor.o <- 30       # Obergrenze des Horizonts
gwl <- 140        # Tiefe des Grundwasserspiegels in cm
bzn <-  0         # Tiefe des Bezugsniveaus
iz  <- -1         # Indikator für Richtung der z-Achse: 
                  # +1 nach oben, -1 nach unten
par(mfrow=c(2,3))
# 
# --- Verlauf der Potentiale, der Wassergehalte und der Leitfähigkeiten ---
#
z   <- -iz*seq(0,gwl,0.5) # Tiefe
gwl <- -iz*gwl
bzn <- -iz*bzn   
psi <- -iz*(z - gwl)
grp <- -iz*(bzn - z)
H   <-  psi + grp
#
Hor.u <- -iz*Hor.u                                          
Hor.o <- -iz*Hor.o                                          
#
n.z <- length(z); theta <- 1:n.z; K_rel <- 1:n.z; K <- 1:n.z
for (i in 1:n.z){
  if((iz*z[i]>=iz*Hor.u)&&(iz*z[i]<=iz*Hor.o)){
    alpha <- alpha.h2; n <- n.h2; theta_s <- theta_s.h2; theta_r <- theta_r.h2
  } else { 
    alpha <- alpha.h1; n <- n.h1; theta_s <- theta_s.h1; theta_r <- theta_r.h1
  }
  theta[i] <- theta_VG(alpha,n,theta_s,theta_r,abs(psi[i]))
  K_rel[i] <- K_VGrel(alpha,n,theta_s,theta_r,abs(psi)[i])
  K[i] <- K_s*K_rel[i]
}
#
if(iz==1) {z.u <- min(z); z.o <- max(z)} else {z.u <- max(z); z.o <- min(z)}
plot(psi,z,     ylim=c(z.u,z.o),type="n",lwd=2,col="blue",main="Matrixpotential")
polygon(c(min(psi)-100,min(psi)-100,max(psi)+100,max(psi)+100),c(Hor.u,Hor.o,Hor.o,Hor.u),
border="gray70",col="gray70")
lines(psi,z,    ylim=c(z.u,z.o),type="l",lwd=2,col="blue")
#
plot(grp,z,     ylim=c(z.u,z.o),type="n",lwd=2,col="orange",main="Gravitationspotential")
polygon(c(min(grp)-100,min(grp)-100,max(grp)+100,max(grp)+100),c(Hor.u,Hor.o,Hor.o,Hor.u),
border="gray70",col="gray70")
lines(grp,z,    ylim=c(z.u,z.o),type="l",lwd=2,col="blue")
#
plot(H,z,       ylim=c(z.u,z.o),type="n",lwd=2,col="black",main="Gesamtpotential")
polygon(c(min(H)-100,min(H)-100,max(H)+100,max(H)+100),c(Hor.u,Hor.o,Hor.o,Hor.u),
border="gray70",col="gray70")
lines(H,z,    ylim=c(z.u,z.o),type="l",lwd=2,col="blue")
#
plot(theta,z,   ylim=c(z.u,z.o),type="n",lwd=2,col="blue",main="Vol. Wassergehalt")
polygon(c(min(theta)-100,min(theta)-100,max(theta)+100,max(theta)+100),c(Hor.u,Hor.o,Hor.o,Hor.u),
border="gray70",col="gray70")
lines(theta,z,    ylim=c(z.u,z.o),type="l",lwd=2,col="blue")
#
plot(K,z,       ylim=c(z.u,z.o),type="n",lwd=2,col="blue",main="Leitfähigkeit")
polygon(c(min(K)-100,min(psi)-100,max(psi)+100,max(psi)+100),c(Hor.u,Hor.o,Hor.o,Hor.u),
border="gray70",col="gray70")
lines(K,z,    ylim=c(z.u,z.o),type="l",lwd=2,col="blue")
#
plot(log10(K),z,ylim=c(z.u,z.o),type="n",lwd=2,col="blue",main="Log10(K)")
polygon(c(min(log10(K))-100,min(log10(K))-100,max(log10(K))+100,max(log10(K))+100),c(Hor.u,Hor.o,Hor.o,Hor.u),
border="gray70",col="gray70")
lines(log10(K),z,    ylim=c(z.u,z.o),type="l",lwd=2,col="blue")
#
stop("++++ hier ++++")
# Skriptende