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
sed_time <- function(d,z) 
{
  sed_time <- 10*(z*18*visc)/(g_acc*(rho_s - rho_w)*d^2)
  # *10 wegen unterschiedl. Einheiten
}
#--------------------------------------------------------------------------
#   Erstellen eines Grids
#--------------------------------------------------------------------------
ndata_d <- 60
ndata_z <- 60
rho_s <- 2.65                                # g/cm^3
rho_w <- 1.00                                # g/cm^3
visc  <- 0.001002                            # kg/(m sec)
g_acc <- 9.81                                # m/sec^2
d_min <- 0.002; d_max <- 0.2                 # mm
z_min <- 5; z_max <- 50                      # cm
z <- seq(z_min,z_max,by=(z_max - z_min)/(ndata_z-1))
d <- seq(d_min,d_max,by=(d_max - d_min)/(ndata_d-1))
pause()
#
# Erzeugen der Datenmatrix mit der Funktion "porosity"
sed_t <- outer(d,z,FUN="sed_time")
#--------------------------------------------------------------------------
# Graphiken
#--------------------------------------------------------------------------
#
#  Contourplot
#
#filled.contour(x,y,por,asp=1,color = topo.colors,nlevels=50
data <- log(sed_t,10)
lgd <- log(d,10)
filled.contour(lgd,z,data,ylim=c(50,5),color = terrain.colors, nlevels=50
, xlab = "lg(diameter (mm))", ylab = "\ndepth (cm)", main="Sedimentation\ntime"
, plot.axes={axis(1);axis(2)
; contour(lgd,z,data,nlevels=20,add=T,lwd=1,col="brown", axes=F)
})
#
pause()
#
#  3D-Plot
#
persp(lgd,z,data, zlim = c(0,max(data)),theta = 30, phi = 20, expand = 0.7, col = "lightblue",
     ltheta = 120, shade = 0.75, ticktype = "detailed",
     xlab = "lg(diameter (mm))", ylab = "\ndepth (cm)", zlab = "lg(time(sec)", main="Sedimentation\ntime",
     r=4,border="lightblue4") -> rotmat
#
#   --- Schluss ---
#