gammaprior_function(l,u,p=0.9)
{

# given a probability interval for the standard deviation
# calculates an appropriate gamma prior for the precision

pp<-(1-p)/2
ppp<-p+pp

f<-function(cc)
	{
	a_cc[1]
	b_cc[2]
	c1<-pgamma(1/(u*u),a,1/b)-pp
	c2<-pgamma(1/(l*l),a,1/b)-ppp
	c1*c1 + c2*c2
	}

sds<-qnorm(ppp)-qnorm(pp)
mean<-0.5*((1/(u*u))+(1/(l*l)))
sd<-((1/(l*l))-(1/(u*u)))/sds
var<-sd*sd
a_mean*mean/var
b_mean/var
print(paste("Normal approx",a,b))

nlmout<-nlm(f,c(a,b),gradtol=1e-10,steptol=1e-10,iterlim=1000)
est<-nlmout$estimate
a<-est[1]
b<-est[2]
print(paste("Numerical solution",a,b))

op_par(mfrow=c(3,1))
vals_rgamma(10000,a,1/b)
hist(vals,30,main=paste("Precision ~ Gamma(",a,",",b,")",sep=""))
hist(1/vals,30,main="Variance")
hist(1/sqrt(vals),30,main="SD")
par(op)

pl<-qgamma(pp,a,1/b)
pu<-qgamma(ppp,a,1/b)
print(paste("Actual",p,"probability interval",sqrt(1/pu),sqrt(1/pl)))

sqerr <- (sqrt(1/pu)-l)^2 + (sqrt(1/pl)-u)^2
if (sqerr > 1) {
	warning(" ***** FAILED TO FIND APPROPRIATE PRIOR ***** ")
}

c(a,b)
}