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 must be real.  X and A must
%   be the same size (or either can be scalar).  The incomplete gamma
%   functions is defined as: 
%
%    gammainc(x,a) = 1 ./ gamma(a) .*
%       integral from 0 to x of t^(a-1) exp(-t) dt
%
%   For any a, 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-97 by The MathWorks, Inc.
%   $Revision: 5.8 $  $Date: 1997/04/08 06:42:55 $
%@@@
realmin=10E-50;
gam = gammaln(a+realmin);
if max(size(x))==1; x = x(ones(size(a))); end
if max(size(a))==1; a = a(ones(size(x))); gam = gam(ones(size(x))); end

b = x;

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

% Series expansion for x < a+1

k = find((a ~= 0) & (x ~= 0) & (x < a+1));
if k~=[]
   ap = a(k);
   sum = 1./ap;
   del = sum;
   while max(abs(del)) >= 10*2.2204e-016*max(abs(sum))
      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 k~=[]
   a0 = ones(size(k));
   a1 = x(k);
   b0 = zeros(size(k));
   b1 = a0;
   fac = 1;
   n = 1;
   g = b1;
   gold = b0;
   while max(abs(g-gold)) >= 10*2.2204e-016*max(abs(g));
      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

end