function b = gammainc(x,a)
%GAMMAINC The incomplete gamma function.
%	Y = GAMMAINC(X,A) evaluates the incomplete gamma function at all
%	the elements of X.  X must be real.  A must be a real scalar.
%	gammainc(x,a) = (integral from 0 to x of t^(a-1) exp(-t) dt)/gamma(a).
%	Note that gammainc(x,a) approaches 1 as x approaches infinity.
%
%	See also GAMMA, GAMMALN.

%	C.R. Denham 6-9-88, C. B. Moler, 10-17-90, 11-4-92.
%	Ref: Abramowitz & Stegun, Handbook of Mathemtical Functions, sec. 6.5.
%	Copyright (c) 1984-93 by The MathWorks, Inc.
%

if length(a) > 1
   error('The second parameter must be a scalar.')
end

b = x;
lngama = gammaln(a);

k = find(x == 0);
if any(k)
   b(k) = 0*k;
end

k = find((x ~= 0) & (x < a+1));
if any(k)
   % Series expansion for x < a+1
   ap = a;
   sum = ones(size(k))/a;
   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*log(x(k)) - lngama);
end

k = find(x >= a+1);
if any(k)
   % Continued fraction for x >= a+1
   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;
      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;
      b(k) = 1 - exp(-x(k) + a*log(x(k)) - lngama) .* g;
   end
end

k = find(x == inf);
if any(k)
   b(k) = ones(size(k));
end
