home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / SPECFUN.DI$ / GAMMAINC.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  1.7 KB  |  71 lines
  1. function b = gammainc(x,a)
  2. %GAMMAINC The incomplete gamma function.
  3. %    Y = GAMMAINC(X,A) evaluates the incomplete gamma function at all
  4. %    the elements of X.  X must be real.  A must be a real scalar.
  5. %    gammainc(x,a) = (integral from 0 to x of t^(a-1) exp(-t) dt)/gamma(a).
  6. %    Note that gammainc(x,a) approaches 1 as x approaches infinity.
  7. %
  8. %    See also GAMMA, GAMMALN.
  9.  
  10. %    C.R. Denham 6-9-88, C. B. Moler, 10-17-90, 11-4-92.
  11. %    Ref: Abramowitz & Stegun, Handbook of Mathemtical Functions, sec. 6.5.
  12. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  13. %
  14.  
  15. if length(a) > 1
  16.    error('The second parameter must be a scalar.')
  17. end
  18.  
  19. b = x;
  20. lngama = gammaln(a);
  21.  
  22. k = find(x == 0);
  23. if any(k)
  24.    b(k) = 0*k;
  25. end
  26.  
  27. k = find((x ~= 0) & (x < a+1));
  28. if any(k)
  29.    % Series expansion for x < a+1
  30.    ap = a;
  31.    sum = ones(size(k))/a;
  32.    del = sum;
  33.    while norm(del,'inf') >= 10*eps*norm(sum,'inf')
  34.       ap = ap + 1;
  35.       del = x(k) .* del/ap;
  36.       sum = sum + del;
  37.    end
  38.    b(k) = sum .* exp(-x(k) + a*log(x(k)) - lngama);
  39. end
  40.  
  41. k = find(x >= a+1);
  42. if any(k)
  43.    % Continued fraction for x >= a+1
  44.    a0 = ones(size(k));
  45.    a1 = x(k);
  46.    b0 = zeros(size(k));
  47.    b1 = a0;
  48.    fac = 1;
  49.    n = 1;
  50.    g = b1;
  51.    gold = b0;
  52.    while norm(g-gold,'inf') >= 10*eps*norm(g,'inf');
  53.       gold = g;
  54.       ana = n - a;
  55.       a0 = (a1 + a0*ana) .* fac;
  56.       b0 = (b1 + b0*ana) .* fac;
  57.       anf = n*fac;
  58.       a1 = x(k) .* a0 + anf .* a1;
  59.       b1 = x(k) .* b0 + anf .* b1;
  60.       fac = 1 ./ a1;
  61.       g = b1 .* fac;
  62.       n = n+1;
  63.       b(k) = 1 - exp(-x(k) + a*log(x(k)) - lngama) .* g;
  64.    end
  65. end
  66.  
  67. k = find(x == inf);
  68. if any(k)
  69.    b(k) = ones(size(k));
  70. end
  71.