rm(list=ls())
#--------------------------------------------------------------------------
#   Einstellung von Graphikparametern
#--------------------------------------------------------------------------
par(mfrow=c(1,1))
#
# eine eigene Pausenfunktion:
pause <- function ()
{
    cat("Pause. Press <Enter> to continue...")
    readline()
    invisible()
}
#
# Funktion zur Berechnung der Porosität
porosity <- function(rho_b,rho_s) 
{
  porosity <- 1-rho_b/rho_s
}
#--------------------------------------------------------------------------
#   Erstellen eines Grids
#--------------------------------------------------------------------------
ndataX <- 40
ndataY <- 40
rhob_min <- 0.5; rhob_max <- 1.8
rhos_min <- 1.3; rhos_max <- 3.0
x <- seq(rhob_min,rhob_max,(rhob_max-rhob_min)/(ndataX-1))
y <- seq(rhos_min,rhos_max,(rhos_max-rhos_min)/(ndataY-1))
pause()
#
# Erzeugen der Datenmatrix mit der Funktion "porosity"
por <- outer(x,y,FUN="porosity")
por[,][por[,] < 0] <- 0 # Elimimieren der negativen Werte
#--------------------------------------------------------------------------
# Graphiken
#--------------------------------------------------------------------------
#
#  Contourplot
#
#filled.contour(x,y,por,asp=1,color = topo.colors,nlevels=50
filled.contour(x,y,por,color = terrain.colors, nlevels=50
, xlab = "bulk density", ylab = "\nparticle density", main="Porosity"
, plot.axes={axis(1);axis(2)
; contour(x,y,por,nlevels=20,add=T,lwd=1,col="brown", axes=F)
})
pause()
#
#  3D-Plot
#
persp(x,y,por, zlim = c(0,1),theta = 30, phi = 20, expand = 0.7, col = "lightblue",
     ltheta = 120, shade = 0.75, ticktype = "detailed",
     xlab = "bulk density", ylab = "\nparticle density", zlab = "porosity", main="Porosity",
     r=4,border="lightblue4") -> rotmat
#
#   --- Schluss ---
#