#--------------------------------------------------------------
#  Tiefenverlauf von Wassergehalt, Matrixpotential 
#                und Gravitationspotential
#  
#  Annahme: steady rate (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) {      
  if(psi >= 0) { theta_VG <- theta_s} else {
    m <- 1-1/n
    theta_VG <- theta_r + (theta_s-theta_r)/(1+(alpha*abs(psi))^n)^m
  }
}
#
# --- Mualemmodell für die hydraulische Leitfähigkeit ---
#
K_VG <- function(alpha,n,theta_s,theta_r,K.sat,psi) {
  if(psi >= 0) { K_VG <- K.sat} else {
    m <- 1-1/n
    K_VG <- K.sat*(1-(alpha*abs(psi))^(n-1)*(1+(alpha*abs(psi))^n)^(-m))^2/
                 (1+(alpha*abs(psi))^n)^(m/2)
  }
}
#
# --- Modellparameter Horizont 1 ---       
#
theta_s.h1 <- 0.55   # Sättigungswassergehalt
theta_r.h1 <- 0.1   # Residualwassergehalt
K_s.h1     <- 20    # Wasserleitfähigkeit in cm/d
alpha.h1   <- 0.05    # emp. Formparameter
n.h1       <- 3   # emp. Formparameter
#
# --- Modellparameter Horizont 2 ---       
#
theta_s.h2 <- 0.60   # Sättigungswassergehalt
theta_r.h2 <- 0.20   # Residualwassergehalt
K_s.h2     <- 2    # Wasserleitfähigkeit in cm/d
alpha.h2   <- 0.01    # emp. Formparameter
n.h2       <- 2    # emp. Formparameter
#
# --- Profil- und Geometriedaten ---
#
Hor.u <- 80       # Untergrenze des zweiten Materials
Hor.o <- 50       # Obergrenze des zweiten Materials
gwl <- 150        # Tiefe des Grundwasserspiegels in cm
#
# Potentiale und Wassergehalte bei vorgegebenem Fluss 
# unter steady-rate-Bedingungen
#
q <- 3
n.z <- 200 # n.z muss gross sein, um Oszillationen zu vermeiden
delta.z <- gwl/n.z
z   <- seq(gwl,0,-delta.z) # Tiefe
n.points <- length(z); 
psix <- 1:n.points; theta <- 1:n.points; K_rel <- 1:n.points; K <- 1:n.points
psix[1] <- 0
theta[1] <- theta_s.h1
#
for (i in 1:(n.points-1)) {
  if((z[i+1] < Hor.u)&&(z[i+1] >= Hor.o)){
    alpha <- alpha.h2; n <- n.h2; theta_s <- theta_s.h2; theta_r <- theta_r.h2
    K_s <- K_s.h2
    cat("+++++ ",i,"  ",K_s,"\n")
  } else { 
    alpha <- alpha.h1; n <- n.h1; theta_s <- theta_s.h1; theta_r <- theta_r.h1
    K_s <- K_s.h1
    cat("*****",i,"  ",K_s,"\n")
  }
  psi.u <- psix[i]
  K.komp <- K_VG(alpha,n,theta_s,theta_r,K_s,psi.u)               
    cat(".....",K.komp,"  ",psi.u,"  \n")
  psi.o <- psi.u + (q/K.komp-1)*delta.z
  psi <- (psi.u+psi.o)/2
  K.komp <- K_VG(alpha,n,theta_s,theta_r,K_s,psi)
  psi.o <- psi.u + (q/K.komp-1)*delta.z
  psix[i+1] <- psi.o
  theta[i+1] <- theta_VG(alpha,n,theta_s,theta_r,psix[i+1])
}
#plot(psix,z,ylim=c(max(z),0),xlim=c(-max(z),100),type="l",lwd=3,col="blue") 
#
par(mfrow=c(1,2))                  
plot(psix,z,ylim=c(max(z),0),type="n",lwd=2,col="blue",main="Matrixpotential")
polygon(c(min(psix)-100,min(psix)-100,max(psix)+100,max(psix)+100),c(Hor.u,Hor.o,Hor.o,Hor.u),
border="gray70",col="gray70")               
lines(psix,z,    ylim=c(max(z),0),type="l",lwd=2,col="blue")
#
plot(theta,z,   ylim=c(max(z),0),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(max(z),0),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