1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93
|
SUBROUTINE bgrat(a,b,x,y,w,eps,ierr)
C-----------------------------------------------------------------------
C ASYMPTOTIC EXPANSION FOR IX(A,B) WHEN A IS LARGER THAN B.
C THE RESULT OF THE EXPANSION IS ADDED TO W. IT IS ASSUMED
C THAT A .GE. 15 AND B .LE. 1. EPS IS THE TOLERANCE USED.
C IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
C-----------------------------------------------------------------------
C .. Scalar Arguments ..
DOUBLE PRECISION a,b,eps,w,x,y
INTEGER ierr
C ..
C .. Local Scalars ..
DOUBLE PRECISION bm1,bp2n,cn,coef,dj,j,l,lnx,n2,nu,p,q,r,s,sum,t,
+ t2,u,v,z
INTEGER i,n,nm1
C ..
C .. Local Arrays ..
DOUBLE PRECISION c(30),d(30)
C ..
C .. External Functions ..
DOUBLE PRECISION algdiv,alnrel,gam1
EXTERNAL algdiv,alnrel,gam1
C ..
C .. External Subroutines ..
EXTERNAL grat1
C ..
C .. Intrinsic Functions ..
INTRINSIC abs,dlog,exp
C ..
C .. Executable Statements ..
C
bm1 = (b-0.5D0) - 0.5D0
nu = a + 0.5D0*bm1
IF (y.GT.0.375D0) GO TO 10
lnx = alnrel(-y)
GO TO 20
10 lnx = dlog(x)
20 z = -nu*lnx
IF (b*z.EQ.0.0D0) GO TO 70
C
C COMPUTATION OF THE EXPANSION
C SET R = EXP(-Z)*Z**B/GAMMA(B)
C
r = b* (1.0D0+gam1(b))*exp(b*dlog(z))
r = r*exp(a*lnx)*exp(0.5D0*bm1*lnx)
u = algdiv(b,a) + b*dlog(nu)
u = r*exp(-u)
IF (u.EQ.0.0D0) GO TO 70
CALL grat1(b,z,r,p,q,eps)
C
v = 0.25D0* (1.0D0/nu)**2
t2 = 0.25D0*lnx*lnx
l = w/u
j = q/r
sum = j
t = 1.0D0
cn = 1.0D0
n2 = 0.0D0
DO 50 n = 1,30
bp2n = b + n2
j = (bp2n* (bp2n+1.0D0)*j+ (z+bp2n+1.0D0)*t)*v
n2 = n2 + 2.0D0
t = t*t2
cn = cn/ (n2* (n2+1.0D0))
c(n) = cn
s = 0.0D0
IF (n.EQ.1) GO TO 40
nm1 = n - 1
coef = b - n
DO 30 i = 1,nm1
s = s + coef*c(i)*d(n-i)
coef = coef + b
30 CONTINUE
40 d(n) = bm1*cn + s/n
dj = d(n)*j
sum = sum + dj
IF (sum.LE.0.0D0) GO TO 70
IF (abs(dj).LE.eps* (sum+l)) GO TO 60
50 CONTINUE
C
C ADD THE RESULTS TO W
C
60 ierr = 0
w = w + u*sum
RETURN
C
C THE EXPANSION CANNOT BE COMPUTED
C
70 ierr = 1
RETURN
END
|