% NORMVALUE  Approximate inverse transform of uniform random numbers.
%    X = NORMVALUE(U) transform the a U(0,1) random value to the corresponding
%    N(0,1)-distributed random value.
%

function X = normvalue(U);

c0 = 2.515517; c1 = 0.802853; c2 = 0.010328;
d1 = 1.432788; d2 = 0.189269; d3 = 0.001308;

rows = size(U,1); cols = size(U,2);
for i = 1:rows
  for j = 1:cols
    if U(i,j) >= 1
      X(i,j) = Inf;
    elseif U(i,j) > 0.5
      t = sqrt(-2*log(1-U(i,j)));
      z = t - sum([c0 + c1*t + c2*t^2])/...
	  (1 + sum([d1*t d2*t^2 d3*t^3]));
      X(i,j) = z;
    elseif U(i,j) == 0.5
      X(i,j) = 0;
    elseif U(i,j) > 0
      t = sqrt(-2*log(U(i,j)));
      z = t - sum([c0 + c1*t + c2*t^2])/...
	  (1 + sum([d1*t d2*t^2 d3*t^3]));
      X(i,j) = -z;
    else
      X(i,j) = -Inf;
    end
  end
end
