# Brooks-Corey (1966)
brookscorey <- function(theta,tred,tmax,ksat,pi){
# Standard values
ths <- seq(0.05,0.4,0.005)
trs <- 0.035
tms <- 0.45                # maximale Feuchte
kss <- 1E-5                # hydraul. Leitfähigkeit in m/s
ips <- 0.2                 # pore size index
# Define values
if (is.null(theta)) {
        theta <- ths       # Default Bodenfeuchte
  }
if (is.null(tred))  {
        tred <- trs        # Default Restfeuchte
  }
if (is.null(tmax))  {
        tmax <- tms        # Default max. Feuchte
  }
if (is.null(ksat))  {
        ksat <- kss        # Default Durchlaessigkeit
  }
if (is.null(pi))    {
        pi <- ips          # Default Porenindex
  }
# Calculate unsaturated hydraulic conductivity
Se <- (theta-tred)/(tmax-tred)  
if ((Se < 0) || (Se > 1)) {
       result <- "Error in the definition of moisture."
       }
else{
  kunsat <- Se^(2+3*pi)
  if (is.null(kunsat)) {
    
                result <- "Error in calculation of unsaturated conductivity."
            }
  else{
    ratio  <- Se^(2+3*pi)
    kunsat <- ratio * ksat
    result <- list(Se,ratio,kunsat)}
  }
return(result)    
}

#Test
unsat <- brookscorey(0.25,0.05,0.4,1E-4,0.1)

# Values
ths <- seq(0.05,0.4,0.005)
trs <- 0.045
tms <- 0.45                # maximale Feuchte
kss <- 1E-5                # hydraul. Leitfähigkeit in m/s
ips <- 0.2                 # pore size index

unsat <- brookscorey(ths,trs,tms,kss,ips)
Se <- unsat[[1]]
ku <- unsat[[3]]

plot(Se,ku)