IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C This file: http://ftp.aset.psu.edu/pub/ger/fortran/hdk/weibull.for
C
C INPUT File: WEIBULL.INP
C OUTPUT File: WEIBULL.OUT
C
C Using the sample data here, this program takes less than 0.2 seconds
C to run on a 200 MHz PC via Lahey Fortran 90 or Digital Visual Fortran.
C
C  F(T;A,B,C)= 1-EXP(-((T-A)/B)**C)
C-------------GIVEN A CENSORED SET OF M DATA POINTS, (X(I),I=1,M),
C                  ORDERED OR UNORDERED, FROM A SAMPLE OF SIZE N,
C                  AND THE VALUE OF THE LOCATION PARAMETER, A, OF THE
C                  WEIBULL DISTRIBUTION WITH DISTRIBUTION FUNCTION...
C                       F(X,A,B,C)= 1-EXP(-((X-A)/B)**C), FOR
C                                   X .GT. A AND B,C .GT. 0
C                  THIS PROGRAM Tries to SOLVES FOR THE MAXIMUM
C                  LIKELIHOOD ESTIMATES FOR THE SCALE PARAMETER B AND THE
C                  SHAPE PARAMETER C. If this computation does not
C                  converge, the program will try to compute the Least
C                  Squares estimates for B and C.
C
C  METHOD... THE ROOT OF D(LOGL)/D(C) WHERE L IS THE LIKELIHOOD FUNCTION
C                  IS COMPUTED VIA A NEWTON-RAPHSON ITERATION WITH THE
C                  STARTING VALUE, C0, BEING READ IN.  IF C0 IS GIVEN AS
C                  ZERO, THE FUNCTION ESTC IS USED TO DERIVE AN INITIAL
C                  GUESS OF C.
C
C
C  FREE FORMAT INPUT...  ITEM 1     -  A,C0,N,TITLE
C                        ITEM 2   X(I) T(I        (I=1, M )
C                        ITEM 3   -1   -1
C                        ITEM 4                   (Blank line)
C
C                  CHARLES E. ANTLE, H. D. KNOBLE - 1972.
C                  S. P. PENNYPACKER, H. D. KNOBLE - 1974.
C                  L. V. MADDEN - 1978/1979.
C                  H. D. Knoble - 1990.
C                  THE PENNSYLVANIA STATE UNIVERSITY
C
C=======================================================================
      CHARACTER*40 TITLE
      CHARACTER*1 FF
      DIMENSION X(100),T(100),TIME(100),DISEAS(100)
      LOGICAL ERROR
      DATA FF/''/
      NDIM=100
      OPEN(UNIT=50,FILE='WEIBULL.INP',STATUS='OLD')
      OPEN(UNIT=60,FILE='WEIBULL.OUT',STATUS='UNKNOWN')
      DO 1000 NSETS=1,99999
      READ(50,*,END=999) A,C0,N,TITLE
      CALL INPUT(X,M,TIME,DISEAS,NPTS)
      IF(M.LE.NDIM .AND. M.LE.N) GO TO 3
      WRITE(60,2) NDIM,M,N
2     FORMAT(' $$$$$ CURRENT DATA ARRAY DIMENSION IS',I5,
     1 ' M=',I6, ' N=',I6,'.  M AND/OR N ARE INCORRECTLY SPECIFIED.')
999   STOP
C-------------TRANSFORM THE DATA TO MAKE THE PROGRAM SCALE INVARIANT.
3     BIGT=-1.D+30
      DO 5 I=1,M
      TI=X(I)-A
      IF(BIGT.GE. TI) GO TO 5
      BIGT=TI
      LOC=I
5     T(I)=TI
      TI=T(M)
      T(M)=T(LOC)
      T(LOC)=TI
      BIGT=1.D0/BIGT
      DO 6 I=1,M
6     T(I)=T(I)*BIGT
C-------------GET AN ESTIMATE OF C.
      C=C0
      IF(C.EQ.0) C=ESTC(N,M,T)
      C0=C
C-------------SOLVE FOR THE ML ESTIMATES OF B AND C.
      CALL NEWTON(B,C,N,M,T,KONVRG)
      B=B/BIGT
      CALL RESID(A,B,C,NPTS,TIME,DISEAS,SUMSQ,RSQ)
      WRITE(60,65) FF,TITLE, SUMSQ,RSQ
65    FORMAT(A1,A40/'MAXIMUM LIKELIHOOD ESTIMATES, SUM OF SQUARES=',
     1 G16.8, 2X,' RSQ= ',G10.3,5X)
      CB=C/B
      WRITE(60,7) A,B,C,M,N,CB
7     FORMAT(/'A=',G16.8,' B=',G16.8,' C=',G16.8,' M=',I5,' N=',I5,2X,
     1' C/B=',G16.8)
      CALL SUMMRY(A,B,C,WMODE,WMEAN,WVAR,ERROR)
      WRITE(60,17) WMODE,WMEAN,WVAR
 17   FORMAT(/'PDF:  MODE= ',G16.8,' MEAN= ',G16.8,' VARIANCE= ',G16.8)
      IF(KONVRG.EQ.1) GO TO 9
      WRITE(60,8) C0
8     FORMAT(' $$$$ SUBROUTINE NEWTON DID NOT ACHIEVE CONVERGENCE WITH I
     1NITIAL GUESS OF C =',G16.8)
9     CALL PLOT(A,B,C,TIME,DISEAS,NPTS)
C-------------SOLVE FOR ITERATIVE LEAST-SQUARES ESTIMATES OF A, B, & C.
      CALL LSFIT(TIME,DISEAS,NPTS,A,B,C,KONVRG)
      IF(KONVRG.EQ.1) GO TO 11
      WRITE(60,10)
10    FORMAT(/' $$$$ SUBROUTINE LSFIT DID NOT CONVERGE.')
11    CALL RESID(A,B,C,NPTS,TIME,DISEAS,SUMSQ,RSQ)
      WRITE(60,12) TITLE,SUMSQ,RSQ
12    FORMAT(///A40/'LEAST SQUARES ESTIMATES, SUM OF SQUARES=',
     1 G16.8,2X,' RSQ= ',G10.3)
      CB=C/B
      WRITE(60,7) A,B,C,M,N,CB
      CALL SUMMRY(A,B,C,WMODE,WMEAN,WVAR,ERROR)
      IF (ERROR) THEN
        WRITE(60,*)'$$ ERROR: GAMMA function error in subroutine SUMMRY'
      ENDIF
      WRITE(60,17) WMODE,WMEAN,WVAR
      CALL PLOT(A,B,C,TIME,DISEAS,NPTS)
1000  CONTINUE
      STOP
      END
      DOUBLE PRECISION FUNCTION ESTC(N,M,X)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C-------------THIS ROUTINE ESTIMATES A STARTING VALUE OF C BY STARTING
C             AT C=1 AND STEPING FORWARD UNITL F=D(LOGL)/D(C)
C             BECOMES NEGATIVE.  THE LAST VALUE OF C WHICH CAUSES F TO
C             BE NEGATIVE IS RETURNED ELSE ESTC=1.  THE MAXIMUM VALUE
C             EVER RETURNED IS DEPENDENT ON THE LOOP PARAMETERS IN LOOP
C             100.
C-----------------------------------------------------------------------
      DIMENSION X(M)
      DO 100 IC=1,50
      C=IC
      SUM=0.
      DIF=N-M
      DO 1 I=1,M
      XMAC=X(I)**C
1     SUM=SUM+XMAC
      B=((SUM+DIF*XMAC)/M)**(1.D0/C)
      SUM=0.
      SUMP=0.
      DO 2 I=1,M
      XMAB=X(I)/B
      XLOG=DLOG(XMAB)
      XMABC=XMAB**C
      SUM=SUM+XLOG
2     SUMP=SUMP+XLOG*XMABC
      F=M/C+SUM-SUMP-DIF*XLOG*XMABC
      ESTC=C
      IF(F.LT.0) RETURN
100   CONTINUE
      ESTC=1.D0
      RETURN
      END
      SUBROUTINE NEWTON(B,C,N,M,X,KONVRG)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C-------------NEWTON USES A NEWTON-RAPHSON ITERATION TO REFINE A
C             M.L. ESTIMATE OF C AND TO CALCULATE A M.L. ESTIMATE
C             OF B.  THAT IS THE ROOT C OF D(LOGL)/D(C) IS FOUND
C             WHERE L IS THE LIKELIHOOD FUNCTION OF THE WEIBULL
C             DISTRIBUTION.
C             KONVRG IS RETURNED EQUAL TO 1 IF CONVERGENCE RESULTED,
C             OTHERWISE KONVRG IS RETURNED EQUAL TO ZERO.
C-----------------------------------------------------------------------
      DIMENSION X(M)
      KONVRG=0
      DO 99 NIT=1,20
      SUML2=0
      SUML=0
      SUM2P=0
      SUM=0
      DIF=N-M
      C1=1.D0/C
      C2=C1*C1
      CM=M**(-C1)
      DO 1 I=1,M
      XIA=X(I)
      XMAC=XIA**C
      XLOG=DLOG(XIA)
      XLOG2=XLOG*XLOG
      XMACL=XMAC*XLOG
      SUML=SUML+XLOG
      SUM=SUM+XMAC
      SUML2=SUML2+XLOG2*XMAC
1     SUM2P=SUM2P+XMACL
      B1=SUM+DIF*XMAC
      B2=B1**C1
C-------------B IS THE ML ESTIMATE OF B.
      B=B2*CM
      IF(KONVRG.EQ.1) RETURN
C-------------DB IS B'.
      DB=DLOG(DBLE(M))*B*C2+(-DLOG(B1)*B2*C2+(SUM2P+XMACL*DIF)
     1*(B2/B1)*C1)*CM
C
      DIFX=DIF*XMAC
      BC=B**(-C)
      BC1=B**(-C-1)
      BLOG=DLOG(B)
      BCD=BC1*C*DB
      DBB=DB*BC1
      BLOGC=BLOG*BC
      BCDB=-BCD-BLOGC
      XLOGMB=(XLOG-BLOG)*DIFX
      BLOGM=BLOG*M
C
      F1=SUML-BLOGM
      F2=-(SUM2P*BC - SUM*BLOGC)
      F3=-XLOGMB*BC+M*C1
C-------------F IS D(LOGL)/D(C)
      F=F1+F2+F3
C
      DF1=-DB*M/B
      DF2=-(SUM2P*BCDB + SUM*(BLOG*BCD+BLOG*BLOGC))
      DF3=SUM*DBB
      DF4=-BC*SUML2+SUM2P*BLOGC
      DF5=-XLOGMB*BCDB
      DF6=DBB*DIFX
      DF7=-XLOG*XLOGMB*BC-M*C2
C-------------DF IS F'.
      DF=DF1+DF2+DF3+DF4+DF5+DF6+DF7
C-------------NEWTON-RAPHSON ITERATION.
      DELTA=F/DF
      C=C-DELTA
      IF(DABS(DELTA/C).GT.1.D-13) GO TO 99
      KONVRG=1
99    CONTINUE
      RETURN
      END
      SUBROUTINE INPUT (U,M,T,X,NPTS)
C-------------INPUT CUMULATIVE DATA.
C             INPUT ITEM II-- TIME(INTEGER), DISEASE(PROPORTION)
C             STOP WITH NEGATIVE VALUE FOR TIME.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION U(100),T(100),X(100)
      INTEGER TT
      NPTS=0
      M=100
55    READ(50,*,END=3) TT,XX
      IF(TT.LT.0) GO TO 3
      NPTS=NPTS+1
      T(NPTS)=TT
      X(NPTS)=XX
      GO TO 55
3     MM=0
      XOLD=X(1)
      DO 4 I=2,NPTS
      IF(X(I).EQ.0) GO TO 4
      K=(X(I)-XOLD)*M + .5D0
      TT=T(I)
      IF(K.LT.1) GO TO 4
      DO 5 J=1,K
      MM=MM+1
        IF(MM.GT.100)MM=100
5     U(MM)=TT
      XOLD=X(I)
4     CONTINUE
      M=MM
      RETURN
      END
      SUBROUTINE PLOT(A,B,C,T,X,NPTS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C-------------PLOT CUMULATIVE DATA AND FIT.
      CHARACTER*1 SYMB(2), LINE(103)
      DIMENSION YMIN(2),SF(2),YP(2),YLAB(11)
      DIMENSION T(NPTS),X(NPTS),F(100)
      DATA SYMB/'O','P'/
      XOLD=0.
      DO 1 I=1,NPTS
      IF(T(I).GE.A) GO TO 4
      F(I)=0.
      GO TO 1
4     F(I)=WDIS(T(I),A,B,C)
      XOLD=F(I)
1     CONTINUE
      XMAX=BIG(X,NPTS)
      FMAX=BIG(F,NPTS)
      XFMAX=DMAX1(XMAX,FMAX)
      XFMAX=1.
      DO 2 I=1,2
      YMIN(I)=0.
2     SF(I)=XFMAX/10.
      CALL SCALEN(2,YMIN,SF)
      DO 6 I=1,11
6     YLAB(I)=(I-1)*SF(1)+YMIN(1)
      WRITE(60,7) YLAB
7     FORMAT(/'TIME  OBS(Y)     PRED(Y)',F7.2,10F10.2/
     1 ' ',28X,10('*.........'),'*'/' ',28X,101('-'))
      DO 3 I=1,NPTS
      YP(1)=X(I)
      YP(2)=F(I)
      CALL GRAPHN(2,LINE,YP,SYMB)
3     WRITE(60,5) T(I),X(I),F(I),LINE
5     FORMAT(' ',F4.0,1X,2G10.4,1X,A1,'I',101A1,'I',A1)
      WRITE(60,8)
8     FORMAT(29X,101('-'))
      RETURN
      END
      DOUBLE PRECISION FUNCTION WDEN(X,A,B,C)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C--------------EVALUATE THE THREE PARAMETER WEIBULL DENSITY FUNCTION.
      XMA=X-A
      D=0.
      CM1=C-1
      ARG1= XMA/B
      ARG2=XMA
      IF(DABS(CM1).LT.5.D-14) GO TO 1
      ARG1=ARG1**C
      ARG2=XMA**CM1
1     IF(ARG1.LT.175.) D=DEXP(-ARG1)
      WDEN=(C*ARG2/B**C)*D
      RETURN
      END
      DOUBLE PRECISION FUNCTION BIG(X,N)
C-------------FIND LARGEST NUMBER.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(N)
      BIG=X(1)
      DO 1 I=1,N
      IF(BIG.LT.X(I)) BIG=X(I)
1     CONTINUE
      RETURN
      END
      DOUBLE PRECISION FUNCTION WDIS(X,A,B,C)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C--------------EVALUATE THE THREE PARAMETER WEIBULL DISTRIBUTION
C              FUNCTION.
      ARG=((X-A)/B)**C
      WDIS=1.D0
      IF(ARG.GT.175.) RETURN
      WDIS=1.D0-DEXP(-ARG)
      RETURN
      END
      SUBROUTINE RESID(A,B,C,N,X,Y,SUMSQ,RSQ)
C-------------CALCULATE ERROR SUMS OF SQUARES AND R-SQUARE.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(N),Y(N)
      SUMSQ=0.
      YSQ=0.0
      YT=0.0
      DO 1 I=1,N
      YT=YT+Y(I)
      YSQ=YSQ+ (Y(I)*Y(I))
1     SUMSQ=SUMSQ+(Y(I)-WDIS(X(I),A,B,C))**2
      SSTO=YSQ- ((YT*YT)/N)
      RSQ=1.0D0-(SUMSQ/SSTO)
      RETURN
      END
      SUBROUTINE LSFIT(X,Y,N,A,BB,C,KONVRG)
C-------------ITERATIVE LEAST-SQUARES PROCEDURE TO SOLVE FOR
C             A, B, & C.  UTILIZES TAYLOR SERIES EXPANSION.
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION GPG(3,4),G(100,3),B(3),DB(3),DBL(3),BN(3),YN(100)
      DIMENSION X(N),Y(N)
      F(Z,B1,B2,B3) = 1-DEXP(-((Z-B1)/B2)**B3)
      F1(Z,B1,B2,B3) = -B3/B2*((Z-B1)/B2)**(B3-1.)*DEXP(-((Z-B1)
     1                 /B2)**B3)
      F2(Z,B1,B2,B3) = -B3/B2*((Z-B1)/B2)**B3*DEXP(-((Z-B1)
     1                 /B2)**B3)
      F3(Z,B1,B2,B3) = ((Z-B1)/B2)**B3*DLOG((Z-B1)/B2)*DEXP
     1                 (-((Z-B1)/B2)**B3)
      D(V,U,B1,B2,B3) = V - F(U, B1, B2, B3)
      NP=3
      KONVRG=0
      B(1)=A
      B(2)=BB
      B(3)=C
      DO 1 I=1,3
1     DBL(I)=1.E-5
      U1 = 0
      KOUNT=0
      DO 6 I = 1, N
    6 U1 = U1 + D(Y(I),X(I),B(1),B(2),B(3))**2
   55 DO 9 I = 1, N
      G(I,1) = F1( X(I), B(1), B(2), B(3))
      G(I,2) = F2(X(I), B(1), B(2), B(3))
      G(I,3) = F3(X(I), B(1), B(2), B(3))
    9 YN(I) = Y(I) - F(X(I),B(1),B(2),B(3))
      FACTOR = 1
      DO 13 I = 1, NP
      DO 12 J = 1, NP
      GPG(I,J) = 0
      DO 12 K = 1, N
   12 GPG (I,J) = GPG(I,J) + G(K,I)*G(K,J)
   13 CONTINUE
      NPGY = NP + 1
      DO 15 I = 1, NP
      GPG(I,NPGY) = 0
      DO 15 K = 1, N
   15 GPG(I,NPGY) = GPG(I,NPGY) + G(K,I)*YN(K)
      DO 22 I = 1, NP
      I1 = I + 1
      Q = 1.D0/GPG(I,I)
      DO 21 J = I1, NPGY
   21 GPG(I,J) = Q*GPG(I,J)
      DO 22 K = 1, NP
      IF (K - I) 23, 22, 23
   23 DO 24 J = I1, NPGY
   24 GPG(K,J) = GPG(K,J) - GPG(I,J)*GPG(K,I)
   22 CONTINUE
   29 DO 25 I = 1, NP
      DB(I) = FACTOR*GPG(I,NPGY)
   25 BN(I) = B(I) + DB(I)
      U2 = 0
      IF (BN(2).LT.0.) BN(2)=.5D0
      IF (BN(1) .GT.X(1)) BN(1)=X(1)-.01
      IF(BN(3).LT.1) BN(3)=1.
      DO 26 I = 1, N
   26 U2 = U2 + D(Y(I),X(I),BN(1),BN(2),BN(3))**2
      IF (U2 .LT. U1) GO TO 30
      FACTOR = FACTOR*.5D0
      IF (FACTOR .GE. 1.0D-13) GO TO 29
      WRITE (60,28) FACTOR
   28 FORMAT (/' .......FACTOR LESS THAN 10**(-13) = ',G16.5)
      GO TO 200
30    DO 35 I = 1, NP
      IF (DABS(DB(I)) .GE. DBL(I)) GO TO 45
   35 CONTINUE
      KONVRG=1
      GO TO 200
   45 KOUNT = KOUNT + 1
      IF(KOUNT.GT.50) GO TO 200
   50 U1 = U2
      DO 51 I = 1, NP
   51 B(I) = BN(I)
      GO TO 55
200   A=B(1)
      BB=B(2)
      C=B(3)
      RETURN
      END
      SUBROUTINE SUMMRY(A,B,C,WMODE,WMEAN,WVAR,ERROR)
C-------------CALCULATE DENSITY FUNCTION MODE, MEAN, & VARIANCE.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL ERROR1, ERROR2, ERROR
      WMODE=A + B*(DABS(1.0D0-(1.0D0/C)))**(1.0/C)
      G1=DGAMMA(1.0D0 +(1.0D0/C),ERROR1)
      G2=DGAMMA(1.0D0 +(2.0D0/C),ERROR2)
      ERROR=ERROR1 .OR. ERROR2
      WMEAN=A +B*G1
      WVAR=B*B*(G2- G1*G1)
      RETURN
      END
      SUBROUTINE GRAPHN(N,CHAR,YC,CH)
      CHARACTER*1 CH(N),CHAR(103),Z
      DOUBLE PRECISION YC(N),SF1(10),Y1(10)
      COMMON /GRAPH/ SF1, Y1
      DATA Z/' '/
      DO 10 I=1,103
 10   CHAR(I)=Z
 20   DO 60 I=1,N
      DI=(YC(I)-Y1(I))*SF1(I) + 2.0
      IF(DI.LT.1.0D0) GO TO 99
      IS=DI
      IF((DI-IS).GE.0.5) IS=IS+1
 30   IF(IS.LE.102)GO TO 59
      IS=103
      GO TO 59
 99   IS=1
 59   CHAR(IS)=CH(I)
 60   CONTINUE
      RETURN
      END
      SUBROUTINE SCALEN(N,Y1,SF)
      DOUBLE PRECISION SF(N), Y1(N), SF1(10), YY1(10)
      COMMON /GRAPH/ SF1, YY1
      DO 1 I=1,N
      YY1(I)=Y1(I)
 1    SF1(I)=10.0D0/SF(I)
      RETURN
      END
      DOUBLE PRECISION FUNCTION DGAMMA(Z,ERROR)
C=== ACM CALGO Algorithm 309, Double Precision Gamma Function.
C    Z is the argument; must be > 0 and < 57 for IBM mainframes, or
C                                       < 171 for IBM PC.
C    MSIZE is the number of decimal digits representable; 15 for both.
C    ERROR is returned LOGICAL; .TRUE. if argument was out of range, or
C                               .FALSE. otherwise (i.e. no error).
C    Translated from Algol 68, March 6, 1989, HDK.
C
      DOUBLE PRECISION X,Y,Z,PI, DELTA,EX,DECIMA,MAXZ,ARGUND,MAXEXP
      DOUBLE PRECISION DLOGAM
      INTEGER MSIZE,PARITY,ENTIER
      LOGICAL ERROR
C---IBM Mainframes:
C     DATA MAXZ/57.D0/, ARGUND/-1.D-74/, MAXEXP/174.D0/, MSIZE/15/
C---IBM PC.
      DATA MAXZ/171.D0/, ARGUND/-1.D-307/, MAXEXP/702.D0/, MSIZE/15/
      DATA PI/3.1415926535897932384626433832795028841971693993751D0/
C     ENTIER(Z)=FLOOR(Z)
      ENTIER(X)=DNINT(X)-DMOD(2.D0+DSIGN(1.D0,X),3.D0)
      PARITY(Z)=MOD(ENTIER(Z),2)
      DECIMA(Z)=DABS(Z-IDINT(Z))
      IF(Z.GT.MAXZ) THEN
        GOTO 99
      ELSE
        IF(Z.EQ.ENTIER(Z)) THEN
          IF(Z.LE.0) THEN
            GOTO 99
          ELSE
            Y=1.D0
            IF(Z.GT.2.D0) THEN
1             Z=Z-1
              Y=Y*Z
              IF(Z.GT.2.D0) GOTO 1
            ENDIF
          ENDIF
        ELSE
          IF(DABS(Z).LT.10**(-MSIZE)) THEN
            Y=1.D0/Z
          ELSE
            IF(Z.LT.0) THEN
              IF(Z.LT.ARGUND) THEN
                Y=0
              ELSE
                DELTA=DECIMA(Z)*PI
                IF(DELTA.LT.10**(-MSIZE*0.5D0)) THEN
                  EX=-DLOG(DECIMA(Z))
                ELSE
                  EX=DLOG(PI/DSIN(DELTA)) - DLOGAM(1.D0-Z, MSIZE)
                ENDIF
                IF(DABS(EX).GT.MAXEXP) THEN
                  Y=0
                ELSE
                  IF(PARITY(Z).EQ.1) THEN
                    Y=DEXP(EX)
                  ELSE
                    Y=-DEXP(EX)
                  ENDIF
                ENDIF
              ENDIF
            ELSE
              Y=DEXP(DLOGAM(Z,MSIZE))
            ENDIF
          ENDIF
        ENDIF
      ENDIF
      ERROR=.FALSE.
      DGAMMA=Y
      RETURN
99    ERROR=.TRUE.
      DGAMMA=0
      RETURN
      END
      DOUBLE PRECISION FUNCTION DLOGAM(T,MSIZE)
C===Log base e of GAMMA(T). See DGAMMA.
      DOUBLE PRECISION F,T,DLGM
      INTEGER TMIN,MSIZE
      TMIN=7
      IF(MSIZE.GT.18) TMIN=MSIZE-10
      IF(T.GT.TMIN) THEN
        DLOGAM=DLGM(T)
      ELSE
        F=T
1       T=T+1
        IF(T.LT.TMIN) THEN
          F=F*T
          GOTO 1
        ENDIF
        DLOGAM=DLGM(T)-DLOG(F)
      ENDIF
      RETURN
      END
      DOUBLE PRECISION FUNCTION DLGM(W)
C===Companion function for DGAMMA and DLOGAM.
      DOUBLE PRECISION C(20), W,W2,PRESUM,CONST,DEN,SUM
      INTEGER I
      C(1)=1.D0/12.D0
      C(2)=-1.D0/360.D0
      C(3)=1.D0/1260.D0
      C(4)=-1.D0/1680.D0
      C(5)=1.D0/1188.D0
      C(6)=-691.D0/360360.D0
      C(7)=1.D0/156.D0
      C(8)=-3617.D0/122400.D0
      C(9)=43867.D0/244188.D0
      C(10)=-174611.D0/125400.D0
      C(11)=77683.D0/5796.D0
      C(12)=-236364091.D0/1506960.D0
      C(13)=657931.D0/300.D0
      C(14)=-3392780147.D0/93960.D0
      C(15)=1723168255201.D0/2492028.D0
      C(16)=-7709321041217.D0/505920.D0
      C(17)=1516286977551.D0/396.D0
      C(18)=-26315271553053477373.D0/2418179400.D0
      C(19)=1542102059911661.D0/444.D0
      C(20)=-261082718496449122051.D0/21106800.D0
      CONST=.91893853320467274178032973640561763986139747363778D0
C---CONST=DLOG(SQRT(2*PI))
      DEN=W
      W2=W*W
      PRESUM=(W-.5D0)*DLOG(W)-W+CONST
      DO 1 I=1,20
        SUM=PRESUM+C(I)/DEN
        IF(SUM.EQ.PRESUM) GOTO 2
        DEN=DEN*W2
1     PRESUM=SUM
2     DLGM=SUM
      RETURN
      END