ComputeDi.MaxBurkhart <- function(hi, d, h, b1, b2, b3, b4, a1, a2){
    x <- hi / h
    x1 <- x - 1
    x2 <- x ^ 2 - 1
    di <- d * sqrt(b1 * x1 + b2 * x2 + b3 * (a1 - x) ^ 2 * ((a1 - x) >= 0.0) + b4 * (a2 - x) ^ 2 * ((a2 - x) >= 0.0))
    return(di)
}