rm(list=ls())
library(DAAG)
#-------------------------------------------------------------
# Dateneilesen und Berechnungen
#-------------------------------------------------------------
daten <- read.table("AB_2_profil.txt",header=T)
attach(daten,warn.conflicts = F)
P <- (1 - BD/SD)*100
LG <- P - WV 
PZ <- P/(100-P)
Zm <- (Zo+Zu)/2 
Delta_Z <- Zu-Zo
#
cat(" Porenraum und Porenziffer:\n",
    "==========================\n")
profile_agg <- data.frame(Zm,P,WV,LG,100-P,PZ)
names(profile_agg) <- c("z ","n ","th","lg ","sg ","pz ")
print(signif(profile_agg,dig=2))
#-------------------------------------------------------------
#  Normalspannung durch trockenen Boden, Wasser; Gesamtgewicht
#  Einfachste Form einer numerischen Integration
#-------------------------------------------------------------
# 
G_acc <- 9.81
Wfc <- rep(0,4); Wdc <- rep(0,4); Wtc <- rep(0,4)
#
W_d <- 10*Delta_Z*BD*G_acc # Faktor 10 überprüfen!!
W_f <- WV*G_acc/100
W_t <- W_d + W_f
#
Wdc[1] <- W_d[1]; Wfc[1] <- W_f[1]; Wtc[1] <- W_t[1] 
for(i in 2:4) {
  W_d[i] <- W_d[i-1] + W_d[i]
  W_f[i] <- W_f[i-1] + W_f[i]
  W_t[i] <- W_t[i-1] + W_t[i]
}
#-------------------------------------------------------------
# Plots: Porenraum, Lagerung
#-------------------------------------------------------------
par(mfrow=c(2,3))
plot(Zm,SD,type="l",lwd=3,col="gray",main="Substanzdichte")
points(Zm,SD,pch=15,cex=1.5,col="blue")
plot(Zm,BD,type="l",lwd=3,col="gray",main="Trockenraumdichte")
points(Zm,BD,pch=15,cex=1.5,col="blue")
plot(Zm,P,type="l",lwd=3,col="gray",main="Porosität")
points(Zm,P,pch=15,cex=1.5,col="blue")
plot(Zm,WV,type="l",lwd=3,col="gray",main="Wassergehalt")
points(Zm,WV,pch=15,cex=1.5,col="blue")
plot(Zm,LG,type="l",lwd=3,col="gray",main="Luftgehalt")
points(Zm,LG,pch=15,cex=1.5,col="blue")
plot(Zm,PZ,type="l",lwd=3,col="gray",main="Porenziffer")
points(Zm,PZ,pch=15,cex=1.5,col="blue")
#-------------------------------------------------------------
# Plots: Auflast
#-------------------------------------------------------------
par(mfrow=c(2,3))
plot(Zu,W_d,type="l",lwd=3,col="gray",main="Auflast\ntrocken")
points(Zu,W_d,pch=15,cex=1.5,col="blue")
plot(Zu,W_f,type="l",lwd=3,col="gray",main="Auflast\nfeucht")
points(Zu,W_f,pch=15,cex=1.5,col="blue")
plot(Zu,W_t,type="l",lwd=3,col="gray",main="Auflast\ngesamt")
points(Zu,W_t,pch=15,cex=1.5,col="blue")
#
detach(daten)