rm(list=ls())
#------------------------------------------------------
#  Auswertung einer bodenhydrologischen Feldbeprobung
#------------------------------------------------------
# Erzeugen von dataframes zur Aufnahme 
# der berechneten Daten:
#
Dbd <- data.frame(cbind(Dbd1=1,Dbd2=1,Dbd3=1,Dbd4=rep(0,10)))
Mbd <- data.frame(cbind(Mbd1=1,Mbd2=1,Mbd3=1,Mbd4=rep(0,10)))
Pco <- data.frame(cbind(Mbd1=1,Mbd2=1,Mbd3=1,Mbd4=rep(0,10)))
Wco <- data.frame(cbind(Wco1=1,Wco2=1,Wco3=1,Wco4=rep(0,10)))
Aco <- data.frame(cbind(Aco1=1,Aco2=1,Aco3=1,Aco4=rep(0,10)))
Sco <- data.frame(cbind(Sco1=1,Sco2=1,Sco3=1,Sco4=rep(0,10)))
#
# Einlesen der Rohdaten
#
werte <- read.table("profile_data.txt"); werte
#
# Durchführung der Berechnungen
#
rho_s <- werte[,2]
#
# Schleife über die 4 Wiederholungen pro Tiefe
# bd in g/cm^3
# Gehalte in %
#
for(i in 1:4) {
  Dbd[,i] <- (werte[,4+(i-1)*3]-werte[,3+(i-1)*3])/100 # dry bulk density
  Mbd[,i] <- (werte[,5+(i-1)*3]-werte[,3+(i-1)*3])/100 # moist bulk density
  Pco[,i] <- (1 - Dbd[,i]/rho_s)*100                   # porosity ompartment
  Wco[,i] <-  werte[,5+(i-1)*3]-werte[,4+(i-1)*3]      # water content compartment
}
Aco <- Pco - Wco # air content in compartment
Sco <- 100 - Pco # solid content in compartment
#
# Berechnung der aggregierten Daten (Mittelwert und Standardabweichung)
#
Pr_agg <- data.frame(z=werte[,1],
               Dbd_mu=apply(Dbd,1,mean),Dbd_sd=apply(Dbd,1,sd),
               Mbd_mu=apply(Mbd,1,mean),Mbd_sd=apply(Mbd,1,sd),
               Pco_mu=apply(Pco,1,mean),Pco_sd=apply(Pco,1,sd),
               Wco_mu=apply(Wco,1,mean),Wco_sd=apply(Wco,1,sd),
               Aco_mu=apply(Aco,1,mean),Aco_sd=apply(Aco,1,sd),
               Sco_mu=apply(Sco,1,mean),Sco_sd=apply(Sco,1,sd))
Pr_agg
#-------------------------------------------------------------
#  Normalspannung durch trockenen Boden, Wasser; Gesamtgewicht
#  Einfachste Form einer numerischen Integration
#  W_d: kumulatives Trockengewicht
#  W_f: kumulatives Gewicht Wasser
#  W_t: kumulatives Gewicht gesamt
#-------------------------------------------------------------
# 
G_acc <- 9.81
W_d <- rep(0,10); W_f <- rep(0,10); W_t <- rep(0,10)
attach(Pr_agg)
#
W_d[1] <- 100*Dbd_mu[1]*G_acc
W_f[1] <- Wco_mu[1]*G_acc
W_t[1] <- W_d[1] + W_f[1]
#
for(i in 2:10) {
  W_d[i] <- W_d[i-1] + 100*Dbd_mu[i]*G_acc
  W_f[i] <- W_f[i-1] +   Wco_mu[i]*G_acc
  W_t[i] <- W_d[i] + W_f[i]
}
NSp <- data.frame(z=werte[,1],W_d=W_d,W_f=W_f,W_t=W_t)
NSp
#
# Ausgabe mit 4 signifikanten Stellen
#
ausgabe <- signif(data.frame(Pr_agg,W_d=W_d,W_f=W_f,W_t=W_t),dig=4);ausgabe
write.table(ausgabe,"profile_erg.txt",sep="\t")
#
detach(Pr_agg)
#
#--- Skriptende ---
#
# Aufgaben:
# ---------
# 1) Analysieren Sie das Skript sehr genau
# 2) Überprüfen Sie die Berechnung der Normalspannungen (Welche Einheiten kommen
#    heraus?
# 3) Schauen Sie sich die zu Beginn des Skripts erzeugten dataframes an
# 4) Warum taucht in der Berechnung der Normalspannungen die 
#    Kompartimentmächtigkeit nicht auf?
#