;$Id: igamma.pro,v 1.15 2001/01/15 22:28:06 scottm Exp $
;
; Copyright (c) 1994-2001, Research Systems, Inc.  All rights reserved.
;       Unauthorized reproduction prohibited.
;+
; NAME:
;       IGAMMA
;
; PURPOSE:
;       This function computes the incomplete gamma function, Px(a).
;
; CATEGORY:
;       Special Functions.
;
; CALLING SEQUENCE:
;       Result = Igamma(a, x)
;
; INPUTS:
;       A:    A positive scalar or array of type integer, float or double that
;             specifies the parametric exponent of the integrand.
;
;       X:    A positive scalar or array of type integer, float or double that
;             specifies the upper limit of integration.
;
; KEYWORD PARAMETERS:
;
;   DOUBLE = Set this keyword to force the computation to be done
;            in double precision.
;
;	EPS = relative accuracy, or tolerance.  The default tolerance
;	      is 3.0e-7 for single precision, and 3.0d-12 for double
;	      precision.
;
;   ITER = Set this keyword equal to a named variable that will contain
;          the actual number of iterations performed.
;
;   ITMAX = Set this keyword to specify the maximum number of iterations.
;           The default value is 100.
;
;   METHOD:  Use this keyword to specify a named variable which returns
;            the method used to compute the incomplete gamma function.
;            A value of 0 indicates that a power series representation
;            was used. A value of 1 indicates that a continued fractions
;            method was used.
;
; EXAMPLE:
;       Compute the incomplete gamma function for the corresponding elements
;       of A and X.
;       Define the parametric exponents.
;         A = [0.10, 0.50, 1.00, 1.10, 6.00, 26.00]
;       Define the the upper limits of integration.
;         X = [0.0316228, 0.0707107, 5.00000, 1.04881, 2.44949, 25.4951]
;       Compute the incomplete gamma functions.
;         result = Igamma(A, X)
;       The result should be:
;         [0.742026, 0.293128, 0.993262, 0.607646, 0.0387318, 0.486387]
;
; PROCEDURE:
;       IGAMMA computes the incomplete gamma function, Px(a), using either
;       a power series representation or a continued fractions method. If X
;       is less than or equal to A+1, the power series representation is used.
;       If X is greater than A+1, the continued fractions method is used.
;
;	The result will be NaN if the algorithm failed to converge
;	after the designated number of iterations.
;
; REFERENCE:
;       Numerical Recipes, The Art of Scientific Computing (Second Edition)
;       Cambridge University Press
;       ISBN 0-521-43108-5
;
; MODIFICATION HISTORY:
;       Written by:  GGS, RSI, September 1994
;                    IGAMMA is based on the routines: gser.c, gcf.c, and
;                    gammln.c described in section 6.2 of Numerical Recipes,
;                    The Art of Scientific Computing (Second Edition), and is
;                    used by permission.
;       Modified:    GGS, RSI, January 1996
;                    Corrected the case of IGAMMA(a, 0) for a > 0.
;		     DMS, Sept, 1999, Added arrays, and more accurate
;			results for double.
;            CT, March 2000, added DOUBLE, ITER keywords.
;-

FUNCTION igamma, a, x, $
	DOUBLE = double, $
	EPS = eps, $
	ITER = iter, $
	ITMAX = itmax, $
	METHOD = method

  on_error, 2

  iter = 0 ; zero iterations so far
  if min(a, /NAN) le 0 then message, 'A must be positive.'
  if min(x, /NAN) lt 0 then message, 'X must be greater than or equal to zero.'

  na = n_elements(a)
  nx = n_elements(x)
  n = na > nx  ; find largest
  dim = 0  ; assume scalar
  ; now find smallest nonscalar dimensions
  IF (SIZE(x,/N_DIMENSION) GT 0) AND (nx LE n) THEN BEGIN
	dim = SIZE(x,/DIMENSIONS)
	n = nx
  ENDIF
  IF (SIZE(a,/N_DIMENSION) GT 0) AND (na LE n) THEN BEGIN
	dim = SIZE(a,/DIMENSIONS)
	n = na
  ENDIF


  IF (N_ELEMENTS(double) LT 1) THEN $
	double = size(x,/type) eq 5 or size(a,/type) eq 5 $
  ELSE $
	double = KEYWORD_SET(double)

  one = double ? 1.0d0 : 1.0
  bad = double ? !values.d_nan : !values.f_nan
  y = n eq 1 ? bad : replicate(bad, n)
  method = n eq 1 ? 0 : intarr(n)
  ap = 0L
  xp = 0L

  if n_elements(eps) eq 0 then eps = double ? 3.0d-12 : 3.0e-7
  fpmin = double ? 1.0d-300 : 1.0e-30 ;Something really small
  if keyword_set(itmax) eq 0 then itmax = 100

  for i=0L, n-1 do begin        ;Each element
      xx = (double) ? DOUBLE(x[xp]) : x[xp]
      aa = (double) ? DOUBLE(a[ap]) : a[ap]
      IF (xx EQ 0) THEN BEGIN  ; avoid underflow
      	k = 0
      	y[i] = 0
      	GOTO, skip
      ENDIF
      if xx le (aa + one) then begin ;Series Representation.
          ap1 = aa
          sum = one / aa
          del = sum
          for k = 1L, itmax do begin
              ap1 = ap1 + one
              del = del * xx / ap1
              sum = sum + del
              if abs(del) lt abs(sum)*eps then begin
                  y[i] = sum * exp(-xx + aa * alog(xx) - lngamma(aa))
                  goto, skip
              endif
          endfor
      endif else begin          ;Continued Fractions.
          method[i] = 1
          b = xx + one - aa
          c = one / fpmin
          d = one / b
          h = d
          for k = 1L, itmax do begin
              an1 = -k * (k - aa)
              b = b + 2
              d = an1 * d + b
              if abs(d) lt fpmin then d = fpmin
              c = b + an1 / c
              if abs(c) lt fpmin then c = fpmin
              d = one / d
              del = d * c
              h = h * del
              if abs(del - 1) lt eps then begin
                  y[i] = one - (exp(-xx + aa * alog(xx) - lngamma(aa)) * h)
                  goto, skip
              endif
          endfor                ;Itmax
      endelse

; if we reach here, then it never converged...
      iter = itmax
      MESSAGE,'Failed to converge within given parameters.'

skip:
      iter = iter > k   ; keep largest # of iterations
      xp = (xp + 1) < (nx - 1)
      ap = (ap + 1) < (na - 1)
  endfor

  IF (dim[0] NE 0) THEN y = REFORM([y], dim, /OVERWRITE)
  return, y
end