% NORMINTERVAL  Generates normally distributed random numbers.
%    X = NORMINTERVAL(MU,SIGMA) uses the approximate inverse transform
%    method to generate a random number from a normal distribution with
%    the mean MU and the standard deviation SIGMA.
%    X = NORMINTERVAL(MU,SIGMA,MIN,MAX) generates a random number larger than 
%    or equal to MIN and less than or equal to MAX.
%
%    [X,XSTAR] = NORMINTERVAL(...) also returns the complementary random
%    number.
%
%    The function uses the uniform random number generator
%    RAND. The sequence of random numbers generated by NORMINTERVAL can
%    be reset by changing the state of the RAND generator.
%

function [randnum,comprandnum] = norminterval(mu,sigma,minval,maxval)

if nargin < 3
  U = rand(1);
  Ustar = 1 - U;
else
  Xmin = (minval-mu)/sigma; Xmax = (maxval-mu)/sigma;
  Umin = phi(Xmin); Umax = phi(Xmax);
  U0 = (Umax-Umin)*rand(1);
  U = U0 + Umin;
  Ustar = Umax - U0;
end
X = normvalue(U);
randnum = X*sigma + mu;
X = normvalue(Ustar);
comprandnum = X*sigma + mu;