1 DOUBLE PRECISION FUNCTION bfrac(a,b,x,y,lambda,eps) 2C----------------------------------------------------------------------- 3C CONTINUED FRACTION EXPANSION FOR IX(A,B) WHEN A,B .GT. 1. 4C IT IS ASSUMED THAT LAMBDA = (A + B)*Y - B. 5C----------------------------------------------------------------------- 6C .. Scalar Arguments .. 7 DOUBLE PRECISION a,b,eps,lambda,x,y 8C .. 9C .. Local Scalars .. 10 DOUBLE PRECISION alpha,an,anp1,beta,bn,bnp1,c,c0,c1,e,n,p,r,r0,s, 11 + t,w,yp1 12C .. 13C .. External Functions .. 14 DOUBLE PRECISION brcomp 15 EXTERNAL brcomp 16C .. 17C .. Intrinsic Functions .. 18 INTRINSIC abs 19C .. 20C .. Executable Statements .. 21C-------------------- 22 bfrac = brcomp(a,b,x,y) 23 IF (bfrac.EQ.0.0D0) RETURN 24C 25 c = 1.0D0 + lambda 26 c0 = b/a 27 c1 = 1.0D0 + 1.0D0/a 28 yp1 = y + 1.0D0 29C 30 n = 0.0D0 31 p = 1.0D0 32 s = a + 1.0D0 33 an = 0.0D0 34 bn = 1.0D0 35 anp1 = 1.0D0 36 bnp1 = c/c1 37 r = c1/c 38C 39C CONTINUED FRACTION CALCULATION 40C 41 10 n = n + 1.0D0 42 t = n/a 43 w = n* (b-n)*x 44 e = a/s 45 alpha = (p* (p+c0)*e*e)* (w*x) 46 e = (1.0D0+t)/ (c1+t+t) 47 beta = n + w/s + e* (c+n*yp1) 48 p = 1.0D0 + t 49 s = s + 2.0D0 50C 51C UPDATE AN, BN, ANP1, AND BNP1 52C 53 t = alpha*an + beta*anp1 54 an = anp1 55 anp1 = t 56 t = alpha*bn + beta*bnp1 57 bn = bnp1 58 bnp1 = t 59C 60 r0 = r 61 r = anp1/bnp1 62 IF (.NOT.(abs(r-r0).GT.eps*r)) GO TO 20 63C 64C RESCALE AN, BN, ANP1, AND BNP1 65C 66 an = an/bnp1 67 bn = bn/bnp1 68 anp1 = r 69 bnp1 = 1.0D0 70 GO TO 10 71C 72C TERMINATION 73C 74 20 bfrac = bfrac*r 75 RETURN 76 77 END 78