function b = gammainc(x,a)
%GAMMAINC Incomplete gamma function.
%   Y = GAMMAINC(X,A) evaluates the incomplete gamma function for
%   corresponding elements of X and A.  X and A must be real and the same
%   size (or either can be a scalar).  A must also be non-negative.
%   The incomplete gamma function is defined as: 
%
%    gammainc(x,a) = 1 ./ gamma(a) .*
%       integral from 0 to x of t^(a-1) exp(-t) dt
%
%   For any a>=0, as x approaches infinity, gammainc(x,a) approaches 1.
%   For small x and a, gammainc(x,a) ~= x^a, so gammainc(0,0) = 1.
%
%   See also GAMMA, GAMMALN.

%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 5.10 $  $Date: 1997/11/21 23:45:35 $

gam = gammaln(a+realmin);
if all(size(x)==1), x = x(ones(size(a))); end
if all(size(a)==1), a = a(ones(size(x))); gam = gam(ones(size(x))); end

b = x;

k = find(x == 0);
if ~isempty(k)
   b(k) = 0;
end
k = find(a == 0);
if ~isempty(k)
   b(k) = 1;
end

% Series expansion for x < a+1

k = find((a ~= 0) & (x ~= 0) & (x < a+1));
if ~isempty(k)
   ap = a(k);
   sum = 1./ap;
   del = sum;
   while norm(del,'inf') >= 10*eps*norm(sum,'inf')
      ap = ap + 1;
      del = x(k) .* del ./ ap;
      sum = sum + del;
   end
   b(k) = sum .* exp(-x(k) + a(k).*log(x(k)) - gam(k));
end

% Continued fraction for x >= a+1

k = find((a ~= 0) & (x ~= 0) & (x >= a+1));
if ~isempty(k)
   a0 = ones(size(k));
   a1 = x(k);
   b0 = zeros(size(k));
   b1 = a0;
   fac = 1;
   n = 1;
   g = b1;
   gold = b0;
   while norm(g-gold,'inf') >= 10*eps*norm(g,'inf');
      gold = g;
      ana = n - a(k);
      a0 = (a1 + a0 .*ana) .* fac;
      b0 = (b1 + b0 .*ana) .* fac;
      anf = n*fac;
      a1 = x(k) .* a0 + anf .* a1;
      b1 = x(k) .* b0 + anf .* b1;
      fac = 1 ./ a1;
      g = b1 .* fac;
      n = n + 1;
   end
   b(k) = 1 - exp(-x(k) + a(k).*log(x(k)) - gam(k)) .* g;
end
