normal.bestfit <- function(mode, p.up, x.up, sd.step = 0.0001){
# This function returns the mean and standard deviation parameters 
# for a normal prior given the mode and upper p percentile.

p.up <- round(p.up, digits = 3)

mean <- mode
sd <- 0

i <- 1
repeat{
sd <- sd + sd.step
percent <- pnorm(x.up,mean,sd)
if(percent <= p.up){
sd.final <- sd
break
}
i <- i + 1
}
mean.final <- mean
final <- c(mean.final, sd.final)
return(final)
}