function bb = betaval(mn,sd,fx)
%BETAVAL(mean, sd, Fx)
% This function calculates a random number 
% from a beta distribution with mean mn, standard deviation
% sd, and cum. distr. function fx.
% This function uses the MATLAB function betainc(x, vv,ww),
% where x is the value of beta, v,w are beta parameters
% that are called a and b in the text. 
if sd == 0; bb = mn;
else
 toler = 0.0001; % this is tolerance of answer: how close
 % the CDF value of the answer must be to the input value (Fx)
 var = sd^2; 
 if var >=(1-mn)*mn disp('sd too high for beta'), pause, end; 
 % this checks that the input mean and st. deviation
 % are possible for a beta. 

 vv = mn*((mn.*(1-mn)/(var))-1); % calculate the beta parameters
 ww = (1-mn).*((mn.*(1-mn)/(var))-1);

 upval = 1; lowval = 0; x = 0.5+ 0.02*rand; 
 % start with a beginning guess x; the use of rand 
 % adds wiggle to the search start to avoid pathologies
 i = betainc(x,vv,ww); % find the CDF value for x
 
 % the following while loop searches for ever better
 % values of x, until the value has a CDF within the 
 % toler of Fx (unless the value
 % is very close to 0 or 1, which will also terminate
 % the search)
  while ( (toler < abs(i-fx))&(x >1e-6)&((1-x)>1e-6) )
	if fx > i 
	   lowval = x; x = (upval+lowval)/2;
        else
	   upval = x; x = (upval+lowval)/2;  
 	end; %if
        i = betainc(x,vv,ww);
  end; % while

 % This makes values of x somewhat random to eliminate
 % pathologies when variance is very small or large.
 % It also truncates values of x, with the 
 % smallest values equal to toler and the biggest
 % equal to 1 - toler.
 bbb = x + toler*0.1*(0.5-rand);
 if bbb < toler; bbb = toler; end;
 if bbb > 1; bbb = 1- toler; end;

 bb=bbb;
end; %else