;$Id: ibeta.pro,v 1.17 2001/01/15 22:28:05 scottm Exp $
;
; Copyright (c) 1994-2001, Research Systems, Inc.  All rights reserved.
;       Unauthorized reproduction prohibited.
;+
; NAME:
;       IBETA
;
; PURPOSE:
;       This function computes the incomplete beta function, Ix(a, b).
;
; CATEGORY:
;       Special Functions.
;
; CALLING SEQUENCE:
;       Result = Ibeta(a, b, x)
;
; INPUTS:
;       A:    A scalar or array of type integer, float or double that
;             specifies the parametric exponent of the integrand.
;
;       B:    A positive scalar or array of type integer, float or double that
;             specifies the parametric exponent of the integrand.
;
;       X:    A scalar, in the interval [0, 1], 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.
;
; EXAMPLE:
;       Compute the incomplete beta function for the corresponding elements
;       of A, B, and X.
;       Define the parametric exponents.
;         A = [0.5, 0.5, 1.0, 5.0, 10.0, 20.0]
;         B = [0.5, 0.5, 0.5, 5.0,  5.0, 10.0]
;       Define the the upper limits of integration.
;         X = [0.01, 0.1, 0.1, 0.5, 1.0, 0.8]
;       Compute the incomplete beta functions.
;         result = Ibeta(A, B, X)
;       The result should be:
;         [0.0637686, 0.204833, 0.0513167, 0.500000, 1.00000, 0.950736]
;
; 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
;                    IBETA is based on the routines: betacf.c, betai.c and
;                    gammln.c described in section 6.2 of Numerical Recipes,
;                    The Art of Scientific Computing (Second Edition), and is
;                    used by permission.
;            CT, March 2000, added DOUBLE, EPS, ITER, ITMAX keywords.
;-

function ibetacf, a, b, x, $
	DOUBLE = double, $
	EPS = eps, $
	ITER = m, $
	ITMAX = maxit
  COMPILE_OPT hidden

  on_error, 2
  m = 0 ; zero iterations so far
  double = KEYWORD_SET(double)
  IF (N_ELEMENTS(eps) LT 1) THEN eps   = (double) ? 3.0d-12 : 3.0e-7
  fpmin = (double) ? 1d-300 : 1.0e-30
  IF (N_ELEMENTS(maxit) LT 1) THEN maxit = 100
  qab = a + b
  qap = a + 1
  qam = a - 1
    c = 1
    d = 1 - qab * x / qap
  if(abs(d) lt fpmin) then d = fpmin
  d = 1 / d
  h = d
  for m = 1, maxit do begin
    m2 = 2 * m
    aa = m * (b - m) * x / ((qam + m2) * (a + m2))
     d = 1 + aa*d
     if(abs(d) lt fpmin) then d = fpmin
     c = 1 + aa / c
     if(abs(c) lt fpmin) then c = fpmin
     d = 1 / d
     h = h * d * c
     aa = -(a + m) *(qab + m) * x/((a + m2) * (qap + m2))
     d = 1 + aa * d
     if(abs(d) lt fpmin) then d = fpmin
     c = 1 + aa / c
     if(abs(c) lt fpmin) then c = fpmin
     d = 1 / d
     del = d * c
     h = h * del
     if(abs(del - 1) lt eps) then return, h
  endfor
  message, 'Failed to converge within given parameters.'
end

function ibeta, a, b, x, $
	DOUBLE = double, $
	EPS = eps, $
	ITER = iter, $
	ITMAX = itmax

  on_error, 2
  iter = 0 ; zero iterations so far
  xmin = min(x, max = xmax, /NAN)
  if xmin lt 0 or xmax gt 1 then message, $
    'X must be in the interval (0, 1).'

  if (min(a, /NAN) le 0) or (min(b, /NAN) le 0) then message, $
    'A and B must be positive'


  na = n_elements(a)
  nb = n_elements(b)
  nx = n_elements(x)
  n = na > nb > 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(b,/N_DIMENSION) GT 0) AND (nb LE n) THEN BEGIN
	dim = SIZE(b,/DIMENSIONS)
	n = nb
  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 BEGIN
	double = (SIZE(x,/TYPE) EQ 5) OR (SIZE(a,/TYPE) EQ 5) $
		OR (SIZE(b,/TYPE) EQ 5)
  ENDIF ELSE $
	double = KEYWORD_SET(double)

  ap = 0L
  bp = 0L
  xp = 0L
  y = (double) ? 0d : 0.0
  if n gt 1 then y = replicate(y, n) ;Make Result

  for i=0L, n-1 do begin         ;Each result
      aa = (double) ? DOUBLE(a[ap]) : a[ap]
      bb = (double) ? DOUBLE(b[bp]) : b[bp]
      xx = (double) ? DOUBLE(x[xp]) : x[xp]
      if (xx lt 1.0 and xx gt 0.0) then $
        bt = exp(lngamma(aa + bb) - lngamma(aa) - lngamma(bb) + $
                 aa * alog(xx) + bb * alog(1.0 - xx)) $
      else bt = (double) ? 0d : 0.0
      if (xx lt (aa + 1.0)/(aa + bb + 2.0)) then $
        y[i] =  bt * ibetacf(aa, bb, xx, $
        	DOUBLE=double, EPS=eps, ITER=m, ITMAX=itmax) / aa $
      else y[i] = 1.0 - bt * ibetacf(bb, aa, 1.0-xx, $
        	DOUBLE=double, EPS=eps, ITER=m, ITMAX=itmax) / bb
      iter = iter > m
      ap = (ap + 1) < (na - 1)
      bp = (bp + 1) < (nb - 1)
      xp = (xp + 1) < (nx - 1)
  endfor
  IF (dim[0] GT 0) THEN y = REFORM([y], dim, /OVERWRITE)
  return, y
end