SUBROUTINE GAMPLT(X,N,GAMMA)
C
C     PURPOSE--THIS SUBROUTINE GENERATES A GAMMA
C              PROBABILITY PLOT
C              (WITH TAIL LENGTH PARAMETER VALUE = GAMMA).
C              THE PROTOTYPE GAMMA DISTRIBUTION USED
C              HEREIN HAS MEAN = GAMMA
C              AND STANDARD DEVIATION = SQRT(GAMMA).
C              THIS DISTRIBUTION IS DEFINED FOR ALL POSITIVE X,
C              AND HAS THE PROBABILITY DENSITY FUNCTION
C              F(X) = (1/CONSTANT) * (X**(GAMMA-1)) * EXP(-X)
C              WHERE THE CONSTANT = THE GAMMA FUNCTION EVALUATED
C              AT THE VALUE GAMMA.
C              AS USED HEREIN, A PROBABILITY PLOT FOR A DISTRIBUTION
C              IS A PLOT OF THE ORDERED OBSERVATIONS VERSUS 
C              THE ORDER STATISTIC MEDIANS FOR THAT DISTRIBUTION.
C              THE GAMMA PROBABILITY PLOT IS USEFUL IN
C              GRAPHICALLY TESTING THE COMPOSITE (THAT IS,
C              LOCATION AND SCALE PARAMETERS NEED NOT BE SPECIFIED)
C              HYPOTHESIS THAT THE UNDERLYING DISTRIBUTION
C              FROM WHICH THE DATA HAVE BEEN RANDOMLY DRAWN 
C              IS THE  GAMMA DISTRIBUTION
C              WITH TAIL LENGTH PARAMETER VALUE = GAMMA.
C              IF THE HYPOTHESIS IS TRUE, THE PROBABILITY PLOT
C              SHOULD BE NEAR-LINEAR.
C              A MEASURE OF SUCH LINEARITY IS GIVEN BY THE
C              CALCULATED PROBABILITY PLOT CORRELATION COEFFICIENT.
C     INPUT  ARGUMENTS--X      = THE SINGLE PRECISION VECTOR OF
C                                (UNSORTED OR SORTED) OBSERVATIONS.
C                     --N      = THE INTEGER NUMBER OF OBSERVATIONS
C                                IN THE VECTOR X. 
C                     --GAMMA  = THE SINGLE PRECISION VALUE OF THE
C                                TAIL LENGTH PARAMETER.
C                                GAMMA SHOULD BE POSITIVE.
C     OUTPUT--A ONE-PAGE GAMMA PROBABILITY PLOT.
C     PRINTING--YES.
C     RESTRICTIONS--THE MAXIMUM ALLOWABLE VALUE OF N
C                   FOR THIS SUBROUTINE IS 7500.
C                 --GAMMA SHOULD BE POSITIVE.
C     OTHER DATAPAC   SUBROUTINES NEEDED--SORT, UNIMED, PLOT.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT, ABS, EXP, DEXP, DLOG. 
C     MODE OF INTERNAL OPERATIONS--SINGLE PRECISION AND DOUBLE PRECISION
C     LANGUAGE--ANSI FORTRAN. 
C     REFERENCES--WILK, GNANADESIKAN, AND HUYETT, 'PROBABILITY
C                 PLOTS FOR THE GAMMA DISTRIBUTION',
C                 TECHNOMETRICS, 1962, PAGES 1-15.
C               --NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS
C                 SERIES 55, 1964, PAGE 257, FORMULA 6.1.41.
C               --FILLIBEN, 'TECHNIQUES FOR TAIL LENGTH ANALYSIS',
C                 PROCEEDINGS OF THE EIGHTEENTH CONFERENCE
C                 ON THE DESIGN OF EXPERIMENTS IN ARMY RESEARCH
C                 DEVELOPMENT AND TESTING (ABERDEEN, MARYLAND,
C                 OCTOBER, 1972), PAGES 425-450.
C               --HAHN AND SHAPIRO, STATISTICAL METHODS IN ENGINEERING,
C                 1967, PAGES 260-308.
C               --JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE
C                 DISTRIBUTIONS--1, 1970, PAGES 166-206.
C     WRITTEN BY--JAMES J. FILLIBEN
C                 STATISTICAL ENGINEERING LABORATORY (205.03)
C                 NATIONAL BUREAU OF STANDARDS
C                 WASHINGTON, D. C. 20234
C                 PHONE:  301-921-2315
C     ORIGINAL VERSION--NOVEMBER  1974. 
C     UPDATED         --SEPTEMBER 1975. 
C     UPDATED         --NOVEMBER  1975. 
C     UPDATED         --FEBRUARY  1976. 
C
C---------------------------------------------------------------------
C
      DOUBLE PRECISION Z,Z2,Z3,Z4,Z5,DEN,A,B,C,D
      DOUBLE PRECISION DEXP,DLOG
      DIMENSION D(10)
      DIMENSION X(1)
      DIMENSION Y(7500),W(7500)
      COMMON /BLOCK2/ WS(15000)
      EQUIVALENCE (Y(1),WS(1)),(W(1),WS(7501))
      DATA C/ .918938533204672741D0/
      DATA D(1),D(2),D(3),D(4),D(5)
     1                 /+.833333333333333333D-1,-.277777777777777778D-2,
     1+.793650793650793651D-3,-.595238095238095238D-3,+.8417508417508417
     151D-3/
      DATA D(6),D(7),D(8),D(9),D(10)
     1     /-.191752691752691753D-2,+.641025641025641025D-2,-.2955065359
     147712418D-1,+.179644372368830573D0,-.139243221690590111D1/
C
      IPR=6
      IUPPER=7500
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      IF(N.LT.1.OR.N.GT.IUPPER)GOTO50
      IF(N.EQ.1)GOTO55
      IF(GAMMA.LE.0.0)GOTO60
      HOLD=X(1)
      DO65I=2,N
      IF(X(I).NE.HOLD)GOTO90
   65 CONTINUE
      WRITE(IPR, 9)HOLD
      RETURN
   50 WRITE(IPR,17)IUPPER
      WRITE(IPR,47)N
      RETURN
   55 WRITE(IPR,18) 
      RETURN
   60 WRITE(IPR,25) 
      WRITE(IPR,46)GAMMA
      RETURN
   90 CONTINUE
    9 FORMAT(1H ,109H***** NON-FATAL DIAGNOSTIC--THE FIRST  INPUT ARGUME
     1NT (A VECTOR) TO THE GAMPLT SUBROUTINE HAS ALL ELEMENTS = ,E15.8,6
     1H *****)
   17 FORMAT(1H , 98H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE
     1 GAMPLT SUBROUTINE IS OUTSIDE THE ALLOWABLE (1,,I6,16H) INTERVAL *
     1****)
   18 FORMAT(1H ,100H***** NON-FATAL DIAGNOSTIC--THE SECOND INPUT ARGUME
     1NT TO THE GAMPLT SUBROUTINE HAS THE VALUE 1 *****)
   25 FORMAT(1H , 91H***** FATAL ERROR--THE THIRD  INPUT ARGUMENT TO THE
     1 GAMPLT SUBROUTINE IS NON-POSITIVE *****)
   46 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,E15.8,6H *****)
   47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8   ,6H *****)
C
C-----START POINT-----------------------------------------------------
C
      AN=N
      DGAMMA=GAMMA
C
C     COMPUTE THE GAMMA FUNCTION USING THE ALGORITHM IN THE 
C     NBS APPLIED MATHEMATICS SERIES REFERENCE.
C     THIS GAMMA FUNCTION NEED BE CALCULATED ONLY ONCE.
C     IT IS USED IN THE CALCULATION OF THE CDF BASED ON
C     THE TENTATIVE VALUE OF THE PPF IN THE ITERATION.
C
      Z=DGAMMA
      DEN=1.0D0
  150 IF(Z.GE.10.0D0)GOTO160
      DEN=DEN*Z
      Z=Z+1.0D0
      GOTO150
  160 Z2=Z*Z
      Z3=Z*Z2
      Z4=Z2*Z2
      Z5=Z2*Z3
      A=(Z-0.5D0)*DLOG(Z)-Z+C 
      B=D(1)/Z+D(2)/Z3+D(3)/Z5+D(4)/(Z2*Z5)+D(5)/(Z4*Z5)+
     1D(6)/(Z*Z5*Z5)+D(7)/(Z3*Z5*Z5)+D(8)/(Z5*Z5*Z5)+D(9)/(Z2*Z5*Z5*Z5)
      G=DEXP(A+B)/DEN
C
C     SORT THE DATA 
C
      CALL SORT(X,N,Y)
C
C     GENERATE UNIFORM ORDER STATISTIC MEDIANS
C
      CALL UNIMED(N,W)
C
C     GENERATE GAMMA DISTRIBUTION ORDER STATISTIC MEDIANS
C
C     DETERMINE LOWER AND UPPER BOUNDS ON THE DESIRED I-TH GAMMA
C     ORDER STATISTIC MEDIAN. 
C     FOR EACH I, A LOWER BOUND IS GIVEN BY
C     (Y(I)*GAMMA*THE GAMMA FUNCTION OF GAMMA)**(1.0/GAMMA) 
C     WHERE Y(I) IS THE CORRESPONDING UNIFORM (0,1) ORDER STATISIC
C     MEDIAN.
C     FOR EACH I EXCEPT I = N, AN UPPER BOUND IS GIVEN BY THE
C     (I+1)-ST GAMMA ORDER STATISTIC MEDIAN (ASSUMEDLY ALREADY
C     CALCULTATED). 
C     FOR I = N, AN UPPER BOUND IS DETERMINED BY COMPUTING
C     MULTIPLES OF THE LOWER BOUND FOR I = N UNTIL A LARGER 
C     VALUE IS OBTAINED.
C     DUE TO THE ABOVE CONSIDERATIONS, THE GAMMA ORDER STATISTIC
C     MEDIANS WILL BE CALCULATED LARGEST TO SMALLEST, THAT IS,
C     IN THE FOLLOWING SEQUENCE:  W(N), W(N-1), ..., W(2), W(1).
C     NOTE ALSO THAT 1) THE CODE IS COMPLICATED SLIGHTLY BY THE
C     FACT THAT PERCENT POINT VALUES INVOLVED IN THE CALCULATION OF
C     THE TAIL LENGTH MEASURE TAU (SEE LABEL 605) ARE GOING ON
C     'SIMULATNEOUSLY'. AND 2) THE VECTOR W WILL AT VARIOUS TIMES
C     IN THE PROGRAM HAVE UNIFORM ORDER STATISTIC MEDIANS AND
C     THEN LATER GRADUALLY FILL UP WITH GAMMA ORDER STATISTIC
C     MEDIANS.
C
      I=N 
      ITAIL=0
  310 IF(ITAIL.EQ.0)U=W(I)
      DP=U
      XMIN0=(U*GAMMA*G)**(1.0/GAMMA)
      XMIN=XMIN0
      IF(I.EQ.N.OR.ITAIL.GE.1)GOTO320
      IP1=I+1
      XMAX=W(IP1)
      GOTO370
  320 ILOOP=1
      ICOUNT=1
  350 ACOUNT=ICOUNT 
      XMAX=ACOUNT*XMIN0
      DX=XMAX
      GOTO1000
  360 IF(PCALC.GE.DP)GOTO370
      XMIN=XMAX
      ICOUNT=ICOUNT+1
      IF(ICOUNT.LE.30000)GOTO350
  370 XMID=(XMIN+XMAX)/2.0
C
C     AT THIS STAGE WE NOW HAVE LOWER AND UPPER LIMITS ON
C     THE DESIRED I-TH GAMMA ORDER STATISITC MEDIAN W(I).
C     NOW ITERATE BY BISECTION UNTIL THE DESIRED ACCURACY IS ACHIEVED 
C     FOR THE I-TH GAMMA ORDER STATISITIC MEDIAN. 
C
      ILOOP=2
      XLOWER=XMIN
      XUPPER=XMAX
      ICOUNT=0
  550 DX=XMID
      GOTO1000
  560 IF(PCALC.EQ.DP)GOTO570
      IF(PCALC.GT.DP)GOTO580
      XLOWER=XMID
      XMID=(XMID+XUPPER)/2.0
      GOTO590
  580 XUPPER=XMID
      XMID=(XMID+XLOWER)/2.0
  590 XDEL=ABS(XMID-XLOWER)
      ICOUNT=ICOUNT+1
      IF(XDEL.LT.0.0000001.OR.ICOUNT.GT.100)GOTO570
      GOTO550
  570 IF(ITAIL.GE.1)GOTO605
      W(I)=XMID
      IF(I.LE.1)GOTO595
      I=I-1
      GOTO310
  595 CONTINUE
C
C     AT THIS POINT, THE GAMMA ORDER STATISTIC MEDIANS ARE ALL COMPUTED.
C     NOW PLOT OUT THE GAMMA PROBABILITY PLOT
C
      CALL PLOT(Y,W,N)
C
C     COMPUTE THE TAIL LENGTH MEASURE OF THE DISTRIBUTION.
C     WRITE OUT THE TAIL LENGTH MEASURE OF THE DISTRIBUTION 
C     AND THE SAMPLE SIZE.
C
  605 IF(ITAIL.EQ.0)GOTO600
      IF(ITAIL.EQ.1)GOTO610
      IF(ITAIL.EQ.2)GOTO620
      IF(ITAIL.EQ.3)GOTO630
      GOTO640
  600 U=.9975
      ITAIL=1
      GOTO310
  610 PP9975=XMID
      U=.0025
      ITAIL=2
      GOTO310
  620 PP0025=XMID
      U=.975
      ITAIL=3
      GOTO310
  630 PP975=XMID
      U=.025
      ITAIL=4
      GOTO310
  640 PP025=XMID
      TAU=(PP9975-PP0025)/(PP975-PP025) 
      WRITE(IPR,655)GAMMA,TAU,N
C
C     COMPUTE THE PROBABILITY PLOT CORRELATION COEFFICIENT. 
C     COMPUTE LOCATION AND SCALE ESTIMATES
C     FROM THE INTERCEPT AND SLOPE OF THE PROBABILITY PLOT. 
C     THEN WRITE THEM OUT.
C
      SUM1=0.0
      SUM2=0.0
      DO660I=1,N
      SUM1=SUM1+Y(I)
      SUM2=SUM2+W(I)
  660 CONTINUE
      YBAR=SUM1/AN
      WBAR=SUM2/AN
      SUM1=0.0
      SUM2=0.0
      SUM3=0.0
      DO670I=1,N
      SUM1=SUM1+(Y(I)-YBAR)*(Y(I)-YBAR) 
      SUM2=SUM2+(Y(I)-YBAR)*(W(I)-WBAR) 
      SUM3=SUM3+(W(I)-WBAR)*(W(I)-WBAR) 
  670 CONTINUE
      CC=SUM2/SQRT(SUM3*SUM1) 
      YSLOPE=SUM2/SUM3
      YINT=YBAR-YSLOPE*WBAR
      WRITE(IPR,675)CC,YINT,YSLOPE
C
  655 FORMAT(1H ,46HGAMMA PROBABILITY PLOT WITH SHAPE PARAMETER = ,
     1E17.10,1X,7H(TAU = ,E15.8,1H),16X,16HSAMPLE SIZE N = ,I7)
  675 FORMAT(1H ,43HPROBABILITY PLOT CORRELATION COEFFICIENT = ,F8.5,5X,
     122HESTIMATED INTERCEPT = ,E15.8,3X,18HESTIMATED SLOPE = ,E15.8) 
C
      RETURN
C
C******************************************************************** 
C     THIS SECTION BELOW IS LOGICALLY SEPARATE FROM THE ABOVE.
C     THIS SECTION COMPUTES A CDF VALUE FOR ANY GIVEN TENTATIVE
C     PERCENT POINT X VALUE AS DEFINED IN EITHER OF THE 2
C     ITERATION LOOPS IN THE ABOVE CODE.
C
C     COMPUTE T-SUB-Q AS DEFINED ON PAGE 4 OF THE WILK, GNANADESIKAN, 
C     AND HUYETT REFERENCE
C
 1000 SUM=1.0/DGAMMA
      TERM=1.0/DGAMMA
      CUT1=DX-DGAMMA
      CUT2=DX*10000000.0
      DO700J=1,1000 
      AJ=J
      TERM=DX*TERM/(DGAMMA+AJ)
      SUM=SUM+TERM
      CUTOFF=CUT1+(CUT2*TERM/SUM)
      IF(AJ.GT.CUTOFF)GOTO750 
  700 CONTINUE
      WRITE(IPR,705)
      WRITE(IPR,707)GAMMA
  705 FORMAT(1H ,48H*****ERROR IN INTERNAL OPERATIONS IN THE GAMPLT , 
     1 53HSUBROUTINE--THE NUMBER OF CDF ITERATIONS EXCEEDS 1000)
  707 FORMAT(1H ,33H     THE INPUT VALUE OF GAMMA IS ,E15.8)
  750 T=SUM
      PCALC=(DX**DGAMMA)*(EXP(-DX))*T/G 
      IF(ILOOP.EQ.1)GOTO360
      GOTO560
C
      END