hpd.gamma <- function(data.df,g.shape,g.rate,target=0.90,steps=300,tol=0.01) {
    if (steps %% 2 == 1) steps <- steps + 1
        g.mode <- sum(data.df$N)/sum(data.df$N*data.df$dur)
        g.range <- seq(qgamma(0.001,g.shape,g.rate), qgamma(0.999,g.shape, g.rate),length=steps)
        g.range <- c(g.range[1:(steps/2)],g.mode,g.range[(steps/2+1):steps])
        g.dens <- dgamma(g.range,g.shape,g.rate)
        g.probs <- pgamma(g.range,g.shape,g.rate)
        for (i in 1:(steps/2)) {
            k.dir <- which(c(g.dens[(steps/2-i)],g.dens[(steps/2+i)]) == max(g.dens[(steps/2-i)],g.dens[(steps/2+i)]))
            k <- c(g.dens[(steps/2-i)],g.dens[(steps/2+i)])[k.dir]
            k.loc <- c((steps/2-i),(steps/2+i))[k.dir]
            if (k.dir == 2) k2.range <- c(1:(steps/2))
            else k2.range <- c((steps/2 + 1):steps)
            k2.min <- which(abs(k-g.dens[k2.range])==min(abs(k-g.dens[k2.range])))
            if (k.dir == 1) k2.min <- k2.min + steps/2
            if (g.probs[k.loc] + (1-g.probs[k2.min]) < 1-target) break
            bounds <- c(g.range[k.loc],g.range[k2.min])
         }
    return(list("cdf.vals"=c(g.probs[k.loc],g.probs[k2.min]), "bounds"=bounds,"k"=k))
}