function log_y = log_normpdf(x,mu,sigma)
%LOG_NORMPDF Logarithm of normal probability density function (pdf).
%   LOG_Y = LOG_NORMPDF(X,MU,SIGMA) Returns the normal pdf with mean, MU, 
%   and standard deviation, SIGMA, at the values in X. 
%
%   The size of LOG_Y is the common size of the input arguments. A scalar input  
%   functions as a constant matrix of the same size as the other inputs.     
%
%   Default values for MU and SIGMA are 0 and 1 respectively.
  
%   References:
%      [1]  M. Abramowitz and I. A. Stegun, "Handbook of Mathematical
%      Functions", Government Printing Office, 1964, 26.1.26.
  
  if nargin < 3, 
    sigma = 1;
  end
  
  if nargin < 2;
    mu = 0;
  end
  
  if nargin < 1, 
    error('Requires at least one input argument.');
  end
  
  [errorcode x mu sigma] = distchck(3,x,mu,sigma);
  
  if errorcode > 0
    error('Requires non-scalar arguments to match in size.');
  end
  
  %   Initialize Y to zero.
  log_y = repmat( -Inf, size(x) );
  
  k = find(sigma > 0);
  if any(k)
    xn = (x(k) - mu(k)) ./ sigma(k);
    log_y(k) = -0.5 * xn .^2 - log(sqrt(2*pi) .* sigma(k));
  end
  
  % Return NaN if SIGMA is negative or zero.
  k1 = find(sigma <= 0);
  if any(k1)
    tmp   = NaN;
    log_y(k1) = tmp(ones(size(k1))); 
  end
