; $Id: beta.pro,v 1.13 2004/01/21 15:54:48 scottm Exp $ ; Copyright (c) 1995-2004, Research Systems, Inc. All rights reserved. ; Unauthorized reproduction prohibited. ;+ ; NAME: ; BETA ; ; PURPOSE: ; Return the Beta function of (possibly complex) Z. ; ; CALLING SEQUENCE: ; Result = BETA(Z) ; ; INPUTS: ; Z: The expression for which the beta function will be evaluated. ; If Z is double-precision, the result is double-precision, ; otherwise the result is floating-point. Z may be complex. ; ; KEYWORD PARAMETERS: ; DOUBLE: Set this keyword to return a double-precision result. ; ; MODIFICATION HISTORY: ; 3 July 1995, AB, RSI. ; AB, 5/4/2001, Switch from using _EXTRA to _STRICT_EXTRA, so that ; incorrect keywords will cause issue proper messages to ; be issued instead of being silently ignored. ; CT, RSI, Dec 2001: Rewrote to use LNGAMMA, which now handles complex. ;- function beta, z, w, DOUBLE=double, _REF_EXTRA=_extra ON_ERROR, 2 tz = SIZE(z,/TYPE) tw = SIZE(w,/TYPE) isComplex = (tz eq 6) or (tz eq 9) or (tw eq 6) or (tw eq 9) doComplex = isComplex ; If either Z or W has a negative number, ; then we must use complex for the LNGAMMA. if (not isComplex) then doComplex = $ (not ARRAY_EQUAL(z,0b,/GE)) or (not ARRAY_EQUAL(w,0b,/GE)) ; Be sure to do the calculation using double precision. z1 = doComplex ? DCOMPLEX(z) : DOUBLE(z) w1 = doComplex ? DCOMPLEX(w) : DOUBLE(w) result = EXP(LNGAMMA(z1, _STRICT_EXTRA=_extra) + $ LNGAMMA(w1, _STRICT_EXTRA=_extra) - $ LNGAMMA(z1+w1, _STRICT_EXTRA=_extra), _STRICT_EXTRA=_extra) ; Return double precision if keyword is set, or either arg is double. doDouble = (N_ELEMENTS(double) gt 0) ? KEYWORD_SET(double) : $ (tz eq 5) or (tz eq 9) or (tw eq 5) or (tw eq 9) ; We may need to convert from DCOMPLEX back to DOUBLE ; because Z or W might have had negative numbers. return, doDouble ? (isComplex ? result : DOUBLE(result)) : $ (isComplex ? COMPLEX(result) : FLOAT(result)) end