| 1 | SUBROUTINE INCOG(A,X,GIP) |
|---|
| 2 | C |
|---|
| 3 | C =================================================== |
|---|
| 4 | C Purpose: Compute the incomplete gamma function |
|---|
| 5 | C r(a,x), â(a,x) and P(a,x) |
|---|
| 6 | C Input : a --- Parameter ( a ó 170 ) |
|---|
| 7 | C x --- Argument |
|---|
| 8 | C Output: GIN --- r(a,x) |
|---|
| 9 | C GIM --- â(a,x) |
|---|
| 10 | C GIP --- P(a,x) |
|---|
| 11 | C Routine called: GAMMA for computing â(x) |
|---|
| 12 | C =================================================== |
|---|
| 13 | C |
|---|
| 14 | IMPLICIT DOUBLE PRECISION (A-H,O-Z) |
|---|
| 15 | XAM=-X+A*DLOG(X) |
|---|
| 16 | IF (XAM.GT.700.0.OR.A.GT.170.0) THEN |
|---|
| 17 | WRITE(*,*)'INCGAMMAFUNC: variables a and/or x too large' |
|---|
| 18 | STOP |
|---|
| 19 | ENDIF |
|---|
| 20 | IF (X.EQ.0.0) THEN |
|---|
| 21 | GIN=0.0 |
|---|
| 22 | CALL GAMMAFUNC(A,GA) |
|---|
| 23 | GIM=GA |
|---|
| 24 | GIP=0.0 |
|---|
| 25 | ELSE IF (X.LE.1.0+A) THEN |
|---|
| 26 | S=1.0D0/A |
|---|
| 27 | R=S |
|---|
| 28 | DO 10 K=1,60 |
|---|
| 29 | R=R*X/(A+K) |
|---|
| 30 | S=S+R |
|---|
| 31 | IF (DABS(R/S).LT.1.0D-15) GO TO 15 |
|---|
| 32 | 10 CONTINUE |
|---|
| 33 | 15 GIN=DEXP(XAM)*S |
|---|
| 34 | CALL GAMMAFUNC(A,GA) |
|---|
| 35 | GIP=GIN/GA |
|---|
| 36 | GIM=GA-GIN |
|---|
| 37 | ELSE IF (X.GT.1.0+A) THEN |
|---|
| 38 | T0=0.0D0 |
|---|
| 39 | DO 20 K=60,1,-1 |
|---|
| 40 | T0=(K-A)/(1.0D0+K/(X+T0)) |
|---|
| 41 | 20 CONTINUE |
|---|
| 42 | GIM=DEXP(XAM)/(X+T0) |
|---|
| 43 | CALL GAMMAFUNC(A,GA) |
|---|
| 44 | GIN=GA-GIM |
|---|
| 45 | GIP=1.0D0-GIM/GA |
|---|
| 46 | ENDIF |
|---|
| 47 | END |
|---|
