function gcf, a, x, gln=gln ;+ ; NAME: ; gcf ; PURPOSE: (one line) ; Returns the incomplete gamma function Q(a; x) ; DESCRIPTION: ; Returns the incomplete gamma function Q(a; x) evaluated by its continued fraction representation ; as gammcf. Also returns ln Gamma(a) as gln. ; CATEGORY: ; Math ; CALLING SEQUENCE: ; y = gcf(a,x, gln=gln) ; INPUTS: ; a - first parameter for Q ; x - second parameter for Q ; OPTIONAL INPUT PARAMETERS: ; none ; KEYWORD INPUT PARAMETERS: ; none ; KEYWORD OUTPUT PARAMETERS: ; gln = ln(Gamma(a)) ; OUTPUTS: ; Q(a; x) ; COMMON BLOCKS: ; None ; SIDE EFFECTS: ; RESTRICTIONS: ; None ; PROCEDURE: ; MODIFICATION HISTORY: ; Written 2003 November, by Leslie Young, SwRI ; based on Numerical Recipes gcf, Section 6.2 ;- itmax = 100 ;Maximum allowed number of iterations. eps = 3.0e-7 ; Relative accuracy. fpmin = 1.0e-30 ;Number near the smallest representable floating-point number. gln=lngamma(a) b=x+1.0-a; Set up for evaluating continued fraction by modi ed Lentz's method (x5.2) with b0 = 0. c=1.0/FPMIN; d=1.0/b; h=d; for i=1, ITMAX do begin ; Iterate to convergence. an = -i*(i-a); b = b + 2.0; d = an*d+b; if (abs(d) lt FPMIN) then d=FPMIN; c = b+an/c; if (abs(c) lt FPMIN) then c=FPMIN; d=1.0/d; del=d*c; h = h * del; if (abs(del-1.0) lt EPS) then goto, jump1; endfor jump1: $ if i gt ITMAX then message, 'gcf - a too large, ITMAX too small'; return, exp(-x+a*alog(x)-(gln))*h; Put factors in front. end